Changeset 5260 for branches/2014/dev_r4650_UKMO10_Tidally_Meaned_Diagnostics/NEMOGCM/NEMO/OPA_SRC/DOM
- Timestamp:
- 2015-05-12T12:37:15+02:00 (9 years ago)
- Location:
- branches/2014/dev_r4650_UKMO10_Tidally_Meaned_Diagnostics/NEMOGCM/NEMO/OPA_SRC/DOM
- Files:
-
- 12 edited
-
daymod.F90 (modified) (1 diff)
-
dom_oce.F90 (modified) (7 diffs)
-
domain.F90 (modified) (2 diffs)
-
domcfg.F90 (modified) (1 diff)
-
domhgr.F90 (modified) (2 diffs)
-
dommsk.F90 (modified) (3 diffs)
-
domvvl.F90 (modified) (24 diffs)
-
domwri.F90 (modified) (4 diffs)
-
domzgr.F90 (modified) (19 diffs)
-
dtatsd.F90 (modified) (2 diffs, 1 prop)
-
istate.F90 (modified) (18 diffs)
-
phycst.F90 (modified) (4 diffs)
Legend:
- Unmodified
- Added
- Removed
-
branches/2014/dev_r4650_UKMO10_Tidally_Meaned_Diagnostics/NEMOGCM/NEMO/OPA_SRC/DOM/daymod.F90
r4162 r5260 123 123 124 124 ! control print 125 IF(lwp) WRITE(numout,'(a,i6,a,i2,a,i2,a,i 6)')' ==============>> 1/2 time step before the start of the run DATE Y/M/D = ', &125 IF(lwp) WRITE(numout,'(a,i6,a,i2,a,i2,a,i8,a,i8)')' =======>> 1/2 time step before the start of the run DATE Y/M/D = ', & 126 126 & nyear, '/', nmonth, '/', nday, ' nsec_day:', nsec_day, ' nsec_week:', nsec_week 127 127 -
branches/2014/dev_r4650_UKMO10_Tidally_Meaned_Diagnostics/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90
r4488 r5260 153 153 INTEGER, PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: nldit , nldjt !: first, last indoor index for each i-domain 154 154 INTEGER, PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: nleit , nlejt !: first, last indoor index for each j-domain 155 INTEGER, PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: nfiimpp, nfipproc, nfilcit 155 156 156 157 !!---------------------------------------------------------------------- … … 161 162 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: gphit, gphiu !: latitude of t-, u-, v- and f-points (degre) 162 163 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: gphiv, gphif !: 163 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, TARGET, DIMENSION(:,:) :: e1t, e2t !: horizontal scale factorsat t-point (m)164 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, TARGET, DIMENSION(:,:) :: e1u, e2u !: horizontal scale factorsat u-point (m)165 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, TARGET, DIMENSION(:,:) :: e1v, e2v !: horizontal scale factorsat v-point (m)166 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, TARGET, DIMENSION(:,:) :: e1f, e2f !: horizontal scale factorsat f-point (m)164 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, TARGET, DIMENSION(:,:) :: e1t, e2t, r1_e1t, r1_e2t !: horizontal scale factors and inverse at t-point (m) 165 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, TARGET, DIMENSION(:,:) :: e1u, e2u, r1_e1u, r1_e2u !: horizontal scale factors and inverse at u-point (m) 166 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, TARGET, DIMENSION(:,:) :: e1v, e2v, r1_e1v, r1_e2v !: horizontal scale factors and inverse at v-point (m) 167 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, TARGET, DIMENSION(:,:) :: e1f, e2f, r1_e1f, r1_e2f !: horizontal scale factors and inverse at f-point (m) 167 168 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: e1e2t !: surface at t-point (m2) 168 169 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ff !: coriolis factor (2.*omega*sin(yphi) ) (s-1) … … 175 176 LOGICAL, PUBLIC :: ln_zps !: z-coordinate - partial step 176 177 LOGICAL, PUBLIC :: ln_sco !: s-coordinate or hybrid z-s coordinate 178 LOGICAL, PUBLIC :: ln_isfcav !: presence of ISF 177 179 178 180 !! All coordinates … … 250 252 INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: mbkt !: vertical index of the bottom last T- ocean level 251 253 INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: mbku, mbkv !: vertical index of the bottom last U- and W- ocean level 252 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: bathy !: ocean depth (meters) 253 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: tmask_i !: interior domain T-point mask 254 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: bmask !: land/ocean mask of barotropic stream function 254 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: bathy !: ocean depth (meters) 255 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: tmask_i, umask_i, vmask_i, fmask_i !: interior domain T-point mask 256 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: bmask !: land/ocean mask of barotropic stream function 257 258 INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: misfdep !: top first ocean level (ISF) 259 INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: mikt, miku, mikv, mikf !: first wet T-, U-, V-, F- ocean level (ISF) 260 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: risfdep !: Iceshelf draft (ISF) 261 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ssmask !: surface domain T-point mask 255 262 256 263 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:), TARGET :: tmask, umask, vmask, fmask !: land/ocean mask at T-, U-, V- and F-pts 264 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:), TARGET :: wmask, wumask, wvmask !: land/ocean mask at WT-, WU- and WV-pts 257 265 258 266 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: tpol, fpol !: north fold mask (jperio= 3 or 4) … … 325 333 INTEGER FUNCTION dom_oce_alloc() 326 334 !!---------------------------------------------------------------------- 327 INTEGER, DIMENSION(1 1) :: ierr335 INTEGER, DIMENSION(12) :: ierr 328 336 !!---------------------------------------------------------------------- 329 337 ierr(:) = 0 330 338 ! 331 ALLOCATE( rdttra(jpk), r2dtra(jpk), mig(jpi), mjg(jpj), STAT=ierr(1) ) 339 ALLOCATE( rdttra(jpk), r2dtra(jpk), mig(jpi), mjg(jpj), nfiimpp(jpni,jpnj), & 340 & nfipproc(jpni,jpnj), nfilcit(jpni,jpnj), STAT=ierr(1) ) 332 341 ! 333 342 ALLOCATE( nimppt(jpnij) , ibonit(jpnij) , nlcit(jpnij) , nlcjt(jpnij) , & … … 337 346 & tpol(jpiglo) , fpol(jpiglo) , STAT=ierr(2) ) 338 347 ! 339 ALLOCATE( glamt(jpi,jpj) , gphit(jpi,jpj) , e1t(jpi,jpj) , e2t(jpi,jpj) , & 340 & glamu(jpi,jpj) , gphiu(jpi,jpj) , e1u(jpi,jpj) , e2u(jpi,jpj) , & 341 & glamv(jpi,jpj) , gphiv(jpi,jpj) , e1v(jpi,jpj) , e2v(jpi,jpj) , e1e2t(jpi,jpj) , & 342 & glamf(jpi,jpj) , gphif(jpi,jpj) , e1f(jpi,jpj) , e2f(jpi,jpj) , ff (jpi,jpj) , STAT=ierr(3) ) 348 ALLOCATE( glamt(jpi,jpj) , gphit(jpi,jpj) , e1t(jpi,jpj) , e2t(jpi,jpj) , r1_e1t(jpi,jpj) , r1_e2t(jpi,jpj) , & 349 & glamu(jpi,jpj) , gphiu(jpi,jpj) , e1u(jpi,jpj) , e2u(jpi,jpj) , r1_e1u(jpi,jpj) , r1_e2u(jpi,jpj) , & 350 & glamv(jpi,jpj) , gphiv(jpi,jpj) , e1v(jpi,jpj) , e2v(jpi,jpj) , r1_e1v(jpi,jpj) , r1_e2v(jpi,jpj) , & 351 & glamf(jpi,jpj) , gphif(jpi,jpj) , e1f(jpi,jpj) , e2f(jpi,jpj) , r1_e1f(jpi,jpj) , r1_e2f(jpi,jpj) , & 352 & e1e2t(jpi,jpj) , ff (jpi,jpj) , STAT=ierr(3) ) 343 353 ! 344 354 ALLOCATE( gdep3w_0(jpi,jpj,jpk) , e3v_0(jpi,jpj,jpk) , e3f_0 (jpi,jpj,jpk) , & … … 379 389 & hift (jpi,jpj) , hifu (jpi,jpj) , rx1 (jpi,jpj) , STAT=ierr(8) ) 380 390 381 ALLOCATE( mbathy(jpi,jpj) , bathy(jpi,jpj) , & 382 & tmask_i(jpi,jpj) , bmask(jpi,jpj) , & 391 ALLOCATE( mbathy(jpi,jpj) , bathy(jpi,jpj) , & 392 & tmask_i(jpi,jpj) , umask_i(jpi,jpj), vmask_i(jpi,jpj), fmask_i(jpi,jpj), & 393 & bmask(jpi,jpj) , & 383 394 & mbkt (jpi,jpj) , mbku (jpi,jpj) , mbkv(jpi,jpj) , STAT=ierr(9) ) 384 395 396 ! (ISF) Allocation of basic array 397 ALLOCATE( misfdep(jpi,jpj) , risfdep(jpi,jpj), & 398 & mikt(jpi,jpj), miku(jpi,jpj), mikv(jpi,jpj) , & 399 & mikf(jpi,jpj), ssmask(jpi,jpj), STAT=ierr(10) ) 400 385 401 ALLOCATE( tmask(jpi,jpj,jpk) , umask(jpi,jpj,jpk), & 386 & vmask(jpi,jpj,jpk) , fmask(jpi,jpj,jpk), STAT=ierr(10) ) 402 & vmask(jpi,jpj,jpk) , fmask(jpi,jpj,jpk), STAT=ierr(11) ) 403 404 ALLOCATE( wmask(jpi,jpj,jpk) , wumask(jpi,jpj,jpk), wvmask(jpi,jpj,jpk) , STAT=ierr(12) ) 387 405 388 406 #if defined key_noslip_accurate 389 ALLOCATE( npcoa(4,jpk), nicoa(2*(jpi+jpj),4,jpk), njcoa(2*(jpi+jpj),4,jpk), STAT=ierr(1 1) )407 ALLOCATE( npcoa(4,jpk), nicoa(2*(jpi+jpj),4,jpk), njcoa(2*(jpi+jpj),4,jpk), STAT=ierr(12) ) 390 408 #endif 391 409 ! -
branches/2014/dev_r4650_UKMO10_Tidally_Meaned_Diagnostics/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90
r4624 r5260 111 111 END DO 112 112 ! ! Inverse of the local depth 113 hur(:,:) = 1._wp / ( hu(:,:) + 1._wp - umask (:,:,1) ) * umask(:,:,1)114 hvr(:,:) = 1._wp / ( hv(:,:) + 1._wp - vmask (:,:,1) ) * vmask(:,:,1)113 hur(:,:) = 1._wp / ( hu(:,:) + 1._wp - umask_i(:,:) ) * umask_i(:,:) 114 hvr(:,:) = 1._wp / ( hv(:,:) + 1._wp - vmask_i(:,:) ) * vmask_i(:,:) 115 115 116 116 CALL dom_stp ! time step … … 365 365 ! 366 366 IF(lk_mpp) THEN 367 CALL mpp_minloc( e1t(:,:), tmask (:,:,1), ze1min, iimi1,ijmi1 )368 CALL mpp_minloc( e2t(:,:), tmask (:,:,1), ze2min, iimi2,ijmi2 )369 CALL mpp_maxloc( e1t(:,:), tmask (:,:,1), ze1max, iima1,ijma1 )370 CALL mpp_maxloc( e2t(:,:), tmask (:,:,1), ze2max, iima2,ijma2 )367 CALL mpp_minloc( e1t(:,:), tmask_i(:,:), ze1min, iimi1,ijmi1 ) 368 CALL mpp_minloc( e2t(:,:), tmask_i(:,:), ze2min, iimi2,ijmi2 ) 369 CALL mpp_maxloc( e1t(:,:), tmask_i(:,:), ze1max, iima1,ijma1 ) 370 CALL mpp_maxloc( e2t(:,:), tmask_i(:,:), ze2max, iima2,ijma2 ) 371 371 ELSE 372 ze1min = MINVAL( e1t(:,:), mask = tmask (:,:,1) == 1._wp )373 ze2min = MINVAL( e2t(:,:), mask = tmask (:,:,1) == 1._wp )374 ze1max = MAXVAL( e1t(:,:), mask = tmask (:,:,1) == 1._wp )375 ze2max = MAXVAL( e2t(:,:), mask = tmask (:,:,1) == 1._wp )376 377 iloc = MINLOC( e1t(:,:), mask = tmask (:,:,1) == 1._wp )372 ze1min = MINVAL( e1t(:,:), mask = tmask_i(:,:) == 1._wp ) 373 ze2min = MINVAL( e2t(:,:), mask = tmask_i(:,:) == 1._wp ) 374 ze1max = MAXVAL( e1t(:,:), mask = tmask_i(:,:) == 1._wp ) 375 ze2max = MAXVAL( e2t(:,:), mask = tmask_i(:,:) == 1._wp ) 376 377 iloc = MINLOC( e1t(:,:), mask = tmask_i(:,:) == 1._wp ) 378 378 iimi1 = iloc(1) + nimpp - 1 379 379 ijmi1 = iloc(2) + njmpp - 1 380 iloc = MINLOC( e2t(:,:), mask = tmask (:,:,1) == 1._wp )380 iloc = MINLOC( e2t(:,:), mask = tmask_i(:,:) == 1._wp ) 381 381 iimi2 = iloc(1) + nimpp - 1 382 382 ijmi2 = iloc(2) + njmpp - 1 383 iloc = MAXLOC( e1t(:,:), mask = tmask (:,:,1) == 1._wp )383 iloc = MAXLOC( e1t(:,:), mask = tmask_i(:,:) == 1._wp ) 384 384 iima1 = iloc(1) + nimpp - 1 385 385 ijma1 = iloc(2) + njmpp - 1 386 iloc = MAXLOC( e2t(:,:), mask = tmask (:,:,1) == 1._wp )386 iloc = MAXLOC( e2t(:,:), mask = tmask_i(:,:) == 1._wp ) 387 387 iima2 = iloc(1) + nimpp - 1 388 388 ijma2 = iloc(2) + njmpp - 1 -
branches/2014/dev_r4650_UKMO10_Tidally_Meaned_Diagnostics/NEMOGCM/NEMO/OPA_SRC/DOM/domcfg.F90
r4245 r5260 82 82 !!---------------------------------------------------------------------- 83 83 ! ! recalculate jpizoom/jpjzoom given lat/lon 84 IF( lk_c1d ) CALL dom_c1d( rn_lat1d, rn_lon1d )84 IF( lk_c1d .AND. ln_c1d_locpt ) CALL dom_c1d( rn_lat1d, rn_lon1d ) 85 85 ! 86 86 ! ! ============== ! -
branches/2014/dev_r4650_UKMO10_Tidally_Meaned_Diagnostics/NEMOGCM/NEMO/OPA_SRC/DOM/domhgr.F90
r4366 r5260 471 471 re2u_e1u(:,:) = e2u(:,:) / e1u(:,:) 472 472 re1v_e2v(:,:) = e1v(:,:) / e2v(:,:) 473 r1_e1t (:,:) = 1._wp / e1t(:,:) 474 r1_e1u (:,:) = 1._wp / e1u(:,:) 475 r1_e1v (:,:) = 1._wp / e1v(:,:) 476 r1_e1f (:,:) = 1._wp / e1f(:,:) 477 r1_e2t (:,:) = 1._wp / e2t(:,:) 478 r1_e2u (:,:) = 1._wp / e2u(:,:) 479 r1_e2v (:,:) = 1._wp / e2v(:,:) 480 r1_e2f (:,:) = 1._wp / e2f(:,:) 473 481 474 482 ! Control printing : Grid informations (if not restart) … … 616 624 CALL iom_open( 'coordinates', inum ) 617 625 618 CALL iom_get( inum, jpdom_data, 'glamt', glamt )619 CALL iom_get( inum, jpdom_data, 'glamu', glamu )620 CALL iom_get( inum, jpdom_data, 'glamv', glamv )621 CALL iom_get( inum, jpdom_data, 'glamf', glamf )622 623 CALL iom_get( inum, jpdom_data, 'gphit', gphit )624 CALL iom_get( inum, jpdom_data, 'gphiu', gphiu )625 CALL iom_get( inum, jpdom_data, 'gphiv', gphiv )626 CALL iom_get( inum, jpdom_data, 'gphif', gphif )627 628 CALL iom_get( inum, jpdom_data, 'e1t', e1t )629 CALL iom_get( inum, jpdom_data, 'e1u', e1u )630 CALL iom_get( inum, jpdom_data, 'e1v', e1v )631 CALL iom_get( inum, jpdom_data, 'e1f', e1f )632 633 CALL iom_get( inum, jpdom_data, 'e2t', e2t )634 CALL iom_get( inum, jpdom_data, 'e2u', e2u )635 CALL iom_get( inum, jpdom_data, 'e2v', e2v )636 CALL iom_get( inum, jpdom_data, 'e2f', e2f )626 CALL iom_get( inum, jpdom_data, 'glamt', glamt, lrowattr=ln_use_jattr ) 627 CALL iom_get( inum, jpdom_data, 'glamu', glamu, lrowattr=ln_use_jattr ) 628 CALL iom_get( inum, jpdom_data, 'glamv', glamv, lrowattr=ln_use_jattr ) 629 CALL iom_get( inum, jpdom_data, 'glamf', glamf, lrowattr=ln_use_jattr ) 630 631 CALL iom_get( inum, jpdom_data, 'gphit', gphit, lrowattr=ln_use_jattr ) 632 CALL iom_get( inum, jpdom_data, 'gphiu', gphiu, lrowattr=ln_use_jattr ) 633 CALL iom_get( inum, jpdom_data, 'gphiv', gphiv, lrowattr=ln_use_jattr ) 634 CALL iom_get( inum, jpdom_data, 'gphif', gphif, lrowattr=ln_use_jattr ) 635 636 CALL iom_get( inum, jpdom_data, 'e1t', e1t, lrowattr=ln_use_jattr ) 637 CALL iom_get( inum, jpdom_data, 'e1u', e1u, lrowattr=ln_use_jattr ) 638 CALL iom_get( inum, jpdom_data, 'e1v', e1v, lrowattr=ln_use_jattr ) 639 CALL iom_get( inum, jpdom_data, 'e1f', e1f, lrowattr=ln_use_jattr ) 640 641 CALL iom_get( inum, jpdom_data, 'e2t', e2t, lrowattr=ln_use_jattr ) 642 CALL iom_get( inum, jpdom_data, 'e2u', e2u, lrowattr=ln_use_jattr ) 643 CALL iom_get( inum, jpdom_data, 'e2v', e2v, lrowattr=ln_use_jattr ) 644 CALL iom_get( inum, jpdom_data, 'e2f', e2f, lrowattr=ln_use_jattr ) 637 645 638 646 CALL iom_close( inum ) 639 647 648 ! need to be define for the extended grid south of -80S 649 ! some point are undefined but you need to have e1 and e2 .NE. 0 650 WHERE (e1t==0.0_wp) 651 e1t=1.0e2 652 END WHERE 653 WHERE (e1v==0.0_wp) 654 e1v=1.0e2 655 END WHERE 656 WHERE (e1u==0.0_wp) 657 e1u=1.0e2 658 END WHERE 659 WHERE (e1f==0.0_wp) 660 e1f=1.0e2 661 END WHERE 662 WHERE (e2t==0.0_wp) 663 e2t=1.0e2 664 END WHERE 665 WHERE (e2v==0.0_wp) 666 e2v=1.0e2 667 END WHERE 668 WHERE (e2u==0.0_wp) 669 e2u=1.0e2 670 END WHERE 671 WHERE (e2f==0.0_wp) 672 e2f=1.0e2 673 END WHERE 674 640 675 END SUBROUTINE hgr_read 641 676 -
branches/2014/dev_r4650_UKMO10_Tidally_Meaned_Diagnostics/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90
r4624 r5260 184 184 END DO 185 185 END DO 186 187 ! (ISF) define barotropic mask and mask the ice shelf point 188 ssmask(:,:)=tmask(:,:,1) ! at this stage ice shelf is not masked 189 190 DO jk = 1, jpk 191 DO jj = 1, jpj 192 DO ji = 1, jpi 193 IF( REAL( misfdep(ji,jj) - jk, wp ) - 0.1_wp >= 0._wp ) THEN 194 tmask(ji,jj,jk) = 0._wp 195 END IF 196 END DO 197 END DO 198 END DO 186 199 187 200 !!gm ???? … … 207 220 ! Interior domain mask (used for global sum) 208 221 ! -------------------- 209 tmask_i(:,:) = tmask(:,:,1)222 tmask_i(:,:) = ssmask(:,:) ! (ISH) tmask_i = 1 even on the ice shelf 210 223 iif = jpreci ! ??? 211 224 iil = nlci - jpreci + 1 … … 250 263 END DO 251 264 END DO 265 ! (ISF) MIN(1,SUM(umask)) is here to check if you have effectively at least 1 wet u point 266 DO jj = 1, jpjm1 267 DO ji = 1, fs_jpim1 ! vector loop 268 umask_i(ji,jj) = ssmask(ji,jj) * ssmask(ji+1,jj ) * MIN(1._wp,SUM(umask(ji,jj,:))) 269 vmask_i(ji,jj) = ssmask(ji,jj) * ssmask(ji ,jj+1) * MIN(1._wp,SUM(vmask(ji,jj,:))) 270 END DO 271 DO ji = 1, jpim1 ! NO vector opt. 272 fmask_i(ji,jj) = ssmask(ji,jj ) * ssmask(ji+1,jj ) & 273 & * ssmask(ji,jj+1) * ssmask(ji+1,jj+1) * MIN(1._wp,SUM(fmask(ji,jj,:))) 274 END DO 275 END DO 252 276 CALL lbc_lnk( umask, 'U', 1._wp ) ! Lateral boundary conditions 253 277 CALL lbc_lnk( vmask, 'V', 1._wp ) 254 278 CALL lbc_lnk( fmask, 'F', 1._wp ) 255 279 CALL lbc_lnk( umask_i, 'U', 1._wp ) ! Lateral boundary conditions 280 CALL lbc_lnk( vmask_i, 'V', 1._wp ) 281 CALL lbc_lnk( fmask_i, 'F', 1._wp ) 282 283 ! 3. Ocean/land mask at wu-, wv- and w points 284 !---------------------------------------------- 285 wmask (:,:,1) = tmask(:,:,1) ! ???????? 286 wumask(:,:,1) = umask(:,:,1) ! ???????? 287 wvmask(:,:,1) = vmask(:,:,1) ! ???????? 288 DO jk=2,jpk 289 wmask (:,:,jk)=tmask(:,:,jk) * tmask(:,:,jk-1) 290 wumask(:,:,jk)=umask(:,:,jk) * umask(:,:,jk-1) 291 wvmask(:,:,jk)=vmask(:,:,jk) * vmask(:,:,jk-1) 292 END DO 256 293 257 294 ! 4. ocean/land mask for the elliptic equation 258 295 ! -------------------------------------------- 259 bmask(:,:) = tmask(:,:,1) ! elliptic equation is written at t-point296 bmask(:,:) = ssmask(:,:) ! elliptic equation is written at t-point 260 297 ! 261 298 ! ! Boundary conditions -
branches/2014/dev_r4650_UKMO10_Tidally_Meaned_Diagnostics/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90
r4624 r5260 8 8 !! 3.3 ! 2011-10 (M. Leclair) totally rewrote domvvl: 9 9 !! vvl option includes z_star and z_tilde coordinates 10 !! 3.6 ! 2014-11 (P. Mathiot) add ice shelf capability 10 11 !!---------------------------------------------------------------------- 11 12 !! 'key_vvl' variable volume … … 44 45 45 46 !!* Namelist nam_vvl 46 LOGICAL , PUBLIC :: ln_vvl_zstar ! zstar vertical coordinate47 LOGICAL , PUBLIC :: ln_vvl_ztilde ! ztilde vertical coordinate48 LOGICAL , PUBLIC :: ln_vvl_layer ! level vertical coordinate49 LOGICAL , PUBLIC :: ln_vvl_ztilde_as_zstar ! ztilde vertical coordinate50 LOGICAL , PUBLIC :: ln_vvl_zstar_at_eqtor ! ztilde vertical coordinate51 LOGICAL , PUBLIC :: ln_vvl_kepe ! kinetic/potential energy transfer52 ! ! conservation: not used yet47 LOGICAL , PUBLIC :: ln_vvl_zstar = .FALSE. ! zstar vertical coordinate 48 LOGICAL , PUBLIC :: ln_vvl_ztilde = .FALSE. ! ztilde vertical coordinate 49 LOGICAL , PUBLIC :: ln_vvl_layer = .FALSE. ! level vertical coordinate 50 LOGICAL , PUBLIC :: ln_vvl_ztilde_as_zstar = .FALSE. ! ztilde vertical coordinate 51 LOGICAL , PUBLIC :: ln_vvl_zstar_at_eqtor = .FALSE. ! ztilde vertical coordinate 52 LOGICAL , PUBLIC :: ln_vvl_kepe = .FALSE. ! kinetic/potential energy transfer 53 ! ! conservation: not used yet 53 54 REAL(wp) :: rn_ahe3 ! thickness diffusion coefficient 54 55 REAL(wp) :: rn_rst_e3t ! ztilde to zstar restoration timescale [days] 55 56 REAL(wp) :: rn_lf_cutoff ! cutoff frequency for low-pass filter [days] 56 57 REAL(wp) :: rn_zdef_max ! maximum fractional e3t deformation 57 LOGICAL , PUBLIC :: ln_vvl_dbg 58 LOGICAL , PUBLIC :: ln_vvl_dbg = .FALSE. ! debug control prints 58 59 59 60 !! * Module variables … … 125 126 INTEGER :: ji,jj,jk 126 127 INTEGER :: ii0, ii1, ij0, ij1 128 REAL(wp):: zcoef 127 129 !!---------------------------------------------------------------------- 128 130 IF( nn_timing == 1 ) CALL timing_start('dom_vvl_init') … … 164 166 ! t- and w- points depth 165 167 ! ---------------------- 168 ! set the isf depth as it is in the initial step 166 169 fsdept_n(:,:,1) = 0.5_wp * fse3w_n(:,:,1) 167 170 fsdepw_n(:,:,1) = 0.0_wp … … 169 172 fsdept_b(:,:,1) = 0.5_wp * fse3w_b(:,:,1) 170 173 fsdepw_b(:,:,1) = 0.0_wp 174 171 175 DO jk = 2, jpk 172 fsdept_n(:,:,jk) = fsdept_n(:,:,jk-1) + fse3w_n(:,:,jk) 173 fsdepw_n(:,:,jk) = fsdepw_n(:,:,jk-1) + fse3t_n(:,:,jk-1) 174 fsde3w_n(:,:,jk) = fsdept_n(:,:,jk ) - sshn (:,:) 175 fsdept_b(:,:,jk) = fsdept_b(:,:,jk-1) + fse3w_b(:,:,jk) 176 fsdepw_b(:,:,jk) = fsdepw_b(:,:,jk-1) + fse3t_b(:,:,jk-1) 176 DO jj = 1,jpj 177 DO ji = 1,jpi 178 ! zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt 179 ! 1 everywhere from mbkt to mikt + 1 or 1 (if no isf) 180 ! 0.5 where jk = mikt 181 zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) 182 fsdepw_n(ji,jj,jk) = fsdepw_n(ji,jj,jk-1) + fse3t_n(ji,jj,jk-1) 183 fsdept_n(ji,jj,jk) = zcoef * ( fsdepw_n(ji,jj,jk ) + 0.5 * fse3w_n(ji,jj,jk)) & 184 & + (1-zcoef) * ( fsdept_n(ji,jj,jk-1) + fse3w_n(ji,jj,jk)) 185 fsde3w_n(ji,jj,jk) = fsdept_n(ji,jj,jk) - sshn(ji,jj) 186 fsdepw_b(ji,jj,jk) = fsdepw_b(ji,jj,jk-1) + fse3t_b(ji,jj,jk-1) 187 fsdept_b(ji,jj,jk) = zcoef * ( fsdepw_b(ji,jj,jk ) + 0.5 * fse3w_b(ji,jj,jk)) & 188 & + (1-zcoef) * ( fsdept_b(ji,jj,jk-1) + fse3w_b(ji,jj,jk)) 189 END DO 190 END DO 177 191 END DO 178 192 … … 185 199 hv_b(:,:) = hv_b(:,:) + fse3v_b(:,:,jk) * vmask(:,:,jk) 186 200 END DO 187 hur_b(:,:) = umask (:,:,1) / ( hu_b(:,:) + 1. - umask(:,:,1) )188 hvr_b(:,:) = vmask (:,:,1) / ( hv_b(:,:) + 1. - vmask(:,:,1) )201 hur_b(:,:) = umask_i(:,:) / ( hu_b(:,:) + 1. - umask_i(:,:) ) 202 hvr_b(:,:) = vmask_i(:,:) / ( hv_b(:,:) + 1. - vmask_i(:,:) ) 189 203 190 204 ! Restoring frequencies for z_tilde coordinate … … 293 307 ! ! --------------------------------------------- ! 294 308 295 z_scale(:,:) = ( ssha(:,:) - sshb(:,:) ) * tmask(:,:,1) / ( ht_0(:,:) + sshn(:,:) + 1. - tmask(:,:,1) )309 z_scale(:,:) = ( ssha(:,:) - sshb(:,:) ) * ssmask(:,:) / ( ht_0(:,:) + sshn(:,:) + 1. - ssmask(:,:) ) 296 310 DO jk = 1, jpkm1 297 311 ! formally this is the same as fse3t_a = e3t_0*(1+ssha/ht_0) … … 313 327 zht (:,:) = zht (:,:) + fse3t_n(:,:,jk) * tmask(:,:,jk) 314 328 END DO 315 zhdiv(:,:) = zhdiv(:,:) / ( zht(:,:) + 1. - tmask (:,:,1) )329 zhdiv(:,:) = zhdiv(:,:) / ( zht(:,:) + 1. - tmask_i(:,:) ) 316 330 317 331 ! 2 - Low frequency baroclinic horizontal divergence (z-tilde case only) … … 338 352 ELSE ! layer case 339 353 DO jk = 1, jpkm1 340 tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - fse3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) 354 tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - fse3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) * tmask(:,:,jk) 341 355 END DO 342 356 END IF … … 386 400 ! d - thickness diffusion transport: boundary conditions 387 401 ! (stored for tracer advction and continuity equation) 388 CALL lbc_lnk( un_td , 'U' , -1. )389 CALL lbc_lnk( vn_td , 'V' , -1. )402 CALL lbc_lnk( un_td , 'U' , -1._wp) 403 CALL lbc_lnk( vn_td , 'V' , -1._wp) 390 404 391 405 ! 4 - Time stepping of baroclinic scale factors … … 398 412 z2dt = 2.0_wp * rdt 399 413 ENDIF 400 CALL lbc_lnk( tilde_e3t_a(:,:,:), 'T', 1. )414 CALL lbc_lnk( tilde_e3t_a(:,:,:), 'T', 1._wp ) 401 415 tilde_e3t_a(:,:,:) = tilde_e3t_b(:,:,:) + z2dt * tmask(:,:,:) * tilde_e3t_a(:,:,:) 402 416 … … 453 467 zht(:,:) = zht(:,:) + tilde_e3t_a(:,:,jk) * tmask(:,:,jk) 454 468 END DO 455 z_scale(:,:) = - zht(:,:) / ( ht_0(:,:) + sshn(:,:) + 1. - tmask(:,:,1) )469 z_scale(:,:) = - zht(:,:) / ( ht_0(:,:) + sshn(:,:) + 1. - ssmask(:,:) ) 456 470 DO jk = 1, jpkm1 457 471 dtilde_e3t_a(:,:,jk) = dtilde_e3t_a(:,:,jk) + fse3t_n(:,:,jk) * z_scale(:,:) * tmask(:,:,jk) … … 463 477 ! ! ---baroclinic part--------- ! 464 478 DO jk = 1, jpkm1 465 fse3t_a(:,:,jk) = fse3t_a(:,:,jk) + dtilde_e3t_a(:,:,jk) 479 fse3t_a(:,:,jk) = fse3t_a(:,:,jk) + dtilde_e3t_a(:,:,jk) * tmask(:,:,jk) 466 480 END DO 467 481 ENDIF … … 531 545 END DO 532 546 ! ! Inverse of the local depth 533 hur_a(:,:) = 1._wp / ( hu_a(:,:) + 1._wp - umask (:,:,1) ) * umask(:,:,1)534 hvr_a(:,:) = 1._wp / ( hv_a(:,:) + 1._wp - vmask (:,:,1) ) * vmask(:,:,1)547 hur_a(:,:) = 1._wp / ( hu_a(:,:) + 1._wp - umask_i(:,:) ) * umask_i(:,:) 548 hvr_a(:,:) = 1._wp / ( hv_a(:,:) + 1._wp - vmask_i(:,:) ) * vmask_i(:,:) 535 549 536 550 CALL wrk_dealloc( jpi, jpj, zht, z_scale, zwu, zwv, zhdiv ) … … 569 583 INTEGER, INTENT( in ) :: kt ! time step 570 584 !! * Local declarations 571 REAL(wp), POINTER, DIMENSION(:,:,:) :: z_e3t_def572 INTEGER :: jk ! dummy loop indices585 INTEGER :: ji,jj,jk ! dummy loop indices 586 REAL(wp) :: zcoef 573 587 !!---------------------------------------------------------------------- 574 588 575 589 IF( nn_timing == 1 ) CALL timing_start('dom_vvl_sf_swp') 576 !577 CALL wrk_alloc( jpi, jpj, jpk, z_e3t_def )578 590 ! 579 591 IF( kt == nit000 ) THEN … … 619 631 ! t- and w- points depth 620 632 ! ---------------------- 633 ! set the isf depth as it is in the initial step 621 634 fsdept_n(:,:,1) = 0.5_wp * fse3w_n(:,:,1) 622 635 fsdepw_n(:,:,1) = 0.0_wp 623 636 fsde3w_n(:,:,1) = fsdept_n(:,:,1) - sshn(:,:) 637 624 638 DO jk = 2, jpk 625 fsdept_n(:,:,jk) = fsdept_n(:,:,jk-1) + fse3w_n(:,:,jk) 626 fsdepw_n(:,:,jk) = fsdepw_n(:,:,jk-1) + fse3t_n(:,:,jk-1) 627 fsde3w_n(:,:,jk) = fsdept_n(:,:,jk ) - sshn (:,:) 639 DO jj = 1,jpj 640 DO ji = 1,jpi 641 ! zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt 642 ! 1 for jk = mikt 643 zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) 644 fsdepw_n(ji,jj,jk) = fsdepw_n(ji,jj,jk-1) + fse3t_n(ji,jj,jk-1) 645 fsdept_n(ji,jj,jk) = zcoef * ( fsdepw_n(ji,jj,jk ) + 0.5 * fse3w_n(ji,jj,jk)) & 646 & + (1-zcoef) * ( fsdept_n(ji,jj,jk-1) + fse3w_n(ji,jj,jk)) 647 fsde3w_n(ji,jj,jk) = fsdept_n(ji,jj,jk) - sshn(ji,jj) 648 END DO 649 END DO 628 650 END DO 651 629 652 ! Local depth and Inverse of the local depth of the water column at u- and v- points 630 653 ! ---------------------------------------------------------------------------------- … … 645 668 ! Write outputs 646 669 ! ============= 647 z_e3t_def(:,:,:) = ( ( fse3t_n(:,:,:) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2 648 CALL iom_put( "cellthc" , fse3t_n (:,:,:) ) 670 CALL iom_put( "e3t" , fse3t_n (:,:,:) ) 671 CALL iom_put( "e3u" , fse3u_n (:,:,:) ) 672 CALL iom_put( "e3v" , fse3v_n (:,:,:) ) 673 CALL iom_put( "e3w" , fse3w_n (:,:,:) ) 649 674 CALL iom_put( "tpt_dep" , fsde3w_n (:,:,:) ) 650 CALL iom_put( "e3tdef" , z_e3t_def(:,:,:) ) 675 IF( iom_use("e3tdef") ) & 676 CALL iom_put( "e3tdef" , ( ( fse3t_n(:,:,:) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2 ) 651 677 652 678 ! write restart file 653 679 ! ================== 654 680 IF( lrst_oce ) CALL dom_vvl_rst( kt, 'WRITE' ) 655 !656 CALL wrk_dealloc( jpi, jpj, jpk, z_e3t_def )657 681 ! 658 682 IF( nn_timing == 1 ) CALL timing_stop('dom_vvl_sf_swp') … … 702 726 IF( l_is_orca ) CALL dom_vvl_orca_fix( pe3_in, pe3_out, pout ) 703 727 ! boundary conditions 704 CALL lbc_lnk( pe3_out(:,:,:), 'U', 1. )728 CALL lbc_lnk( pe3_out(:,:,:), 'U', 1._wp ) 705 729 pe3_out(:,:,:) = pe3_out(:,:,:) + e3u_0(:,:,:) 706 730 ! ! ------------------------------------- ! … … 720 744 IF( l_is_orca ) CALL dom_vvl_orca_fix( pe3_in, pe3_out, pout ) 721 745 ! boundary conditions 722 CALL lbc_lnk( pe3_out(:,:,:), 'V', 1. )746 CALL lbc_lnk( pe3_out(:,:,:), 'V', 1._wp ) 723 747 pe3_out(:,:,:) = pe3_out(:,:,:) + e3v_0(:,:,:) 724 748 ! ! ------------------------------------- ! … … 738 762 IF( l_is_orca ) CALL dom_vvl_orca_fix( pe3_in, pe3_out, pout ) 739 763 ! boundary conditions 740 CALL lbc_lnk( pe3_out(:,:,:), 'F', 1. )764 CALL lbc_lnk( pe3_out(:,:,:), 'F', 1._wp ) 741 765 pe3_out(:,:,:) = pe3_out(:,:,:) + e3f_0(:,:,:) 742 766 ! ! ------------------------------------- ! … … 808 832 id3 = iom_varid( numror, 'tilde_e3t_b', ldstop = .FALSE. ) 809 833 id4 = iom_varid( numror, 'tilde_e3t_n', ldstop = .FALSE. ) 810 id5 = iom_varid( numror, 'hdi f_lf', ldstop = .FALSE. )834 id5 = iom_varid( numror, 'hdiv_lf', ldstop = .FALSE. ) 811 835 ! ! --------- ! 812 836 ! ! all cases ! … … 815 839 CALL iom_get( numror, jpdom_autoglo, 'fse3t_b', fse3t_b(:,:,:) ) 816 840 CALL iom_get( numror, jpdom_autoglo, 'fse3t_n', fse3t_n(:,:,:) ) 841 ! needed to restart if land processor not computed 842 IF(lwp) write(numout,*) 'dom_vvl_rst : fse3t_b and fse3t_n found in restart files' 843 WHERE ( tmask(:,:,:) == 0.0_wp ) 844 fse3t_n(:,:,:) = e3t_0(:,:,:) 845 fse3t_b(:,:,:) = e3t_0(:,:,:) 846 END WHERE 817 847 IF( neuler == 0 ) THEN 818 848 fse3t_b(:,:,:) = fse3t_n(:,:,:) 819 849 ENDIF 820 850 ELSE IF( id1 > 0 ) THEN 851 IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : fse3t_n not found in restart files' 852 IF(lwp) write(numout,*) 'fse3t_n set equal to fse3t_b.' 853 IF(lwp) write(numout,*) 'neuler is forced to 0' 854 CALL iom_get( numror, jpdom_autoglo, 'fse3t_b', fse3t_b(:,:,:) ) 855 fse3t_n(:,:,:) = fse3t_b(:,:,:) 856 neuler = 0 857 ELSE IF( id2 > 0 ) THEN 821 858 IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : fse3t_b not found in restart files' 822 859 IF(lwp) write(numout,*) 'fse3t_b set equal to fse3t_n.' 823 860 IF(lwp) write(numout,*) 'neuler is forced to 0' 861 CALL iom_get( numror, jpdom_autoglo, 'fse3t_n', fse3t_n(:,:,:) ) 824 862 fse3t_b(:,:,:) = fse3t_n(:,:,:) 825 863 neuler = 0 … … 830 868 DO jk=1,jpk 831 869 fse3t_n(:,:,jk) = e3t_0(:,:,jk) * ( ht_0(:,:) + sshn(:,:) ) & 832 & / ( ht_0(:,:) + 1._wp - tmask(:,:,1) ) * tmask(:,:,jk) &870 & / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk) & 833 871 & + e3t_0(:,:,jk) * (1._wp -tmask(:,:,jk)) 834 872 END DO … … 964 1002 965 1003 IF( ioptio /= 1 ) CALL ctl_stop( 'Choose ONE vertical coordinate in namelist nam_vvl' ) 1004 IF( .NOT. ln_vvl_zstar .AND. nn_isf .NE. 0) CALL ctl_stop( 'Only vvl_zstar has been tested with ice shelf cavity' ) 966 1005 967 1006 IF(lwp) THEN ! Print the choice -
branches/2014/dev_r4650_UKMO10_Tidally_Meaned_Diagnostics/NEMOGCM/NEMO/OPA_SRC/DOM/domwri.F90
r4292 r5260 132 132 133 133 CALL dom_uniq( zprw, 'T' ) 134 zprt = tmask(:,:,1) * zprw ! ! unique point mask 134 DO jj = 1, jpj 135 DO ji = 1, jpi 136 jk=mikt(ji,jj) 137 zprt(ji,jj) = tmask(ji,jj,jk) * zprw(ji,jj) ! ! unique point mask 138 END DO 139 END DO ! ! unique point mask 135 140 CALL iom_rstput( 0, 0, inum2, 'tmaskutil', zprt, ktype = jp_i1 ) 136 141 CALL dom_uniq( zprw, 'U' ) 137 zprt = umask(:,:,1) * zprw 142 DO jj = 1, jpj 143 DO ji = 1, jpi 144 jk=miku(ji,jj) 145 zprt(ji,jj) = umask(ji,jj,jk) * zprw(ji,jj) ! ! unique point mask 146 END DO 147 END DO 138 148 CALL iom_rstput( 0, 0, inum2, 'umaskutil', zprt, ktype = jp_i1 ) 139 149 CALL dom_uniq( zprw, 'V' ) 140 zprt = vmask(:,:,1) * zprw 150 DO jj = 1, jpj 151 DO ji = 1, jpi 152 jk=mikv(ji,jj) 153 zprt(ji,jj) = vmask(ji,jj,jk) * zprw(ji,jj) ! ! unique point mask 154 END DO 155 END DO 141 156 CALL iom_rstput( 0, 0, inum2, 'vmaskutil', zprt, ktype = jp_i1 ) 142 157 CALL dom_uniq( zprw, 'F' ) 143 zprt = fmask(:,:,1) * zprw 158 DO jj = 1, jpj 159 DO ji = 1, jpi 160 jk=mikf(ji,jj) 161 zprt(ji,jj) = fmask(ji,jj,jk) * zprw(ji,jj) ! ! unique point mask 162 END DO 163 END DO 144 164 CALL iom_rstput( 0, 0, inum2, 'fmaskutil', zprt, ktype = jp_i1 ) 145 165 … … 168 188 169 189 ! note that mbkt is set to 1 over land ==> use surface tmask 170 zprt(:,:) = tmask(:,:,1) * REAL( mbkt(:,:) , wp )190 zprt(:,:) = ssmask(:,:) * REAL( mbkt(:,:) , wp ) 171 191 CALL iom_rstput( 0, 0, inum4, 'mbathy', zprt, ktype = jp_i2 ) ! ! nb of ocean T-points 192 zprt(:,:) = ssmask(:,:) * REAL( mikt(:,:) , wp ) 193 CALL iom_rstput( 0, 0, inum4, 'misf', zprt, ktype = jp_i2 ) ! ! nb of ocean T-points 194 zprt(:,:) = ssmask(:,:) * REAL( risfdep(:,:) , wp ) 195 CALL iom_rstput( 0, 0, inum4, 'isfdraft', zprt, ktype = jp_r4 ) ! ! nb of ocean T-points 172 196 173 197 IF( ln_sco ) THEN ! s-coordinate … … 203 227 DO jj = 1,jpj 204 228 DO ji = 1,jpi 205 e3tp(ji,jj) = e3t_0(ji,jj,mbkt(ji,jj)) * tmask(ji,jj,1)206 e3wp(ji,jj) = e3w_0(ji,jj,mbkt(ji,jj)) * tmask(ji,jj,1)229 e3tp(ji,jj) = e3t_0(ji,jj,mbkt(ji,jj)) * ssmask(ji,jj) 230 e3wp(ji,jj) = e3w_0(ji,jj,mbkt(ji,jj)) * ssmask(ji,jj) 207 231 END DO 208 232 END DO … … 228 252 DO jj = 1,jpj 229 253 DO ji = 1,jpi 230 zprt(ji,jj) = gdept_0(ji,jj,mbkt(ji,jj) ) * tmask(ji,jj,1)231 zprw(ji,jj) = gdepw_0(ji,jj,mbkt(ji,jj)+1) * tmask(ji,jj,1)254 zprt(ji,jj) = gdept_0(ji,jj,mbkt(ji,jj) ) * ssmask(ji,jj) 255 zprw(ji,jj) = gdepw_0(ji,jj,mbkt(ji,jj)+1) * ssmask(ji,jj) 232 256 END DO 233 257 END DO -
branches/2014/dev_r4650_UKMO10_Tidally_Meaned_Diagnostics/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90
r4624 r5260 17 17 !! 3.4 ! 2012-08 (J. Siddorn) added Siddorn and Furner stretching function 18 18 !! 3.4 ! 2012-12 (R. Bourdalle-Badie and G. Reffray) modify C1D case 19 !! 3.6 ! 2014-11 (P. Mathiot and C. Harris) add ice shelf capabilitye 19 20 !!---------------------------------------------------------------------- 20 21 … … 101 102 INTEGER :: ios 102 103 ! 103 NAMELIST/namzgr/ ln_zco, ln_zps, ln_sco 104 NAMELIST/namzgr/ ln_zco, ln_zps, ln_sco, ln_isfcav 104 105 !!---------------------------------------------------------------------- 105 106 ! … … 120 121 WRITE(numout,*) '~~~~~~~' 121 122 WRITE(numout,*) ' Namelist namzgr : set vertical coordinate' 122 WRITE(numout,*) ' z-coordinate - full steps ln_zco = ', ln_zco 123 WRITE(numout,*) ' z-coordinate - partial steps ln_zps = ', ln_zps 124 WRITE(numout,*) ' s- or hybrid z-s-coordinate ln_sco = ', ln_sco 123 WRITE(numout,*) ' z-coordinate - full steps ln_zco = ', ln_zco 124 WRITE(numout,*) ' z-coordinate - partial steps ln_zps = ', ln_zps 125 WRITE(numout,*) ' s- or hybrid z-s-coordinate ln_sco = ', ln_sco 126 WRITE(numout,*) ' ice shelf cavities ln_isfcav = ', ln_isfcav 125 127 ENDIF 126 128 … … 145 147 IF( .NOT.lk_c1d ) CALL zgr_bat_ctl ! check bathymetry (mbathy) and suppress isolated ocean points 146 148 CALL zgr_bot_level ! deepest ocean level for t-, u- and v-points 149 CALL zgr_top_level ! shallowest ocean level for T-, U-, V- points 147 150 ! 148 151 IF( lk_c1d ) THEN ! 1D config.: same mbathy value over the 3x3 domain … … 295 298 ENDIF 296 299 300 IF ( ln_isfcav ) THEN 301 ! need to be like this to compute the pressure gradient with ISF. If not, level beneath the ISF are not aligned (sum(e3t) /= depth) 302 ! define e3t_0 and e3w_0 as the differences between gdept and gdepw respectively 303 DO jk = 1, jpkm1 304 e3t_1d(jk) = gdepw_1d(jk+1)-gdepw_1d(jk) 305 END DO 306 e3t_1d(jpk) = e3t_1d(jpk-1) ! we don't care because this level is masked in NEMO 307 308 DO jk = 2, jpk 309 e3w_1d(jk) = gdept_1d(jk) - gdept_1d(jk-1) 310 END DO 311 e3w_1d(1 ) = 2._wp * (gdept_1d(1) - gdepw_1d(1)) 312 END IF 313 297 314 !!gm BUG in s-coordinate this does not work! 298 315 ! deepest/shallowest W level Above/Below ~10m … … 350 367 INTEGER :: ji, jj, jl, jk ! dummy loop indices 351 368 INTEGER :: inum ! temporary logical unit 369 INTEGER :: ierror ! error flag 352 370 INTEGER :: ii_bump, ij_bump, ih ! bump center position 353 371 INTEGER :: ii0, ii1, ij0, ij1, ik ! local indices 354 372 REAL(wp) :: r_bump , h_bump , h_oce ! bump characteristics 355 373 REAL(wp) :: zi, zj, zh, zhmin ! local scalars 356 INTEGER , POINTER, DIMENSION(:,:) :: idta ! global domain integer data357 REAL(wp), POINTER, DIMENSION(:,:) :: zdta ! global domain scalar data374 INTEGER , ALLOCATABLE, DIMENSION(:,:) :: idta ! global domain integer data 375 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdta ! global domain scalar data 358 376 !!---------------------------------------------------------------------- 359 377 ! 360 378 IF( nn_timing == 1 ) CALL timing_start('zgr_bat') 361 !362 CALL wrk_alloc( jpidta, jpjdta, idta )363 CALL wrk_alloc( jpidta, jpjdta, zdta )364 379 ! 365 380 IF(lwp) WRITE(numout,*) … … 370 385 ! ! ================== ! 371 386 ! ! global domain level and meter bathymetry (idta,zdta) 387 ! 388 ALLOCATE( idta(jpidta,jpjdta), STAT=ierror ) 389 IF( ierror > 0 ) CALL ctl_stop( 'STOP', 'zgr_bat: unable to allocate idta array' ) 390 ALLOCATE( zdta(jpidta,jpjdta), STAT=ierror ) 391 IF( ierror > 0 ) CALL ctl_stop( 'STOP', 'zgr_bat: unable to allocate zdta array' ) 372 392 ! 373 393 IF( ntopo == 0 ) THEN ! flat basin … … 450 470 END DO 451 471 END DO 472 risfdep(:,:)=0.e0 473 misfdep(:,:)=1 474 ! 475 ! (ISF) TODO build ice draft netcdf file for isomip and build the corresponding part of code 476 IF( cp_cfg == "isomip" .AND. ln_isfcav ) THEN 477 risfdep(:,:)=200.e0 478 misfdep(:,:)=1 479 ij0 = 1 ; ij1 = 40 480 DO jj = mj0(ij0), mj1(ij1) 481 risfdep(:,jj)=700.0_wp-(gphit(:,jj)+80.0_wp)*125.0_wp 482 END DO 483 WHERE( bathy(:,:) <= 0._wp ) risfdep(:,:) = 0._wp 484 ! 485 ELSEIF ( cp_cfg == "isomip2" .AND. ln_isfcav ) THEN 486 ! 487 risfdep(:,:)=0.e0 488 misfdep(:,:)=1 489 ij0 = 1 ; ij1 = 40 490 DO jj = mj0(ij0), mj1(ij1) 491 risfdep(:,jj)=700.0_wp-(gphit(:,jj)+80.0_wp)*125.0_wp 492 END DO 493 WHERE( bathy(:,:) <= 0._wp ) risfdep(:,:) = 0._wp 494 END IF 495 ! 496 DEALLOCATE( idta, zdta ) 452 497 ! 453 498 ! ! ================ ! … … 490 535 IF( ln_zps .OR. ln_sco ) THEN ! zps or sco : read meter bathymetry 491 536 CALL iom_open ( 'bathy_meter.nc', inum ) 492 CALL iom_get ( inum, jpdom_data, 'Bathymetry', bathy ) 537 IF ( ln_isfcav ) THEN 538 CALL iom_get ( inum, jpdom_data, 'Bathymetry_isf', bathy, lrowattr=.false. ) 539 ELSE 540 CALL iom_get ( inum, jpdom_data, 'Bathymetry' , bathy, lrowattr=ln_use_jattr ) 541 END IF 493 542 CALL iom_close( inum ) 494 543 ! 544 risfdep(:,:)=0._wp 545 misfdep(:,:)=1 546 IF ( ln_isfcav ) THEN 547 CALL iom_open ( 'isf_draft_meter.nc', inum ) 548 CALL iom_get ( inum, jpdom_data, 'isf_draft', risfdep ) 549 CALL iom_close( inum ) 550 WHERE( bathy(:,:) <= 0._wp ) risfdep(:,:) = 0._wp 551 END IF 552 ! 495 553 IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN ! ORCA R2 configuration 496 554 ! … … 530 588 ! 531 589 IF ( .not. ln_sco ) THEN !== set a minimum depth ==! 590 ! patch to avoid case bathy = ice shelf draft and bathy between 0 and zhmin 591 IF ( ln_isfcav ) THEN 592 WHERE (bathy == risfdep) 593 bathy = 0.0_wp ; risfdep = 0.0_wp 594 END WHERE 595 END IF 596 ! end patch 532 597 IF( rn_hmin < 0._wp ) THEN ; ik = - INT( rn_hmin ) ! from a nb of level 533 598 ELSE ; ik = MINLOC( gdepw_1d, mask = gdepw_1d > rn_hmin, dim = 1 ) ! from a depth … … 539 604 IF(lwp) write(numout,*) 'Minimum ocean depth: ', zhmin, ' minimum number of ocean levels : ', ik 540 605 ENDIF 541 !542 CALL wrk_dealloc( jpidta, jpjdta, idta )543 CALL wrk_dealloc( jpidta, jpjdta, zdta )544 606 ! 545 607 IF( nn_timing == 1 ) CALL timing_stop('zgr_bat') … … 783 845 END SUBROUTINE zgr_bot_level 784 846 847 SUBROUTINE zgr_top_level 848 !!---------------------------------------------------------------------- 849 !! *** ROUTINE zgr_bot_level *** 850 !! 851 !! ** Purpose : defines the vertical index of ocean top (mik. arrays) 852 !! 853 !! ** Method : computes from misfdep with a minimum value of 1 854 !! 855 !! ** Action : mikt, miku, mikv : vertical indices of the shallowest 856 !! ocean level at t-, u- & v-points 857 !! (min value = 1) 858 !!---------------------------------------------------------------------- 859 !! 860 INTEGER :: ji, jj ! dummy loop indices 861 REAL(wp), POINTER, DIMENSION(:,:) :: zmik 862 !!---------------------------------------------------------------------- 863 ! 864 IF( nn_timing == 1 ) CALL timing_start('zgr_top_level') 865 ! 866 CALL wrk_alloc( jpi, jpj, zmik ) 867 ! 868 IF(lwp) WRITE(numout,*) 869 IF(lwp) WRITE(numout,*) ' zgr_top_level : ocean top k-index of T-, U-, V- and W-levels ' 870 IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~~~' 871 ! 872 mikt(:,:) = MAX( misfdep(:,:) , 1 ) ! top k-index of T-level (=1) 873 ! ! top k-index of W-level (=mikt) 874 DO jj = 1, jpjm1 ! top k-index of U- (U-) level 875 DO ji = 1, jpim1 876 miku(ji,jj) = MAX( mikt(ji+1,jj ) , mikt(ji,jj) ) 877 mikv(ji,jj) = MAX( mikt(ji ,jj+1) , mikt(ji,jj) ) 878 mikf(ji,jj) = MAX( mikt(ji ,jj+1) , mikt(ji,jj), mikt(ji+1,jj ), mikt(ji+1,jj+1) ) 879 END DO 880 END DO 881 882 ! converte into REAL to use lbc_lnk ; impose a min value of 1 as a zero can be set in lbclnk 883 zmik(:,:) = REAL( miku(:,:), wp ) ; CALL lbc_lnk(zmik,'U',1.) ; miku (:,:) = MAX( INT( zmik(:,:) ), 1 ) 884 zmik(:,:) = REAL( mikv(:,:), wp ) ; CALL lbc_lnk(zmik,'V',1.) ; mikv (:,:) = MAX( INT( zmik(:,:) ), 1 ) 885 zmik(:,:) = REAL( mikf(:,:), wp ) ; CALL lbc_lnk(zmik,'F',1.) ; mikf (:,:) = MAX( INT( zmik(:,:) ), 1 ) 886 ! 887 CALL wrk_dealloc( jpi, jpj, zmik ) 888 ! 889 IF( nn_timing == 1 ) CALL timing_stop('zgr_top_level') 890 ! 891 END SUBROUTINE zgr_top_level 785 892 786 893 SUBROUTINE zgr_zco … … 906 1013 WHERE( 0._wp < bathy(:,:) .AND. bathy(:,:) <= zdepth ) mbathy(:,:) = jk-1 907 1014 END DO 1015 1016 IF ( ln_isfcav ) CALL zgr_isf 908 1017 909 1018 ! Scale factors and depth at T- and W-points … … 938 1047 !gm Bug? check the gdepw_1d 939 1048 ! ... on ik 940 gdept_0(ji,jj,ik) = gdepw_1d(ik) + ( gdepw_0 941 & * ((gdept_1d(ik ) - gdepw_1d(ik) ) &942 & / ( gdepw_1d(ik+1) - gdepw_1d(ik) ))943 e3t_0 (ji,jj,ik) = e3t_1d(ik) * ( gdepw_0 (ji,jj,ik+1) - gdepw_1d(ik) ) &944 & / ( gdepw_1d( ik+1) - gdepw_1d(ik) )1049 gdept_0(ji,jj,ik) = gdepw_1d(ik) + ( gdepw_0(ji,jj,ik+1) - gdepw_1d(ik) ) & 1050 & * ((gdept_1d( ik ) - gdepw_1d(ik) ) & 1051 & / ( gdepw_1d( ik+1) - gdepw_1d(ik) )) 1052 e3t_0 (ji,jj,ik) = e3t_1d (ik) * ( gdepw_0 (ji,jj,ik+1) - gdepw_1d(ik) ) & 1053 & / ( gdepw_1d( ik+1) - gdepw_1d(ik) ) 945 1054 e3w_0(ji,jj,ik) = 0.5_wp * ( gdepw_0(ji,jj,ik+1) + gdepw_1d(ik+1) - 2._wp * gdepw_1d(ik) ) & 946 1055 & * ( e3w_1d(ik) / ( gdepw_1d(ik+1) - gdepw_1d(ik) ) ) … … 973 1082 END DO 974 1083 END DO 1084 ! 1085 IF ( ln_isfcav ) THEN 1086 ! (ISF) Definition of e3t, u, v, w for ISF case 1087 DO jj = 1, jpj 1088 DO ji = 1, jpi 1089 ik = misfdep(ji,jj) 1090 IF( ik > 1 ) THEN ! ice shelf point only 1091 IF( risfdep(ji,jj) < gdepw_1d(ik) ) risfdep(ji,jj)= gdepw_1d(ik) 1092 gdepw_0(ji,jj,ik) = risfdep(ji,jj) 1093 !gm Bug? check the gdepw_0 1094 ! ... on ik 1095 gdept_0(ji,jj,ik) = gdepw_1d(ik+1) - ( gdepw_1d(ik+1) - gdepw_0(ji,jj,ik) ) & 1096 & * ( gdepw_1d(ik+1) - gdept_1d(ik) ) & 1097 & / ( gdepw_1d(ik+1) - gdepw_1d(ik) ) 1098 e3t_0 (ji,jj,ik ) = gdepw_1d(ik+1) - gdepw_0(ji,jj,ik) 1099 e3w_0 (ji,jj,ik+1) = gdept_1d(ik+1) - gdept_0(ji,jj,ik) 1100 1101 IF( ik + 1 == mbathy(ji,jj) ) THEN ! ice shelf point only (2 cell water column) 1102 e3w_0 (ji,jj,ik+1) = gdept_0(ji,jj,ik+1) - gdept_0(ji,jj,ik) 1103 ENDIF 1104 ! ... on ik / ik-1 1105 e3w_0 (ji,jj,ik ) = 2._wp * (gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik)) 1106 e3t_0 (ji,jj,ik-1) = gdepw_0(ji,jj,ik) - gdepw_1d(ik-1) 1107 ! The next line isn't required and doesn't affect results - included for consistency with bathymetry code 1108 gdept_0(ji,jj,ik-1) = gdept_1d(ik-1) 1109 ENDIF 1110 END DO 1111 END DO 1112 ! 1113 it = 0 1114 DO jj = 1, jpj 1115 DO ji = 1, jpi 1116 ik = misfdep(ji,jj) 1117 IF( ik > 1 ) THEN ! ice shelf point only 1118 e3tp (ji,jj) = e3t_0(ji,jj,ik ) 1119 e3wp (ji,jj) = e3w_0(ji,jj,ik+1 ) 1120 ! test 1121 zdiff= gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik ) 1122 IF( zdiff <= 0. .AND. lwp ) THEN 1123 it = it + 1 1124 WRITE(numout,*) ' it = ', it, ' ik = ', ik, ' (i,j) = ', ji, jj 1125 WRITE(numout,*) ' risfdep = ', risfdep(ji,jj) 1126 WRITE(numout,*) ' gdept = ', gdept_0(ji,jj,ik), ' gdepw = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff 1127 WRITE(numout,*) ' e3tp = ', e3tp(ji,jj), ' e3wp = ', e3wp(ji,jj) 1128 ENDIF 1129 ENDIF 1130 END DO 1131 END DO 1132 END IF 1133 ! END (ISF) 975 1134 976 1135 ! Scale factors and depth at U-, V-, UW and VW-points … … 991 1150 END DO 992 1151 END DO 1152 IF ( ln_isfcav ) THEN 1153 ! (ISF) define e3uw (adapted for 2 cells in the water column) 1154 ! Need to test if the modification of only mikt and mbkt levels is enough 1155 DO jk = 2,jpk 1156 DO jj = 1, jpjm1 1157 DO ji = 1, fs_jpim1 ! vector opt. 1158 e3uw_0(ji,jj,jk) = MIN( gdept_0(ji,jj,jk), gdept_0(ji+1,jj ,jk) ) & 1159 & - MAX( gdept_0(ji,jj,jk-1), gdept_0(ji+1,jj ,jk-1) ) 1160 e3vw_0(ji,jj,jk) = MIN( gdept_0(ji,jj,jk), gdept_0(ji ,jj+1,jk) ) & 1161 & - MAX( gdept_0(ji,jj,jk-1), gdept_0(ji ,jj+1,jk-1) ) 1162 END DO 1163 END DO 1164 END DO 1165 END IF 1166 993 1167 CALL lbc_lnk( e3u_0 , 'U', 1._wp ) ; CALL lbc_lnk( e3uw_0, 'U', 1._wp ) ! lateral boundary conditions 994 1168 CALL lbc_lnk( e3v_0 , 'V', 1._wp ) ; CALL lbc_lnk( e3vw_0, 'V', 1._wp ) … … 1032 1206 1033 1207 ! Compute gdep3w_0 (vertical sum of e3w) 1034 gdep3w_0(:,:,1) = 0.5_wp * e3w_0(:,:,1) 1035 DO jk = 2, jpk 1036 gdep3w_0(:,:,jk) = gdep3w_0(:,:,jk-1) + e3w_0(:,:,jk) 1037 END DO 1038 1208 IF ( ln_isfcav ) THEN ! if cavity 1209 WHERE (misfdep == 0) misfdep = 1 1210 DO jj = 1,jpj 1211 DO ji = 1,jpi 1212 gdep3w_0(ji,jj,1) = 0.5_wp * e3w_0(ji,jj,1) 1213 DO jk = 2, misfdep(ji,jj) 1214 gdep3w_0(ji,jj,jk) = gdep3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk) 1215 END DO 1216 IF (misfdep(ji,jj) .GE. 2) gdep3w_0(ji,jj,misfdep(ji,jj)) = risfdep(ji,jj) + 0.5_wp * e3w_0(ji,jj,misfdep(ji,jj)) 1217 DO jk = misfdep(ji,jj) + 1, jpk 1218 gdep3w_0(ji,jj,jk) = gdep3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk) 1219 END DO 1220 END DO 1221 END DO 1222 ELSE ! no cavity 1223 gdep3w_0(:,:,1) = 0.5_wp * e3w_0(:,:,1) 1224 DO jk = 2, jpk 1225 gdep3w_0(:,:,jk) = gdep3w_0(:,:,jk-1) + e3w_0(:,:,jk) 1226 END DO 1227 END IF 1039 1228 ! ! ================= ! 1040 1229 IF(lwp .AND. ll_print) THEN ! Control print ! … … 1070 1259 ! 1071 1260 END SUBROUTINE zgr_zps 1261 1262 SUBROUTINE zgr_isf 1263 !!---------------------------------------------------------------------- 1264 !! *** ROUTINE zgr_isf *** 1265 !! 1266 !! ** Purpose : check the bathymetry in levels 1267 !! 1268 !! ** Method : THe water column have to contained at least 2 cells 1269 !! Bathymetry and isfdraft are modified (dig/close) to respect 1270 !! this criterion. 1271 !! 1272 !! 1273 !! ** Action : - test compatibility between isfdraft and bathy 1274 !! - bathy and isfdraft are modified 1275 !!---------------------------------------------------------------------- 1276 !! 1277 INTEGER :: ji, jj, jk, jl ! dummy loop indices 1278 INTEGER :: ik, it ! temporary integers 1279 INTEGER :: id, jd, nprocd 1280 INTEGER :: icompt, ibtest, ibtestim1, ibtestip1, ibtestjm1, ibtestjp1 ! (ISF) 1281 LOGICAL :: ll_print ! Allow control print for debugging 1282 REAL(wp) :: ze3tp , ze3wp ! Last ocean level thickness at T- and W-points 1283 REAL(wp) :: zdepwp, zdepth ! Ajusted ocean depth to avoid too small e3t 1284 REAL(wp) :: zmax, zmin ! Maximum and minimum depth 1285 REAL(wp) :: zdiff ! temporary scalar 1286 REAL(wp) :: zrefdep ! temporary scalar 1287 REAL(wp) :: zbathydiff, zrisfdepdiff ! isf temporary scalar 1288 REAL(wp), POINTER, DIMENSION(:,:) :: zrisfdep, zbathy, zmask ! 2D workspace (ISH) 1289 INTEGER , POINTER, DIMENSION(:,:) :: zmbathy, zmisfdep ! 2D workspace (ISH) 1290 !!--------------------------------------------------------------------- 1291 ! 1292 IF( nn_timing == 1 ) CALL timing_start('zgr_isf') 1293 ! 1294 CALL wrk_alloc( jpi, jpj, zbathy, zmask, zrisfdep) 1295 CALL wrk_alloc( jpi, jpj, zmisfdep, zmbathy ) 1296 1297 1298 ! (ISF) compute misfdep 1299 WHERE( risfdep(:,:) == 0._wp .AND. bathy(:,:) .NE. 0) ; misfdep(:,:) = 1 ! open water : set misfdep to 1 1300 ELSEWHERE ; misfdep(:,:) = 2 ! iceshelf : initialize misfdep to second level 1301 END WHERE 1302 1303 ! Compute misfdep for ocean points (i.e. first wet level) 1304 ! find the first ocean level such that the first level thickness 1305 ! is larger than the bot_level of e3zps_min and e3zps_rat * e3t_0 (where 1306 ! e3t_0 is the reference level thickness 1307 DO jk = 2, jpkm1 1308 zdepth = gdepw_1d(jk+1) - MIN( e3zps_min, e3t_1d(jk)*e3zps_rat ) 1309 WHERE( 0._wp < risfdep(:,:) .AND. risfdep(:,:) >= zdepth ) misfdep(:,:) = jk+1 1310 END DO 1311 WHERE (risfdep(:,:) <= e3t_1d(1) .AND. risfdep(:,:) .GT. 0._wp) 1312 risfdep(:,:) = 0. ; misfdep(:,:) = 1 1313 END WHERE 1314 1315 ! basic check for the compatibility of bathy and risfdep. I think it should be offline because it is not perfect and cannot solved all the situation 1316 icompt = 0 1317 ! run the bathy check 10 times to be sure all the modif in the bathy or iceshelf draft are compatible together 1318 DO jl = 1, 10 1319 WHERE (bathy(:,:) .EQ. risfdep(:,:) ) 1320 misfdep(:,:) = 0 ; risfdep(:,:) = 0._wp 1321 mbathy (:,:) = 0 ; bathy (:,:) = 0._wp 1322 END WHERE 1323 WHERE (mbathy(:,:) <= 0) 1324 misfdep(:,:) = 0; risfdep(:,:) = 0._wp 1325 mbathy (:,:) = 0; bathy (:,:) = 0._wp 1326 ENDWHERE 1327 IF( lk_mpp ) THEN 1328 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1329 CALL lbc_lnk( zbathy, 'T', 1. ) 1330 misfdep(:,:) = INT( zbathy(:,:) ) 1331 CALL lbc_lnk( risfdep, 'T', 1. ) 1332 CALL lbc_lnk( bathy, 'T', 1. ) 1333 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1334 CALL lbc_lnk( zbathy, 'T', 1. ) 1335 mbathy(:,:) = INT( zbathy(:,:) ) 1336 ENDIF 1337 IF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 ) THEN 1338 misfdep( 1 ,:) = misfdep(jpim1,:) ! local domain is cyclic east-west 1339 misfdep(jpi,:) = misfdep( 2 ,:) 1340 ENDIF 1341 1342 IF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 ) THEN 1343 mbathy( 1 ,:) = mbathy(jpim1,:) ! local domain is cyclic east-west 1344 mbathy(jpi,:) = mbathy( 2 ,:) 1345 ENDIF 1346 1347 ! split last cell if possible (only where water column is 2 cell or less) 1348 DO jk = jpkm1, 1, -1 1349 zmax = gdepw_1d(jk) + MIN( e3zps_min, e3t_1d(jk)*e3zps_rat ) 1350 WHERE( gdepw_1d(jk) < bathy(:,:) .AND. bathy(:,:) <= zmax .AND. misfdep + 1 >= mbathy) 1351 mbathy(:,:) = jk 1352 bathy(:,:) = zmax 1353 END WHERE 1354 END DO 1355 1356 ! split top cell if possible (only where water column is 2 cell or less) 1357 DO jk = 2, jpkm1 1358 zmax = gdepw_1d(jk+1) - MIN( e3zps_min, e3t_1d(jk)*e3zps_rat ) 1359 WHERE( gdepw_1d(jk+1) > risfdep(:,:) .AND. risfdep(:,:) >= zmax .AND. misfdep + 1 >= mbathy) 1360 misfdep(:,:) = jk 1361 risfdep(:,:) = zmax 1362 END WHERE 1363 END DO 1364 1365 1366 ! Case where bathy and risfdep compatible but not the level variable mbathy/misfdep because of partial cell condition 1367 DO jj = 1, jpj 1368 DO ji = 1, jpi 1369 ! find the minimum change option: 1370 ! test bathy 1371 IF (risfdep(ji,jj) .GT. 1) THEN 1372 zbathydiff =ABS(bathy(ji,jj) - (gdepw_1d(mbathy (ji,jj)+1) & 1373 & + MIN( e3zps_min, e3t_1d(mbathy (ji,jj)+1)*e3zps_rat ))) 1374 zrisfdepdiff=ABS(risfdep(ji,jj) - (gdepw_1d(misfdep(ji,jj) ) & 1375 & - MIN( e3zps_min, e3t_1d(misfdep(ji,jj)-1)*e3zps_rat ))) 1376 1377 IF (bathy(ji,jj) .GT. risfdep(ji,jj) .AND. mbathy(ji,jj) .LT. misfdep(ji,jj)) THEN 1378 IF (zbathydiff .LE. zrisfdepdiff) THEN 1379 bathy(ji,jj) = gdepw_1d(mbathy(ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj)+1)*e3zps_rat ) 1380 mbathy(ji,jj)= mbathy(ji,jj) + 1 1381 ELSE 1382 risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj)) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj)-1)*e3zps_rat ) 1383 misfdep(ji,jj) = misfdep(ji,jj) - 1 1384 END IF 1385 END IF 1386 END IF 1387 END DO 1388 END DO 1389 1390 ! At least 2 levels for water thickness at T, U, and V point. 1391 DO jj = 1, jpj 1392 DO ji = 1, jpi 1393 ! find the minimum change option: 1394 ! test bathy 1395 IF( misfdep(ji,jj) == mbathy(ji,jj) .AND. mbathy(ji,jj) .GT. 1) THEN 1396 zbathydiff =ABS(bathy(ji,jj) - (gdepw_1d(mbathy (ji,jj)+1)& 1397 & + MIN( e3zps_min,e3t_1d(mbathy (ji,jj)+1)*e3zps_rat ))) 1398 zrisfdepdiff=ABS(risfdep(ji,jj) - (gdepw_1d(misfdep(ji,jj) ) & 1399 & - MIN( e3zps_min,e3t_1d(misfdep(ji,jj)-1)*e3zps_rat ))) 1400 IF (zbathydiff .LE. zrisfdepdiff) THEN 1401 mbathy(ji,jj) = mbathy(ji,jj) + 1 1402 bathy(ji,jj) = gdepw_1d(mbathy (ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj) +1)*e3zps_rat ) 1403 ELSE 1404 misfdep(ji,jj)= misfdep(ji,jj) - 1 1405 risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj))*e3zps_rat ) 1406 END IF 1407 ENDIF 1408 END DO 1409 END DO 1410 1411 ! point V mbathy(ji,jj) EQ misfdep(ji,jj+1) 1412 DO jj = 1, jpjm1 1413 DO ji = 1, jpim1 1414 IF( misfdep(ji,jj+1) == mbathy(ji,jj) .AND. mbathy(ji,jj) .GT. 1) THEN 1415 zbathydiff =ABS(bathy(ji,jj ) - (gdepw_1d(mbathy (ji,jj)+1) & 1416 & + MIN( e3zps_min, e3t_1d(mbathy (ji,jj )+1)*e3zps_rat ))) 1417 zrisfdepdiff=ABS(risfdep(ji,jj+1) - (gdepw_1d(misfdep(ji,jj+1)) & 1418 & - MIN( e3zps_min, e3t_1d(misfdep(ji,jj+1)-1)*e3zps_rat ))) 1419 IF (zbathydiff .LE. zrisfdepdiff) THEN 1420 mbathy(ji,jj) = mbathy(ji,jj) + 1 1421 bathy(ji,jj) = gdepw_1d(mbathy (ji,jj )) & 1422 & + MIN( e3zps_min, e3t_1d(mbathy(ji,jj )+1)*e3zps_rat ) 1423 ELSE 1424 misfdep(ji,jj+1) = misfdep(ji,jj+1) - 1 1425 risfdep (ji,jj+1) = gdepw_1d(misfdep(ji,jj+1)+1) & 1426 & - MIN( e3zps_min, e3t_1d(misfdep(ji,jj+1))*e3zps_rat ) 1427 END IF 1428 ENDIF 1429 END DO 1430 END DO 1431 1432 IF( lk_mpp ) THEN 1433 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1434 CALL lbc_lnk( zbathy, 'T', 1. ) 1435 misfdep(:,:) = INT( zbathy(:,:) ) 1436 CALL lbc_lnk( risfdep, 'T', 1. ) 1437 CALL lbc_lnk( bathy, 'T', 1. ) 1438 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1439 CALL lbc_lnk( zbathy, 'T', 1. ) 1440 mbathy(:,:) = INT( zbathy(:,:) ) 1441 ENDIF 1442 ! point V misdep(ji,jj) EQ mbathy(ji,jj+1) 1443 DO jj = 1, jpjm1 1444 DO ji = 1, jpim1 1445 IF( misfdep(ji,jj) == mbathy(ji,jj+1) .AND. mbathy(ji,jj) .GT. 1) THEN 1446 zbathydiff =ABS( bathy(ji,jj+1) - (gdepw_1d(mbathy (ji,jj+1)+1) & 1447 & + MIN( e3zps_min, e3t_1d(mbathy (ji,jj+1)+1)*e3zps_rat ))) 1448 zrisfdepdiff=ABS(risfdep(ji,jj ) - (gdepw_1d(misfdep(ji,jj ) ) & 1449 & - MIN( e3zps_min, e3t_1d(misfdep(ji,jj )-1)*e3zps_rat ))) 1450 IF (zbathydiff .LE. zrisfdepdiff) THEN 1451 mbathy (ji,jj+1) = mbathy(ji,jj+1) + 1 1452 bathy (ji,jj+1) = gdepw_1d(mbathy (ji,jj+1) ) + MIN( e3zps_min, e3t_1d(mbathy (ji,jj+1)+1)*e3zps_rat ) 1453 ELSE 1454 misfdep(ji,jj) = misfdep(ji,jj) - 1 1455 risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj )+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj ) )*e3zps_rat ) 1456 END IF 1457 ENDIF 1458 END DO 1459 END DO 1460 1461 1462 IF( lk_mpp ) THEN 1463 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1464 CALL lbc_lnk( zbathy, 'T', 1. ) 1465 misfdep(:,:) = INT( zbathy(:,:) ) 1466 CALL lbc_lnk( risfdep, 'T', 1. ) 1467 CALL lbc_lnk( bathy, 'T', 1. ) 1468 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1469 CALL lbc_lnk( zbathy, 'T', 1. ) 1470 mbathy(:,:) = INT( zbathy(:,:) ) 1471 ENDIF 1472 1473 ! point U mbathy(ji,jj) EQ misfdep(ji,jj+1) 1474 DO jj = 1, jpjm1 1475 DO ji = 1, jpim1 1476 IF( misfdep(ji+1,jj) == mbathy(ji,jj) .AND. mbathy(ji,jj) .GT. 1) THEN 1477 zbathydiff =ABS( bathy(ji ,jj) - (gdepw_1d(mbathy (ji,jj)+1) & 1478 & + MIN( e3zps_min, e3t_1d(mbathy (ji ,jj)+1)*e3zps_rat ))) 1479 zrisfdepdiff=ABS(risfdep(ji+1,jj) - (gdepw_1d(misfdep(ji+1,jj)) & 1480 & - MIN( e3zps_min, e3t_1d(misfdep(ji+1,jj)-1)*e3zps_rat ))) 1481 IF (zbathydiff .LE. zrisfdepdiff) THEN 1482 mbathy(ji,jj) = mbathy(ji,jj) + 1 1483 bathy(ji,jj) = gdepw_1d(mbathy (ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj) +1)*e3zps_rat ) 1484 ELSE 1485 misfdep(ji+1,jj)= misfdep(ji+1,jj) - 1 1486 risfdep(ji+1,jj) = gdepw_1d(misfdep(ji+1,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji+1,jj))*e3zps_rat ) 1487 END IF 1488 ENDIF 1489 ENDDO 1490 ENDDO 1491 1492 IF( lk_mpp ) THEN 1493 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1494 CALL lbc_lnk( zbathy, 'T', 1. ) 1495 misfdep(:,:) = INT( zbathy(:,:) ) 1496 CALL lbc_lnk( risfdep, 'T', 1. ) 1497 CALL lbc_lnk( bathy, 'T', 1. ) 1498 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1499 CALL lbc_lnk( zbathy, 'T', 1. ) 1500 mbathy(:,:) = INT( zbathy(:,:) ) 1501 ENDIF 1502 1503 ! point U misfdep(ji,jj) EQ bathy(ji,jj+1) 1504 DO jj = 1, jpjm1 1505 DO ji = 1, jpim1 1506 IF( misfdep(ji,jj) == mbathy(ji+1,jj) .AND. mbathy(ji,jj) .GT. 1) THEN 1507 zbathydiff =ABS( bathy(ji+1,jj) - (gdepw_1d(mbathy (ji+1,jj)+1) & 1508 & + MIN( e3zps_min, e3t_1d(mbathy (ji+1,jj)+1)*e3zps_rat ))) 1509 zrisfdepdiff=ABS(risfdep(ji ,jj) - (gdepw_1d(misfdep(ji ,jj) ) & 1510 & - MIN( e3zps_min, e3t_1d(misfdep(ji ,jj)-1)*e3zps_rat ))) 1511 IF (zbathydiff .LE. zrisfdepdiff) THEN 1512 mbathy(ji+1,jj) = mbathy (ji+1,jj) + 1 1513 bathy (ji+1,jj) = gdepw_1d(mbathy (ji+1,jj) ) & 1514 & + MIN( e3zps_min, e3t_1d(mbathy (ji+1,jj) +1)*e3zps_rat ) 1515 ELSE 1516 misfdep(ji,jj) = misfdep(ji ,jj) - 1 1517 risfdep(ji,jj) = gdepw_1d(misfdep(ji ,jj)+1) & 1518 & - MIN( e3zps_min, e3t_1d(misfdep(ji ,jj) )*e3zps_rat ) 1519 END IF 1520 ENDIF 1521 ENDDO 1522 ENDDO 1523 1524 IF( lk_mpp ) THEN 1525 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1526 CALL lbc_lnk( zbathy, 'T', 1. ) 1527 misfdep(:,:) = INT( zbathy(:,:) ) 1528 CALL lbc_lnk( risfdep, 'T', 1. ) 1529 CALL lbc_lnk( bathy, 'T', 1. ) 1530 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1531 CALL lbc_lnk( zbathy, 'T', 1. ) 1532 mbathy(:,:) = INT( zbathy(:,:) ) 1533 ENDIF 1534 END DO 1535 ! end dig bathy/ice shelf to be compatible 1536 ! now fill single point in "coastline" of ice shelf, bathy, hole, and test again one cell tickness 1537 DO jl = 1,20 1538 1539 ! remove single point "bay" on isf coast line in the ice shelf draft' 1540 DO jk = 1, jpk 1541 WHERE (misfdep==0) misfdep=jpk 1542 zmask=0 1543 WHERE (misfdep .LE. jk) zmask=1 1544 DO jj = 2, jpjm1 1545 DO ji = 2, jpim1 1546 IF (misfdep(ji,jj) .EQ. jk) THEN 1547 ibtest = zmask(ji-1,jj) + zmask(ji+1,jj) + zmask(ji,jj-1) + zmask(ji,jj+1) 1548 IF (ibtest .LE. 1) THEN 1549 risfdep(ji,jj)=gdepw_1d(jk+1) ; misfdep(ji,jj)=jk+1 1550 IF (misfdep(ji,jj) .GT. mbathy(ji,jj)) misfdep(ji,jj) = jpk 1551 END IF 1552 END IF 1553 END DO 1554 END DO 1555 END DO 1556 WHERE (misfdep==jpk) 1557 misfdep=0 ; risfdep=0. ; mbathy=0 ; bathy=0. 1558 END WHERE 1559 IF( lk_mpp ) THEN 1560 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1561 CALL lbc_lnk( zbathy, 'T', 1. ) 1562 misfdep(:,:) = INT( zbathy(:,:) ) 1563 CALL lbc_lnk( risfdep, 'T', 1. ) 1564 CALL lbc_lnk( bathy, 'T', 1. ) 1565 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1566 CALL lbc_lnk( zbathy, 'T', 1. ) 1567 mbathy(:,:) = INT( zbathy(:,:) ) 1568 ENDIF 1569 1570 ! remove single point "bay" on bathy coast line beneath an ice shelf' 1571 DO jk = jpk,1,-1 1572 zmask=0 1573 WHERE (mbathy .GE. jk ) zmask=1 1574 DO jj = 2, jpjm1 1575 DO ji = 2, jpim1 1576 IF (mbathy(ji,jj) .EQ. jk .AND. misfdep(ji,jj) .GE. 2) THEN 1577 ibtest = zmask(ji-1,jj) + zmask(ji+1,jj) + zmask(ji,jj-1) + zmask(ji,jj+1) 1578 IF (ibtest .LE. 1) THEN 1579 bathy(ji,jj)=gdepw_1d(jk) ; mbathy(ji,jj)=jk-1 1580 IF (misfdep(ji,jj) .GT. mbathy(ji,jj)) mbathy(ji,jj) = 0 1581 END IF 1582 END IF 1583 END DO 1584 END DO 1585 END DO 1586 WHERE (mbathy==0) 1587 misfdep=0 ; risfdep=0. ; mbathy=0 ; bathy=0. 1588 END WHERE 1589 IF( lk_mpp ) THEN 1590 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1591 CALL lbc_lnk( zbathy, 'T', 1. ) 1592 misfdep(:,:) = INT( zbathy(:,:) ) 1593 CALL lbc_lnk( risfdep, 'T', 1. ) 1594 CALL lbc_lnk( bathy, 'T', 1. ) 1595 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1596 CALL lbc_lnk( zbathy, 'T', 1. ) 1597 mbathy(:,:) = INT( zbathy(:,:) ) 1598 ENDIF 1599 1600 ! fill hole in ice shelf 1601 zmisfdep = misfdep 1602 zrisfdep = risfdep 1603 WHERE (zmisfdep .LE. 1) zmisfdep=jpk 1604 DO jj = 2, jpjm1 1605 DO ji = 2, jpim1 1606 ibtestim1 = zmisfdep(ji-1,jj ) ; ibtestip1 = zmisfdep(ji+1,jj ) 1607 ibtestjm1 = zmisfdep(ji ,jj-1) ; ibtestjp1 = zmisfdep(ji ,jj+1) 1608 IF( zmisfdep(ji,jj) .GE. mbathy(ji-1,jj ) ) ibtestim1 = jpk!MAX(0, mbathy(ji-1,jj ) - 1) 1609 IF( zmisfdep(ji,jj) .GE. mbathy(ji+1,jj ) ) ibtestip1 = jpk!MAX(0, mbathy(ji+1,jj ) - 1) 1610 IF( zmisfdep(ji,jj) .GE. mbathy(ji ,jj-1) ) ibtestjm1 = jpk!MAX(0, mbathy(ji ,jj-1) - 1) 1611 IF( zmisfdep(ji,jj) .GE. mbathy(ji ,jj+1) ) ibtestjp1 = jpk!MAX(0, mbathy(ji ,jj+1) - 1) 1612 ibtest=MIN(ibtestim1, ibtestip1, ibtestjm1, ibtestjp1) 1613 IF( ibtest == jpk .AND. misfdep(ji,jj) .GE. 2) THEN 1614 mbathy(ji,jj) = 0 ; bathy(ji,jj) = 0.0_wp ; misfdep(ji,jj) = 0 ; risfdep(ji,jj) = 0.0_wp 1615 END IF 1616 IF( zmisfdep(ji,jj) < ibtest .AND. misfdep(ji,jj) .GE. 2) THEN 1617 misfdep(ji,jj) = ibtest 1618 risfdep(ji,jj) = gdepw_1d(ibtest) 1619 ENDIF 1620 ENDDO 1621 ENDDO 1622 1623 IF( lk_mpp ) THEN 1624 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1625 CALL lbc_lnk( zbathy, 'T', 1. ) 1626 misfdep(:,:) = INT( zbathy(:,:) ) 1627 CALL lbc_lnk( risfdep, 'T', 1. ) 1628 CALL lbc_lnk( bathy, 'T', 1. ) 1629 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1630 CALL lbc_lnk( zbathy, 'T', 1. ) 1631 mbathy(:,:) = INT( zbathy(:,:) ) 1632 ENDIF 1633 ! 1634 !! fill hole in bathymetry 1635 zmbathy (:,:)=mbathy (:,:) 1636 DO jj = 2, jpjm1 1637 DO ji = 2, jpim1 1638 ibtestim1 = zmbathy(ji-1,jj ) ; ibtestip1 = zmbathy(ji+1,jj ) 1639 ibtestjm1 = zmbathy(ji ,jj-1) ; ibtestjp1 = zmbathy(ji ,jj+1) 1640 IF( zmbathy(ji,jj) .LT. misfdep(ji-1,jj ) ) ibtestim1 = 0!MIN(jpk-1, misfdep(ji-1,jj ) + 1) 1641 IF( zmbathy(ji,jj) .LT. misfdep(ji+1,jj ) ) ibtestip1 = 0 1642 IF( zmbathy(ji,jj) .LT. misfdep(ji ,jj-1) ) ibtestjm1 = 0 1643 IF( zmbathy(ji,jj) .LT. misfdep(ji ,jj+1) ) ibtestjp1 = 0 1644 ibtest=MAX(ibtestim1, ibtestip1, ibtestjm1, ibtestjp1) 1645 IF( ibtest == 0 .AND. misfdep(ji,jj) .GE. 2) THEN 1646 mbathy(ji,jj) = 0 ; bathy(ji,jj) = 0.0_wp ; misfdep(ji,jj) = 0 ; risfdep(ji,jj) = 0.0_wp ; 1647 END IF 1648 IF( ibtest < zmbathy(ji,jj) .AND. misfdep(ji,jj) .GE. 2) THEN 1649 mbathy(ji,jj) = ibtest 1650 bathy(ji,jj) = gdepw_1d(ibtest+1) 1651 ENDIF 1652 END DO 1653 END DO 1654 IF( lk_mpp ) THEN 1655 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1656 CALL lbc_lnk( zbathy, 'T', 1. ) 1657 misfdep(:,:) = INT( zbathy(:,:) ) 1658 CALL lbc_lnk( risfdep, 'T', 1. ) 1659 CALL lbc_lnk( bathy, 'T', 1. ) 1660 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1661 CALL lbc_lnk( zbathy, 'T', 1. ) 1662 mbathy(:,:) = INT( zbathy(:,:) ) 1663 ENDIF 1664 ! if not compatible after all check (ie U point water column less than 2 cells), mask U 1665 DO jj = 1, jpjm1 1666 DO ji = 1, jpim1 1667 IF (mbathy(ji,jj) == misfdep(ji+1,jj) .AND. mbathy(ji,jj) .GE. 1 .AND. mbathy(ji+1,jj) .GE. 1) THEN 1668 mbathy(ji,jj) = mbathy(ji,jj) - 1 ; bathy(ji,jj) = gdepw_1d(mbathy(ji,jj)+1) ; 1669 END IF 1670 END DO 1671 END DO 1672 IF( lk_mpp ) THEN 1673 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1674 CALL lbc_lnk( zbathy, 'T', 1. ) 1675 misfdep(:,:) = INT( zbathy(:,:) ) 1676 CALL lbc_lnk( risfdep, 'T', 1. ) 1677 CALL lbc_lnk( bathy, 'T', 1. ) 1678 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1679 CALL lbc_lnk( zbathy, 'T', 1. ) 1680 mbathy(:,:) = INT( zbathy(:,:) ) 1681 ENDIF 1682 ! if not compatible after all check (ie U point water column less than 2 cells), mask U 1683 DO jj = 1, jpjm1 1684 DO ji = 1, jpim1 1685 IF (misfdep(ji,jj) == mbathy(ji+1,jj) .AND. mbathy(ji,jj) .GE. 1 .AND. mbathy(ji+1,jj) .GE. 1) THEN 1686 mbathy(ji+1,jj) = mbathy(ji+1,jj) - 1; bathy(ji+1,jj) = gdepw_1d(mbathy(ji+1,jj)+1) ; 1687 END IF 1688 END DO 1689 END DO 1690 IF( lk_mpp ) THEN 1691 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1692 CALL lbc_lnk( zbathy, 'T', 1. ) 1693 misfdep(:,:) = INT( zbathy(:,:) ) 1694 CALL lbc_lnk( risfdep, 'T', 1. ) 1695 CALL lbc_lnk( bathy, 'T', 1. ) 1696 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1697 CALL lbc_lnk( zbathy, 'T', 1. ) 1698 mbathy(:,:) = INT( zbathy(:,:) ) 1699 ENDIF 1700 ! if not compatible after all check (ie V point water column less than 2 cells), mask V 1701 DO jj = 1, jpjm1 1702 DO ji = 1, jpi 1703 IF (mbathy(ji,jj) == misfdep(ji,jj+1) .AND. mbathy(ji,jj) .GE. 1 .AND. mbathy(ji,jj+1) .GE. 1) THEN 1704 mbathy(ji,jj) = mbathy(ji,jj) - 1 ; bathy(ji,jj) = gdepw_1d(mbathy(ji,jj)+1) ; 1705 END IF 1706 END DO 1707 END DO 1708 IF( lk_mpp ) THEN 1709 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1710 CALL lbc_lnk( zbathy, 'T', 1. ) 1711 misfdep(:,:) = INT( zbathy(:,:) ) 1712 CALL lbc_lnk( risfdep, 'T', 1. ) 1713 CALL lbc_lnk( bathy, 'T', 1. ) 1714 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1715 CALL lbc_lnk( zbathy, 'T', 1. ) 1716 mbathy(:,:) = INT( zbathy(:,:) ) 1717 ENDIF 1718 ! if not compatible after all check (ie V point water column less than 2 cells), mask V 1719 DO jj = 1, jpjm1 1720 DO ji = 1, jpi 1721 IF (misfdep(ji,jj) == mbathy(ji,jj+1) .AND. mbathy(ji,jj) .GE. 1 .AND. mbathy(ji,jj+1) .GE. 1) THEN 1722 mbathy(ji,jj+1) = mbathy(ji,jj+1) - 1 ; bathy(ji,jj+1) = gdepw_1d(mbathy(ji,jj+1)+1) ; 1723 END IF 1724 END DO 1725 END DO 1726 IF( lk_mpp ) THEN 1727 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1728 CALL lbc_lnk( zbathy, 'T', 1. ) 1729 misfdep(:,:) = INT( zbathy(:,:) ) 1730 CALL lbc_lnk( risfdep, 'T', 1. ) 1731 CALL lbc_lnk( bathy, 'T', 1. ) 1732 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1733 CALL lbc_lnk( zbathy, 'T', 1. ) 1734 mbathy(:,:) = INT( zbathy(:,:) ) 1735 ENDIF 1736 ! if not compatible after all check, mask T 1737 DO jj = 1, jpj 1738 DO ji = 1, jpi 1739 IF (mbathy(ji,jj) <= misfdep(ji,jj)) THEN 1740 misfdep(ji,jj) = 0 ; risfdep(ji,jj) = 0._wp ; mbathy(ji,jj) = 0 ; bathy(ji,jj) = 0._wp ; 1741 END IF 1742 END DO 1743 END DO 1744 1745 WHERE (mbathy(:,:) == 1) 1746 mbathy = 0; bathy = 0.0_wp ; misfdep = 0 ; risfdep = 0.0_wp 1747 END WHERE 1748 END DO 1749 ! end check compatibility ice shelf/bathy 1750 ! remove very shallow ice shelf (less than ~ 10m if 75L) 1751 WHERE (misfdep(:,:) <= 5) 1752 misfdep = 1; risfdep = 0.0_wp; 1753 END WHERE 1754 1755 IF( icompt == 0 ) THEN 1756 IF(lwp) WRITE(numout,*)' no points with ice shelf too close to bathymetry' 1757 ELSE 1758 IF(lwp) WRITE(numout,*)' ',icompt,' ocean grid points with ice shelf thickness reduced to avoid bathymetry' 1759 ENDIF 1760 1761 CALL wrk_dealloc( jpi, jpj, zmask, zbathy, zrisfdep ) 1762 CALL wrk_dealloc( jpi, jpj, zmisfdep, zmbathy ) 1763 1764 IF( nn_timing == 1 ) CALL timing_stop('zgr_isf') 1765 1766 END SUBROUTINE 1072 1767 1073 1768 SUBROUTINE zgr_sco … … 1445 2140 DO jk = 1, jpkm1 1446 2141 IF( scobot(ji,jj) >= fsdept(ji,jj,jk) ) mbathy(ji,jj) = MAX( 2, jk ) 1447 IF( scobot(ji,jj) == 0._wp ) mbathy(ji,jj) = 01448 END DO2142 END DO 2143 IF( scobot(ji,jj) == 0._wp ) mbathy(ji,jj) = 0 1449 2144 END DO 1450 2145 END DO -
branches/2014/dev_r4650_UKMO10_Tidally_Meaned_Diagnostics/NEMOGCM/NEMO/OPA_SRC/DOM/dtatsd.F90
- Property svn:keywords set to Id
r4624 r5260 39 39 !!---------------------------------------------------------------------- 40 40 !! NEMO/OPA 3.3 , NEMO Consortium (2010) 41 !! $Id : dtatem.F90 2392 2010-11-15 21:20:05Z gm$41 !! $Id$ 42 42 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 43 43 !!---------------------------------------------------------------------- … … 263 263 ptsd(ji,jj,ik,jp_sal) = (1.-zl) * ptsd(ji,jj,ik,jp_sal) + zl * ptsd(ji,jj,ik-1,jp_sal) 264 264 ENDIF 265 ik = mikt(ji,jj) 266 IF( ik > 1 ) THEN 267 zl = ( gdept_0(ji,jj,ik) - gdept_1d(ik) ) / ( gdept_1d(ik+1) - gdept_1d(ik) ) 268 ptsd(ji,jj,ik,jp_tem) = (1.-zl) * ptsd(ji,jj,ik,jp_tem) + zl * ptsd(ji,jj,ik+1,jp_tem) 269 ptsd(ji,jj,ik,jp_sal) = (1.-zl) * ptsd(ji,jj,ik,jp_sal) + zl * ptsd(ji,jj,ik+1,jp_sal) 270 END IF 265 271 END DO 266 272 END DO -
branches/2014/dev_r4650_UKMO10_Tidally_Meaned_Diagnostics/NEMOGCM/NEMO/OPA_SRC/DOM/istate.F90
r4370 r5260 34 34 USE dtatsd ! data temperature and salinity (dta_tsd routine) 35 35 USE dtauvd ! data: U & V current (dta_uvd routine) 36 USE in_out_manager ! I/O manager37 USE iom ! I/O library38 36 USE zpshde ! partial step: hor. derivative (zps_hde routine) 39 37 USE eosbn2 ! equation of state (eos bn2 routine) … … 42 40 USE dynspg_flt ! filtered free surface 43 41 USE sol_oce ! ocean solver variables 42 ! 43 USE in_out_manager ! I/O manager 44 USE iom ! I/O library 44 45 USE lib_mpp ! MPP library 45 46 USE restart ! restart … … 56 57 # include "vectopt_loop_substitute.h90" 57 58 !!---------------------------------------------------------------------- 58 !! NEMO/OPA 3. 3 , NEMO Consortium (2010)59 !! NEMO/OPA 3.7 , NEMO Consortium (2014) 59 60 !! $Id$ 60 61 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) … … 68 69 !! ** Purpose : Initialization of the dynamics and tracer fields. 69 70 !!---------------------------------------------------------------------- 70 ! - ML - needed for initialization of e3t_b 71 INTEGER :: ji,jj,jk ! dummy loop indices 72 REAL(wp), POINTER, DIMENSION(:,:,:,:) :: zuvd ! U & V data workspace 73 !!---------------------------------------------------------------------- 74 ! 75 IF( nn_timing == 1 ) CALL timing_start('istate_init') 76 ! 77 78 IF(lwp) WRITE(numout,*) 71 INTEGER :: ji, jj, jk ! dummy loop indices 72 REAL(wp), POINTER, DIMENSION(:,:,:,:) :: zuvd ! U & V data workspace 73 !!---------------------------------------------------------------------- 74 ! 75 IF( nn_timing == 1 ) CALL timing_start('istate_init') 76 ! 77 78 IF(lwp) WRITE(numout,*) ' ' 79 79 IF(lwp) WRITE(numout,*) 'istate_ini : Initialization of the dynamics and tracers' 80 80 IF(lwp) WRITE(numout,*) '~~~~~~~~~~' … … 83 83 IF( lk_c1d ) CALL dta_uvd_init ! Initialization of U & V input data 84 84 85 rhd (:,:,: ) = 0. e086 r hop (:,:,: ) = 0.e087 rn2 (:,:,: ) = 0.e088 tsa (:,:,:,:) = 0.e085 rhd (:,:,: ) = 0._wp ; rhop (:,:,: ) = 0._wp ! set one for all to 0 at level jpk 86 rn2b (:,:,: ) = 0._wp ; rn2 (:,:,: ) = 0._wp ! set one for all to 0 at levels 1 and jpk 87 tsa (:,:,:,:) = 0._wp ! set one for all to 0 at level jpk 88 rab_b(:,:,:,:) = 0._wp ; rab_n(:,:,:,:) = 0._wp ! set one for all to 0 at level jpk 89 89 90 90 IF( ln_rstart ) THEN ! Restart from a file … … 110 110 ELSEIF( cp_cfg == 'gyre' ) THEN 111 111 CALL istate_gyre ! GYRE configuration : start from pre-defined T-S fields 112 ELSEIF( cp_cfg == 'isomip' .OR. cp_cfg == 'isomip2') THEN 113 IF(lwp) WRITE(numout,*) 'Initialization of T+S for ISOMIP domain' 114 tsn(:,:,:,jp_tem)=-1.9*tmask(:,:,:) ! ISOMIP configuration : start from constant T+S fields 115 tsn(:,:,:,jp_sal)=34.4*tmask(:,:,:) 116 tsb(:,:,:,:)=tsn(:,:,:,:) 112 117 ELSE ! Initial T-S, U-V fields read in files 113 118 IF ( ln_tsd_init ) THEN ! read 3D T and S data at nit000 … … 129 134 CALL eos( tsb, rhd, rhop, gdept_0(:,:,:) ) ! before potential and in situ densities 130 135 #if ! defined key_c1d 131 IF( ln_zps ) CALL zps_hde( nit000, jpts, tsb, gtsu, gtsv, & ! zps: before hor. gradient 132 & rhd, gru , grv ) ! of t,s,rd at ocean bottom 136 IF( ln_zps .AND. .NOT. ln_isfcav) & 137 & CALL zps_hde ( nit000, jpts, tsb, gtsu, gtsv, & ! Partial steps: before horizontal gradient 138 & rhd, gru , grv ) ! of t, s, rd at the last ocean level 139 IF( ln_zps .AND. ln_isfcav) & 140 & CALL zps_hde_isf( nit000, jpts, tsb, gtsu, gtsv, & ! Partial steps for top cell (ISF) 141 & rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv , & 142 & gtui, gtvi, grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi ) ! of t, s, rd at the last ocean level 133 143 #endif 134 144 ! … … 162 172 ! 163 173 DO jk = 1, jpkm1 164 #if defined key_vectopt_loop165 DO jj = 1, 1 !Vector opt. => forced unrolling166 DO ji = 1, jpij167 #else168 174 DO jj = 1, jpj 169 175 DO ji = 1, jpi 170 #endif171 176 un_b(ji,jj) = un_b(ji,jj) + fse3u_n(ji,jj,jk) * un(ji,jj,jk) * umask(ji,jj,jk) 172 177 vn_b(ji,jj) = vn_b(ji,jj) + fse3v_n(ji,jj,jk) * vn(ji,jj,jk) * vmask(ji,jj,jk) … … 185 190 ! 186 191 ! 187 IF( nn_timing == 1 ) CALL timing_stop('istate_init')192 IF( nn_timing == 1 ) CALL timing_stop('istate_init') 188 193 ! 189 194 END SUBROUTINE istate_init 195 190 196 191 197 SUBROUTINE istate_t_s … … 219 225 END SUBROUTINE istate_t_s 220 226 227 221 228 SUBROUTINE istate_eel 222 229 !!---------------------------------------------------------------------- … … 233 240 USE divcur ! hor. divergence & rel. vorticity (div_cur routine) 234 241 USE iom 235 242 ! 236 243 INTEGER :: inum ! temporary logical unit 237 244 INTEGER :: ji, jj, jk ! dummy loop indices … … 244 251 REAL(wp), DIMENSION(jpiglo,jpjglo) :: zssh ! initial ssh over the global domain 245 252 !!---------------------------------------------------------------------- 246 253 ! 247 254 SELECT CASE ( jp_cfg ) 248 255 ! ! ==================== … … 375 382 INTEGER, PARAMETER :: ntsinit = 0 ! (0/1) (analytical/input data files) T&S initialization 376 383 !!---------------------------------------------------------------------- 377 384 ! 378 385 SELECT CASE ( ntsinit) 379 386 ! 380 387 CASE ( 0 ) ! analytical T/S profil deduced from LEVITUS 381 388 IF(lwp) WRITE(numout,*) 382 389 IF(lwp) WRITE(numout,*) 'istate_gyre : initial analytical T and S profil deduced from LEVITUS ' 383 390 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 384 391 ! 385 392 DO jk = 1, jpk 386 393 DO jj = 1, jpj … … 407 414 END DO 408 415 END DO 409 416 ! 410 417 CASE ( 1 ) ! T/S data fields read in dta_tem.nc/data_sal.nc files 411 418 IF(lwp) WRITE(numout,*) … … 431 438 tsn(:,:,:,jp_sal) = tsn(:,:,:,jp_sal) * tmask(:,:,:) 432 439 tsb(:,:,:,jp_sal) = tsn(:,:,:,jp_sal) 433 440 ! 434 441 END SELECT 435 442 ! 436 443 IF(lwp) THEN 437 444 WRITE(numout,*) … … 440 447 WRITE(numout, "(10x, i4, 3f10.2)" ) ( jk, gdept_1d(jk), tsn(2,2,jk,jp_tem), tsn(2,2,jk,jp_sal), jk = 1, jpk ) 441 448 ENDIF 442 449 ! 443 450 END SUBROUTINE istate_gyre 451 444 452 445 453 SUBROUTINE istate_uvg … … 457 465 USE divcur ! hor. divergence & rel. vorticity (div_cur routine) 458 466 USE lbclnk ! ocean lateral boundary condition (or mpp link) 459 467 ! 460 468 INTEGER :: ji, jj, jk ! dummy loop indices 461 469 INTEGER :: indic ! ??? … … 567 575 !!===================================================================== 568 576 END MODULE istate 569 -
branches/2014/dev_r4650_UKMO10_Tidally_Meaned_Diagnostics/NEMOGCM/NEMO/OPA_SRC/DOM/phycst.F90
r3625 r5260 41 41 REAL(wp), PUBLIC :: rt0 = 273.15_wp !: freezing point of fresh water [Kelvin] 42 42 #if defined key_lim3 43 REAL(wp), PUBLIC :: rt0_snow = 273.1 6_wp !: melting point of snow [Kelvin]44 REAL(wp), PUBLIC :: rt0_ice = 273.1 6_wp !: melting point of ice [Kelvin]43 REAL(wp), PUBLIC :: rt0_snow = 273.15_wp !: melting point of snow [Kelvin] 44 REAL(wp), PUBLIC :: rt0_ice = 273.15_wp !: melting point of ice [Kelvin] 45 45 #else 46 46 REAL(wp), PUBLIC :: rt0_snow = 273.15_wp !: melting point of snow [Kelvin] 47 47 REAL(wp), PUBLIC :: rt0_ice = 273.05_wp !: melting point of ice [Kelvin] 48 48 #endif 49 #if defined key_cice 50 REAL(wp), PUBLIC :: rau0 = 1026._wp !: volumic mass of reference [kg/m3] 51 #else 52 REAL(wp), PUBLIC :: rau0 = 1035._wp !: volumic mass of reference [kg/m3] 53 #endif 49 REAL(wp), PUBLIC :: rau0 !: volumic mass of reference [kg/m3] 54 50 REAL(wp), PUBLIC :: r1_rau0 !: = 1. / rau0 [m3/kg] 55 REAL(wp), PUBLIC :: rauw = 1000._wp !: volumic mass of pure water [m3/kg] 56 REAL(wp), PUBLIC :: rcp = 4.e3_wp !: ocean specific heat [J/Kelvin] 51 REAL(wp), PUBLIC :: rcp !: ocean specific heat [J/Kelvin] 57 52 REAL(wp), PUBLIC :: r1_rcp !: = 1. / rcp [Kelvin/J] 53 REAL(wp), PUBLIC :: rau0_rcp !: = rau0 * rcp 58 54 REAL(wp), PUBLIC :: r1_rau0_rcp !: = 1. / ( rau0 * rcp ) 59 55 … … 87 83 REAL(wp), PUBLIC :: xlic = 300.33e+6_wp !: volumetric latent heat fusion of ice [J/m3] 88 84 REAL(wp), PUBLIC :: xsn = 2.8e+6_wp !: volumetric latent heat of sublimation of snow [J/m3] 85 #endif 86 #if defined key_lim3 87 REAL(wp), PUBLIC :: r1_rhoic !: 1 / rhoic 88 REAL(wp), PUBLIC :: r1_rhosn !: 1 / rhosn 89 89 #endif 90 90 !!---------------------------------------------------------------------- … … 163 163 IF(lwp) WRITE(numout,*) ' melting point of ice rt0_ice = ', rt0_ice , ' K' 164 164 165 r1_rau0 = 1._wp / rau0 166 r1_rcp = 1._wp / rcp 167 r1_rau0_rcp = 1._wp / ( rau0 * rcp ) 168 IF(lwp) WRITE(numout,*) 169 IF(lwp) WRITE(numout,*) ' volumic mass of pure water rauw = ', rauw , ' kg/m^3' 170 IF(lwp) WRITE(numout,*) ' volumic mass of reference rau0 = ', rau0 , ' kg/m^3' 171 IF(lwp) WRITE(numout,*) ' 1. / rau0 r1_rau0 = ', r1_rau0, ' m^3/kg' 172 IF(lwp) WRITE(numout,*) ' ocean specific heat rcp = ', rcp , ' J/Kelvin' 173 IF(lwp) WRITE(numout,*) ' 1. / ( rau0 * rcp ) r1_rau0_rcp = ', r1_rau0_rcp 174 175 165 IF(lwp) WRITE(numout,*) ' reference density and heat capacity now defined in eosbn2.f90' 166 176 167 #if defined key_lim3 || defined key_cice 177 168 xlsn = lfus * rhosn ! volumetric latent heat fusion of snow [J/m3] … … 180 171 lfus = xlsn / rhosn ! latent heat of fusion of fresh ice 181 172 #endif 182 173 #if defined key_lim3 174 r1_rhoic = 1._wp / rhoic 175 r1_rhosn = 1._wp / rhosn 176 #endif 183 177 IF(lwp) THEN 184 178 WRITE(numout,*)
Note: See TracChangeset
for help on using the changeset viewer.
