Changeset 4990 for trunk/NEMOGCM/NEMO/OPA_SRC/DOM
- Timestamp:
- 2014-12-15T17:42:49+01:00 (9 years ago)
- Location:
- trunk/NEMOGCM/NEMO/OPA_SRC/DOM
- Files:
-
- 10 edited
-
dom_oce.F90 (modified) (3 diffs)
-
domain.F90 (modified) (2 diffs)
-
domhgr.F90 (modified) (1 diff)
-
dommsk.F90 (modified) (3 diffs)
-
domvvl.F90 (modified) (18 diffs)
-
domwri.F90 (modified) (4 diffs)
-
domzgr.F90 (modified) (18 diffs)
-
dtatsd.F90 (modified) (1 diff)
-
istate.F90 (modified) (18 diffs)
-
phycst.F90 (modified) (3 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90
r4671 r4990 176 176 LOGICAL, PUBLIC :: ln_zps !: z-coordinate - partial step 177 177 LOGICAL, PUBLIC :: ln_sco !: s-coordinate or hybrid z-s coordinate 178 LOGICAL, PUBLIC :: ln_isfcav !: presence of ISF 178 179 179 180 !! All coordinates … … 251 252 INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: mbkt !: vertical index of the bottom last T- ocean level 252 253 INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: mbku, mbkv !: vertical index of the bottom last U- and W- ocean level 253 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: bathy !: ocean depth (meters) 254 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: tmask_i !: interior domain T-point mask 255 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 256 262 257 263 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:), TARGET :: tmask, umask, vmask, fmask !: land/ocean mask at T-, U-, V- and F-pts … … 381 387 & hift (jpi,jpj) , hifu (jpi,jpj) , rx1 (jpi,jpj) , STAT=ierr(8) ) 382 388 383 ALLOCATE( mbathy(jpi,jpj) , bathy(jpi,jpj) , & 384 & tmask_i(jpi,jpj) , bmask(jpi,jpj) , & 389 ALLOCATE( mbathy(jpi,jpj) , bathy(jpi,jpj) , & 390 & tmask_i(jpi,jpj) , umask_i(jpi,jpj), vmask_i(jpi,jpj), fmask_i(jpi,jpj), & 391 & bmask(jpi,jpj) , & 385 392 & mbkt (jpi,jpj) , mbku (jpi,jpj) , mbkv(jpi,jpj) , STAT=ierr(9) ) 386 393 394 ! (ISF) Allocation of basic array 395 ALLOCATE( misfdep(jpi,jpj) , risfdep(jpi,jpj), & 396 & mikt(jpi,jpj), miku(jpi,jpj), mikv(jpi,jpj) , & 397 & mikf(jpi,jpj), ssmask(jpi,jpj), STAT=ierr(10) ) 398 387 399 ALLOCATE( tmask(jpi,jpj,jpk) , umask(jpi,jpj,jpk), & 388 & vmask(jpi,jpj,jpk) , fmask(jpi,jpj,jpk), STAT=ierr(1 0) )400 & vmask(jpi,jpj,jpk) , fmask(jpi,jpj,jpk), STAT=ierr(11) ) 389 401 390 402 #if defined key_noslip_accurate -
trunk/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90
r4624 r4990 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 -
trunk/NEMOGCM/NEMO/OPA_SRC/DOM/domhgr.F90
r4366 r4990 638 638 CALL iom_close( inum ) 639 639 640 ! need to be define for the extended grid south of -80S 641 ! some point are undefined but you need to have e1 and e2 .NE. 0 642 WHERE (e1t==0.0_wp) 643 e1t=1.0e2 644 END WHERE 645 WHERE (e1v==0.0_wp) 646 e1v=1.0e2 647 END WHERE 648 WHERE (e1u==0.0_wp) 649 e1u=1.0e2 650 END WHERE 651 WHERE (e1f==0.0_wp) 652 e1f=1.0e2 653 END WHERE 654 WHERE (e2t==0.0_wp) 655 e2t=1.0e2 656 END WHERE 657 WHERE (e2v==0.0_wp) 658 e2v=1.0e2 659 END WHERE 660 WHERE (e2u==0.0_wp) 661 e2u=1.0e2 662 END WHERE 663 WHERE (e2f==0.0_wp) 664 e2f=1.0e2 665 END WHERE 666 640 667 END SUBROUTINE hgr_read 641 668 -
trunk/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90
r4624 r4990 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 ) 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 ) 255 282 256 283 257 284 ! 4. ocean/land mask for the elliptic equation 258 285 ! -------------------------------------------- 259 bmask(:,:) = tmask(:,:,1) ! elliptic equation is written at t-point286 bmask(:,:) = ssmask(:,:) ! elliptic equation is written at t-point 260 287 ! 261 288 ! ! Boundary conditions -
trunk/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90
r4795 r4990 169 169 fsdept_b(:,:,1) = 0.5_wp * fse3w_b(:,:,1) 170 170 fsdepw_b(:,:,1) = 0.0_wp 171 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) 171 DO jj = 1,jpj 172 DO ji = 1,jpi 173 DO jk = 2,mikt(ji,jj)-1 174 fsdept_n(ji,jj,jk) = gdept_0(ji,jj,jk) 175 fsdepw_n(ji,jj,jk) = gdepw_0(ji,jj,jk) 176 fsde3w_n(ji,jj,jk) = gdept_0(ji,jj,jk) - sshn(ji,jj) 177 fsdept_b(ji,jj,jk) = gdept_0(ji,jj,jk) 178 fsdepw_b(ji,jj,jk) = gdepw_0(ji,jj,jk) 179 END DO 180 IF (mikt(ji,jj) .GT. 1) THEN 181 jk = mikt(ji,jj) 182 fsdept_n(ji,jj,jk) = gdepw_0(ji,jj,jk) + 0.5_wp * fse3w_n(ji,jj,jk) 183 fsdepw_n(ji,jj,jk) = gdepw_0(ji,jj,jk) 184 fsde3w_n(ji,jj,jk) = fsdept_n(ji,jj,jk ) - sshn (ji,jj) 185 fsdept_b(ji,jj,jk) = gdepw_0(ji,jj,jk) + 0.5_wp * fse3w_b(ji,jj,jk) 186 fsdepw_b(ji,jj,jk) = gdepw_0(ji,jj,jk) 187 END IF 188 DO jk = mikt(ji,jj)+1, jpk 189 fsdept_n(ji,jj,jk) = fsdept_n(ji,jj,jk-1) + fse3w_n(ji,jj,jk) 190 fsdepw_n(ji,jj,jk) = fsdepw_n(ji,jj,jk-1) + fse3t_n(ji,jj,jk-1) 191 fsde3w_n(ji,jj,jk) = fsdept_n(ji,jj,jk ) - sshn (ji,jj) 192 fsdept_b(ji,jj,jk) = fsdept_b(ji,jj,jk-1) + fse3w_b(ji,jj,jk) 193 fsdepw_b(ji,jj,jk) = fsdepw_b(ji,jj,jk-1) + fse3t_b(ji,jj,jk-1) 194 END DO 195 END DO 177 196 END DO 178 197 … … 185 204 hv_b(:,:) = hv_b(:,:) + fse3v_b(:,:,jk) * vmask(:,:,jk) 186 205 END DO 187 hur_b(:,:) = umask (:,:,1) / ( hu_b(:,:) + 1. - umask(:,:,1) )188 hvr_b(:,:) = vmask (:,:,1) / ( hv_b(:,:) + 1. - vmask(:,:,1) )206 hur_b(:,:) = umask_i(:,:) / ( hu_b(:,:) + 1. - umask_i(:,:) ) 207 hvr_b(:,:) = vmask_i(:,:) / ( hv_b(:,:) + 1. - vmask_i(:,:) ) 189 208 190 209 ! Restoring frequencies for z_tilde coordinate … … 293 312 ! ! --------------------------------------------- ! 294 313 295 z_scale(:,:) = ( ssha(:,:) - sshb(:,:) ) * tmask(:,:,1) / ( ht_0(:,:) + sshn(:,:) + 1. - tmask(:,:,1) )314 z_scale(:,:) = ( ssha(:,:) - sshb(:,:) ) * ssmask(:,:) / ( ht_0(:,:) + sshn(:,:) + 1. - ssmask(:,:) ) 296 315 DO jk = 1, jpkm1 297 316 ! formally this is the same as fse3t_a = e3t_0*(1+ssha/ht_0) … … 313 332 zht (:,:) = zht (:,:) + fse3t_n(:,:,jk) * tmask(:,:,jk) 314 333 END DO 315 zhdiv(:,:) = zhdiv(:,:) / ( zht(:,:) + 1. - tmask (:,:,1) )334 zhdiv(:,:) = zhdiv(:,:) / ( zht(:,:) + 1. - tmask_i(:,:) ) 316 335 317 336 ! 2 - Low frequency baroclinic horizontal divergence (z-tilde case only) … … 338 357 ELSE ! layer case 339 358 DO jk = 1, jpkm1 340 tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - fse3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) 359 tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - fse3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) * tmask(:,:,jk) 341 360 END DO 342 361 END IF … … 386 405 ! d - thickness diffusion transport: boundary conditions 387 406 ! (stored for tracer advction and continuity equation) 388 CALL lbc_lnk( un_td , 'U' , -1. )389 CALL lbc_lnk( vn_td , 'V' , -1. )407 CALL lbc_lnk( un_td , 'U' , -1._wp) 408 CALL lbc_lnk( vn_td , 'V' , -1._wp) 390 409 391 410 ! 4 - Time stepping of baroclinic scale factors … … 398 417 z2dt = 2.0_wp * rdt 399 418 ENDIF 400 CALL lbc_lnk( tilde_e3t_a(:,:,:), 'T', 1. )419 CALL lbc_lnk( tilde_e3t_a(:,:,:), 'T', 1._wp ) 401 420 tilde_e3t_a(:,:,:) = tilde_e3t_b(:,:,:) + z2dt * tmask(:,:,:) * tilde_e3t_a(:,:,:) 402 421 … … 453 472 zht(:,:) = zht(:,:) + tilde_e3t_a(:,:,jk) * tmask(:,:,jk) 454 473 END DO 455 z_scale(:,:) = - zht(:,:) / ( ht_0(:,:) + sshn(:,:) + 1. - tmask(:,:,1) )474 z_scale(:,:) = - zht(:,:) / ( ht_0(:,:) + sshn(:,:) + 1. - ssmask(:,:) ) 456 475 DO jk = 1, jpkm1 457 476 dtilde_e3t_a(:,:,jk) = dtilde_e3t_a(:,:,jk) + fse3t_n(:,:,jk) * z_scale(:,:) * tmask(:,:,jk) … … 463 482 ! ! ---baroclinic part--------- ! 464 483 DO jk = 1, jpkm1 465 fse3t_a(:,:,jk) = fse3t_a(:,:,jk) + dtilde_e3t_a(:,:,jk) 484 fse3t_a(:,:,jk) = fse3t_a(:,:,jk) + dtilde_e3t_a(:,:,jk) * tmask(:,:,jk) 466 485 END DO 467 486 ENDIF … … 531 550 END DO 532 551 ! ! 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)552 hur_a(:,:) = 1._wp / ( hu_a(:,:) + 1._wp - umask_i(:,:) ) * umask_i(:,:) 553 hvr_a(:,:) = 1._wp / ( hv_a(:,:) + 1._wp - vmask_i(:,:) ) * vmask_i(:,:) 535 554 536 555 CALL wrk_dealloc( jpi, jpj, zht, z_scale, zwu, zwv, zhdiv ) … … 570 589 !! * Local declarations 571 590 REAL(wp), POINTER, DIMENSION(:,:,:) :: z_e3t_def 572 INTEGER :: j k ! dummy loop indices591 INTEGER :: ji,jj,jk ! dummy loop indices 573 592 !!---------------------------------------------------------------------- 574 593 … … 622 641 fsdepw_n(:,:,1) = 0.0_wp 623 642 fsde3w_n(:,:,1) = fsdept_n(:,:,1) - sshn(:,:) 624 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 (:,:) 643 DO jj = 1,jpj 644 DO ji = 1,jpi 645 DO jk = 2,mikt(ji,jj)-1 646 fsdept_n(ji,jj,jk) = gdept_0(ji,jj,jk) 647 fsdepw_n(ji,jj,jk) = gdepw_0(ji,jj,jk) 648 fsde3w_n(ji,jj,jk) = gdept_0(ji,jj,jk) - sshn(ji,jj) 649 END DO 650 IF (mikt(ji,jj) .GT. 1) THEN 651 jk = mikt(ji,jj) 652 fsdept_n(ji,jj,jk) = gdepw_0(ji,jj,jk) + 0.5_wp * fse3w_n(ji,jj,jk) 653 fsdepw_n(ji,jj,jk) = gdepw_0(ji,jj,jk) 654 fsde3w_n(ji,jj,jk) = fsdept_n(ji,jj,jk ) - sshn (ji,jj) 655 END IF 656 DO jk = mikt(ji,jj)+1, jpk 657 fsdept_n(ji,jj,jk) = fsdept_n(ji,jj,jk-1) + fse3w_n(ji,jj,jk) 658 fsdepw_n(ji,jj,jk) = fsdepw_n(ji,jj,jk-1) + fse3t_n(ji,jj,jk-1) 659 fsde3w_n(ji,jj,jk) = fsdept_n(ji,jj,jk ) - sshn (ji,jj) 660 END DO 661 END DO 628 662 END DO 629 663 ! Local depth and Inverse of the local depth of the water column at u- and v- points … … 702 736 IF( l_is_orca ) CALL dom_vvl_orca_fix( pe3_in, pe3_out, pout ) 703 737 ! boundary conditions 704 CALL lbc_lnk( pe3_out(:,:,:), 'U', 1. )738 CALL lbc_lnk( pe3_out(:,:,:), 'U', 1._wp ) 705 739 pe3_out(:,:,:) = pe3_out(:,:,:) + e3u_0(:,:,:) 706 740 ! ! ------------------------------------- ! … … 720 754 IF( l_is_orca ) CALL dom_vvl_orca_fix( pe3_in, pe3_out, pout ) 721 755 ! boundary conditions 722 CALL lbc_lnk( pe3_out(:,:,:), 'V', 1. )756 CALL lbc_lnk( pe3_out(:,:,:), 'V', 1._wp ) 723 757 pe3_out(:,:,:) = pe3_out(:,:,:) + e3v_0(:,:,:) 724 758 ! ! ------------------------------------- ! … … 738 772 IF( l_is_orca ) CALL dom_vvl_orca_fix( pe3_in, pe3_out, pout ) 739 773 ! boundary conditions 740 CALL lbc_lnk( pe3_out(:,:,:), 'F', 1. )774 CALL lbc_lnk( pe3_out(:,:,:), 'F', 1._wp ) 741 775 pe3_out(:,:,:) = pe3_out(:,:,:) + e3f_0(:,:,:) 742 776 ! ! ------------------------------------- ! … … 815 849 CALL iom_get( numror, jpdom_autoglo, 'fse3t_b', fse3t_b(:,:,:) ) 816 850 CALL iom_get( numror, jpdom_autoglo, 'fse3t_n', fse3t_n(:,:,:) ) 851 ! needed to restart if land processor not computed 852 IF(lwp) write(numout,*) 'dom_vvl_rst : fse3t_b and fse3t_n found in restart files' 853 WHERE ( tmask(:,:,:) == 0.0_wp ) 854 fse3t_n(:,:,:) = e3t_0(:,:,:) 855 fse3t_b(:,:,:) = e3t_0(:,:,:) 856 END WHERE 817 857 IF( neuler == 0 ) THEN 818 858 fse3t_b(:,:,:) = fse3t_n(:,:,:) 819 859 ENDIF 820 860 ELSE IF( id1 > 0 ) THEN 861 IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : fse3t_n not found in restart files' 862 IF(lwp) write(numout,*) 'fse3t_n set equal to fse3t_b.' 863 IF(lwp) write(numout,*) 'neuler is forced to 0' 864 CALL iom_get( numror, jpdom_autoglo, 'fse3t_b', fse3t_b(:,:,:) ) 865 fse3t_n(:,:,:) = fse3t_b(:,:,:) 866 neuler = 0 867 ELSE IF( id2 > 0 ) THEN 821 868 IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : fse3t_b not found in restart files' 822 869 IF(lwp) write(numout,*) 'fse3t_b set equal to fse3t_n.' 823 870 IF(lwp) write(numout,*) 'neuler is forced to 0' 871 CALL iom_get( numror, jpdom_autoglo, 'fse3t_n', fse3t_n(:,:,:) ) 824 872 fse3t_b(:,:,:) = fse3t_n(:,:,:) 825 873 neuler = 0 … … 830 878 DO jk=1,jpk 831 879 fse3t_n(:,:,jk) = e3t_0(:,:,jk) * ( ht_0(:,:) + sshn(:,:) ) & 832 & / ( ht_0(:,:) + 1._wp - tmask(:,:,1) ) * tmask(:,:,jk) &880 & / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk) & 833 881 & + e3t_0(:,:,jk) * (1._wp -tmask(:,:,jk)) 834 882 END DO … … 964 1012 965 1013 IF( ioptio /= 1 ) CALL ctl_stop( 'Choose ONE vertical coordinate in namelist nam_vvl' ) 1014 IF( .NOT. ln_vvl_zstar .AND. nn_isf .NE. 0) CALL ctl_stop( 'Only vvl_zstar has been tested with ice shelf cavity' ) 966 1015 967 1016 IF(lwp) THEN ! Print the choice -
trunk/NEMOGCM/NEMO/OPA_SRC/DOM/domwri.F90
r4292 r4990 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 -
trunk/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90
r4687 r4990 35 35 USE oce ! ocean variables 36 36 USE dom_oce ! ocean domain 37 USE sbc_oce ! surface variable (isf) 37 38 USE closea ! closed seas 38 39 USE c1d ! 1D vertical configuration … … 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 … … 294 297 gdepw_1d(1) = 0._wp ! force first w-level to be exactly at zero 295 298 ENDIF 299 300 ! need to be like this to compute the pressure gradient with ISF. If not, level beneath the ISF are not aligned (sum(e3t) /= depth) 301 ! define e3t_0 and e3w_0 as the differences between gdept and gdepw respectively 302 DO jk = 1, jpkm1 303 e3t_1d(jk) = gdepw_1d(jk+1)-gdepw_1d(jk) 304 END DO 305 e3t_1d(jpk) = e3t_1d(jpk-1) ! we don't care because this level is masked in NEMO 306 307 DO jk = 2, jpk 308 e3w_1d(jk) = gdept_1d(jk) - gdept_1d(jk-1) 309 END DO 310 e3w_1d(1 ) = 2._wp * (gdept_1d(1) - gdepw_1d(1)) 296 311 297 312 !!gm BUG in s-coordinate this does not work! … … 451 466 END DO 452 467 ! 468 ! (ISF) TODO build ice draft netcdf file for isomip and build the corresponding part of code 469 IF( cp_cfg == "isomip" ) THEN 470 ! 471 risfdep(:,:)=200.e0 472 misfdep(:,:)=1 473 ij0 = 1 ; ij1 = 40 474 DO jj = mj0(ij0), mj1(ij1) 475 risfdep(:,jj)=700.0_wp-(gphit(:,jj)+80.0_wp)*125.0_wp 476 END DO 477 WHERE( bathy(:,:) <= 0._wp ) risfdep(:,:) = 0._wp 478 ! 479 ELSEIF ( cp_cfg == "isomip2" ) THEN 480 ! 481 risfdep(:,:)=0.e0 482 misfdep(:,:)=1 483 ij0 = 1 ; ij1 = 40 484 DO jj = mj0(ij0), mj1(ij1) 485 risfdep(:,jj)=700.0_wp-(gphit(:,jj)+80.0_wp)*125.0_wp 486 END DO 487 WHERE( bathy(:,:) <= 0._wp ) risfdep(:,:) = 0._wp 488 END IF 489 ! 453 490 ! ! ================ ! 454 491 ELSEIF( ntopo == 1 ) THEN ! read in file ! (over the local domain) … … 492 529 CALL iom_get ( inum, jpdom_data, 'Bathymetry', bathy ) 493 530 CALL iom_close( inum ) 494 ! 531 ! 532 risfdep(:,:)=0._wp 533 misfdep(:,:)=1 534 IF ( ln_isfcav ) THEN 535 CALL iom_open ( 'isf_draft_meter.nc', inum ) 536 CALL iom_get ( inum, jpdom_data, 'isf_draft', risfdep ) 537 CALL iom_close( inum ) 538 WHERE( bathy(:,:) <= 0._wp ) risfdep(:,:) = 0._wp 539 END IF 540 ! 495 541 IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN ! ORCA R2 configuration 496 542 ! … … 530 576 ! 531 577 IF ( .not. ln_sco ) THEN !== set a minimum depth ==! 578 ! patch to avoid case bathy = ice shelf draft and bathy between 0 and zhmin 579 WHERE (bathy == risfdep) 580 bathy = 0.0_wp ; risfdep = 0.0_wp 581 END WHERE 582 ! end patch 532 583 IF( rn_hmin < 0._wp ) THEN ; ik = - INT( rn_hmin ) ! from a nb of level 533 584 ELSE ; ik = MINLOC( gdepw_1d, mask = gdepw_1d > rn_hmin, dim = 1 ) ! from a depth … … 783 834 END SUBROUTINE zgr_bot_level 784 835 836 SUBROUTINE zgr_top_level 837 !!---------------------------------------------------------------------- 838 !! *** ROUTINE zgr_bot_level *** 839 !! 840 !! ** Purpose : defines the vertical index of ocean top (mik. arrays) 841 !! 842 !! ** Method : computes from misfdep with a minimum value of 1 843 !! 844 !! ** Action : mikt, miku, mikv : vertical indices of the shallowest 845 !! ocean level at t-, u- & v-points 846 !! (min value = 1) 847 !!---------------------------------------------------------------------- 848 !! 849 INTEGER :: ji, jj ! dummy loop indices 850 REAL(wp), POINTER, DIMENSION(:,:) :: zmik 851 !!---------------------------------------------------------------------- 852 ! 853 IF( nn_timing == 1 ) CALL timing_start('zgr_top_level') 854 ! 855 CALL wrk_alloc( jpi, jpj, zmik ) 856 ! 857 IF(lwp) WRITE(numout,*) 858 IF(lwp) WRITE(numout,*) ' zgr_top_level : ocean top k-index of T-, U-, V- and W-levels ' 859 IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~~~' 860 ! 861 mikt(:,:) = MAX( misfdep(:,:) , 1 ) ! top k-index of T-level (=1) 862 ! ! top k-index of W-level (=mikt) 863 DO jj = 1, jpjm1 ! top k-index of U- (U-) level 864 DO ji = 1, jpim1 865 miku(ji,jj) = MAX( mikt(ji+1,jj ) , mikt(ji,jj) ) 866 mikv(ji,jj) = MAX( mikt(ji ,jj+1) , mikt(ji,jj) ) 867 mikf(ji,jj) = MAX( mikt(ji ,jj+1) , mikt(ji,jj), mikt(ji+1,jj ), mikt(ji+1,jj+1) ) 868 END DO 869 END DO 870 871 ! converte into REAL to use lbc_lnk ; impose a min value of 1 as a zero can be set in lbclnk 872 zmik(:,:) = REAL( miku(:,:), wp ) ; CALL lbc_lnk(zmik,'U',1.) ; miku (:,:) = MAX( INT( zmik(:,:) ), 1 ) 873 zmik(:,:) = REAL( mikv(:,:), wp ) ; CALL lbc_lnk(zmik,'V',1.) ; mikv (:,:) = MAX( INT( zmik(:,:) ), 1 ) 874 zmik(:,:) = REAL( mikf(:,:), wp ) ; CALL lbc_lnk(zmik,'F',1.) ; mikf (:,:) = MAX( INT( zmik(:,:) ), 1 ) 875 ! 876 CALL wrk_dealloc( jpi, jpj, zmik ) 877 ! 878 IF( nn_timing == 1 ) CALL timing_stop('zgr_top_level') 879 ! 880 END SUBROUTINE zgr_top_level 785 881 786 882 SUBROUTINE zgr_zco … … 861 957 !!---------------------------------------------------------------------- 862 958 !! 863 INTEGER :: ji, jj, jk 959 INTEGER :: ji, jj, jk, jl ! dummy loop indices 864 960 INTEGER :: ik, it ! temporary integers 961 INTEGER :: id, jd, nprocd 962 INTEGER :: icompt, ibtest, ibtestim1, ibtestip1, ibtestjm1, ibtestjp1 ! (ISF) 865 963 LOGICAL :: ll_print ! Allow control print for debugging 866 964 REAL(wp) :: ze3tp , ze3wp ! Last ocean level thickness at T- and W-points 867 965 REAL(wp) :: zdepwp, zdepth ! Ajusted ocean depth to avoid too small e3t 868 REAL(wp) :: zmax ! Maximum depth966 REAL(wp) :: zmax, zmin ! Maximum and minimum depth 869 967 REAL(wp) :: zdiff ! temporary scalar 870 968 REAL(wp) :: zrefdep ! temporary scalar 969 REAL(wp) :: zbathydiff, zrisfdepdiff 970 REAL(wp), POINTER, DIMENSION(:,:) :: zrisfdep, zbathy, zmask ! 3D workspace (ISH) 971 INTEGER , POINTER, DIMENSION(:,:) :: zmbathy, zmisfdep ! 3D workspace (ISH) 871 972 REAL(wp), POINTER, DIMENSION(:,:,:) :: zprt 872 973 !!--------------------------------------------------------------------- … … 875 976 ! 876 977 CALL wrk_alloc( jpi, jpj, jpk, zprt ) 978 CALL wrk_alloc( jpi, jpj, zbathy, zmask, zrisfdep) 979 CALL wrk_alloc( jpi, jpj, zmbathy, zmisfdep) 877 980 ! 878 981 IF(lwp) WRITE(numout,*) … … 889 992 ENDIF 890 993 891 892 994 ! bathymetry in level (from bathy_meter) 893 995 ! =================== … … 906 1008 WHERE( 0._wp < bathy(:,:) .AND. bathy(:,:) <= zdepth ) mbathy(:,:) = jk-1 907 1009 END DO 1010 ! (ISF) compute misfdep 1011 WHERE( risfdep(:,:) == 0._wp .AND. bathy(:,:) .NE. 0) ; misfdep(:,:) = 1 ! open water : set misfdep to 1 1012 ELSEWHERE ; misfdep(:,:) = 2 ! iceshelf : initialize misfdep to second level 1013 END WHERE 1014 1015 ! Compute misfdep for ocean points (i.e. first wet level) 1016 ! find the first ocean level such that the first level thickness 1017 ! is larger than the bot_level of e3zps_min and e3zps_rat * e3t_0 (where 1018 ! e3t_0 is the reference level thickness 1019 DO jk = 2, jpkm1 1020 zdepth = gdepw_1d(jk+1) - MIN( e3zps_min, e3t_1d(jk)*e3zps_rat ) 1021 WHERE( 0._wp < risfdep(:,:) .AND. risfdep(:,:) >= zdepth ) misfdep(:,:) = jk+1 1022 END DO 1023 WHERE (risfdep(:,:) <= e3t_1d(1) .AND. risfdep(:,:) .GT. 0._wp) 1024 risfdep(:,:) = 0. ; misfdep(:,:) = 1 1025 END WHERE 1026 1027 ! 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 1028 icompt = 0 1029 ! run the bathy check 10 times to be sure all the modif in the bathy or iceshelf draft are compatible together 1030 DO jl = 1, 10 1031 WHERE (bathy(:,:) .EQ. risfdep(:,:) ) 1032 misfdep(:,:) = 0 ; risfdep(:,:) = 0._wp 1033 mbathy (:,:) = 0 ; bathy (:,:) = 0._wp 1034 END WHERE 1035 WHERE (mbathy(:,:) <= 0) 1036 misfdep(:,:) = 0; risfdep(:,:) = 0._wp 1037 mbathy (:,:) = 0; bathy (:,:) = 0._wp 1038 ENDWHERE 1039 IF( lk_mpp ) THEN 1040 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1041 CALL lbc_lnk( zbathy, 'T', 1. ) 1042 misfdep(:,:) = INT( zbathy(:,:) ) 1043 CALL lbc_lnk( risfdep, 'T', 1. ) 1044 CALL lbc_lnk( bathy, 'T', 1. ) 1045 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1046 CALL lbc_lnk( zbathy, 'T', 1. ) 1047 mbathy(:,:) = INT( zbathy(:,:) ) 1048 ENDIF 1049 IF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 ) THEN 1050 misfdep( 1 ,:) = misfdep(jpim1,:) ! local domain is cyclic east-west 1051 misfdep(jpi,:) = misfdep( 2 ,:) 1052 ENDIF 1053 1054 IF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 ) THEN 1055 mbathy( 1 ,:) = mbathy(jpim1,:) ! local domain is cyclic east-west 1056 mbathy(jpi,:) = mbathy( 2 ,:) 1057 ENDIF 1058 1059 ! split last cell if possible (only where water column is 2 cell or less) 1060 DO jk = jpkm1, 1, -1 1061 zmax = gdepw_1d(jk) + MIN( e3zps_min, e3t_1d(jk)*e3zps_rat ) 1062 WHERE( gdepw_1d(jk) < bathy(:,:) .AND. bathy(:,:) <= zmax .AND. misfdep + 1 >= mbathy) 1063 mbathy(:,:) = jk 1064 bathy(:,:) = zmax 1065 END WHERE 1066 END DO 1067 1068 ! split top cell if possible (only where water column is 2 cell or less) 1069 DO jk = 2, jpkm1 1070 zmax = gdepw_1d(jk+1) - MIN( e3zps_min, e3t_1d(jk)*e3zps_rat ) 1071 WHERE( gdepw_1d(jk+1) > risfdep(:,:) .AND. risfdep(:,:) >= zmax .AND. misfdep + 1 >= mbathy) 1072 misfdep(:,:) = jk 1073 risfdep(:,:) = zmax 1074 END WHERE 1075 END DO 1076 1077 1078 ! Case where bathy and risfdep compatible but not the level variable mbathy/misfdep because of partial cell condition 1079 DO jj = 1, jpj 1080 DO ji = 1, jpi 1081 ! find the minimum change option: 1082 ! test bathy 1083 IF (risfdep(ji,jj) .GT. 1) THEN 1084 zbathydiff =ABS(bathy(ji,jj) - (gdepw_1d(mbathy (ji,jj)+1) & 1085 & + MIN( e3zps_min, e3t_1d(mbathy (ji,jj)+1)*e3zps_rat ))) 1086 zrisfdepdiff=ABS(risfdep(ji,jj) - (gdepw_1d(misfdep(ji,jj) ) & 1087 & - MIN( e3zps_min, e3t_1d(misfdep(ji,jj)-1)*e3zps_rat ))) 1088 1089 IF (bathy(ji,jj) .GT. risfdep(ji,jj) .AND. mbathy(ji,jj) .LT. misfdep(ji,jj)) THEN 1090 IF (zbathydiff .LE. zrisfdepdiff) THEN 1091 bathy(ji,jj) = gdepw_1d(mbathy(ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj)+1)*e3zps_rat ) 1092 mbathy(ji,jj)= mbathy(ji,jj) + 1 1093 ELSE 1094 risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj)) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj)-1)*e3zps_rat ) 1095 misfdep(ji,jj) = misfdep(ji,jj) - 1 1096 END IF 1097 END IF 1098 END IF 1099 END DO 1100 END DO 1101 1102 ! At least 2 levels for water thickness at T, U, and V point. 1103 DO jj = 1, jpj 1104 DO ji = 1, jpi 1105 ! find the minimum change option: 1106 ! test bathy 1107 IF( misfdep(ji,jj) == mbathy(ji,jj) .AND. mbathy(ji,jj) .GT. 1) THEN 1108 zbathydiff =ABS(bathy(ji,jj) - (gdepw_1d(mbathy (ji,jj)+1)& 1109 & + MIN( e3zps_min,e3t_1d(mbathy (ji,jj)+1)*e3zps_rat ))) 1110 zrisfdepdiff=ABS(risfdep(ji,jj) - (gdepw_1d(misfdep(ji,jj) ) & 1111 & - MIN( e3zps_min,e3t_1d(misfdep(ji,jj)-1)*e3zps_rat ))) 1112 IF (zbathydiff .LE. zrisfdepdiff) THEN 1113 mbathy(ji,jj) = mbathy(ji,jj) + 1 1114 bathy(ji,jj) = gdepw_1d(mbathy (ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj) +1)*e3zps_rat ) 1115 ELSE 1116 misfdep(ji,jj)= misfdep(ji,jj) - 1 1117 risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj))*e3zps_rat ) 1118 END IF 1119 ENDIF 1120 END DO 1121 END DO 1122 1123 ! point V mbathy(ji,jj) EQ misfdep(ji,jj+1) 1124 DO jj = 1, jpjm1 1125 DO ji = 1, jpim1 1126 IF( misfdep(ji,jj+1) == mbathy(ji,jj) .AND. mbathy(ji,jj) .GT. 1) THEN 1127 zbathydiff =ABS(bathy(ji,jj ) - (gdepw_1d(mbathy (ji,jj)+1) & 1128 & + MIN( e3zps_min, e3t_1d(mbathy (ji,jj )+1)*e3zps_rat ))) 1129 zrisfdepdiff=ABS(risfdep(ji,jj+1) - (gdepw_1d(misfdep(ji,jj+1)) & 1130 & - MIN( e3zps_min, e3t_1d(misfdep(ji,jj+1)-1)*e3zps_rat ))) 1131 IF (zbathydiff .LE. zrisfdepdiff) THEN 1132 mbathy(ji,jj) = mbathy(ji,jj) + 1 1133 bathy(ji,jj) = gdepw_1d(mbathy (ji,jj )) & 1134 & + MIN( e3zps_min, e3t_1d(mbathy(ji,jj )+1)*e3zps_rat ) 1135 ELSE 1136 misfdep(ji,jj+1) = misfdep(ji,jj+1) - 1 1137 risfdep (ji,jj+1) = gdepw_1d(misfdep(ji,jj+1)+1) & 1138 & - MIN( e3zps_min, e3t_1d(misfdep(ji,jj+1))*e3zps_rat ) 1139 END IF 1140 ENDIF 1141 END DO 1142 END DO 1143 1144 IF( lk_mpp ) THEN 1145 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1146 CALL lbc_lnk( zbathy, 'T', 1. ) 1147 misfdep(:,:) = INT( zbathy(:,:) ) 1148 CALL lbc_lnk( risfdep, 'T', 1. ) 1149 CALL lbc_lnk( bathy, 'T', 1. ) 1150 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1151 CALL lbc_lnk( zbathy, 'T', 1. ) 1152 mbathy(:,:) = INT( zbathy(:,:) ) 1153 ENDIF 1154 ! point V misdep(ji,jj) EQ mbathy(ji,jj+1) 1155 DO jj = 1, jpjm1 1156 DO ji = 1, jpim1 1157 IF( misfdep(ji,jj) == mbathy(ji,jj+1) .AND. mbathy(ji,jj) .GT. 1) THEN 1158 zbathydiff =ABS( bathy(ji,jj+1) - (gdepw_1d(mbathy (ji,jj+1)+1) & 1159 & + MIN( e3zps_min, e3t_1d(mbathy (ji,jj+1)+1)*e3zps_rat ))) 1160 zrisfdepdiff=ABS(risfdep(ji,jj ) - (gdepw_1d(misfdep(ji,jj ) ) & 1161 & - MIN( e3zps_min, e3t_1d(misfdep(ji,jj )-1)*e3zps_rat ))) 1162 IF (zbathydiff .LE. zrisfdepdiff) THEN 1163 mbathy (ji,jj+1) = mbathy(ji,jj+1) + 1 1164 bathy (ji,jj+1) = gdepw_1d(mbathy (ji,jj+1) ) + MIN( e3zps_min, e3t_1d(mbathy (ji,jj+1)+1)*e3zps_rat ) 1165 ELSE 1166 misfdep(ji,jj) = misfdep(ji,jj) - 1 1167 risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj )+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj ) )*e3zps_rat ) 1168 END IF 1169 ENDIF 1170 END DO 1171 END DO 1172 1173 1174 IF( lk_mpp ) THEN 1175 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1176 CALL lbc_lnk( zbathy, 'T', 1. ) 1177 misfdep(:,:) = INT( zbathy(:,:) ) 1178 CALL lbc_lnk( risfdep, 'T', 1. ) 1179 CALL lbc_lnk( bathy, 'T', 1. ) 1180 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1181 CALL lbc_lnk( zbathy, 'T', 1. ) 1182 mbathy(:,:) = INT( zbathy(:,:) ) 1183 ENDIF 1184 1185 ! point U mbathy(ji,jj) EQ misfdep(ji,jj+1) 1186 DO jj = 1, jpjm1 1187 DO ji = 1, jpim1 1188 IF( misfdep(ji+1,jj) == mbathy(ji,jj) .AND. mbathy(ji,jj) .GT. 1) THEN 1189 zbathydiff =ABS( bathy(ji ,jj) - (gdepw_1d(mbathy (ji,jj)+1) & 1190 & + MIN( e3zps_min, e3t_1d(mbathy (ji ,jj)+1)*e3zps_rat ))) 1191 zrisfdepdiff=ABS(risfdep(ji+1,jj) - (gdepw_1d(misfdep(ji+1,jj)) & 1192 & - MIN( e3zps_min, e3t_1d(misfdep(ji+1,jj)-1)*e3zps_rat ))) 1193 IF (zbathydiff .LE. zrisfdepdiff) THEN 1194 mbathy(ji,jj) = mbathy(ji,jj) + 1 1195 bathy(ji,jj) = gdepw_1d(mbathy (ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj) +1)*e3zps_rat ) 1196 ELSE 1197 misfdep(ji+1,jj)= misfdep(ji+1,jj) - 1 1198 risfdep(ji+1,jj) = gdepw_1d(misfdep(ji+1,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji+1,jj))*e3zps_rat ) 1199 END IF 1200 ENDIF 1201 ENDDO 1202 ENDDO 1203 1204 IF( lk_mpp ) THEN 1205 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1206 CALL lbc_lnk( zbathy, 'T', 1. ) 1207 misfdep(:,:) = INT( zbathy(:,:) ) 1208 CALL lbc_lnk( risfdep, 'T', 1. ) 1209 CALL lbc_lnk( bathy, 'T', 1. ) 1210 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1211 CALL lbc_lnk( zbathy, 'T', 1. ) 1212 mbathy(:,:) = INT( zbathy(:,:) ) 1213 ENDIF 1214 1215 ! point U misfdep(ji,jj) EQ bathy(ji,jj+1) 1216 DO jj = 1, jpjm1 1217 DO ji = 1, jpim1 1218 IF( misfdep(ji,jj) == mbathy(ji+1,jj) .AND. mbathy(ji,jj) .GT. 1) THEN 1219 zbathydiff =ABS( bathy(ji+1,jj) - (gdepw_1d(mbathy (ji+1,jj)+1) & 1220 & + MIN( e3zps_min, e3t_1d(mbathy (ji+1,jj)+1)*e3zps_rat ))) 1221 zrisfdepdiff=ABS(risfdep(ji ,jj) - (gdepw_1d(misfdep(ji ,jj) ) & 1222 & - MIN( e3zps_min, e3t_1d(misfdep(ji ,jj)-1)*e3zps_rat ))) 1223 IF (zbathydiff .LE. zrisfdepdiff) THEN 1224 mbathy(ji+1,jj) = mbathy (ji+1,jj) + 1 1225 bathy (ji+1,jj) = gdepw_1d(mbathy (ji+1,jj) ) & 1226 & + MIN( e3zps_min, e3t_1d(mbathy (ji+1,jj) +1)*e3zps_rat ) 1227 ELSE 1228 misfdep(ji,jj) = misfdep(ji ,jj) - 1 1229 risfdep(ji,jj) = gdepw_1d(misfdep(ji ,jj)+1) & 1230 & - MIN( e3zps_min, e3t_1d(misfdep(ji ,jj) )*e3zps_rat ) 1231 END IF 1232 ENDIF 1233 ENDDO 1234 ENDDO 1235 1236 IF( lk_mpp ) THEN 1237 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1238 CALL lbc_lnk( zbathy, 'T', 1. ) 1239 misfdep(:,:) = INT( zbathy(:,:) ) 1240 CALL lbc_lnk( risfdep, 'T', 1. ) 1241 CALL lbc_lnk( bathy, 'T', 1. ) 1242 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1243 CALL lbc_lnk( zbathy, 'T', 1. ) 1244 mbathy(:,:) = INT( zbathy(:,:) ) 1245 ENDIF 1246 END DO 1247 ! end dig bathy/ice shelf to be compatible 1248 ! now fill single point in "coastline" of ice shelf, bathy, hole, and test again one cell tickness 1249 DO jl = 1,20 1250 1251 ! remove single point "bay" on isf coast line in the ice shelf draft' 1252 DO jk = 1, jpk 1253 WHERE (misfdep==0) misfdep=jpk 1254 zmask=0 1255 WHERE (misfdep .LE. jk) zmask=1 1256 DO jj = 2, jpjm1 1257 DO ji = 2, jpim1 1258 IF (misfdep(ji,jj) .EQ. jk) THEN 1259 ibtest = zmask(ji-1,jj) + zmask(ji+1,jj) + zmask(ji,jj-1) + zmask(ji,jj+1) 1260 IF (ibtest .LE. 1) THEN 1261 risfdep(ji,jj)=gdepw_1d(jk+1) ; misfdep(ji,jj)=jk+1 1262 IF (misfdep(ji,jj) .GT. mbathy(ji,jj)) misfdep(ji,jj) = jpk 1263 END IF 1264 END IF 1265 END DO 1266 END DO 1267 END DO 1268 WHERE (misfdep==jpk) 1269 misfdep=0 ; risfdep=0. ; mbathy=0 ; bathy=0. 1270 END WHERE 1271 IF( lk_mpp ) THEN 1272 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1273 CALL lbc_lnk( zbathy, 'T', 1. ) 1274 misfdep(:,:) = INT( zbathy(:,:) ) 1275 CALL lbc_lnk( risfdep, 'T', 1. ) 1276 CALL lbc_lnk( bathy, 'T', 1. ) 1277 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1278 CALL lbc_lnk( zbathy, 'T', 1. ) 1279 mbathy(:,:) = INT( zbathy(:,:) ) 1280 ENDIF 1281 1282 ! remove single point "bay" on bathy coast line beneath an ice shelf' 1283 DO jk = jpk,1,-1 1284 zmask=0 1285 WHERE (mbathy .GE. jk ) zmask=1 1286 DO jj = 2, jpjm1 1287 DO ji = 2, jpim1 1288 IF (mbathy(ji,jj) .EQ. jk .AND. misfdep(ji,jj) .GE. 2) THEN 1289 ibtest = zmask(ji-1,jj) + zmask(ji+1,jj) + zmask(ji,jj-1) + zmask(ji,jj+1) 1290 IF (ibtest .LE. 1) THEN 1291 bathy(ji,jj)=gdepw_1d(jk) ; mbathy(ji,jj)=jk-1 1292 IF (misfdep(ji,jj) .GT. mbathy(ji,jj)) mbathy(ji,jj) = 0 1293 END IF 1294 END IF 1295 END DO 1296 END DO 1297 END DO 1298 WHERE (mbathy==0) 1299 misfdep=0 ; risfdep=0. ; mbathy=0 ; bathy=0. 1300 END WHERE 1301 IF( lk_mpp ) THEN 1302 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1303 CALL lbc_lnk( zbathy, 'T', 1. ) 1304 misfdep(:,:) = INT( zbathy(:,:) ) 1305 CALL lbc_lnk( risfdep, 'T', 1. ) 1306 CALL lbc_lnk( bathy, 'T', 1. ) 1307 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1308 CALL lbc_lnk( zbathy, 'T', 1. ) 1309 mbathy(:,:) = INT( zbathy(:,:) ) 1310 ENDIF 1311 1312 ! fill hole in ice shelf 1313 zmisfdep = misfdep 1314 zrisfdep = risfdep 1315 WHERE (zmisfdep .LE. 1) zmisfdep=jpk 1316 DO jj = 2, jpjm1 1317 DO ji = 2, jpim1 1318 ibtestim1 = zmisfdep(ji-1,jj ) ; ibtestip1 = zmisfdep(ji+1,jj ) 1319 ibtestjm1 = zmisfdep(ji ,jj-1) ; ibtestjp1 = zmisfdep(ji ,jj+1) 1320 IF( zmisfdep(ji,jj) .GE. mbathy(ji-1,jj ) ) ibtestim1 = jpk!MAX(0, mbathy(ji-1,jj ) - 1) 1321 IF( zmisfdep(ji,jj) .GE. mbathy(ji+1,jj ) ) ibtestip1 = jpk!MAX(0, mbathy(ji+1,jj ) - 1) 1322 IF( zmisfdep(ji,jj) .GE. mbathy(ji ,jj-1) ) ibtestjm1 = jpk!MAX(0, mbathy(ji ,jj-1) - 1) 1323 IF( zmisfdep(ji,jj) .GE. mbathy(ji ,jj+1) ) ibtestjp1 = jpk!MAX(0, mbathy(ji ,jj+1) - 1) 1324 ibtest=MIN(ibtestim1, ibtestip1, ibtestjm1, ibtestjp1) 1325 IF( ibtest == jpk .AND. misfdep(ji,jj) .GE. 2) THEN 1326 mbathy(ji,jj) = 0 ; bathy(ji,jj) = 0.0_wp ; misfdep(ji,jj) = 0 ; risfdep(ji,jj) = 0.0_wp 1327 END IF 1328 IF( zmisfdep(ji,jj) < ibtest .AND. misfdep(ji,jj) .GE. 2) THEN 1329 misfdep(ji,jj) = ibtest 1330 risfdep(ji,jj) = gdepw_1d(ibtest) 1331 ENDIF 1332 ENDDO 1333 ENDDO 1334 1335 IF( lk_mpp ) THEN 1336 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1337 CALL lbc_lnk( zbathy, 'T', 1. ) 1338 misfdep(:,:) = INT( zbathy(:,:) ) 1339 CALL lbc_lnk( risfdep, 'T', 1. ) 1340 CALL lbc_lnk( bathy, 'T', 1. ) 1341 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1342 CALL lbc_lnk( zbathy, 'T', 1. ) 1343 mbathy(:,:) = INT( zbathy(:,:) ) 1344 ENDIF 1345 ! 1346 !! fill hole in bathymetry 1347 zmbathy (:,:)=mbathy (:,:) 1348 DO jj = 2, jpjm1 1349 DO ji = 2, jpim1 1350 ibtestim1 = zmbathy(ji-1,jj ) ; ibtestip1 = zmbathy(ji+1,jj ) 1351 ibtestjm1 = zmbathy(ji ,jj-1) ; ibtestjp1 = zmbathy(ji ,jj+1) 1352 IF( zmbathy(ji,jj) .LT. misfdep(ji-1,jj ) ) ibtestim1 = 0!MIN(jpk-1, misfdep(ji-1,jj ) + 1) 1353 IF( zmbathy(ji,jj) .LT. misfdep(ji+1,jj ) ) ibtestip1 = 0 1354 IF( zmbathy(ji,jj) .LT. misfdep(ji ,jj-1) ) ibtestjm1 = 0 1355 IF( zmbathy(ji,jj) .LT. misfdep(ji ,jj+1) ) ibtestjp1 = 0 1356 ibtest=MAX(ibtestim1, ibtestip1, ibtestjm1, ibtestjp1) 1357 IF( ibtest == 0 ) THEN 1358 mbathy(ji,jj) = 0 ; bathy(ji,jj) = 0.0_wp ; misfdep(ji,jj) = 0 ; risfdep(ji,jj) = 0.0_wp ; 1359 END IF 1360 IF( ibtest < zmbathy(ji,jj) .AND. misfdep(ji,jj) .GE. 2) THEN 1361 mbathy(ji,jj) = ibtest 1362 bathy(ji,jj) = gdepw_1d(ibtest+1) 1363 ENDIF 1364 END DO 1365 END DO 1366 IF( lk_mpp ) THEN 1367 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1368 CALL lbc_lnk( zbathy, 'T', 1. ) 1369 misfdep(:,:) = INT( zbathy(:,:) ) 1370 CALL lbc_lnk( risfdep, 'T', 1. ) 1371 CALL lbc_lnk( bathy, 'T', 1. ) 1372 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1373 CALL lbc_lnk( zbathy, 'T', 1. ) 1374 mbathy(:,:) = INT( zbathy(:,:) ) 1375 ENDIF 1376 ! if not compatible after all check (ie U point water column less than 2 cells), mask U 1377 DO jj = 1, jpjm1 1378 DO ji = 1, jpim1 1379 IF (mbathy(ji,jj) == misfdep(ji+1,jj) .AND. mbathy(ji,jj) .GE. 1 .AND. mbathy(ji+1,jj) .GE. 1) THEN 1380 mbathy(ji,jj) = mbathy(ji,jj) - 1 ; bathy(ji,jj) = gdepw_1d(mbathy(ji,jj)+1) ; 1381 END IF 1382 END DO 1383 END DO 1384 IF( lk_mpp ) THEN 1385 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1386 CALL lbc_lnk( zbathy, 'T', 1. ) 1387 misfdep(:,:) = INT( zbathy(:,:) ) 1388 CALL lbc_lnk( risfdep, 'T', 1. ) 1389 CALL lbc_lnk( bathy, 'T', 1. ) 1390 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1391 CALL lbc_lnk( zbathy, 'T', 1. ) 1392 mbathy(:,:) = INT( zbathy(:,:) ) 1393 ENDIF 1394 ! if not compatible after all check (ie U point water column less than 2 cells), mask U 1395 DO jj = 1, jpjm1 1396 DO ji = 1, jpim1 1397 IF (misfdep(ji,jj) == mbathy(ji+1,jj) .AND. mbathy(ji,jj) .GE. 1 .AND. mbathy(ji+1,jj) .GE. 1) THEN 1398 mbathy(ji+1,jj) = mbathy(ji+1,jj) - 1; bathy(ji+1,jj) = gdepw_1d(mbathy(ji+1,jj)+1) ; 1399 END IF 1400 END DO 1401 END DO 1402 IF( lk_mpp ) THEN 1403 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1404 CALL lbc_lnk( zbathy, 'T', 1. ) 1405 misfdep(:,:) = INT( zbathy(:,:) ) 1406 CALL lbc_lnk( risfdep, 'T', 1. ) 1407 CALL lbc_lnk( bathy, 'T', 1. ) 1408 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1409 CALL lbc_lnk( zbathy, 'T', 1. ) 1410 mbathy(:,:) = INT( zbathy(:,:) ) 1411 ENDIF 1412 ! if not compatible after all check (ie V point water column less than 2 cells), mask V 1413 DO jj = 1, jpjm1 1414 DO ji = 1, jpi 1415 IF (mbathy(ji,jj) == misfdep(ji,jj+1) .AND. mbathy(ji,jj) .GE. 1 .AND. mbathy(ji,jj+1) .GE. 1) THEN 1416 mbathy(ji,jj) = mbathy(ji,jj) - 1 ; bathy(ji,jj) = gdepw_1d(mbathy(ji,jj)+1) ; 1417 END IF 1418 END DO 1419 END DO 1420 IF( lk_mpp ) THEN 1421 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1422 CALL lbc_lnk( zbathy, 'T', 1. ) 1423 misfdep(:,:) = INT( zbathy(:,:) ) 1424 CALL lbc_lnk( risfdep, 'T', 1. ) 1425 CALL lbc_lnk( bathy, 'T', 1. ) 1426 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1427 CALL lbc_lnk( zbathy, 'T', 1. ) 1428 mbathy(:,:) = INT( zbathy(:,:) ) 1429 ENDIF 1430 ! if not compatible after all check (ie V point water column less than 2 cells), mask V 1431 DO jj = 1, jpjm1 1432 DO ji = 1, jpi 1433 IF (misfdep(ji,jj) == mbathy(ji,jj+1) .AND. mbathy(ji,jj) .GE. 1 .AND. mbathy(ji,jj+1) .GE. 1) THEN 1434 mbathy(ji,jj+1) = mbathy(ji,jj+1) - 1 ; bathy(ji,jj+1) = gdepw_1d(mbathy(ji,jj+1)+1) ; 1435 END IF 1436 END DO 1437 END DO 1438 IF( lk_mpp ) THEN 1439 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1440 CALL lbc_lnk( zbathy, 'T', 1. ) 1441 misfdep(:,:) = INT( zbathy(:,:) ) 1442 CALL lbc_lnk( risfdep, 'T', 1. ) 1443 CALL lbc_lnk( bathy, 'T', 1. ) 1444 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1445 CALL lbc_lnk( zbathy, 'T', 1. ) 1446 mbathy(:,:) = INT( zbathy(:,:) ) 1447 ENDIF 1448 ! if not compatible after all check, mask T 1449 DO jj = 1, jpj 1450 DO ji = 1, jpi 1451 IF (mbathy(ji,jj) <= misfdep(ji,jj)) THEN 1452 misfdep(ji,jj) = 0 ; risfdep(ji,jj) = 0._wp ; mbathy(ji,jj) = 0 ; bathy(ji,jj) = 0._wp ; 1453 END IF 1454 END DO 1455 END DO 1456 1457 WHERE (mbathy(:,:) == 1) 1458 mbathy = 0; bathy = 0.0_wp ; misfdep = 0 ; risfdep = 0.0_wp 1459 END WHERE 1460 END DO 1461 ! end check compatibility ice shelf/bathy 1462 ! remove very shallow ice shelf (less than ~ 10m if 75L) 1463 WHERE (misfdep(:,:) <= 5) 1464 misfdep = 1; risfdep = 0.0_wp; 1465 END WHERE 1466 1467 IF( icompt == 0 ) THEN 1468 IF(lwp) WRITE(numout,*)' no points with ice shelf too close to bathymetry' 1469 ELSE 1470 IF(lwp) WRITE(numout,*)' ',icompt,' ocean grid points with ice shelf thickness reduced to avoid bathymetry' 1471 ENDIF 908 1472 909 1473 ! Scale factors and depth at T- and W-points … … 938 1502 !gm Bug? check the gdepw_1d 939 1503 ! ... 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) ))1504 gdept_0(ji,jj,ik) = gdepw_1d(ik) + ( gdepw_0(ji,jj,ik+1) - gdepw_1d(ik) ) & 1505 & * ((gdept_1d( ik ) - gdepw_1d(ik) ) & 1506 & / ( gdepw_1d( ik+1) - gdepw_1d(ik) )) 943 1507 e3t_0(ji,jj,ik) = e3t_1d (ik) * ( gdepw_0 (ji,jj,ik+1) - gdepw_1d(ik) ) & 944 1508 & / ( gdepw_1d( ik+1) - gdepw_1d(ik) ) … … 973 1537 END DO 974 1538 END DO 1539 ! 1540 ! (ISF) Definition of e3t, u, v, w for ISF case 1541 DO jj = 1, jpj 1542 DO ji = 1, jpi 1543 ik = misfdep(ji,jj) 1544 IF( ik > 1 ) THEN ! ice shelf point only 1545 IF( risfdep(ji,jj) < gdepw_1d(ik) ) risfdep(ji,jj)= gdepw_1d(ik) 1546 gdepw_0(ji,jj,ik) = risfdep(ji,jj) 1547 !gm Bug? check the gdepw_0 1548 ! ... on ik 1549 gdept_0(ji,jj,ik) = gdepw_1d(ik+1) - ( gdepw_1d(ik+1) - gdepw_0(ji,jj,ik) ) & 1550 & * ( gdepw_1d(ik+1) - gdept_1d(ik) ) & 1551 & / ( gdepw_1d(ik+1) - gdepw_1d(ik) ) 1552 e3t_0 (ji,jj,ik ) = gdepw_1d(ik+1) - gdepw_0(ji,jj,ik) 1553 e3w_0 (ji,jj,ik+1) = gdept_1d(ik+1) - gdept_0(ji,jj,ik) 1554 1555 IF( ik + 1 == mbathy(ji,jj) ) THEN ! ice shelf point only (2 cell water column) 1556 e3w_0 (ji,jj,ik+1) = gdept_0(ji,jj,ik+1) - gdept_0(ji,jj,ik) 1557 ENDIF 1558 ! ... on ik / ik-1 1559 e3w_0 (ji,jj,ik ) = 2._wp * (gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik)) 1560 e3t_0 (ji,jj,ik-1) = gdepw_0(ji,jj,ik) - gdepw_1d(ik-1) 1561 ! The next line isn't required and doesn't affect results - included for consistency with bathymetry code 1562 gdept_0(ji,jj,ik-1) = gdept_1d(ik-1) 1563 ENDIF 1564 END DO 1565 END DO 1566 ! 1567 it = 0 1568 DO jj = 1, jpj 1569 DO ji = 1, jpi 1570 ik = misfdep(ji,jj) 1571 IF( ik > 1 ) THEN ! ice shelf point only 1572 e3tp (ji,jj) = e3t_0(ji,jj,ik ) 1573 e3wp (ji,jj) = e3w_0(ji,jj,ik+1 ) 1574 ! test 1575 zdiff= gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik ) 1576 IF( zdiff <= 0. .AND. lwp ) THEN 1577 it = it + 1 1578 WRITE(numout,*) ' it = ', it, ' ik = ', ik, ' (i,j) = ', ji, jj 1579 WRITE(numout,*) ' risfdep = ', risfdep(ji,jj) 1580 WRITE(numout,*) ' gdept = ', gdept_0(ji,jj,ik), ' gdepw = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff 1581 WRITE(numout,*) ' e3tp = ', e3tp(ji,jj), ' e3wp = ', e3wp(ji,jj) 1582 ENDIF 1583 ENDIF 1584 END DO 1585 END DO 1586 ! END (ISF) 975 1587 976 1588 ! Scale factors and depth at U-, V-, UW and VW-points … … 991 1603 END DO 992 1604 END DO 1605 ! (ISF) define e3uw 1606 DO jk = 2,jpk 1607 DO jj = 1, jpjm1 1608 DO ji = 1, fs_jpim1 ! vector opt. 1609 e3uw_0(ji,jj,jk) = MIN( gdept_0(ji,jj,jk), gdept_0(ji+1,jj ,jk) ) & 1610 & - MAX( gdept_0(ji,jj,jk-1), gdept_0(ji+1,jj ,jk-1) ) 1611 e3vw_0(ji,jj,jk) = MIN( gdept_0(ji,jj,jk), gdept_0(ji ,jj+1,jk) ) & 1612 & - MAX( gdept_0(ji,jj,jk-1), gdept_0(ji ,jj+1,jk-1) ) 1613 END DO 1614 END DO 1615 END DO 1616 !End (ISF) 1617 993 1618 CALL lbc_lnk( e3u_0 , 'U', 1._wp ) ; CALL lbc_lnk( e3uw_0, 'U', 1._wp ) ! lateral boundary conditions 994 1619 CALL lbc_lnk( e3v_0 , 'V', 1._wp ) ; CALL lbc_lnk( e3vw_0, 'V', 1._wp ) … … 1032 1657 1033 1658 ! 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 1659 WHERE (misfdep == 0) misfdep = 1 1660 DO jj = 1,jpj 1661 DO ji = 1,jpi 1662 gdep3w_0(ji,jj,1) = 0.5_wp * e3w_0(ji,jj,1) 1663 DO jk = 2, misfdep(ji,jj) 1664 gdep3w_0(ji,jj,jk) = gdep3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk) 1665 END DO 1666 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)) 1667 DO jk = misfdep(ji,jj) + 1, jpk 1668 gdep3w_0(ji,jj,jk) = gdep3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk) 1669 END DO 1670 END DO 1671 END DO 1039 1672 ! ! ================= ! 1040 1673 IF(lwp .AND. ll_print) THEN ! Control print ! … … 1066 1699 ! 1067 1700 CALL wrk_dealloc( jpi, jpj, jpk, zprt ) 1701 CALL wrk_dealloc( jpi, jpj, zmask, zbathy, zrisfdep ) 1702 CALL wrk_dealloc( jpi, jpj, zmisfdep, zmbathy ) 1068 1703 ! 1069 1704 IF( nn_timing == 1 ) CALL timing_stop('zgr_zps') -
trunk/NEMOGCM/NEMO/OPA_SRC/DOM/dtatsd.F90
r4624 r4990 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 -
trunk/NEMOGCM/NEMO/OPA_SRC/DOM/istate.F90
r4370 r4990 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) … … 73 74 !!---------------------------------------------------------------------- 74 75 ! 75 IF( nn_timing == 1 ) CALL timing_start('istate_init')76 ! 77 78 IF(lwp) WRITE(numout,*) 76 IF( nn_timing == 1 ) CALL timing_start('istate_init') 77 ! 78 79 IF(lwp) WRITE(numout,*) ' ' 79 80 IF(lwp) WRITE(numout,*) 'istate_ini : Initialization of the dynamics and tracers' 80 81 IF(lwp) WRITE(numout,*) '~~~~~~~~~~' … … 83 84 IF( lk_c1d ) CALL dta_uvd_init ! Initialization of U & V input data 84 85 85 rhd (:,:,: ) = 0.e0 86 rhop (:,:,: ) = 0.e0 87 rn2 (:,:,: ) = 0.e0 88 tsa (:,:,:,:) = 0.e0 86 rhd (:,:,: ) = 0._wp 87 rhop (:,:,: ) = 0._wp 88 rn2 (:,:,: ) = 0._wp 89 tsa (:,:,:,:) = 0._wp 90 rab_b(:,:,:,:) = 0._wp 91 rab_n(:,:,:,:) = 0._wp 89 92 90 93 IF( ln_rstart ) THEN ! Restart from a file … … 110 113 ELSEIF( cp_cfg == 'gyre' ) THEN 111 114 CALL istate_gyre ! GYRE configuration : start from pre-defined T-S fields 115 ELSEIF( cp_cfg == 'isomip' .OR. cp_cfg == 'isomip2') THEN 116 IF(lwp) WRITE(numout,*) 'Initialization of T+S for ISOMIP domain' 117 tsn(:,:,:,jp_tem)=-1.9*tmask(:,:,:) ! ISOMIP configuration : start from constant T+S fields 118 tsn(:,:,:,jp_sal)=34.4*tmask(:,:,:) 119 tsb(:,:,:,:)=tsn(:,:,:,:) 112 120 ELSE ! Initial T-S, U-V fields read in files 113 121 IF ( ln_tsd_init ) THEN ! read 3D T and S data at nit000 … … 129 137 CALL eos( tsb, rhd, rhop, gdept_0(:,:,:) ) ! before potential and in situ densities 130 138 #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 139 IF( ln_zps ) CALL zps_hde( nit000, jpts, tsb, gtsu, gtsv, & ! Partial steps: before horizontal gradient 140 & rhd, gru , grv, aru, arv, gzu, gzv, ge3ru, ge3rv, & ! 141 & gtui, gtvi, grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi ) ! of t, s, rd at the last ocean level 133 142 #endif 134 143 ! … … 162 171 ! 163 172 DO jk = 1, jpkm1 164 #if defined key_vectopt_loop165 DO jj = 1, 1 !Vector opt. => forced unrolling166 DO ji = 1, jpij167 #else168 173 DO jj = 1, jpj 169 174 DO ji = 1, jpi 170 #endif171 175 un_b(ji,jj) = un_b(ji,jj) + fse3u_n(ji,jj,jk) * un(ji,jj,jk) * umask(ji,jj,jk) 172 176 vn_b(ji,jj) = vn_b(ji,jj) + fse3v_n(ji,jj,jk) * vn(ji,jj,jk) * vmask(ji,jj,jk) … … 185 189 ! 186 190 ! 187 IF( nn_timing == 1 ) CALL timing_stop('istate_init')191 IF( nn_timing == 1 ) CALL timing_stop('istate_init') 188 192 ! 189 193 END SUBROUTINE istate_init 194 190 195 191 196 SUBROUTINE istate_t_s … … 219 224 END SUBROUTINE istate_t_s 220 225 226 221 227 SUBROUTINE istate_eel 222 228 !!---------------------------------------------------------------------- … … 233 239 USE divcur ! hor. divergence & rel. vorticity (div_cur routine) 234 240 USE iom 235 241 ! 236 242 INTEGER :: inum ! temporary logical unit 237 243 INTEGER :: ji, jj, jk ! dummy loop indices … … 244 250 REAL(wp), DIMENSION(jpiglo,jpjglo) :: zssh ! initial ssh over the global domain 245 251 !!---------------------------------------------------------------------- 246 252 ! 247 253 SELECT CASE ( jp_cfg ) 248 254 ! ! ==================== … … 375 381 INTEGER, PARAMETER :: ntsinit = 0 ! (0/1) (analytical/input data files) T&S initialization 376 382 !!---------------------------------------------------------------------- 377 383 ! 378 384 SELECT CASE ( ntsinit) 379 385 ! 380 386 CASE ( 0 ) ! analytical T/S profil deduced from LEVITUS 381 387 IF(lwp) WRITE(numout,*) 382 388 IF(lwp) WRITE(numout,*) 'istate_gyre : initial analytical T and S profil deduced from LEVITUS ' 383 389 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 384 390 ! 385 391 DO jk = 1, jpk 386 392 DO jj = 1, jpj … … 407 413 END DO 408 414 END DO 409 415 ! 410 416 CASE ( 1 ) ! T/S data fields read in dta_tem.nc/data_sal.nc files 411 417 IF(lwp) WRITE(numout,*) … … 431 437 tsn(:,:,:,jp_sal) = tsn(:,:,:,jp_sal) * tmask(:,:,:) 432 438 tsb(:,:,:,jp_sal) = tsn(:,:,:,jp_sal) 433 439 ! 434 440 END SELECT 435 441 ! 436 442 IF(lwp) THEN 437 443 WRITE(numout,*) … … 440 446 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 447 ENDIF 442 448 ! 443 449 END SUBROUTINE istate_gyre 450 444 451 445 452 SUBROUTINE istate_uvg … … 457 464 USE divcur ! hor. divergence & rel. vorticity (div_cur routine) 458 465 USE lbclnk ! ocean lateral boundary condition (or mpp link) 459 466 ! 460 467 INTEGER :: ji, jj, jk ! dummy loop indices 461 468 INTEGER :: indic ! ??? … … 567 574 !!===================================================================== 568 575 END MODULE istate 569 -
trunk/NEMOGCM/NEMO/OPA_SRC/DOM/phycst.F90
r4689 r4990 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/kg/K] 57 REAL(wp), PUBLIC :: r1_rcp !: = 1. / rcp [kg.K/J] 51 REAL(wp), PUBLIC :: rcp !: ocean specific heat [J/Kelvin] 52 REAL(wp), PUBLIC :: r1_rcp !: = 1. / rcp [Kelvin/J] 58 53 REAL(wp), PUBLIC :: r1_rau0_rcp !: = 1. / ( rau0 * rcp ) 59 54 … … 69 64 #if defined key_lim3 || defined key_cice 70 65 REAL(wp), PUBLIC :: rhoic = 917._wp !: volumic mass of sea ice [kg/m3] 71 REAL(wp), PUBLIC :: rcdic = 2.034396_wp !: thermal conductivity of fresh ice [W/m/K]72 REAL(wp), PUBLIC :: rcdsn = 0.31_wp !: thermal conductivity of snow [W/m/K]73 REAL(wp), PUBLIC :: cpic = 2067.0_wp !: specific heat for ice [J/kg/K]66 REAL(wp), PUBLIC :: rcdic = 2.034396_wp !: thermal conductivity of fresh ice 67 REAL(wp), PUBLIC :: rcdsn = 0.31_wp !: thermal conductivity of snow 68 REAL(wp), PUBLIC :: cpic = 2067.0_wp !: specific heat for ice 74 69 REAL(wp), PUBLIC :: lsub = 2.834e+6_wp !: pure ice latent heat of sublimation [J/kg] 75 70 REAL(wp), PUBLIC :: lfus = 0.334e+6_wp !: latent heat of fusion of fresh ice [J/kg] 76 REAL(wp), PUBLIC :: tmut = 0.054_wp !: decrease of seawater meltpoint with salinity [degC/ppt]71 REAL(wp), PUBLIC :: tmut = 0.054_wp !: decrease of seawater meltpoint with salinity 77 72 REAL(wp), PUBLIC :: xlsn !: = lfus*rhosn (volumetric latent heat fusion of snow) [J/m3] 78 73 #else … … 163 158 IF(lwp) WRITE(numout,*) ' melting point of ice rt0_ice = ', rt0_ice , ' K' 164 159 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 160 IF(lwp) WRITE(numout,*) ' reference density and heat capacity now defined in eosbn2.f90' 161 176 162 #if defined key_lim3 || defined key_cice 177 163 xlsn = lfus * rhosn ! volumetric latent heat fusion of snow [J/m3]
Note: See TracChangeset
for help on using the changeset viewer.
