00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025 module ice_grid
00026
00027
00028
00029 use ice_kinds_mod
00030 use ice_boundary
00031 use ice_communicate, only: my_task, master_task
00032 use ice_constants
00033 use ice_blocks
00034 use ice_domain_size
00035 use ice_domain
00036 use ice_fileunits
00037 use ice_gather_scatter
00038 use ice_read_write
00039 use ice_timers
00040 use ice_probability
00041 use ice_exit
00042
00043
00044
00045 implicit none
00046
00047
00048 character (len=char_len_long), save ::
00049 grid_format ,
00050 grid_file ,
00051 kmt_file ,
00052 grid_type
00053
00054
00055 real (kind=dbl_kind), dimension (nx_block,ny_block,max_blocks), save ::
00056 dxt ,
00057 dyt ,
00058 dxu ,
00059 dyu ,
00060 HTE ,
00061 HTN ,
00062 tarea ,
00063 uarea ,
00064 tarear ,
00065 uarear ,
00066 tinyarea,
00067 tarean ,
00068 tareas ,
00069 ULON ,
00070 ULAT ,
00071 TLON ,
00072 TLAT ,
00073 ANGLE ,
00074 ANGLET
00075
00076 real (kind=dbl_kind), dimension (nx_block,ny_block,max_blocks), save ::
00077 cyp ,
00078 cxp ,
00079 cym ,
00080 cxm ,
00081 dxhy ,
00082 dyhx
00083
00084
00085 real (kind=dbl_kind), dimension (4,nx_block,ny_block,max_blocks), save ::
00086 lont_bounds,
00087 latt_bounds,
00088 lonu_bounds,
00089 latu_bounds
00090
00091
00092 real (kind=dbl_kind), dimension (nx_block,ny_block,max_blocks), save ::
00093 xav ,
00094 yav ,
00095 xxav ,
00096 xyav ,
00097 yyav ,
00098 xxxav,
00099 xxyav,
00100 xyyav,
00101 yyyav
00102
00103 real (kind=dbl_kind),
00104 dimension (2,2,nx_block,ny_block,max_blocks), save ::
00105 mne,
00106 mnw,
00107 mse,
00108 msw
00109
00110
00111 real (kind=dbl_kind), dimension (nx_block,ny_block,max_blocks), save ::
00112 hm ,
00113 uvm
00114
00115 real (kind=dbl_kind), dimension (nx_block,ny_block,max_blocks), save ::
00116 ocn_gridcell_frac
00117
00118
00119 logical (kind=log_kind),
00120 dimension (nx_block,ny_block,max_blocks), save ::
00121 tmask ,
00122 umask ,
00123 lmask_n,
00124 lmask_s
00125
00126
00127 real (kind=dbl_kind), parameter ::
00128 dxrect = 30.e5_dbl_kind ,
00129 dyrect = 30.e5_dbl_kind
00130
00131 real (kind=dbl_kind), dimension (nx_block,ny_block,max_blocks), save ::
00132 rndex_global
00133
00134 real (kind=dbl_kind), private,
00135 dimension(nx_block,ny_block,max_blocks) ::
00136 work1
00137
00138
00139
00140 contains
00141
00142
00143
00144
00145
00146
00147
00148
00149 subroutine init_grid1
00150
00151
00152
00153
00154
00155
00156
00157
00158
00159
00160
00161
00162 use ice_broadcast
00163 use ice_work, only: work_g1, work_g2
00164
00165
00166
00167
00168
00169
00170 integer (kind=int_kind) ::
00171 i, j, iblk,
00172 fid_grid,
00173 fid_kmt
00174
00175 real (dbl_kind), allocatable, dimension(:,:) :: bStats
00176
00177 integer(kind=int_kind), allocatable :: WorkPerBlock(:),blockType(:)
00178 real(kind=dbl_kind), allocatable :: ProbPerBlock(:)
00179
00180 character (char_len) ::
00181 fieldname
00182
00183
00184
00185
00186
00187 allocate(work_g1(nx_global,ny_global))
00188 allocate(work_g2(nx_global,ny_global))
00189
00190 if (trim(grid_type) == 'displaced_pole' .or. &
00191 trim(grid_type) == 'tripole' ) then
00192
00193 if (trim(grid_format) == 'nc') then
00194
00195 call ice_open_nc(grid_file,fid_grid)
00196 call ice_open_nc(kmt_file,fid_kmt)
00197
00198 fieldname='ulat'
00199 call ice_read_global_nc(fid_grid,1,fieldname,work_g1,.true.)
00200 fieldname='kmt'
00201 call ice_read_global_nc(fid_kmt,1,fieldname,work_g2,.true.)
00202
00203 if (my_task == master_task) then
00204 call ice_close_nc(fid_grid)
00205 call ice_close_nc(fid_kmt)
00206 endif
00207
00208 else
00209
00210 call ice_open(nu_grid,grid_file,64)
00211 call ice_open(nu_kmt, kmt_file, 32)
00212
00213 call ice_read_global(nu_grid,1,work_g1,'rda8',.true.)
00214 call ice_read_global(nu_kmt, 1,work_g2,'ida4',.true.)
00215
00216 if (my_task == master_task) then
00217 close (nu_grid)
00218 close (nu_kmt)
00219 endif
00220
00221 endif
00222
00223 elseif (trim(grid_type) == 'panarctic') then
00224
00225 call ice_open(nu_grid,grid_file,64)
00226
00227 call ice_read_global(nu_grid,1,work_g2,'ida8',.true.)
00228 call ice_read_global(nu_grid,2,work_g1,'rda8',.true.)
00229
00230
00231
00232 if (my_task == master_task) close (nu_grid)
00233
00234 else
00235
00236 work_g1(:,:) = 75._dbl_kind/rad_to_deg
00237 work_g2(:,:) = c1
00238
00239 endif
00240
00241 call broadcast_array(work_g1, master_task)
00242 call broadcast_array(work_g2, master_task)
00243
00244 allocate(WorkPerBlock(nblocks_tot),ProbPerBlock(nblocks_tot),blockType(nblocks_tot))
00245 allocate(bStats(numCoeff,nblocks_tot))
00246 call CalcWorkPerBlock(distribution_wght, work_g2,work_g1, &
00247 WorkPerBlock,ProbPerBlock,blockType,bStats)
00248
00249
00250
00251
00252
00253
00254
00255
00256
00257
00258
00259 call init_domain_distribution(work_g2, work_g1, &
00260 WorkPerBlock,ProbPerBlock,blockType,bStats)
00261
00262
00263 deallocate(bStats)
00264 deallocate(WorkPerBlock,ProbPerBlock,blockType)
00265
00266 deallocate(work_g1)
00267 deallocate(work_g2)
00268
00269
00270
00271
00272
00273 if (my_task == master_task) then
00274 write(nu_diag,'(a26,i6)') ' Block size: nx_block = ',nx_block
00275 write(nu_diag,'(a26,i6)') ' ny_block = ',ny_block
00276 endif
00277
00278 end subroutine init_grid1
00279
00280
00281
00282
00283
00284
00285
00286
00287 subroutine init_grid2
00288
00289
00290
00291
00292
00293
00294
00295
00296
00297
00298
00299
00300
00301
00302
00303
00304
00305
00306
00307 use ice_work, only: work_g1
00308 use ice_exit
00309 use ice_blocks, only: get_block, block
00310
00311
00312
00313
00314
00315 integer (kind=int_kind) ::
00316 i, j, iblk,
00317 ilo,ihi,jlo,jhi
00318
00319 real (kind=dbl_kind) ::
00320 angle_0, angle_w, angle_s, angle_sw
00321
00322 logical (kind=log_kind), dimension(nx_block,ny_block,max_blocks)::
00323 out_of_range
00324
00325 type (block) ::
00326 this_block
00327
00328
00329
00330
00331
00332 if (trim(grid_type) == 'displaced_pole' .or. &
00333 trim(grid_type) == 'tripole' ) then
00334 if (trim(grid_format) == 'nc') then
00335 call popgrid_nc
00336 else
00337 call popgrid
00338 endif
00339 elseif (trim(grid_type) == 'panarctic') then
00340 call panarctic_grid
00341 #ifdef CCSMCOUPLED
00342 elseif (trim(grid_type) == 'latlon') then
00343 call latlongrid
00344 return
00345 #endif
00346 else
00347 call rectgrid
00348 endif
00349
00350
00351
00352
00353
00354
00355 do iblk = 1, nblocks
00356 this_block = get_block(blocks_ice(iblk),iblk)
00357 ilo = this_block%ilo
00358 ihi = this_block%ihi
00359 jlo = this_block%jlo
00360 jhi = this_block%jhi
00361
00362 do j = jlo, jhi
00363 do i = ilo, ihi
00364 tarea(i,j,iblk) = dxt(i,j,iblk)*dyt(i,j,iblk)
00365 uarea(i,j,iblk) = dxu(i,j,iblk)*dyu(i,j,iblk)
00366 if (tarea(i,j,iblk) > c0) then
00367 tarear(i,j,iblk) = c1/tarea(i,j,iblk)
00368 else
00369 tarear(i,j,iblk) = c0
00370 endif
00371 if (uarea(i,j,iblk) > c0) then
00372 uarear(i,j,iblk) = c1/uarea(i,j,iblk)
00373 else
00374 uarear(i,j,iblk) = c0
00375 endif
00376 tinyarea(i,j,iblk) = puny*tarea(i,j,iblk)
00377
00378 dxhy(i,j,iblk) = p5*(HTE(i,j,iblk) - HTE(i-1,j,iblk))
00379 dyhx(i,j,iblk) = p5*(HTN(i,j,iblk) - HTN(i,j-1,iblk))
00380 enddo
00381 enddo
00382
00383 do j = jlo, jhi+1
00384 do i = ilo, ihi+1
00385 cyp(i,j,iblk) = (c1p5*HTE(i,j,iblk) - p5*HTE(i-1,j,iblk))
00386 cxp(i,j,iblk) = (c1p5*HTN(i,j,iblk) - p5*HTN(i,j-1,iblk))
00387
00388 cym(i,j,iblk) = -(c1p5*HTE(i-1,j,iblk) - p5*HTE(i,j,iblk))
00389 cxm(i,j,iblk) = -(c1p5*HTN(i,j-1,iblk) - p5*HTN(i,j,iblk))
00390 enddo
00391 enddo
00392
00393 enddo
00394
00395
00396
00397
00398
00399
00400
00401
00402
00403
00404
00405
00406 call ice_timer_start(timer_bound)
00407 call ice_HaloUpdate (tarea, halo_info, &
00408 field_loc_center, field_type_scalar, &
00409 fillValue=c1)
00410 call ice_HaloUpdate (uarea, halo_info, &
00411 field_loc_NEcorner, field_type_scalar, &
00412 fillValue=c1)
00413 call ice_HaloUpdate (tarear, halo_info, &
00414 field_loc_center, field_type_scalar, &
00415 fillValue=c1)
00416 call ice_HaloUpdate (uarear, halo_info, &
00417 field_loc_NEcorner, field_type_scalar, &
00418 fillValue=c1)
00419 call ice_HaloUpdate (tinyarea, halo_info, &
00420 field_loc_center, field_type_scalar, &
00421 fillValue=c1)
00422 call ice_HaloUpdate (dxhy, halo_info, &
00423 field_loc_center, field_type_vector, &
00424 fillValue=c1)
00425 call ice_HaloUpdate (dyhx, halo_info, &
00426 field_loc_center, field_type_vector, &
00427 fillValue=c1)
00428 call ice_timer_stop(timer_bound)
00429
00430
00431
00432
00433
00434
00435 out_of_range = .false.
00436 where (ANGLE < -pi .or. ANGLE > pi) out_of_range = .true.
00437 if (count(out_of_range) > 0) then
00438 call abort_ice ('ice: init_grid: ANGLE out of expected range')
00439 endif
00440
00441
00442
00443
00444 ANGLET = c0
00445
00446
00447
00448 do iblk = 1, nblocks
00449 this_block = get_block(blocks_ice(iblk),iblk)
00450 ilo = this_block%ilo
00451 ihi = this_block%ihi
00452 jlo = this_block%jlo
00453 jhi = this_block%jhi
00454
00455 do j = jlo, jhi
00456 do i = ilo, ihi
00457 angle_0 = ANGLE(i ,j ,iblk)
00458 angle_w = ANGLE(i-1,j ,iblk)
00459 angle_s = ANGLE(i, j-1,iblk)
00460 angle_sw = ANGLE(i-1,j-1,iblk)
00461
00462 if ( angle_0 < c0 ) then
00463 if ( abs(angle_w - angle_0) > pi) &
00464 angle_w = angle_w - pi2
00465 if ( abs(angle_s - angle_0) > pi) &
00466 angle_s = angle_s - pi2
00467 if ( abs(angle_sw - angle_0) > pi) &
00468 angle_sw = angle_sw - pi2
00469 endif
00470
00471 ANGLET(i,j,iblk) = angle_0 * p25 + angle_w * p25 &
00472 + angle_s * p25 + angle_sw* p25
00473 enddo
00474 enddo
00475 enddo
00476
00477
00478 call ice_timer_start(timer_bound)
00479 call ice_HaloUpdate (ANGLET, halo_info, &
00480 field_loc_center, field_type_angle, &
00481 fillValue=c1)
00482 call ice_timer_stop(timer_bound)
00483
00484 call makemask
00485
00486 call Tlatlon
00487
00488
00489
00490
00491
00492 call gridbox_corners
00493
00494
00495
00496
00497
00498 if (my_task==master_task) then
00499 allocate(work_g1(nx_global,ny_global))
00500 do j=1,ny_global
00501 do i=1,nx_global
00502 work_g1(i,j) = real((j-1)*nx_global + i,kind=dbl_kind)
00503 enddo
00504 enddo
00505 else
00506 allocate(work_g1(1,1))
00507 endif
00508
00509 call scatter_global(rndex_global, work_g1, &
00510 master_task, distrb_info, &
00511 field_loc_center, field_type_scalar)
00512
00513 deallocate(work_g1)
00514
00515 end subroutine init_grid2
00516
00517
00518
00519
00520
00521
00522
00523
00524
00525 subroutine popgrid
00526
00527
00528
00529
00530
00531
00532
00533
00534
00535
00536
00537
00538
00539
00540
00541
00542
00543
00544
00545
00546
00547 use ice_work, only: work_g1
00548
00549
00550
00551
00552
00553
00554 integer (kind=int_kind) ::
00555 i, j, iblk,
00556 ilo,ihi,jlo,jhi
00557
00558 logical (kind=log_kind) :: diag
00559
00560 type (block) ::
00561 this_block
00562
00563 call ice_open(nu_grid,grid_file,64)
00564 call ice_open(nu_kmt,kmt_file,32)
00565
00566 diag = .true.
00567
00568
00569
00570
00571
00572 call ice_read(nu_kmt,1,work1,'ida4',diag, &
00573 field_loc=field_loc_center, &
00574 field_type=field_type_scalar)
00575
00576 hm(:,:,:) = c0
00577
00578 do iblk = 1, nblocks
00579 this_block = get_block(blocks_ice(iblk),iblk)
00580 ilo = this_block%ilo
00581 ihi = this_block%ihi
00582 jlo = this_block%jlo
00583 jhi = this_block%jhi
00584
00585 do j = jlo, jhi
00586 do i = ilo, ihi
00587 hm(i,j,iblk) = work1(i,j,iblk)
00588 if (hm(i,j,iblk) >= c1) hm(i,j,iblk) = c1
00589 enddo
00590 enddo
00591 enddo
00592
00593
00594
00595
00596
00597
00598 if (my_task == master_task) then
00599 allocate(work_g1(nx_global,ny_global))
00600 else
00601 allocate(work_g1(1,1))
00602 endif
00603
00604 call ice_read_global(nu_grid,1,work_g1,'rda8',.true.)
00605 call gridbox_verts(work_g1,latt_bounds)
00606 call scatter_global(ULAT, work_g1, master_task, distrb_info, &
00607 field_loc_NEcorner, field_type_scalar)
00608 call ice_HaloExtrapolate(ULAT, distrb_info, &
00609 ew_boundary_type, ns_boundary_type)
00610
00611 call ice_read_global(nu_grid,2,work_g1,'rda8',.true.)
00612 call gridbox_verts(work_g1,lont_bounds)
00613 call scatter_global(ULON, work_g1, master_task, distrb_info, &
00614 field_loc_NEcorner, field_type_scalar)
00615 call ice_HaloExtrapolate(ULON, distrb_info, &
00616 ew_boundary_type, ns_boundary_type)
00617
00618 call ice_read_global(nu_grid,7,work_g1,'rda8',.true.)
00619 call scatter_global(ANGLE, work_g1, master_task, distrb_info, &
00620 field_loc_NEcorner, field_type_angle)
00621
00622
00623
00624
00625
00626
00627
00628 call ice_read_global(nu_grid,3,work_g1,'rda8',.true.)
00629 call primary_grid_lengths_HTN(work_g1)
00630
00631 call ice_read_global(nu_grid,4,work_g1,'rda8',.true.)
00632 call primary_grid_lengths_HTE(work_g1)
00633
00634 deallocate(work_g1)
00635
00636 if (my_task == master_task) then
00637 close (nu_grid)
00638 close (nu_kmt)
00639 endif
00640
00641 end subroutine popgrid
00642
00643
00644
00645
00646
00647
00648
00649
00650
00651 subroutine popgrid_nc
00652
00653 #ifdef ncdf
00654
00655
00656
00657
00658
00659
00660
00661
00662
00663
00664
00665
00666
00667
00668
00669
00670
00671
00672
00673
00674
00675
00676 use ice_work, only: work_g1
00677
00678
00679
00680
00681
00682
00683 integer (kind=int_kind) ::
00684 i, j, iblk,
00685 ilo,ihi,jlo,jhi,
00686 fid_grid,
00687 fid_kmt
00688
00689 logical (kind=log_kind) :: diag
00690
00691 character (char_len) ::
00692 fieldname
00693
00694 type (block) ::
00695 this_block
00696
00697 call ice_open_nc(grid_file,fid_grid)
00698 call ice_open_nc(kmt_file,fid_kmt)
00699
00700 diag = .true.
00701
00702
00703
00704
00705
00706 fieldname='kmt'
00707 call ice_read_nc(fid_kmt,1,fieldname,work1,diag, &
00708 field_loc=field_loc_center, &
00709 field_type=field_type_scalar)
00710
00711 hm(:,:,:) = c0
00712
00713 do iblk = 1, nblocks
00714 this_block = get_block(blocks_ice(iblk),iblk)
00715 ilo = this_block%ilo
00716 ihi = this_block%ihi
00717 jlo = this_block%jlo
00718 jhi = this_block%jhi
00719
00720 do j = jlo, jhi
00721 do i = ilo, ihi
00722 hm(i,j,iblk) = work1(i,j,iblk)
00723 if (hm(i,j,iblk) >= c1) hm(i,j,iblk) = c1
00724 enddo
00725 enddo
00726 enddo
00727
00728
00729
00730
00731
00732
00733 if (my_task == master_task) then
00734 allocate(work_g1(nx_global,ny_global))
00735 else
00736 allocate(work_g1(1,1))
00737 endif
00738
00739 fieldname='ulat'
00740 call ice_read_global_nc(fid_grid,1,fieldname,work_g1,diag)
00741 call gridbox_verts(work_g1,latt_bounds)
00742 call scatter_global(ULAT, work_g1, master_task, distrb_info, &
00743 field_loc_NEcorner, field_type_scalar)
00744 call ice_HaloExtrapolate(ULAT, distrb_info, &
00745 ew_boundary_type, ns_boundary_type)
00746
00747 fieldname='ulon'
00748 call ice_read_global_nc(fid_grid,2,fieldname,work_g1,diag)
00749 call gridbox_verts(work_g1,lont_bounds)
00750 call scatter_global(ULON, work_g1, master_task, distrb_info, &
00751 field_loc_NEcorner, field_type_scalar)
00752 call ice_HaloExtrapolate(ULON, distrb_info, &
00753 ew_boundary_type, ns_boundary_type)
00754
00755 fieldname='angle'
00756 call ice_read_global_nc(fid_grid,7,fieldname,work_g1,diag)
00757 call scatter_global(ANGLE, work_g1, master_task, distrb_info, &
00758 field_loc_NEcorner, field_type_angle)
00759
00760
00761 where (ANGLE > pi) ANGLE = pi
00762 where (ANGLE < -pi) ANGLE = -pi
00763
00764
00765
00766
00767
00768
00769
00770 fieldname='htn'
00771 call ice_read_global_nc(fid_grid,3,fieldname,work_g1,diag)
00772 call primary_grid_lengths_HTN(work_g1)
00773
00774 fieldname='hte'
00775 call ice_read_global_nc(fid_grid,4,fieldname,work_g1,diag)
00776 call primary_grid_lengths_HTE(work_g1)
00777
00778 deallocate(work_g1)
00779
00780 if (my_task == master_task) then
00781 call ice_close_nc(fid_grid)
00782 call ice_close_nc(fid_kmt)
00783 endif
00784
00785 #endif
00786 end subroutine popgrid_nc
00787
00788
00789
00790
00791
00792
00793
00794
00795 subroutine panarctic_grid
00796
00797
00798
00799
00800
00801
00802
00803
00804
00805
00806
00807
00808 use ice_domain_size
00809 use ice_work, only: work_g1
00810
00811
00812
00813
00814
00815
00816
00817
00818
00819
00820
00821
00822
00823
00824
00825
00826
00827
00828
00829
00830
00831
00832
00833
00834 integer (kind=int_kind) ::
00835 i, j, iblk,
00836 ilo,ihi,jlo,jhi
00837
00838 logical (kind=log_kind) :: diag
00839
00840 type (block) ::
00841 this_block
00842
00843 call ice_open(nu_grid,grid_file,64)
00844
00845 diag = .true.
00846
00847 if (my_task == master_task) &
00848 write (nu_diag,*) '** Reading pan-Arctic grid **'
00849
00850
00851
00852
00853
00854 call ice_read(nu_grid,1,work1,'ida8',diag, &
00855 field_loc=field_loc_center, &
00856 field_type=field_type_scalar)
00857
00858 hm(:,:,:) = c0
00859
00860 do iblk = 1, nblocks
00861 this_block = get_block(blocks_ice(iblk),iblk)
00862 ilo = this_block%ilo
00863 ihi = this_block%ihi
00864 jlo = this_block%jlo
00865 jhi = this_block%jhi
00866
00867 do j = jlo, jhi
00868 do i = ilo, ihi
00869 hm(i,j,iblk) = work1(i,j,iblk)
00870 if (hm(i,j,iblk) >= c1) hm(i,j,iblk) = c1
00871 enddo
00872 enddo
00873 enddo
00874
00875
00876
00877
00878
00879
00880 if (my_task == master_task) then
00881 allocate(work_g1(nx_global,ny_global))
00882 else
00883 allocate(work_g1(1,1))
00884 endif
00885
00886 call ice_read_global(nu_grid,2,work_g1,'rda8',.true.)
00887 call gridbox_verts(work_g1,latt_bounds)
00888 call scatter_global(ULAT, work_g1, master_task, distrb_info, &
00889 field_loc_NEcorner, field_type_scalar)
00890 call ice_HaloExtrapolate(ULAT, distrb_info, &
00891 ew_boundary_type, ns_boundary_type)
00892
00893 call ice_read_global(nu_grid,3,work_g1,'rda8',.true.)
00894 call gridbox_verts(work_g1,lont_bounds)
00895 call scatter_global(ULON, work_g1, master_task, distrb_info, &
00896 field_loc_NEcorner, field_type_scalar)
00897 call ice_HaloExtrapolate(ULON, distrb_info, &
00898 ew_boundary_type, ns_boundary_type)
00899
00900 call ice_read_global(nu_grid,8,work_g1,'rda8',.true.)
00901 call scatter_global(ANGLE, work_g1, master_task, distrb_info, &
00902 field_loc_NEcorner, field_type_angle)
00903
00904
00905
00906
00907
00908
00909
00910 call ice_read_global(nu_grid,4,work_g1,'rda8',.true.)
00911 call primary_grid_lengths_HTN(work_g1)
00912
00913 call ice_read_global(nu_grid,5,work_g1,'rda8',.true.)
00914 call primary_grid_lengths_HTE(work_g1)
00915
00916 deallocate(work_g1)
00917
00918 if (my_task == master_task) close (nu_grid)
00919
00920 end subroutine panarctic_grid
00921 #ifdef CCSMCOUPLED
00922
00923
00924
00925
00926
00927
00928
00929 subroutine latlongrid
00930
00931
00932
00933
00934
00935
00936
00937
00938
00939
00940
00941
00942
00943
00944 use ice_domain_size
00945 use ice_scam, only : scmlat, scmlon, single_column
00946 use netcdf
00947
00948
00949 include "netcdf.inc"
00950
00951
00952
00953 integer (kind=int_kind) ::
00954 i, j, iblk
00955
00956 integer (kind=int_kind) ::
00957 ni, nj, ncid, dimid, varid, ier
00958
00959 character (len=char_len) ::
00960 subname='latlongrid'
00961
00962 type (block) ::
00963 this_block
00964
00965 integer (kind=int_kind) ::
00966 ilo,ihi,jlo,jhi
00967
00968 real (kind=dbl_kind) ::
00969 closelat,
00970 closelon,
00971 closelatidx,
00972 closelonidx
00973
00974 integer (kind=int_kind) ::
00975 start(2),
00976 count(2)
00977
00978 integer (kind=int_kind) ::
00979 start3(3),
00980 count3(3)
00981
00982 integer (kind=int_kind) ::
00983 status
00984
00985 real (kind=dbl_kind), allocatable ::
00986 lats(:),lons(:),pos_lons(:), glob_grid(:,:)
00987
00988 real (kind=dbl_kind) ::
00989 pos_scmlon,
00990 scamdata
00991
00992
00993
00994
00995
00996
00997
00998
00999
01000
01001 if (my_task == master_task) then
01002 call ice_open_nc(kmt_file, ncid)
01003
01004 status = nf90_inq_dimid (ncid, 'ni', dimid)
01005 status = nf90_inquire_dimension(ncid, dimid, len=ni)
01006 status = nf90_inq_dimid (ncid, 'nj', dimid)
01007 status = nf90_inquire_dimension(ncid, dimid, len=nj)
01008 end if
01009
01010
01011
01012
01013 if (single_column) then
01014
01015 if (my_task == master_task) then
01016 if ((nx_global /= 1).or. (ny_global /= 1)) then
01017 write(nu_diag,*) 'Because you have selected the column model flag'
01018 write(nu_diag,*) 'Please set nx_global=ny_global=1 in file'
01019 write(nu_diag,*) 'ice_domain_size.F and recompile'
01020 call abort_ice ('latlongrid: check nx_global, ny_global')
01021 endif
01022 end if
01023
01024
01025 allocate(lats(nj))
01026 allocate(lons(ni))
01027 allocate(pos_lons(ni))
01028 allocate(glob_grid(ni,nj))
01029
01030 start3=(/1,1,1/)
01031 count3=(/ni,nj,1/)
01032 call check_ret(nf_inq_varid(ncid, 'xc' , varid), subname)
01033 call check_ret(nf_get_vara_double(ncid, varid,start3, count3, glob_grid), subname)
01034 do i = 1,ni
01035 lons(i) = glob_grid(i,1)
01036 end do
01037
01038 call check_ret(nf_inq_varid(ncid, 'yc' , varid), subname)
01039 call check_ret(nf_get_vara_double(ncid, varid,start3, count3, glob_grid), subname)
01040 do j = 1,nj
01041 lats(j) = glob_grid(1,j)
01042 end do
01043
01044
01045
01046
01047 pos_lons(:)= mod(lons(:) + 360._dbl_kind,360._dbl_kind)
01048 pos_scmlon = mod(scmlon + 360._dbl_kind,360._dbl_kind)
01049 start(1) = (MINLOC(abs(pos_lons-pos_scmlon),dim=1))
01050 start(2) = (MINLOC(abs(lats -scmlat ),dim=1))
01051 count(1) = 1
01052 count(2) = 1
01053
01054 deallocate(lats)
01055 deallocate(lons)
01056 deallocate(pos_lons)
01057 deallocate(glob_grid)
01058
01059 call check_ret(nf_inq_varid(ncid, 'xc' , varid), subname)
01060 call check_ret(nf_get_vara_double(ncid, varid, start, count, scamdata), subname)
01061 TLON = scamdata
01062 call check_ret(nf_inq_varid(ncid, 'yc' , varid), subname)
01063 call check_ret(nf_get_vara_double(ncid, varid, start, count, scamdata), subname)
01064 TLAT = scamdata
01065 call check_ret(nf_inq_varid(ncid, 'area' , varid), subname)
01066 call check_ret(nf_get_vara_double(ncid, varid, start, count, scamdata), subname)
01067 tarea = scamdata
01068 call check_ret(nf_inq_varid(ncid, 'mask' , varid), subname)
01069 call check_ret(nf_get_vara_double(ncid, varid, start, count, scamdata), subname)
01070 hm = scamdata
01071 call check_ret(nf_inq_varid(ncid, 'frac' , varid), subname)
01072 call check_ret(nf_get_vara_double(ncid, varid, start, count, scamdata), subname)
01073 ocn_gridcell_frac = scamdata
01074 else
01075
01076 if (my_task == master_task) then
01077 if (nx_global /= ni .and. ny_global /= nj) then
01078 call abort_ice ('latlongrid: ni,ny not equal to nx_global,ny_global')
01079 end if
01080 end if
01081
01082
01083 call ice_read_nc(ncid, 1, 'xc' , TLON , diag=.true.)
01084 call ice_read_nc(ncid, 1, 'yc' , TLAT , diag=.true.)
01085 call ice_read_nc(ncid, 1, 'area', tarea , diag=.true., &
01086 field_loc=field_loc_center,field_type=field_type_scalar)
01087 call ice_read_nc(ncid, 1, 'mask', hm , diag=.true.)
01088 call ice_read_nc(ncid, 1, 'frac', ocn_gridcell_frac, diag=.true.)
01089 end if
01090
01091 if (my_task == master_task) then
01092 call ice_close_nc(ncid)
01093 end if
01094
01095
01096 do iblk = 1,nblocks
01097 this_block = get_block(blocks_ice(iblk),iblk)
01098 ilo = this_block%ilo
01099 ihi = this_block%ihi
01100 jlo = this_block%jlo
01101 jhi = this_block%jhi
01102
01103 do j = jlo, jhi
01104 do i = ilo, ihi
01105
01106 TLON(i,j,iblk) = pi*TLON(i,j,iblk)/180._dbl_kind
01107
01108
01109 TLAT(i,j,iblk) = pi*TLAT(i,j,iblk)/180._dbl_kind
01110
01111
01112
01113 tarea(i,j,iblk) = tarea(i,j,iblk) * (radius*radius)
01114 end do
01115 end do
01116 end do
01117
01118
01119
01120
01121
01122
01123
01124
01125
01126
01127
01128
01129
01130
01131
01132
01133 do iblk = 1, nblocks
01134 this_block = get_block(blocks_ice(iblk),iblk)
01135 ilo = this_block%ilo
01136 ihi = this_block%ihi
01137 jlo = this_block%jlo
01138 jhi = this_block%jhi
01139
01140 do j = jlo, jhi
01141 do i = ilo, ihi
01142
01143 uarea(i,j,iblk) = p25* &
01144 (tarea(i,j, iblk) + tarea(i+1,j, iblk) &
01145 + tarea(i,j+1,iblk) + tarea(i+1,j+1,iblk))
01146 tarear(i,j,iblk) = c1/tarea(i,j,iblk)
01147 uarear(i,j,iblk) = c1/uarea(i,j,iblk)
01148 tinyarea(i,j,iblk) = puny*tarea(i,j,iblk)
01149
01150 if (single_column) then
01151 ULAT (i,j,iblk) = TLAT(i,j,iblk)+(pi/nj)
01152 else
01153 ULAT (i,j,iblk) = TLAT(i,j,iblk)+(pi/ny_global)
01154 endif
01155 ULON (i,j,iblk) = c0
01156 ANGLE (i,j,iblk) = c0
01157
01158 ANGLET(i,j,iblk) = c0
01159 HTN (i,j,iblk) = 1.e36_dbl_kind
01160 HTE (i,j,iblk) = 1.e36_dbl_kind
01161 dxt (i,j,iblk) = 1.e36_dbl_kind
01162 dyt (i,j,iblk) = 1.e36_dbl_kind
01163 dxu (i,j,iblk) = 1.e36_dbl_kind
01164 dyu (i,j,iblk) = 1.e36_dbl_kind
01165 dxhy (i,j,iblk) = 1.e36_dbl_kind
01166 dyhx (i,j,iblk) = 1.e36_dbl_kind
01167 cyp (i,j,iblk) = 1.e36_dbl_kind
01168 cxp (i,j,iblk) = 1.e36_dbl_kind
01169 cym (i,j,iblk) = 1.e36_dbl_kind
01170 cxm (i,j,iblk) = 1.e36_dbl_kind
01171 enddo
01172 enddo
01173 enddo
01174
01175
01176 call makemask
01177
01178 end subroutine latlongrid
01179
01180
01181
01182
01183
01184
01185
01186 subroutine check_ret(ret, calling)
01187
01188
01189
01190
01191
01192 use netcdf
01193
01194 implicit none
01195 integer, intent(in) :: ret
01196 character(len=*) :: calling
01197
01198
01199
01200
01201
01202
01203
01204 if (ret /= NF90_NOERR) then
01205 write(nu_diag,*)'netcdf error from ',trim(calling)
01206 write(nu_diag,*)'netcdf strerror = ',trim(NF90_STRERROR(ret))
01207 call abort_ice('ice ice_grid: netcdf check_ret error')
01208 end if
01209
01210 end subroutine check_ret
01211
01212 #endif
01213
01214
01215
01216
01217
01218
01219
01220 subroutine rectgrid
01221
01222
01223
01224
01225
01226
01227
01228
01229
01230
01231
01232 use ice_domain_size
01233 use ice_work, only: work_g1
01234
01235
01236
01237
01238
01239 integer (kind=int_kind) ::
01240 i, j, iblk,
01241 imid, jmid
01242
01243 real (kind=dbl_kind) :: length
01244
01245
01246
01247
01248
01249
01250 do iblk = 1, nblocks
01251 do j = 1, ny_block
01252 do i = 1, nx_block
01253 ANGLE(i,j,iblk) = c0
01254 enddo
01255 enddo
01256 enddo
01257
01258
01259 if (my_task == master_task) then
01260 allocate(work_g1(nx_global,ny_global))
01261 else
01262 allocate(work_g1(1,1))
01263 endif
01264
01265
01266
01267
01268 if (my_task == master_task) then
01269 work_g1 = c0
01270 length = dxrect*cm_to_m/radius*rad_to_deg
01271 work_g1(1,:) = -55._dbl_kind
01272 do j = 1, ny_global
01273 do i = 2, nx_global
01274 work_g1(i,j) = work_g1(i-1,j) + length
01275 enddo
01276 enddo
01277 work_g1(:,:) = work_g1(:,:) / rad_to_deg
01278 endif
01279 call scatter_global(ULON, work_g1, master_task, distrb_info, &
01280 field_loc_NEcorner, field_type_scalar)
01281
01282 if (my_task == master_task) then
01283 work_g1 = c0
01284 length = dyrect*cm_to_m/radius*rad_to_deg
01285 work_g1(:,1) = -75._dbl_kind
01286 do i = 1, nx_global
01287 do j = 2, ny_global
01288 work_g1(i,j) = work_g1(i,j-1) + length
01289 enddo
01290 enddo
01291 work_g1(:,:) = work_g1(:,:) / rad_to_deg
01292 endif
01293 call scatter_global(ULAT, work_g1, master_task, distrb_info, &
01294 field_loc_NEcorner, field_type_scalar)
01295
01296 if (my_task == master_task) then
01297 do j = 1, ny_global
01298 do i = 1, nx_global
01299 work_g1(i,j) = dxrect
01300 enddo
01301 enddo
01302 endif
01303 call primary_grid_lengths_HTN(work_g1)
01304
01305 if (my_task == master_task) then
01306 do j = 1, ny_global
01307 do i = 1, nx_global
01308 work_g1(i,j) = dyrect
01309 enddo
01310 enddo
01311 endif
01312 call primary_grid_lengths_HTE(work_g1)
01313
01314
01315
01316
01317
01318
01319 if (my_task == master_task) then
01320 work_g1(:,:) = c0
01321
01322 if (trim(ew_boundary_type) == 'cyclic') then
01323
01324 do j = 3,ny_global-2
01325 do i = 1,nx_global
01326 work_g1(i,j) = c1
01327 enddo
01328 enddo
01329
01330 elseif (trim(ew_boundary_type) == 'closed') then
01331
01332 do j = 3,ny_global-2
01333 do i = 3,nx_global-2
01334 work_g1(i,j) = c1
01335 enddo
01336 enddo
01337
01338 elseif (trim(ew_boundary_type) == 'open') then
01339
01340
01341
01342 imid = aint(real(nx_global)/c2,kind=int_kind)
01343 jmid = aint(real(ny_global)/c2,kind=int_kind)
01344
01345 do j = 3,ny_global-2
01346 do i = 3,nx_global-2
01347 work_g1(i,j) = c1
01348 enddo
01349 enddo
01350
01351 do j = 1, jmid+2
01352 do i = 1, imid+2
01353 work_g1(i,j) = c1
01354 enddo
01355 enddo
01356
01357 do j = jmid-2, ny_global
01358 do i = imid-2, nx_global
01359 work_g1(i,j) = c1
01360 enddo
01361 enddo
01362
01363 endif
01364 endif
01365
01366 call scatter_global(hm, work_g1, master_task, distrb_info, &
01367 field_loc_center, field_type_scalar)
01368
01369 deallocate(work_g1)
01370
01371 end subroutine rectgrid
01372
01373
01374
01375
01376
01377
01378
01379
01380 subroutine primary_grid_lengths_HTN(work_g)
01381
01382
01383
01384
01385
01386
01387
01388
01389
01390
01391
01392
01393
01394 use ice_work, only: work_g2
01395
01396
01397
01398
01399
01400
01401
01402 real (kind=dbl_kind), dimension(:,:) :: work_g
01403
01404 integer (kind=int_kind) ::
01405 i, j,
01406 ip1
01407
01408 if (my_task == master_task) then
01409 allocate(work_g2(nx_global,ny_global))
01410 else
01411 allocate(work_g2(1,1))
01412 endif
01413
01414 if (my_task == master_task) then
01415 do j = 1, ny_global
01416 do i = 1, nx_global
01417 work_g(i,j) = work_g(i,j) * cm_to_m
01418 enddo
01419 enddo
01420 do j = 1, ny_global
01421 do i = 1, nx_global
01422
01423 ip1 = i+1
01424 if (i == nx_global) ip1 = 1
01425 work_g2(i,j) = p5*(work_g(i,j) + work_g(ip1,j))
01426 enddo
01427 enddo
01428 endif
01429 call scatter_global(HTN, work_g, master_task, distrb_info, &
01430 field_loc_Nface, field_type_scalar)
01431 call scatter_global(dxu, work_g2, master_task, distrb_info, &
01432 field_loc_NEcorner, field_type_scalar)
01433
01434 if (my_task == master_task) then
01435 do j = 2, ny_global
01436 do i = 1, nx_global
01437 work_g2(i,j) = p5*(work_g(i,j) + work_g(i,j-1))
01438 enddo
01439 enddo
01440
01441 do i = 1, nx_global
01442 work_g2(i,1) = c2*work_g(i,2) - work_g(i,3)
01443 enddo
01444 endif
01445 call scatter_global(dxt, work_g2, master_task, distrb_info, &
01446 field_loc_center, field_type_scalar)
01447
01448 deallocate(work_g2)
01449
01450 end subroutine primary_grid_lengths_HTN
01451
01452
01453
01454
01455
01456
01457
01458
01459 subroutine primary_grid_lengths_HTE(work_g)
01460
01461
01462
01463
01464
01465
01466
01467
01468
01469
01470
01471
01472
01473 use ice_work, only: work_g2
01474
01475
01476
01477
01478
01479
01480
01481 real (kind=dbl_kind), dimension(:,:) :: work_g
01482
01483 integer (kind=int_kind) ::
01484 i, j, nyg,
01485 im1
01486
01487 if (my_task == master_task) then
01488 allocate(work_g2(nx_global,ny_global))
01489 else
01490 allocate(work_g2(1,1))
01491 endif
01492
01493 if (my_task == master_task) then
01494 do j = 1, ny_global
01495 do i = 1, nx_global
01496 work_g(i,j) = work_g(i,j) * cm_to_m
01497 enddo
01498 enddo
01499 do j = 1, ny_global-1
01500 do i = 1, nx_global
01501 work_g2(i,j) = p5*(work_g(i,j) + work_g(i,j+1))
01502 enddo
01503 enddo
01504
01505
01506 nyg = ny_global
01507 do i = 1, nx_global
01508 work_g2(i,nyg) = c2*work_g(i,nyg-1) &
01509 - work_g(i,nyg-2)
01510 enddo
01511 endif
01512 call scatter_global(HTE, work_g, master_task, distrb_info, &
01513 field_loc_Eface, field_type_scalar)
01514 call scatter_global(dyu, work_g2, master_task, distrb_info, &
01515 field_loc_NEcorner, field_type_scalar)
01516
01517 if (my_task == master_task) then
01518 do j = 1, ny_global
01519 do i = 1, nx_global
01520
01521 im1 = i-1
01522 if (i == 1) im1 = nx_global
01523 work_g2(i,j) = p5*(work_g(i,j) + work_g(im1,j))
01524 enddo
01525 enddo
01526 endif
01527 call scatter_global(dyt, work_g2, master_task, distrb_info, &
01528 field_loc_center, field_type_scalar)
01529
01530 deallocate(work_g2)
01531
01532 end subroutine primary_grid_lengths_HTE
01533
01534
01535
01536
01537
01538
01539
01540
01541 subroutine makemask
01542
01543
01544
01545
01546
01547
01548
01549
01550
01551
01552
01553
01554
01555
01556
01557
01558
01559 integer (kind=int_kind) ::
01560 i, j, iblk,
01561 ilo,ihi,jlo,jhi
01562
01563 type (block) ::
01564 this_block
01565
01566 call ice_timer_start(timer_bound)
01567 call ice_HaloUpdate (hm, halo_info, &
01568 field_loc_center, field_type_scalar)
01569 call ice_timer_stop(timer_bound)
01570
01571
01572
01573
01574
01575
01576 do iblk = 1, nblocks
01577 this_block = get_block(blocks_ice(iblk),iblk)
01578 ilo = this_block%ilo
01579 ihi = this_block%ihi
01580 jlo = this_block%jlo
01581 jhi = this_block%jhi
01582
01583 do j = jlo, jhi
01584 do i = ilo, ihi
01585 uvm(i,j,iblk) = min (hm(i,j, iblk), hm(i+1,j, iblk), &
01586 hm(i,j+1,iblk), hm(i+1,j+1,iblk))
01587 enddo
01588 enddo
01589 enddo
01590
01591
01592 call ice_timer_start(timer_bound)
01593 call ice_HaloUpdate (uvm, halo_info, &
01594 field_loc_NEcorner, field_type_scalar)
01595 call ice_timer_stop(timer_bound)
01596
01597
01598 do iblk = 1, nblocks
01599 do j = 1, ny_block
01600 do i = 1, nx_block
01601 tmask(i,j,iblk) = .false.
01602 umask(i,j,iblk) = .false.
01603 if ( hm(i,j,iblk) > p5) tmask(i,j,iblk) = .true.
01604 if (uvm(i,j,iblk) > p5) umask(i,j,iblk) = .true.
01605 enddo
01606 enddo
01607
01608
01609
01610
01611
01612 lmask_n(:,:,iblk) = .false.
01613 lmask_s(:,:,iblk) = .false.
01614
01615 tarean(:,:,iblk) = c0
01616 tareas(:,:,iblk) = c0
01617
01618 do j = 1, ny_block
01619 do i = 1, nx_block
01620
01621 if (ULAT(i,j,iblk) >= -puny) lmask_n(i,j,iblk) = .true.
01622 if (ULAT(i,j,iblk) < -puny) lmask_s(i,j,iblk) = .true.
01623
01624
01625 if (lmask_n(i,j,iblk)) tarean(i,j,iblk) = tarea(i,j,iblk) &
01626 * hm(i,j,iblk)
01627
01628
01629 if (lmask_s(i,j,iblk)) tareas(i,j,iblk) = tarea(i,j,iblk) &
01630 * hm(i,j,iblk)
01631
01632 enddo
01633 enddo
01634
01635 enddo
01636
01637
01638 end subroutine makemask
01639
01640
01641
01642
01643
01644
01645
01646
01647 subroutine Tlatlon
01648
01649
01650
01651
01652
01653
01654
01655
01656
01657
01658
01659
01660 use ice_domain_size
01661 use ice_global_reductions, only: global_minval, global_maxval
01662
01663
01664
01665
01666
01667 save
01668
01669 integer (kind=int_kind) ::
01670 i, j, iblk ,
01671 ig, jg ,
01672 im1 ,
01673 ilo,ihi,jlo,jhi
01674
01675 real (kind=dbl_kind) ::
01676 z1,x1,y1,z2,x2,y2,z3,x3,y3,z4,x4,y4,tx,ty,tz,da
01677
01678 type (block) ::
01679 this_block
01680
01681 TLAT(:,:,:) = c0
01682 TLON(:,:,:) = c0
01683
01684
01685
01686
01687 do iblk = 1, nblocks
01688 this_block = get_block(blocks_ice(iblk),iblk)
01689 ilo = this_block%ilo
01690 ihi = this_block%ihi
01691 jlo = this_block%jlo
01692 jhi = this_block%jhi
01693
01694 do j = jlo, jhi
01695 do i = ilo, ihi
01696
01697 z1 = cos(ULAT(i-1,j-1,iblk))
01698 x1 = cos(ULON(i-1,j-1,iblk))*z1
01699 y1 = sin(ULON(i-1,j-1,iblk))*z1
01700 z1 = sin(ULAT(i-1,j-1,iblk))
01701
01702 z2 = cos(ULAT(i,j-1,iblk))
01703 x2 = cos(ULON(i,j-1,iblk))*z2
01704 y2 = sin(ULON(i,j-1,iblk))*z2
01705 z2 = sin(ULAT(i,j-1,iblk))
01706
01707 z3 = cos(ULAT(i-1,j,iblk))
01708 x3 = cos(ULON(i-1,j,iblk))*z3
01709 y3 = sin(ULON(i-1,j,iblk))*z3
01710 z3 = sin(ULAT(i-1,j,iblk))
01711
01712 z4 = cos(ULAT(i,j,iblk))
01713 x4 = cos(ULON(i,j,iblk))*z4
01714 y4 = sin(ULON(i,j,iblk))*z4
01715 z4 = sin(ULAT(i,j,iblk))
01716
01717 tx = (x1+x2+x3+x4)/c4
01718 ty = (y1+y2+y3+y4)/c4
01719 tz = (z1+z2+z3+z4)/c4
01720 da = sqrt(tx**2+ty**2+tz**2)
01721
01722 tz = tz/da
01723
01724
01725 TLON(i,j,iblk) = c0
01726 if (tx /= c0 .or. ty /= c0) TLON(i,j,iblk) = atan2(ty,tx)
01727
01728
01729 TLAT(i,j,iblk) = asin(tz)
01730
01731 enddo
01732 enddo
01733 enddo
01734
01735
01736 call ice_timer_start(timer_bound)
01737 call ice_HaloUpdate (TLON, halo_info, &
01738 field_loc_center, field_type_scalar, &
01739 fillValue=c1)
01740 call ice_HaloUpdate (TLAT, halo_info, &
01741 field_loc_center, field_type_scalar, &
01742 fillValue=c1)
01743 call ice_HaloExtrapolate(TLON, distrb_info, &
01744 ew_boundary_type, ns_boundary_type)
01745 call ice_HaloExtrapolate(TLAT, distrb_info, &
01746 ew_boundary_type, ns_boundary_type)
01747 call ice_timer_stop(timer_bound)
01748
01749 x1 = global_minval(TLON, distrb_info, tmask)
01750 x2 = global_maxval(TLON, distrb_info, tmask)
01751 x3 = global_minval(TLAT, distrb_info, tmask)
01752 x4 = global_maxval(TLAT, distrb_info, tmask)
01753
01754 y1 = global_minval(ULON, distrb_info, umask)
01755 y2 = global_maxval(ULON, distrb_info, umask)
01756 y3 = global_minval(ULAT, distrb_info, umask)
01757 y4 = global_maxval(ULAT, distrb_info, umask)
01758
01759 if (my_task==master_task) then
01760 write(nu_diag,*) ' '
01761 write(nu_diag,*) 'min/max ULON:', y1*rad_to_deg, y2*rad_to_deg
01762 write(nu_diag,*) 'min/max TLON:', x1*rad_to_deg, x2*rad_to_deg
01763 write(nu_diag,*) 'min/max ULAT:', y3*rad_to_deg, y4*rad_to_deg
01764 write(nu_diag,*) 'min/max TLAT:', x3*rad_to_deg, x4*rad_to_deg
01765 endif
01766
01767 end subroutine Tlatlon
01768
01769
01770
01771
01772
01773
01774
01775
01776 subroutine t2ugrid_vector (work)
01777
01778
01779
01780
01781
01782
01783
01784
01785
01786
01787
01788
01789
01790
01791 real (kind=dbl_kind), dimension(nx_block,ny_block,max_blocks),
01792 intent(inout) ::
01793 work
01794
01795 integer (int_kind) :: iblk
01796
01797
01798
01799 work1(:,:,:) = work(:,:,:)
01800
01801
01802
01803
01804 call ice_timer_start(timer_bound)
01805 call ice_HaloUpdate (work1, halo_info, &
01806 field_loc_center, field_type_vector)
01807 call ice_timer_stop(timer_bound)
01808
01809 call to_ugrid(work1,work)
01810
01811 end subroutine t2ugrid_vector
01812
01813
01814
01815
01816
01817
01818
01819
01820 subroutine to_ugrid(work1,work2)
01821
01822
01823
01824
01825
01826
01827
01828
01829
01830
01831
01832
01833
01834
01835
01836
01837 real (kind=dbl_kind), intent(in) ::
01838 work1(nx_block,ny_block,max_blocks)
01839
01840 real (kind=dbl_kind), intent(out) ::
01841 work2(nx_block,ny_block,max_blocks)
01842
01843 type (block) ::
01844 this_block
01845
01846
01847
01848 integer (kind=int_kind) ::
01849 i, j, iblk,
01850 ilo,ihi,jlo,jhi
01851
01852 work2(:,:,:) = c0
01853
01854
01855 do iblk = 1, nblocks
01856 this_block = get_block(blocks_ice(iblk),iblk)
01857 ilo = this_block%ilo
01858 ihi = this_block%ihi
01859 jlo = this_block%jlo
01860 jhi = this_block%jhi
01861
01862 do j = jlo, jhi
01863 do i = ilo, ihi
01864 work2(i,j,iblk) = p25 * &
01865 (work1(i, j, iblk)*tarea(i, j, iblk) &
01866 + work1(i+1,j, iblk)*tarea(i+1,j, iblk) &
01867 + work1(i, j+1,iblk)*tarea(i, j+1,iblk) &
01868 + work1(i+1,j+1,iblk)*tarea(i+1,j+1,iblk)) &
01869 / uarea(i, j, iblk)
01870 enddo
01871 enddo
01872 enddo
01873
01874
01875 end subroutine to_ugrid
01876
01877
01878
01879
01880
01881
01882
01883
01884 subroutine u2tgrid_vector (work)
01885
01886
01887
01888
01889
01890
01891
01892
01893
01894
01895
01896
01897
01898
01899
01900
01901 real (kind=dbl_kind), dimension (nx_block,ny_block,max_blocks),
01902 intent(inout) ::
01903 work
01904
01905
01906
01907 work1(:,:,:) = work(:,:,:)
01908
01909 call ice_timer_start(timer_bound)
01910 call ice_HaloUpdate (work1, halo_info, &
01911 field_loc_NEcorner, field_type_vector)
01912 call ice_timer_stop(timer_bound)
01913
01914 call to_tgrid(work1,work)
01915
01916 end subroutine u2tgrid_vector
01917
01918
01919
01920
01921
01922
01923
01924
01925 subroutine to_tgrid(work1, work2)
01926
01927
01928
01929
01930
01931
01932
01933
01934
01935
01936
01937
01938
01939
01940
01941
01942 real (kind=dbl_kind) :: work1(nx_block,ny_block,max_blocks),
01943 work2(nx_block,ny_block,max_blocks)
01944
01945
01946
01947 integer (kind=int_kind) ::
01948 i, j, iblk,
01949 ilo,ihi,jlo,jhi
01950
01951 type (block) ::
01952 this_block
01953
01954
01955 do iblk = 1, nblocks
01956 this_block = get_block(blocks_ice(iblk),iblk)
01957 ilo = this_block%ilo
01958 ihi = this_block%ihi
01959 jlo = this_block%jlo
01960 jhi = this_block%jhi
01961
01962 do j = jlo, jhi
01963 do i = ilo, ihi
01964 work2(i,j,iblk) = p25 * &
01965 (work1(i, j ,iblk) * uarea(i, j, iblk) &
01966 + work1(i-1,j ,iblk) * uarea(i-1,j, iblk) &
01967 + work1(i, j-1,iblk) * uarea(i, j-1,iblk) &
01968 + work1(i-1,j-1,iblk) * uarea(i-1,j-1,iblk)) &
01969 / tarea(i, j, iblk)
01970 enddo
01971 enddo
01972 enddo
01973
01974
01975 end subroutine to_tgrid
01976
01977
01978
01979
01980
01981
01982
01983
01984 subroutine sinlat(a, k)
01985
01986
01987
01988
01989
01990
01991
01992
01993
01994
01995
01996
01997
01998
01999
02000 use ice_exit
02001
02002
02003
02004
02005
02006 real (kind=dbl_kind), dimension(ny_global), intent(out) ::
02007 a
02008
02009 integer (kind=int_kind), intent(in) ::
02010 k
02011
02012
02013
02014 real (kind=dbl_kind), dimension(k) ::
02015 sinlats
02016
02017 real (kind=dbl_kind) ::
02018 eps ,
02019 c ,
02020 fk ,
02021 xz ,
02022 pkm1 ,
02023 pkm2 ,
02024 pkmrk ,
02025 pk ,
02026 sp ,
02027 avsp ,
02028 fn ,
02029 avsp_prev,
02030 testdiff
02031
02032 real (kind=dbl_kind), parameter ::
02033 eps27 = 1.e-27_dbl_kind
02034
02035 integer (kind=int_kind) ::
02036 n, l ,
02037 iter ,
02038 is ,
02039 kk
02040
02041
02042
02043
02044
02045 c = (c1-(c2/pi)**2)*p25
02046 fk = k
02047 kk = k/2
02048 call bsslzr(sinlats,kk)
02049 do is=1,kk
02050 xz = cos(sinlats(is)/sqrt((fk+p5)**2+c))
02051
02052
02053
02054 iter = 0
02055 avsp = c100
02056 10 continue
02057 avsp_prev = avsp
02058 pkm2 = c1
02059 pkm1 = xz
02060 iter = iter + 1
02061 if (iter > 100) then
02062 call abort_ice('SINLAT: no convergence in 100 iterations')
02063 end if
02064
02065
02066
02067 do n=2,k
02068 fn = n
02069 pk = ((c2*fn-1._dbl_kind)*xz*pkm1-(fn-c1)*pkm2)/fn
02070 pkm2 = pkm1
02071 pkm1 = pk
02072 enddo
02073 pkm1 = pkm2
02074 pkmrk = (fk*(pkm1-xz*pk))/(c1-xz**2)
02075 sp = pk/pkmrk
02076 xz = xz - sp
02077 avsp = abs(sp)
02078 testdiff = avsp_prev - avsp
02079 if (testdiff > eps27) go to 10
02080 sinlats(is) = xz
02081 end do
02082
02083
02084
02085
02086 do n=1,kk
02087 l = k + 1 - n
02088 a(n) = sinlats(n)
02089 a(l) = -sinlats(n)
02090 end do
02091
02092 end subroutine sinlat
02093
02094
02095
02096
02097
02098
02099
02100
02101 subroutine bsslzr(bes, n)
02102
02103
02104
02105
02106
02107
02108
02109
02110
02111
02112
02113
02114
02115
02116
02117
02118
02119
02120
02121
02122
02123
02124
02125
02126 integer (kind=int_kind), intent(in) ::
02127 n
02128
02129 real (kind=dbl_kind), dimension(n), intent(inout) ::
02130 bes
02131
02132
02133
02134
02135
02136
02137
02138 integer (kind=int_kind) ::
02139 nn, j
02140
02141 real (kind=dbl_kind), dimension(50) :: bz
02142
02143 save bz
02144
02145
02146
02147
02148
02149 data bz/ &
02150 2.4048255577_dbl_kind, 5.5200781103_dbl_kind, 8.6537279129_dbl_kind, &
02151 11.7915344391_dbl_kind, 14.9309177086_dbl_kind, 18.0710639679_dbl_kind, &
02152 21.2116366299_dbl_kind, 24.3524715308_dbl_kind, 27.4934791320_dbl_kind, &
02153 30.6346064684_dbl_kind, 33.7758202136_dbl_kind, 36.9170983537_dbl_kind, &
02154 40.0584257646_dbl_kind, 43.1997917132_dbl_kind, 46.3411883717_dbl_kind, &
02155 49.4826098974_dbl_kind, 52.6240518411_dbl_kind, 55.7655107550_dbl_kind, &
02156 58.9069839261_dbl_kind, 62.0484691902_dbl_kind, 65.1899648002_dbl_kind, &
02157 68.3314693299_dbl_kind, 71.4729816036_dbl_kind, 74.6145006437_dbl_kind, &
02158 77.7560256304_dbl_kind, 80.8975558711_dbl_kind, 84.0390907769_dbl_kind, &
02159 87.1806298436_dbl_kind, 90.3221726372_dbl_kind, 93.4637187819_dbl_kind, &
02160 96.6052679510_dbl_kind, 99.7468198587_dbl_kind, 102.8883742542_dbl_kind, &
02161 106.0299309165_dbl_kind, 109.1714896498_dbl_kind, 112.3130502805_dbl_kind, &
02162 115.4546126537_dbl_kind, 118.5961766309_dbl_kind, 121.7377420880_dbl_kind, &
02163 124.8793089132_dbl_kind, 128.0208770059_dbl_kind, 131.1624462752_dbl_kind, &
02164 134.3040166383_dbl_kind, 137.4455880203_dbl_kind, 140.5871603528_dbl_kind, &
02165 143.7287335737_dbl_kind, 146.8703076258_dbl_kind, 150.0118824570_dbl_kind, &
02166 153.1534580192_dbl_kind, 156.2950342685_dbl_kind/
02167
02168 nn = n
02169 if (n > 50) then
02170 bes(50) = bz(50)
02171 do j = 51, n
02172 bes(j) = bes(j-1) + pi
02173 end do
02174 nn = 49
02175 endif
02176 bes(:nn) = bz(:nn)
02177
02178 end subroutine bsslzr
02179
02180
02181
02182
02183
02184
02185
02186
02187
02188
02189
02190
02191 subroutine gridbox_corners
02192
02193
02194
02195
02196
02197
02198
02199
02200
02201
02202
02203
02204
02205 use ice_work, only: work_g2
02206
02207
02208
02209
02210
02211 integer (kind=int_kind) ::
02212 i,j,iblk,icorner,
02213 ilo,ihi,jlo,jhi
02214
02215 type (block) ::
02216 this_block
02217
02218
02219
02220
02221
02222
02223
02224 do iblk = 1, nblocks
02225 this_block = get_block(blocks_ice(iblk),iblk)
02226 ilo = this_block%ilo
02227 ihi = this_block%ihi
02228 jlo = this_block%jlo
02229 jhi = this_block%jhi
02230
02231 do j = jlo, jhi
02232 do i = ilo, ihi
02233
02234 latu_bounds(1,i,j,iblk)=TLAT(i ,j ,iblk)*rad_to_deg
02235 latu_bounds(2,i,j,iblk)=TLAT(i+1,j ,iblk)*rad_to_deg
02236 latu_bounds(3,i,j,iblk)=TLAT(i+1,j+1,iblk)*rad_to_deg
02237 latu_bounds(4,i,j,iblk)=TLAT(i ,j+1,iblk)*rad_to_deg
02238
02239 lonu_bounds(1,i,j,iblk)=TLON(i ,j ,iblk)*rad_to_deg
02240 lonu_bounds(2,i,j,iblk)=TLON(i+1,j ,iblk)*rad_to_deg
02241 lonu_bounds(3,i,j,iblk)=TLON(i+1,j+1,iblk)*rad_to_deg
02242 lonu_bounds(4,i,j,iblk)=TLON(i ,j+1,iblk)*rad_to_deg
02243
02244 enddo
02245 enddo
02246 enddo
02247
02248
02249
02250
02251
02252
02253 if (my_task == master_task) then
02254 allocate(work_g2(nx_global,ny_global))
02255 else
02256 allocate(work_g2(1,1))
02257 endif
02258
02259 work1(:,:,:) = latu_bounds(2,:,:,:)
02260 call gather_global(work_g2, work1, master_task, distrb_info)
02261 if (my_task == master_task) then
02262 do j = 1, ny_global
02263 work_g2(nx_global,j) = c2*work_g2(nx_global-1,j) &
02264 - work_g2(nx_global-2,j)
02265 enddo
02266 endif
02267 call scatter_global(work1, work_g2, &
02268 master_task, distrb_info, &
02269 field_loc_NEcorner, field_type_scalar)
02270 latu_bounds(2,:,:,:) = work1(:,:,:)
02271
02272 work1(:,:,:) = latu_bounds(3,:,:,:)
02273 call gather_global(work_g2, work1, master_task, distrb_info)
02274 if (my_task == master_task) then
02275 do i = 1, nx_global
02276 work_g2(i,ny_global) = c2*work_g2(i,ny_global-1) &
02277 - work_g2(i,ny_global-2)
02278 enddo
02279 do j = 1, ny_global
02280 work_g2(nx_global,j) = c2*work_g2(nx_global-1,j) &
02281 - work_g2(nx_global-2,j)
02282 enddo
02283 endif
02284 call scatter_global(work1, work_g2, &
02285 master_task, distrb_info, &
02286 field_loc_NEcorner, field_type_scalar)
02287 latu_bounds(3,:,:,:) = work1(:,:,:)
02288
02289 work1(:,:,:) = latu_bounds(4,:,:,:)
02290 call gather_global(work_g2, work1, master_task, distrb_info)
02291 if (my_task == master_task) then
02292 do i = 1, nx_global
02293 work_g2(i,ny_global) = c2*work_g2(i,ny_global-1) &
02294 - work_g2(i,ny_global-2)
02295 enddo
02296 endif
02297 call scatter_global(work1, work_g2, &
02298 master_task, distrb_info, &
02299 field_loc_NEcorner, field_type_scalar)
02300 latu_bounds(4,:,:,:) = work1(:,:,:)
02301
02302 work1(:,:,:) = lonu_bounds(2,:,:,:)
02303 call gather_global(work_g2, work1, master_task, distrb_info)
02304 if (my_task == master_task) then
02305 do j = 1, ny_global
02306 work_g2(nx_global,j) = c2*work_g2(nx_global-1,j) &
02307 - work_g2(nx_global-2,j)
02308 enddo
02309 endif
02310 call scatter_global(work1, work_g2, &
02311 master_task, distrb_info, &
02312 field_loc_NEcorner, field_type_scalar)
02313 lonu_bounds(2,:,:,:) = work1(:,:,:)
02314
02315 work1(:,:,:) = lonu_bounds(3,:,:,:)
02316 call gather_global(work_g2, work1, master_task, distrb_info)
02317 if (my_task == master_task) then
02318 do i = 1, nx_global
02319 work_g2(i,ny_global) = c2*work_g2(i,ny_global-1) &
02320 - work_g2(i,ny_global-2)
02321 enddo
02322 do j = 1, ny_global
02323 work_g2(nx_global,j) = c2*work_g2(nx_global-1,j) &
02324 - work_g2(nx_global-2,j)
02325 enddo
02326 endif
02327 call scatter_global(work1, work_g2, &
02328 master_task, distrb_info, &
02329 field_loc_NEcorner, field_type_scalar)
02330 lonu_bounds(3,:,:,:) = work1(:,:,:)
02331
02332 work1(:,:,:) = lonu_bounds(4,:,:,:)
02333 call gather_global(work_g2, work1, master_task, distrb_info)
02334 if (my_task == master_task) then
02335 do i = 1, nx_global
02336 work_g2(i,ny_global) = c2*work_g2(i,ny_global-1) &
02337 - work_g2(i,ny_global-2)
02338 enddo
02339 endif
02340 call scatter_global(work1, work_g2, &
02341 master_task, distrb_info, &
02342 field_loc_NEcorner, field_type_scalar)
02343 lonu_bounds(4,:,:,:) = work1(:,:,:)
02344
02345 deallocate(work_g2)
02346
02347
02348
02349
02350
02351 allocate(work_g2(nx_block,ny_block))
02352
02353
02354
02355 do iblk = 1, nblocks
02356 do icorner = 1, 4
02357 work_g2(:,:) = lont_bounds(icorner,:,:,iblk) + c360
02358 where (work_g2 > c360) work_g2 = work_g2 - c360
02359 where (work_g2 < c0 ) work_g2 = work_g2 + c360
02360 lont_bounds(icorner,:,:,iblk) = work_g2(:,:)
02361
02362 work_g2(:,:) = lonu_bounds(icorner,:,:,iblk) + c360
02363 where (work_g2 > c360) work_g2 = work_g2 - c360
02364 where (work_g2 < c0 ) work_g2 = work_g2 + c360
02365 lonu_bounds(icorner,:,:,iblk) = work_g2(:,:)
02366 enddo
02367 enddo
02368
02369 deallocate(work_g2)
02370
02371 end subroutine gridbox_corners
02372
02373
02374
02375
02376
02377
02378
02379
02380
02381 subroutine gridbox_verts(work_g,vbounds)
02382
02383
02384
02385
02386
02387
02388
02389
02390
02391
02392
02393
02394
02395 use ice_work, only: work_g2
02396
02397
02398
02399
02400
02401 real (kind=dbl_kind), dimension(:,:), intent(in) :: work_g
02402
02403 real (kind=dbl_kind),
02404 dimension(4,nx_block,ny_block,max_blocks),
02405 intent(out) :: vbounds
02406
02407 integer (kind=int_kind) ::
02408 i,j,
02409 ilo,ihi,jlo,jhi
02410
02411 type (block) ::
02412 this_block
02413
02414 if (my_task == master_task) then
02415 allocate(work_g2(nx_global,ny_global))
02416 else
02417 allocate(work_g2(1,1))
02418 endif
02419
02420
02421
02422
02423
02424
02425 work_g2(:,:) = c0
02426 if (my_task == master_task) then
02427 do j = 2, ny_global
02428 do i = 2, nx_global
02429 work_g2(i,j) = work_g(i-1,j-1) * rad_to_deg
02430 enddo
02431 enddo
02432
02433 do j = 1, ny_global
02434 work_g2(1,j) = c2*work_g2(2,j) - work_g2(3,j)
02435 enddo
02436 do i = 1, nx_global
02437 work_g2(i,1) = c2*work_g2(i,2) - work_g2(i,3)
02438 enddo
02439 endif
02440 call scatter_global(work1, work_g2, &
02441 master_task, distrb_info, &
02442 field_loc_NEcorner, field_type_scalar)
02443 vbounds(1,:,:,:) = work1(:,:,:)
02444
02445 work_g2(:,:) = c0
02446 if (my_task == master_task) then
02447 do j = 2, ny_global
02448 do i = 1, nx_global
02449 work_g2(i,j) = work_g(i,j-1) * rad_to_deg
02450 enddo
02451 enddo
02452
02453 do i = 1, nx_global
02454 work_g2(i,1) = (c2*work_g2(i,2) - work_g2(i,3))
02455 enddo
02456 endif
02457 call scatter_global(work1, work_g2, &
02458 master_task, distrb_info, &
02459 field_loc_NEcorner, field_type_scalar)
02460 vbounds(2,:,:,:) = work1(:,:,:)
02461
02462 work_g2(:,:) = c0
02463 if (my_task == master_task) then
02464 do j = 1, ny_global
02465 do i = 1, nx_global
02466 work_g2(i,j) = work_g(i,j) * rad_to_deg
02467 enddo
02468 enddo
02469 endif
02470 call scatter_global(work1, work_g2, &
02471 master_task, distrb_info, &
02472 field_loc_NEcorner, field_type_scalar)
02473 vbounds(3,:,:,:) = work1(:,:,:)
02474
02475 work_g2(:,:) = c0
02476 if (my_task == master_task) then
02477 do j = 1, ny_global
02478 do i = 2, nx_global
02479 work_g2(i,j) = work_g(i-1,j ) * rad_to_deg
02480 enddo
02481 enddo
02482
02483 do j = 1, ny_global
02484 work_g2(1,j) = c2*work_g2(2,j) - work_g2(3,j)
02485 enddo
02486 endif
02487 call scatter_global(work1, work_g2, &
02488 master_task, distrb_info, &
02489 field_loc_NEcorner, field_type_scalar)
02490 vbounds(4,:,:,:) = work1(:,:,:)
02491
02492 deallocate (work_g2)
02493
02494 end subroutine gridbox_verts
02495
02496
02497
02498 end module ice_grid
02499
02500