00001 module ice_probability
00002
00003 use ice_kinds_mod
00004 use ice_communicate, only: my_task, master_task, get_num_procs
00005 use ice_constants, only: pi, puny, c1
00006 use ice_fileunits, only: nu_timing
00007 use ice_domain_size, only: nx_global, ny_global, max_blocks
00008 use ice_domain, only: nblocks, blocks_ice, distrb_info, distribution_wght_file
00009 use ice_blocks, only: nx_block, ny_block, nblocks_tot, block, get_block
00010 use ice_read_write, only: ice_open_nc, ice_read_global_nc, ice_close_nc
00011 use ice_broadcast, only: broadcast_array
00012 use ice_global_reductions, only: sum_vector_dbl
00013 use ice_timers, only:
00014 use ice_work, only: work_gr
00015 use ice_probability_tools
00016
00017 implicit none
00018
00019 integer (int_kind), public, parameter ::
00020 lndType = 0,
00021 icefreeType = 1,
00022 iceType = 2
00023
00024
00025 integer (int_kind), public :: dynCnt
00026 real (dbl_kind), allocatable :: lnumIceCells(:)
00027
00028 public :: ReadProbabilityFile
00029
00030 public :: CalcWorkPerBlock
00031
00032
00033
00034 public :: init_numIceCells, &
00035 accum_numIceCells, &
00036 accum_numIceCells2, &
00037 print_numIceCells
00038
00039 public :: set_numIceCells,write_numIceCells
00040
00041
00042 contains
00043
00044 subroutine init_numIceCells()
00045
00046 dynCnt = 0
00047
00048
00049 allocate(lnumIceCells(nblocks_tot))
00050 lnumIceCells = 0.0
00051
00052 end subroutine init_numIceCells
00053
00054 subroutine accum_numIceCells(iblk,icells)
00055
00056 integer (int_kind) :: iblk,icells
00057
00058 lnumIceCells(iblk) = lnumIceCells(iblk) + real(icells,kind=dbl_kind)
00059
00060 end subroutine accum_numIceCells
00061
00062 subroutine accum_numIceCells2(aice)
00063
00064 real (dbl_kind) :: aice(nx_block,ny_block,max_blocks)
00065
00066 integer (int_kind) :: igblk,iblk
00067
00068 real (dbl_kind) :: tmp
00069 type (block) :: this_block
00070 integer (int_kind) :: i,j,ihi,ilo,jhi,jlo
00071
00072
00073 do iblk = 1,nblocks
00074 igblk = blocks_ice(iblk)
00075 this_block = get_block(igblk,iblk)
00076 ilo = this_block%ilo
00077 ihi = this_block%ihi
00078 jlo = this_block%jlo
00079 jhi = this_block%jhi
00080 tmp = 0.0
00081 do j = jlo,jhi
00082 do i = ilo,ihi
00083 if(aice(i,j,iblk) > puny) then
00084 tmp = tmp + c1
00085 endif
00086 enddo
00087 enddo
00088 lnumIceCells(igblk) = lnumIceCells(igblk) + tmp
00089 enddo
00090
00091 end subroutine accum_numIceCells2
00092
00093
00094
00095 subroutine set_numIceCells(iblk,ncells)
00096 integer (int_kind) :: iblk
00097 real (dbl_kind) :: ncells
00098
00099 lnumIceCells(iblk) = ncells
00100
00101 end subroutine set_numIceCells
00102
00103 subroutine write_numIceCells
00104 real (dbl_kind), allocatable :: gnumIceCells(:)
00105 integer :: n
00106
00107 allocate(gnumIceCells(nblocks_tot))
00108
00109 call sum_vector_dbl(lnumIceCells,gnumIceCells,distrb_info)
00110
00111 gnumIceCells=gnumIceCells/real(dynCnt,kind=dbl_kind)
00112 if(my_task == master_task) then
00113 open(nu_timing,file='numCells2.bin',recl=8*nblocks_tot, &
00114 form = 'unformatted', access = 'direct', status = 'unknown')
00115 write(nu_timing,rec=1) gnumIceCells
00116 close(nu_timing)
00117 endif
00118
00119
00120
00121 end subroutine write_numIceCells
00122
00123 subroutine print_numIceCells
00124
00125 real (dbl_kind), allocatable :: gnumIceCells(:)
00126 real (dbl_kind), allocatable :: gnumIceCells2(:)
00127 integer :: ii,n
00128
00129 allocate(gnumIceCells(nblocks_tot))
00130 allocate(gnumIceCells2(nblocks_tot))
00131
00132 call sum_vector_dbl(lnumIceCells, gnumIceCells, distrb_info)
00133
00134 gnumIceCells = gnumIceCells/real(dynCnt,kind=dbl_kind)
00135 if(my_task == master_task) then
00136
00137 ii =0
00138 do n=1,nblocks_tot
00139 if(nocn(n) > 0) then
00140 ii = ii+1
00141 gnumIceCells2(ii) = gnumIceCells(n)
00142 endif
00143 enddo
00144 open(nu_timing,file='numCells.bin',recl=8*ii, &
00145 form = 'unformatted', access = 'direct', status = 'unknown')
00146 write(nu_timing,rec=1) gnumIceCells2(1:ii)
00147 close(nu_timing)
00148 print *,'numCells: ',gnumIceCells2(1:ii)
00149 endif
00150
00151 deallocate(gnumIceCells,gnumIceCells2)
00152
00153
00154 end subroutine print_numIceCells
00155
00156 subroutine ReadProbabilityFile(distribution_wght_file,Prob)
00157
00158 character(char_len_long), intent(in) :: distribution_wght_file
00159 real(real_kind), intent(inout) :: Prob(:,:)
00160
00161 type(block) :: this_block
00162 integer(int_kind) :: ilo,ihi
00163 integer(int_kind) :: jlo,jhi
00164
00165 integer(int_kind) :: fid_prob
00166 integer(int_kind) :: amode,ncid,iostat,varid
00167 integer(int_kind) :: i,j,n,ncnt,ierr,ig,jg
00168 real(dbl_kind) :: val,sum,avg
00169
00170 character (char_len) ::
00171 fieldname
00172
00173
00174 call ice_open_nc(distribution_wght_file,fid_prob)
00175
00176
00177 fieldname = 'ice_present'
00178 call ice_read_global_nc(fid_prob,1,fieldname,Prob,.true.)
00179
00180 if(my_task == master_task) then
00181 call ice_close_nc(fid_prob)
00182 print *,'MAXVAL(Prob): ',MAXVAL(Prob)
00183 print *,'MINVAL(Prob): ',MINVAL(Prob)
00184 print *,'COUNT(Prob > 0.5): ',COUNT(Prob > 0.5)
00185 endif
00186
00187
00188
00189
00190 end subroutine ReadProbabilityFile
00191
00192
00193
00194
00195 subroutine CalcWorkPerBlock(distribution_wght, KMTG,ULATG,work_per_block, prob_per_block,blockType,bStats)
00196
00197 character (char_len), intent(in) :: distribution_wght
00198
00199 real (dbl_kind), dimension(nx_global,ny_global), intent(in) ::
00200 KMTG ,
00201 ULATG
00202
00203 integer (int_kind), intent(inout) ::
00204 work_per_block(nblocks_tot)
00205
00206 real (dbl_kind), intent(inout) ::
00207 prob_per_block(nblocks_tot)
00208
00209 integer (int_kind), intent(inout) ::
00210 blockType(nblocks_tot)
00211
00212
00213
00214 real (dbl_kind), intent(inout) :: bStats(:,:)
00215
00216
00217 type (block) ::
00218 this_block
00219
00220 integer (int_kind), parameter ::
00221 max_work_unit=10
00222
00223 integer(int_kind) :: jg,j,n,ig,i,work_unit
00224 integer(int_kind) :: maxkmt,mkmt,ilo,jlo,ihi,jhi
00225 integer(int_kind) :: activePts
00226
00227 integer(int_kind) :: bid(max_blocks)
00228 integer(int_kind) :: bsize_x,bsize_y
00229
00230 integer(int_kind), allocatable :: iLocation(:)
00231
00232 real(dbl_kind), allocatable, dimension(:) :: prob,work
00233
00234 real(dbl_kind) :: lat
00235 real(dbl_kind), parameter :: w0 = 1.,
00236 w1 = 10.
00237
00238 real(dbl_kind) :: maxlat,maxalat
00239 real(dbl_kind) :: shlatT,nhlatT
00240 real(real_kind) :: tmp
00241
00242 integer(int_kind) :: ActiveBlocks,totWork
00243
00244 logical, parameter :: Debug = .FALSE.
00245
00246 integer(int_kind) :: numLnd,numIcefree,numIce
00247
00248 integer(int_kind) :: ierr
00249
00250
00251
00252
00253
00254
00255
00256
00257
00258
00259 allocate(work_gr(nx_global,ny_global))
00260 allocate(prob(nblocks_tot),work(nblocks_tot))
00261 allocate(nocn(nblocks_tot))
00262 allocate(nice005(nblocks_tot),nice010(nblocks_tot),nice050(nblocks_tot), &
00263 nice100(nblocks_tot),nice250(nblocks_tot),nice500(nblocks_tot))
00264
00265 nocn=0
00266 nice005=0
00267 nice010=0
00268 nice050=0
00269 nice100=0
00270 nice250=0
00271 nice500=0
00272
00273 ActiveBlocks = 0
00274 select case (distribution_wght)
00275 case('file')
00276 call ReadProbabilityFile(distribution_wght_file,work_gr)
00277 case('erfc')
00278 work_gr = ErfcProbability(ULATG)
00279 case default
00280 work_gr = ErfcProbability(ULATG)
00281 end select
00282
00283
00284
00285
00286
00287
00288
00289
00290
00291
00292
00293
00294
00295
00296
00297
00298 if(my_task == master_task) then
00299 work_per_block=0
00300 do n=1,nblocks_tot
00301
00302 this_block=get_block(n,n)
00303 ilo = this_block%ilo;ihi = this_block%ihi
00304 jlo = this_block%jlo;jhi = this_block%jhi
00305
00306
00307
00308 tmp = 0.0
00309 do j=jlo,jhi
00310 jg=this_block%j_glob(j)
00311 if(jg>0) then
00312 do i=ilo,ihi
00313 ig = this_block%i_glob(i)
00314 if(ig>0) then
00315 if(KMTG(ig,jg)>puny) then
00316 nocn(n) = nocn(n) + 1
00317 if(work_gr(ig,jg) > 0.005) nice005(n) = nice005(n) + 1
00318 if(work_gr(ig,jg) > 0.010) nice010(n) = nice010(n) + 1
00319 if(work_gr(ig,jg) > 0.050) nice050(n) = nice050(n) + 1
00320 if(work_gr(ig,jg) > 0.100) nice100(n) = nice100(n) + 1
00321 if(work_gr(ig,jg) > 0.250) nice250(n) = nice250(n) + 1
00322 if(work_gr(ig,jg) > 0.500) nice500(n) = nice500(n) + 1
00323 tmp = tmp + work_gr(ig,jg)
00324 endif
00325 endif
00326 enddo
00327 endif
00328 enddo
00329 if(nocn(n) > 0) then
00330
00331 prob_per_block(n) = real(tmp,kind=dbl_kind)/real(nocn(n),kind=dbl_kind)
00332
00333
00334
00335
00336 work_per_block(n) = &
00337 CEILING(w0 + (real(nocn(n),kind=dbl_kind)/real((ihi-ilo+1)*(jhi-jlo+1),kind=dbl_kind))
00338 *prob_per_block(n)*w1,kind=int_kind)
00339 else
00340 prob_per_block(n) = 0.0d0
00341 endif
00342
00343
00344
00345
00346 if(prob_per_block(n) >0.005 .and. nocn(n) > 0) then
00347 blockType(n) = iceType
00348 elseif (nocn(n) > 0) then
00349 blockType(n) = icefreeType
00350 elseif (nocn(n) == 0) then
00351 blockType(n) = lndType
00352 endif
00353
00354 enddo
00355
00356 numIceFree = COUNT(blockType .eq. iceFreeType)
00357 numLnd = COUNT(blockType .eq. lndType)
00358 numIce = COUNT(blockTYpe .eq. iceType)
00359
00360
00361
00362 write(*,23) nblocks_tot,numIce,numIceFree,numLnd
00363
00364 23 format('CalcWorkPerBlock: Total blocks: ',i5,' Ice blocks: ',i5,' IceFree blocks: ',i5,' Land blocks: ',i5)
00365
00366 endif
00367
00368
00369
00370 call broadcast_array(work_per_block,master_task)
00371 call broadcast_array(prob_per_block,master_task)
00372 call broadcast_array(blockType,master_task)
00373 call broadcast_array(nocn,master_task)
00374 call broadcast_array(nice005,master_task)
00375 call broadcast_array(nice010,master_task)
00376 call broadcast_array(nice050,master_task)
00377 call broadcast_array(nice100,master_task)
00378 call broadcast_array(nice250,master_task)
00379 call broadcast_array(nice500,master_task)
00380
00381 allocate(iLocation(nblocks_tot))
00382 do i=1,nblocks_tot
00383 iLocation(i) = i
00384 enddo
00385
00386 call BuildProbabilityStats2(iLocation,bStats)
00387
00388 deallocate(work_gr)
00389 deallocate(iLocation)
00390
00391
00392 end subroutine CalcWorkPerBlock
00393
00394 function ErfcProbability(lat) result(prob)
00395
00396 real (kind=dbl_kind), intent(in),
00397 dimension (nx_global,ny_global) :: lat
00398
00399 real (kind=real_kind),
00400 dimension (nx_global,ny_global) :: prob
00401
00402
00403 real (kind=dbl_kind) :: ltmp,thetai,sigmai,arg
00404 integer :: i,j
00405
00406
00407 do j=1,ny_global
00408 do i=1,nx_global
00409 ltmp = lat(i,j)
00410 if(ltmp > 0.0d0) then
00411
00412
00413
00414
00415
00416 thetai = (55.0_dbl_kind/180.0_dbl_kind)*pi
00417 sigmai = (5.0_dbl_kind/180.0_dbl_kind)*pi
00418 arg = (thetai-ltmp)/sigmai
00419 else
00420
00421
00422
00423
00424
00425 thetai = (60.0_dbl_kind/180.0_dbl_kind)*pi
00426 sigmai = (5.0_dbl_kind/180.0_dbl_kind)*pi
00427 arg = (thetai+ltmp)/sigmai
00428 endif
00429 prob(i,j) = real(erfc06(arg),kind=real_kind)/2.0_real_kind
00430 enddo
00431 enddo
00432
00433 end function ErfcProbability
00434
00435 function erfc06(x) result(erfc)
00436
00437 implicit none
00438 real(dbl_kind), intent(in) :: x
00439 real(dbl_kind) :: erfc
00440 real(dbl_kind) :: t,u
00441 real(dbl_kind), parameter :: pa = 3.97886080735226000e+00_dbl_kind
00442 real(dbl_kind), parameter :: p00 = 2.75374741597376782e-01_dbl_kind
00443 real(dbl_kind), parameter :: p01 = 4.90165080585318424e-01_dbl_kind
00444 real(dbl_kind), parameter :: p02 = 7.74368199119538609e-01_dbl_kind
00445 real(dbl_kind), parameter :: p03 = 1.07925515155856677e+00_dbl_kind
00446 real(dbl_kind), parameter :: p04 = 1.31314653831023098e+00_dbl_kind
00447 real(dbl_kind), parameter :: p05 = 1.37040217682338167e+00_dbl_kind
00448 real(dbl_kind), parameter :: p06 = 1.18902982909273333e+00_dbl_kind
00449 real(dbl_kind), parameter :: p07 = 8.05276408752910567e-01_dbl_kind
00450 real(dbl_kind), parameter :: p08 = 3.57524274449531043e-01_dbl_kind
00451 real(dbl_kind), parameter :: p09 = 1.66207924969367356e-02_dbl_kind
00452 real(dbl_kind), parameter :: p10 = -1.19463959964325415e-01_dbl_kind
00453 real(dbl_kind), parameter :: p11 = -8.38864557023001992e-02_dbl_kind
00454 real(dbl_kind), parameter :: p12 = 2.49367200053503304e-03_dbl_kind
00455 real(dbl_kind), parameter :: p13 = 3.90976845588484035e-02_dbl_kind
00456 real(dbl_kind), parameter :: p14 = 1.61315329733252248e-02_dbl_kind
00457 real(dbl_kind), parameter :: p15 = -1.33823644533460069e-02_dbl_kind
00458 real(dbl_kind), parameter :: p16 = -1.27223813782122755e-02_dbl_kind
00459 real(dbl_kind), parameter :: p17 = 3.83335126264887303e-03_dbl_kind
00460 real(dbl_kind), parameter :: p18 = 7.73672528313526668e-03_dbl_kind
00461 real(dbl_kind), parameter :: p19 = -8.70779635317295828e-04_dbl_kind
00462 real(dbl_kind), parameter :: p20 = -3.96385097360513500e-03_dbl_kind
00463 real(dbl_kind), parameter :: p21 = 1.19314022838340944e-04_dbl_kind
00464 real(dbl_kind), parameter :: p22 = 1.27109764952614092e-03_dbl_kind
00465
00466
00467
00468
00469 t = pa/(pa + abs(x))
00470 u = t - 0.5_dbl_kind
00471 erfc = ((((((((((((((((((((((p22 * u + p21) * u + p20) &
00472 * u + p19) * u + p18) &
00473 * u + p17) * u + p16) &
00474 * u + p15) * u + p14) &
00475 * u + p13) * u + p12) &
00476 * u + p11) * u + p10) &
00477 * u + p09) * u + p08) &
00478 * u + p07) * u + p06) &
00479 * u + p05) * u + p04) &
00480 * u + p03) * u + p02) &
00481 * u + p01) * u + p00) * t * exp(-x**2)
00482 if (x < 0.0_dbl_kind) erfc = 2.0_dbl_kind - erfc
00483
00484 end function erfc06
00485
00486
00487
00488 end module ice_probability