- Timestamp:
- 2015-05-29T11:46:03+02:00 (9 years ago)
- Location:
- branches/2014/dev_r4650_UKMO11_restart_functionality/NEMOGCM/NEMO/OPA_SRC/DOM
- Files:
-
- 7 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2014/dev_r4650_UKMO11_restart_functionality/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90
r5312 r5313 162 162 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: gphit, gphiu !: latitude of t-, u-, v- and f-points (degre) 163 163 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: gphiv, gphif !: 164 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, TARGET, DIMENSION(:,:) :: e1t, e2t !: horizontal scale factorsat t-point (m)165 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, TARGET, DIMENSION(:,:) :: e1u, e2u !: horizontal scale factorsat u-point (m)166 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, TARGET, DIMENSION(:,:) :: e1v, e2v !: horizontal scale factorsat v-point (m)167 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) 168 168 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: e1e2t !: surface at t-point (m2) 169 169 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ff !: coriolis factor (2.*omega*sin(yphi) ) (s-1) … … 262 262 263 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 264 265 265 266 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: tpol, fpol !: north fold mask (jperio= 3 or 4) … … 332 333 INTEGER FUNCTION dom_oce_alloc() 333 334 !!---------------------------------------------------------------------- 334 INTEGER, DIMENSION(1 1) :: ierr335 INTEGER, DIMENSION(12) :: ierr 335 336 !!---------------------------------------------------------------------- 336 337 ierr(:) = 0 … … 345 346 & tpol(jpiglo) , fpol(jpiglo) , STAT=ierr(2) ) 346 347 ! 347 ALLOCATE( glamt(jpi,jpj) , gphit(jpi,jpj) , e1t(jpi,jpj) , e2t(jpi,jpj) , & 348 & glamu(jpi,jpj) , gphiu(jpi,jpj) , e1u(jpi,jpj) , e2u(jpi,jpj) , & 349 & glamv(jpi,jpj) , gphiv(jpi,jpj) , e1v(jpi,jpj) , e2v(jpi,jpj) , e1e2t(jpi,jpj) , & 350 & 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) ) 351 353 ! 352 354 ALLOCATE( gdep3w_0(jpi,jpj,jpk) , e3v_0(jpi,jpj,jpk) , e3f_0 (jpi,jpj,jpk) , & … … 400 402 & vmask(jpi,jpj,jpk) , fmask(jpi,jpj,jpk), STAT=ierr(11) ) 401 403 404 ALLOCATE( wmask(jpi,jpj,jpk) , wumask(jpi,jpj,jpk), wvmask(jpi,jpj,jpk) , STAT=ierr(12) ) 405 402 406 #if defined key_noslip_accurate 403 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) ) 404 408 #endif 405 409 ! -
branches/2014/dev_r4650_UKMO11_restart_functionality/NEMOGCM/NEMO/OPA_SRC/DOM/domhgr.F90
r5312 r5313 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 ) -
branches/2014/dev_r4650_UKMO11_restart_functionality/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90
r5312 r5313 281 281 CALL lbc_lnk( fmask_i, 'F', 1._wp ) 282 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 283 293 284 294 ! 4. ocean/land mask for the elliptic equation -
branches/2014/dev_r4650_UKMO11_restart_functionality/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90
r5312 r5313 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 … … 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 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) 174 175 DO jk = 2, jpk 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)) 190 182 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) 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) 193 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)) 194 189 END DO 195 190 END DO … … 589 584 !! * Local declarations 590 585 INTEGER :: ji,jj,jk ! dummy loop indices 586 REAL(wp) :: zcoef 591 587 !!---------------------------------------------------------------------- 592 588 … … 635 631 ! t- and w- points depth 636 632 ! ---------------------- 633 ! set the isf depth as it is in the initial step 637 634 fsdept_n(:,:,1) = 0.5_wp * fse3w_n(:,:,1) 638 635 fsdepw_n(:,:,1) = 0.0_wp 639 636 fsde3w_n(:,:,1) = fsdept_n(:,:,1) - sshn(:,:) 640 DO jj = 1,jpj 641 DO ji = 1,jpi 642 DO jk = 2,mikt(ji,jj)-1 643 fsdept_n(ji,jj,jk) = gdept_0(ji,jj,jk) 644 fsdepw_n(ji,jj,jk) = gdepw_0(ji,jj,jk) 645 fsde3w_n(ji,jj,jk) = gdept_0(ji,jj,jk) - sshn(ji,jj) 646 END DO 647 IF (mikt(ji,jj) .GT. 1) THEN 648 jk = mikt(ji,jj) 649 fsdept_n(ji,jj,jk) = gdepw_0(ji,jj,jk) + 0.5_wp * fse3w_n(ji,jj,jk) 650 fsdepw_n(ji,jj,jk) = gdepw_0(ji,jj,jk) 651 fsde3w_n(ji,jj,jk) = fsdept_n(ji,jj,jk ) - sshn (ji,jj) 652 END IF 653 DO jk = mikt(ji,jj)+1, jpk 654 fsdept_n(ji,jj,jk) = fsdept_n(ji,jj,jk-1) + fse3w_n(ji,jj,jk) 637 638 DO jk = 2, jpk 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)) 655 644 fsdepw_n(ji,jj,jk) = fsdepw_n(ji,jj,jk-1) + fse3t_n(ji,jj,jk-1) 656 fsde3w_n(ji,jj,jk) = fsdept_n(ji,jj,jk ) - sshn (ji,jj) 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) 657 648 END DO 658 649 END DO 659 650 END DO 651 660 652 ! Local depth and Inverse of the local depth of the water column at u- and v- points 661 653 ! ---------------------------------------------------------------------------------- -
branches/2014/dev_r4650_UKMO11_restart_functionality/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90
r5312 r5313 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 … … 35 36 USE oce ! ocean variables 36 37 USE dom_oce ! ocean domain 37 USE sbc_oce ! surface variable (isf)38 38 USE closea ! closed seas 39 39 USE c1d ! 1D vertical configuration … … 298 298 ENDIF 299 299 300 IF ( ln_isfcav ) THEN 300 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) 301 302 ! 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)) 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 311 313 312 314 !!gm BUG in s-coordinate this does not work! … … 472 474 ! 473 475 ! (ISF) TODO build ice draft netcdf file for isomip and build the corresponding part of code 474 IF( cp_cfg == "isomip" ) THEN 475 ! 476 risfdep(:,:)=200.e0 477 misfdep(:,:)=1 478 ij0 = 1 ; ij1 = 40 479 DO jj = mj0(ij0), mj1(ij1) 480 risfdep(:,jj)=700.0_wp-(gphit(:,jj)+80.0_wp)*125.0_wp 481 END DO 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 482 483 WHERE( bathy(:,:) <= 0._wp ) risfdep(:,:) = 0._wp 483 484 ELSEIF ( cp_cfg == "isomip2" ) THEN484 ! 485 ELSEIF ( cp_cfg == "isomip2" .AND. ln_isfcav ) THEN 485 486 ! 486 487 risfdep(:,:)=0.e0 … … 534 535 IF( ln_zps .OR. ln_sco ) THEN ! zps or sco : read meter bathymetry 535 536 CALL iom_open ( 'bathy_meter.nc', inum ) 536 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 537 542 CALL iom_close( inum ) 538 ! 543 ! 539 544 risfdep(:,:)=0._wp 540 545 misfdep(:,:)=1 … … 584 589 IF ( .not. ln_sco ) THEN !== set a minimum depth ==! 585 590 ! patch to avoid case bathy = ice shelf draft and bathy between 0 and zhmin 586 WHERE (bathy == risfdep) 587 bathy = 0.0_wp ; risfdep = 0.0_wp 588 END WHERE 591 IF ( ln_isfcav ) THEN 592 WHERE (bathy == risfdep) 593 bathy = 0.0_wp ; risfdep = 0.0_wp 594 END WHERE 595 END IF 589 596 ! end patch 590 597 IF( rn_hmin < 0._wp ) THEN ; ik = - INT( rn_hmin ) ! from a nb of level … … 961 968 !!---------------------------------------------------------------------- 962 969 !! 970 INTEGER :: ji, jj, jk ! dummy loop indices 971 INTEGER :: ik, it ! temporary integers 972 LOGICAL :: ll_print ! Allow control print for debugging 973 REAL(wp) :: ze3tp , ze3wp ! Last ocean level thickness at T- and W-points 974 REAL(wp) :: zdepwp, zdepth ! Ajusted ocean depth to avoid too small e3t 975 REAL(wp) :: zmax ! Maximum depth 976 REAL(wp) :: zdiff ! temporary scalar 977 REAL(wp) :: zrefdep ! temporary scalar 978 REAL(wp), POINTER, DIMENSION(:,:,:) :: zprt 979 !!--------------------------------------------------------------------- 980 ! 981 IF( nn_timing == 1 ) CALL timing_start('zgr_zps') 982 ! 983 CALL wrk_alloc( jpi, jpj, jpk, zprt ) 984 ! 985 IF(lwp) WRITE(numout,*) 986 IF(lwp) WRITE(numout,*) ' zgr_zps : z-coordinate with partial steps' 987 IF(lwp) WRITE(numout,*) ' ~~~~~~~ ' 988 IF(lwp) WRITE(numout,*) ' mbathy is recomputed : bathy_level file is NOT used' 989 990 ll_print = .FALSE. ! Local variable for debugging 991 992 IF(lwp .AND. ll_print) THEN ! control print of the ocean depth 993 WRITE(numout,*) 994 WRITE(numout,*) 'dom_zgr_zps: bathy (in hundred of meters)' 995 CALL prihre( bathy, jpi, jpj, 1,jpi, 1, 1, jpj, 1, 1.e-2, numout ) 996 ENDIF 997 998 999 ! bathymetry in level (from bathy_meter) 1000 ! =================== 1001 zmax = gdepw_1d(jpk) + e3t_1d(jpk) ! maximum depth (i.e. the last ocean level thickness <= 2*e3t_1d(jpkm1) ) 1002 bathy(:,:) = MIN( zmax , bathy(:,:) ) ! bounded value of bathy (min already set at the end of zgr_bat) 1003 WHERE( bathy(:,:) == 0._wp ) ; mbathy(:,:) = 0 ! land : set mbathy to 0 1004 ELSE WHERE ; mbathy(:,:) = jpkm1 ! ocean : initialize mbathy to the max ocean level 1005 END WHERE 1006 1007 ! Compute mbathy for ocean points (i.e. the number of ocean levels) 1008 ! find the number of ocean levels such that the last level thickness 1009 ! is larger than the minimum of e3zps_min and e3zps_rat * e3t_1d (where 1010 ! e3t_1d is the reference level thickness 1011 DO jk = jpkm1, 1, -1 1012 zdepth = gdepw_1d(jk) + MIN( e3zps_min, e3t_1d(jk)*e3zps_rat ) 1013 WHERE( 0._wp < bathy(:,:) .AND. bathy(:,:) <= zdepth ) mbathy(:,:) = jk-1 1014 END DO 1015 1016 IF ( ln_isfcav ) CALL zgr_isf 1017 1018 ! Scale factors and depth at T- and W-points 1019 DO jk = 1, jpk ! intitialization to the reference z-coordinate 1020 gdept_0(:,:,jk) = gdept_1d(jk) 1021 gdepw_0(:,:,jk) = gdepw_1d(jk) 1022 e3t_0 (:,:,jk) = e3t_1d (jk) 1023 e3w_0 (:,:,jk) = e3w_1d (jk) 1024 END DO 1025 ! 1026 DO jj = 1, jpj 1027 DO ji = 1, jpi 1028 ik = mbathy(ji,jj) 1029 IF( ik > 0 ) THEN ! ocean point only 1030 ! max ocean level case 1031 IF( ik == jpkm1 ) THEN 1032 zdepwp = bathy(ji,jj) 1033 ze3tp = bathy(ji,jj) - gdepw_1d(ik) 1034 ze3wp = 0.5_wp * e3w_1d(ik) * ( 1._wp + ( ze3tp/e3t_1d(ik) ) ) 1035 e3t_0(ji,jj,ik ) = ze3tp 1036 e3t_0(ji,jj,ik+1) = ze3tp 1037 e3w_0(ji,jj,ik ) = ze3wp 1038 e3w_0(ji,jj,ik+1) = ze3tp 1039 gdepw_0(ji,jj,ik+1) = zdepwp 1040 gdept_0(ji,jj,ik ) = gdept_1d(ik-1) + ze3wp 1041 gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + ze3tp 1042 ! 1043 ELSE ! standard case 1044 IF( bathy(ji,jj) <= gdepw_1d(ik+1) ) THEN ; gdepw_0(ji,jj,ik+1) = bathy(ji,jj) 1045 ELSE ; gdepw_0(ji,jj,ik+1) = gdepw_1d(ik+1) 1046 ENDIF 1047 !gm Bug? check the gdepw_1d 1048 ! ... on 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) ) 1054 e3w_0(ji,jj,ik) = 0.5_wp * ( gdepw_0(ji,jj,ik+1) + gdepw_1d(ik+1) - 2._wp * gdepw_1d(ik) ) & 1055 & * ( e3w_1d(ik) / ( gdepw_1d(ik+1) - gdepw_1d(ik) ) ) 1056 ! ... on ik+1 1057 e3w_0 (ji,jj,ik+1) = e3t_0 (ji,jj,ik) 1058 e3t_0 (ji,jj,ik+1) = e3t_0 (ji,jj,ik) 1059 gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + e3t_0(ji,jj,ik) 1060 ENDIF 1061 ENDIF 1062 END DO 1063 END DO 1064 ! 1065 it = 0 1066 DO jj = 1, jpj 1067 DO ji = 1, jpi 1068 ik = mbathy(ji,jj) 1069 IF( ik > 0 ) THEN ! ocean point only 1070 e3tp (ji,jj) = e3t_0(ji,jj,ik) 1071 e3wp (ji,jj) = e3w_0(ji,jj,ik) 1072 ! test 1073 zdiff= gdepw_0(ji,jj,ik+1) - gdept_0(ji,jj,ik ) 1074 IF( zdiff <= 0._wp .AND. lwp ) THEN 1075 it = it + 1 1076 WRITE(numout,*) ' it = ', it, ' ik = ', ik, ' (i,j) = ', ji, jj 1077 WRITE(numout,*) ' bathy = ', bathy(ji,jj) 1078 WRITE(numout,*) ' gdept_0 = ', gdept_0(ji,jj,ik), ' gdepw_0 = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff 1079 WRITE(numout,*) ' e3tp = ', e3t_0 (ji,jj,ik), ' e3wp = ', e3w_0 (ji,jj,ik ) 1080 ENDIF 1081 ENDIF 1082 END DO 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) 1134 1135 ! Scale factors and depth at U-, V-, UW and VW-points 1136 DO jk = 1, jpk ! initialisation to z-scale factors 1137 e3u_0 (:,:,jk) = e3t_1d(jk) 1138 e3v_0 (:,:,jk) = e3t_1d(jk) 1139 e3uw_0(:,:,jk) = e3w_1d(jk) 1140 e3vw_0(:,:,jk) = e3w_1d(jk) 1141 END DO 1142 DO jk = 1,jpk ! Computed as the minimum of neighbooring scale factors 1143 DO jj = 1, jpjm1 1144 DO ji = 1, fs_jpim1 ! vector opt. 1145 e3u_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji+1,jj,jk) ) 1146 e3v_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji,jj+1,jk) ) 1147 e3uw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji+1,jj,jk) ) 1148 e3vw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji,jj+1,jk) ) 1149 END DO 1150 END DO 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 1167 CALL lbc_lnk( e3u_0 , 'U', 1._wp ) ; CALL lbc_lnk( e3uw_0, 'U', 1._wp ) ! lateral boundary conditions 1168 CALL lbc_lnk( e3v_0 , 'V', 1._wp ) ; CALL lbc_lnk( e3vw_0, 'V', 1._wp ) 1169 ! 1170 DO jk = 1, jpk ! set to z-scale factor if zero (i.e. along closed boundaries) 1171 WHERE( e3u_0 (:,:,jk) == 0._wp ) e3u_0 (:,:,jk) = e3t_1d(jk) 1172 WHERE( e3v_0 (:,:,jk) == 0._wp ) e3v_0 (:,:,jk) = e3t_1d(jk) 1173 WHERE( e3uw_0(:,:,jk) == 0._wp ) e3uw_0(:,:,jk) = e3w_1d(jk) 1174 WHERE( e3vw_0(:,:,jk) == 0._wp ) e3vw_0(:,:,jk) = e3w_1d(jk) 1175 END DO 1176 1177 ! Scale factor at F-point 1178 DO jk = 1, jpk ! initialisation to z-scale factors 1179 e3f_0(:,:,jk) = e3t_1d(jk) 1180 END DO 1181 DO jk = 1, jpk ! Computed as the minimum of neighbooring V-scale factors 1182 DO jj = 1, jpjm1 1183 DO ji = 1, fs_jpim1 ! vector opt. 1184 e3f_0(ji,jj,jk) = MIN( e3v_0(ji,jj,jk), e3v_0(ji+1,jj,jk) ) 1185 END DO 1186 END DO 1187 END DO 1188 CALL lbc_lnk( e3f_0, 'F', 1._wp ) ! Lateral boundary conditions 1189 ! 1190 DO jk = 1, jpk ! set to z-scale factor if zero (i.e. along closed boundaries) 1191 WHERE( e3f_0(:,:,jk) == 0._wp ) e3f_0(:,:,jk) = e3t_1d(jk) 1192 END DO 1193 !!gm bug ? : must be a do loop with mj0,mj1 1194 ! 1195 e3t_0(:,mj0(1),:) = e3t_0(:,mj0(2),:) ! we duplicate factor scales for jj = 1 and jj = 2 1196 e3w_0(:,mj0(1),:) = e3w_0(:,mj0(2),:) 1197 e3u_0(:,mj0(1),:) = e3u_0(:,mj0(2),:) 1198 e3v_0(:,mj0(1),:) = e3v_0(:,mj0(2),:) 1199 e3f_0(:,mj0(1),:) = e3f_0(:,mj0(2),:) 1200 1201 ! Control of the sign 1202 IF( MINVAL( e3t_0 (:,:,:) ) <= 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r e3t_0 <= 0' ) 1203 IF( MINVAL( e3w_0 (:,:,:) ) <= 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r e3w_0 <= 0' ) 1204 IF( MINVAL( gdept_0(:,:,:) ) < 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r gdept_0 < 0' ) 1205 IF( MINVAL( gdepw_0(:,:,:) ) < 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r gdepw_0 < 0' ) 1206 1207 ! Compute gdep3w_0 (vertical sum of e3w) 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 1228 ! ! ================= ! 1229 IF(lwp .AND. ll_print) THEN ! Control print ! 1230 ! ! ================= ! 1231 DO jj = 1,jpj 1232 DO ji = 1, jpi 1233 ik = MAX( mbathy(ji,jj), 1 ) 1234 zprt(ji,jj,1) = e3t_0 (ji,jj,ik) 1235 zprt(ji,jj,2) = e3w_0 (ji,jj,ik) 1236 zprt(ji,jj,3) = e3u_0 (ji,jj,ik) 1237 zprt(ji,jj,4) = e3v_0 (ji,jj,ik) 1238 zprt(ji,jj,5) = e3f_0 (ji,jj,ik) 1239 zprt(ji,jj,6) = gdep3w_0(ji,jj,ik) 1240 END DO 1241 END DO 1242 WRITE(numout,*) 1243 WRITE(numout,*) 'domzgr e3t(mbathy)' ; CALL prihre(zprt(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 1244 WRITE(numout,*) 1245 WRITE(numout,*) 'domzgr e3w(mbathy)' ; CALL prihre(zprt(:,:,2),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 1246 WRITE(numout,*) 1247 WRITE(numout,*) 'domzgr e3u(mbathy)' ; CALL prihre(zprt(:,:,3),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 1248 WRITE(numout,*) 1249 WRITE(numout,*) 'domzgr e3v(mbathy)' ; CALL prihre(zprt(:,:,4),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 1250 WRITE(numout,*) 1251 WRITE(numout,*) 'domzgr e3f(mbathy)' ; CALL prihre(zprt(:,:,5),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 1252 WRITE(numout,*) 1253 WRITE(numout,*) 'domzgr gdep3w(mbathy)' ; CALL prihre(zprt(:,:,6),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 1254 ENDIF 1255 ! 1256 CALL wrk_dealloc( jpi, jpj, jpk, zprt ) 1257 ! 1258 IF( nn_timing == 1 ) CALL timing_stop('zgr_zps') 1259 ! 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 !! 963 1277 INTEGER :: ji, jj, jk, jl ! dummy loop indices 964 1278 INTEGER :: ik, it ! temporary integers … … 971 1285 REAL(wp) :: zdiff ! temporary scalar 972 1286 REAL(wp) :: zrefdep ! temporary scalar 973 REAL(wp) :: zbathydiff, zrisfdepdiff 974 REAL(wp), POINTER, DIMENSION(:,:) :: zrisfdep, zbathy, zmask ! 3D workspace (ISH) 975 INTEGER , POINTER, DIMENSION(:,:) :: zmbathy, zmisfdep ! 3D workspace (ISH) 976 REAL(wp), POINTER, DIMENSION(:,:,:) :: zprt 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) 977 1290 !!--------------------------------------------------------------------- 978 1291 ! 979 IF( nn_timing == 1 ) CALL timing_start('zgr_zps') 980 ! 981 CALL wrk_alloc( jpi, jpj, jpk, zprt ) 1292 IF( nn_timing == 1 ) CALL timing_start('zgr_isf') 1293 ! 982 1294 CALL wrk_alloc( jpi, jpj, zbathy, zmask, zrisfdep) 983 CALL wrk_alloc( jpi, jpj, zmbathy, zmisfdep) 984 ! 985 IF(lwp) WRITE(numout,*) 986 IF(lwp) WRITE(numout,*) ' zgr_zps : z-coordinate with partial steps' 987 IF(lwp) WRITE(numout,*) ' ~~~~~~~ ' 988 IF(lwp) WRITE(numout,*) ' mbathy is recomputed : bathy_level file is NOT used' 989 990 ll_print = .FALSE. ! Local variable for debugging 991 992 IF(lwp .AND. ll_print) THEN ! control print of the ocean depth 993 WRITE(numout,*) 994 WRITE(numout,*) 'dom_zgr_zps: bathy (in hundred of meters)' 995 CALL prihre( bathy, jpi, jpj, 1,jpi, 1, 1, jpj, 1, 1.e-2, numout ) 996 ENDIF 997 998 ! bathymetry in level (from bathy_meter) 999 ! =================== 1000 zmax = gdepw_1d(jpk) + e3t_1d(jpk) ! maximum depth (i.e. the last ocean level thickness <= 2*e3t_1d(jpkm1) ) 1001 bathy(:,:) = MIN( zmax , bathy(:,:) ) ! bounded value of bathy (min already set at the end of zgr_bat) 1002 WHERE( bathy(:,:) == 0._wp ) ; mbathy(:,:) = 0 ! land : set mbathy to 0 1003 ELSE WHERE ; mbathy(:,:) = jpkm1 ! ocean : initialize mbathy to the max ocean level 1004 END WHERE 1005 1006 ! Compute mbathy for ocean points (i.e. the number of ocean levels) 1007 ! find the number of ocean levels such that the last level thickness 1008 ! is larger than the minimum of e3zps_min and e3zps_rat * e3t_1d (where 1009 ! e3t_1d is the reference level thickness 1010 DO jk = jpkm1, 1, -1 1011 zdepth = gdepw_1d(jk) + MIN( e3zps_min, e3t_1d(jk)*e3zps_rat ) 1012 WHERE( 0._wp < bathy(:,:) .AND. bathy(:,:) <= zdepth ) mbathy(:,:) = jk-1 1013 END DO 1295 CALL wrk_alloc( jpi, jpj, zmisfdep, zmbathy ) 1296 1297 1014 1298 ! (ISF) compute misfdep 1015 1299 WHERE( risfdep(:,:) == 0._wp .AND. bathy(:,:) .NE. 0) ; misfdep(:,:) = 1 ! open water : set misfdep to 1 … … 1055 1339 misfdep(jpi,:) = misfdep( 2 ,:) 1056 1340 ENDIF 1057 1341 1058 1342 IF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 ) THEN 1059 1343 mbathy( 1 ,:) = mbathy(jpim1,:) ! local domain is cyclic east-west 1060 1344 mbathy(jpi,:) = mbathy( 2 ,:) 1061 1345 ENDIF 1062 1346 1063 1347 ! split last cell if possible (only where water column is 2 cell or less) 1064 1348 DO jk = jpkm1, 1, -1 … … 1078 1362 END WHERE 1079 1363 END DO 1080 1364 1081 1365 1082 1366 ! Case where bathy and risfdep compatible but not the level variable mbathy/misfdep because of partial cell condition … … 1359 1643 IF( zmbathy(ji,jj) .LT. misfdep(ji ,jj+1) ) ibtestjp1 = 0 1360 1644 ibtest=MAX(ibtestim1, ibtestip1, ibtestjm1, ibtestjp1) 1361 IF( ibtest == 0 ) THEN1645 IF( ibtest == 0 .AND. misfdep(ji,jj) .GE. 2) THEN 1362 1646 mbathy(ji,jj) = 0 ; bathy(ji,jj) = 0.0_wp ; misfdep(ji,jj) = 0 ; risfdep(ji,jj) = 0.0_wp ; 1363 1647 END IF … … 1475 1759 ENDIF 1476 1760 1477 ! Scale factors and depth at T- and W-points1478 DO jk = 1, jpk ! intitialization to the reference z-coordinate1479 gdept_0(:,:,jk) = gdept_1d(jk)1480 gdepw_0(:,:,jk) = gdepw_1d(jk)1481 e3t_0 (:,:,jk) = e3t_1d (jk)1482 e3w_0 (:,:,jk) = e3w_1d (jk)1483 END DO1484 !1485 DO jj = 1, jpj1486 DO ji = 1, jpi1487 ik = mbathy(ji,jj)1488 IF( ik > 0 ) THEN ! ocean point only1489 ! max ocean level case1490 IF( ik == jpkm1 ) THEN1491 zdepwp = bathy(ji,jj)1492 ze3tp = bathy(ji,jj) - gdepw_1d(ik)1493 ze3wp = 0.5_wp * e3w_1d(ik) * ( 1._wp + ( ze3tp/e3t_1d(ik) ) )1494 e3t_0(ji,jj,ik ) = ze3tp1495 e3t_0(ji,jj,ik+1) = ze3tp1496 e3w_0(ji,jj,ik ) = ze3wp1497 e3w_0(ji,jj,ik+1) = ze3tp1498 gdepw_0(ji,jj,ik+1) = zdepwp1499 gdept_0(ji,jj,ik ) = gdept_1d(ik-1) + ze3wp1500 gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + ze3tp1501 !1502 ELSE ! standard case1503 IF( bathy(ji,jj) <= gdepw_1d(ik+1) ) THEN ; gdepw_0(ji,jj,ik+1) = bathy(ji,jj)1504 ELSE ; gdepw_0(ji,jj,ik+1) = gdepw_1d(ik+1)1505 ENDIF1506 !gm Bug? check the gdepw_1d1507 ! ... on ik1508 gdept_0(ji,jj,ik) = gdepw_1d(ik) + ( gdepw_0(ji,jj,ik+1) - gdepw_1d(ik) ) &1509 & * ((gdept_1d( ik ) - gdepw_1d(ik) ) &1510 & / ( gdepw_1d( ik+1) - gdepw_1d(ik) ))1511 e3t_0(ji,jj,ik) = e3t_1d (ik) * ( gdepw_0 (ji,jj,ik+1) - gdepw_1d(ik) ) &1512 & / ( gdepw_1d( ik+1) - gdepw_1d(ik) )1513 e3w_0(ji,jj,ik) = 0.5_wp * ( gdepw_0(ji,jj,ik+1) + gdepw_1d(ik+1) - 2._wp * gdepw_1d(ik) ) &1514 & * ( e3w_1d(ik) / ( gdepw_1d(ik+1) - gdepw_1d(ik) ) )1515 ! ... on ik+11516 e3w_0 (ji,jj,ik+1) = e3t_0 (ji,jj,ik)1517 e3t_0 (ji,jj,ik+1) = e3t_0 (ji,jj,ik)1518 gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + e3t_0(ji,jj,ik)1519 ENDIF1520 ENDIF1521 END DO1522 END DO1523 !1524 it = 01525 DO jj = 1, jpj1526 DO ji = 1, jpi1527 ik = mbathy(ji,jj)1528 IF( ik > 0 ) THEN ! ocean point only1529 e3tp (ji,jj) = e3t_0(ji,jj,ik)1530 e3wp (ji,jj) = e3w_0(ji,jj,ik)1531 ! test1532 zdiff= gdepw_0(ji,jj,ik+1) - gdept_0(ji,jj,ik )1533 IF( zdiff <= 0._wp .AND. lwp ) THEN1534 it = it + 11535 WRITE(numout,*) ' it = ', it, ' ik = ', ik, ' (i,j) = ', ji, jj1536 WRITE(numout,*) ' bathy = ', bathy(ji,jj)1537 WRITE(numout,*) ' gdept_0 = ', gdept_0(ji,jj,ik), ' gdepw_0 = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff1538 WRITE(numout,*) ' e3tp = ', e3t_0 (ji,jj,ik), ' e3wp = ', e3w_0 (ji,jj,ik )1539 ENDIF1540 ENDIF1541 END DO1542 END DO1543 !1544 ! (ISF) Definition of e3t, u, v, w for ISF case1545 DO jj = 1, jpj1546 DO ji = 1, jpi1547 ik = misfdep(ji,jj)1548 IF( ik > 1 ) THEN ! ice shelf point only1549 IF( risfdep(ji,jj) < gdepw_1d(ik) ) risfdep(ji,jj)= gdepw_1d(ik)1550 gdepw_0(ji,jj,ik) = risfdep(ji,jj)1551 !gm Bug? check the gdepw_01552 ! ... on ik1553 gdept_0(ji,jj,ik) = gdepw_1d(ik+1) - ( gdepw_1d(ik+1) - gdepw_0(ji,jj,ik) ) &1554 & * ( gdepw_1d(ik+1) - gdept_1d(ik) ) &1555 & / ( gdepw_1d(ik+1) - gdepw_1d(ik) )1556 e3t_0 (ji,jj,ik ) = gdepw_1d(ik+1) - gdepw_0(ji,jj,ik)1557 e3w_0 (ji,jj,ik+1) = gdept_1d(ik+1) - gdept_0(ji,jj,ik)1558 1559 IF( ik + 1 == mbathy(ji,jj) ) THEN ! ice shelf point only (2 cell water column)1560 e3w_0 (ji,jj,ik+1) = gdept_0(ji,jj,ik+1) - gdept_0(ji,jj,ik)1561 ENDIF1562 ! ... on ik / ik-11563 e3w_0 (ji,jj,ik ) = 2._wp * (gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik))1564 e3t_0 (ji,jj,ik-1) = gdepw_0(ji,jj,ik) - gdepw_1d(ik-1)1565 ! The next line isn't required and doesn't affect results - included for consistency with bathymetry code1566 gdept_0(ji,jj,ik-1) = gdept_1d(ik-1)1567 ENDIF1568 END DO1569 END DO1570 !1571 it = 01572 DO jj = 1, jpj1573 DO ji = 1, jpi1574 ik = misfdep(ji,jj)1575 IF( ik > 1 ) THEN ! ice shelf point only1576 e3tp (ji,jj) = e3t_0(ji,jj,ik )1577 e3wp (ji,jj) = e3w_0(ji,jj,ik+1 )1578 ! test1579 zdiff= gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik )1580 IF( zdiff <= 0. .AND. lwp ) THEN1581 it = it + 11582 WRITE(numout,*) ' it = ', it, ' ik = ', ik, ' (i,j) = ', ji, jj1583 WRITE(numout,*) ' risfdep = ', risfdep(ji,jj)1584 WRITE(numout,*) ' gdept = ', gdept_0(ji,jj,ik), ' gdepw = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff1585 WRITE(numout,*) ' e3tp = ', e3tp(ji,jj), ' e3wp = ', e3wp(ji,jj)1586 ENDIF1587 ENDIF1588 END DO1589 END DO1590 ! END (ISF)1591 1592 ! Scale factors and depth at U-, V-, UW and VW-points1593 DO jk = 1, jpk ! initialisation to z-scale factors1594 e3u_0 (:,:,jk) = e3t_1d(jk)1595 e3v_0 (:,:,jk) = e3t_1d(jk)1596 e3uw_0(:,:,jk) = e3w_1d(jk)1597 e3vw_0(:,:,jk) = e3w_1d(jk)1598 END DO1599 DO jk = 1,jpk ! Computed as the minimum of neighbooring scale factors1600 DO jj = 1, jpjm11601 DO ji = 1, fs_jpim1 ! vector opt.1602 e3u_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji+1,jj,jk) )1603 e3v_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji,jj+1,jk) )1604 e3uw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji+1,jj,jk) )1605 e3vw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji,jj+1,jk) )1606 END DO1607 END DO1608 END DO1609 ! (ISF) define e3uw1610 DO jk = 2,jpk1611 DO jj = 1, jpjm11612 DO ji = 1, fs_jpim1 ! vector opt.1613 e3uw_0(ji,jj,jk) = MIN( gdept_0(ji,jj,jk), gdept_0(ji+1,jj ,jk) ) &1614 & - MAX( gdept_0(ji,jj,jk-1), gdept_0(ji+1,jj ,jk-1) )1615 e3vw_0(ji,jj,jk) = MIN( gdept_0(ji,jj,jk), gdept_0(ji ,jj+1,jk) ) &1616 & - MAX( gdept_0(ji,jj,jk-1), gdept_0(ji ,jj+1,jk-1) )1617 END DO1618 END DO1619 END DO1620 !End (ISF)1621 1622 CALL lbc_lnk( e3u_0 , 'U', 1._wp ) ; CALL lbc_lnk( e3uw_0, 'U', 1._wp ) ! lateral boundary conditions1623 CALL lbc_lnk( e3v_0 , 'V', 1._wp ) ; CALL lbc_lnk( e3vw_0, 'V', 1._wp )1624 !1625 DO jk = 1, jpk ! set to z-scale factor if zero (i.e. along closed boundaries)1626 WHERE( e3u_0 (:,:,jk) == 0._wp ) e3u_0 (:,:,jk) = e3t_1d(jk)1627 WHERE( e3v_0 (:,:,jk) == 0._wp ) e3v_0 (:,:,jk) = e3t_1d(jk)1628 WHERE( e3uw_0(:,:,jk) == 0._wp ) e3uw_0(:,:,jk) = e3w_1d(jk)1629 WHERE( e3vw_0(:,:,jk) == 0._wp ) e3vw_0(:,:,jk) = e3w_1d(jk)1630 END DO1631 1632 ! Scale factor at F-point1633 DO jk = 1, jpk ! initialisation to z-scale factors1634 e3f_0(:,:,jk) = e3t_1d(jk)1635 END DO1636 DO jk = 1, jpk ! Computed as the minimum of neighbooring V-scale factors1637 DO jj = 1, jpjm11638 DO ji = 1, fs_jpim1 ! vector opt.1639 e3f_0(ji,jj,jk) = MIN( e3v_0(ji,jj,jk), e3v_0(ji+1,jj,jk) )1640 END DO1641 END DO1642 END DO1643 CALL lbc_lnk( e3f_0, 'F', 1._wp ) ! Lateral boundary conditions1644 !1645 DO jk = 1, jpk ! set to z-scale factor if zero (i.e. along closed boundaries)1646 WHERE( e3f_0(:,:,jk) == 0._wp ) e3f_0(:,:,jk) = e3t_1d(jk)1647 END DO1648 !!gm bug ? : must be a do loop with mj0,mj11649 !1650 e3t_0(:,mj0(1),:) = e3t_0(:,mj0(2),:) ! we duplicate factor scales for jj = 1 and jj = 21651 e3w_0(:,mj0(1),:) = e3w_0(:,mj0(2),:)1652 e3u_0(:,mj0(1),:) = e3u_0(:,mj0(2),:)1653 e3v_0(:,mj0(1),:) = e3v_0(:,mj0(2),:)1654 e3f_0(:,mj0(1),:) = e3f_0(:,mj0(2),:)1655 1656 ! Control of the sign1657 IF( MINVAL( e3t_0 (:,:,:) ) <= 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r e3t_0 <= 0' )1658 IF( MINVAL( e3w_0 (:,:,:) ) <= 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r e3w_0 <= 0' )1659 IF( MINVAL( gdept_0(:,:,:) ) < 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r gdept_0 < 0' )1660 IF( MINVAL( gdepw_0(:,:,:) ) < 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r gdepw_0 < 0' )1661 1662 ! Compute gdep3w_0 (vertical sum of e3w)1663 WHERE (misfdep == 0) misfdep = 11664 DO jj = 1,jpj1665 DO ji = 1,jpi1666 gdep3w_0(ji,jj,1) = 0.5_wp * e3w_0(ji,jj,1)1667 DO jk = 2, misfdep(ji,jj)1668 gdep3w_0(ji,jj,jk) = gdep3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk)1669 END DO1670 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))1671 DO jk = misfdep(ji,jj) + 1, jpk1672 gdep3w_0(ji,jj,jk) = gdep3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk)1673 END DO1674 END DO1675 END DO1676 ! ! ================= !1677 IF(lwp .AND. ll_print) THEN ! Control print !1678 ! ! ================= !1679 DO jj = 1,jpj1680 DO ji = 1, jpi1681 ik = MAX( mbathy(ji,jj), 1 )1682 zprt(ji,jj,1) = e3t_0 (ji,jj,ik)1683 zprt(ji,jj,2) = e3w_0 (ji,jj,ik)1684 zprt(ji,jj,3) = e3u_0 (ji,jj,ik)1685 zprt(ji,jj,4) = e3v_0 (ji,jj,ik)1686 zprt(ji,jj,5) = e3f_0 (ji,jj,ik)1687 zprt(ji,jj,6) = gdep3w_0(ji,jj,ik)1688 END DO1689 END DO1690 WRITE(numout,*)1691 WRITE(numout,*) 'domzgr e3t(mbathy)' ; CALL prihre(zprt(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)1692 WRITE(numout,*)1693 WRITE(numout,*) 'domzgr e3w(mbathy)' ; CALL prihre(zprt(:,:,2),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)1694 WRITE(numout,*)1695 WRITE(numout,*) 'domzgr e3u(mbathy)' ; CALL prihre(zprt(:,:,3),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)1696 WRITE(numout,*)1697 WRITE(numout,*) 'domzgr e3v(mbathy)' ; CALL prihre(zprt(:,:,4),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)1698 WRITE(numout,*)1699 WRITE(numout,*) 'domzgr e3f(mbathy)' ; CALL prihre(zprt(:,:,5),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)1700 WRITE(numout,*)1701 WRITE(numout,*) 'domzgr gdep3w(mbathy)' ; CALL prihre(zprt(:,:,6),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)1702 ENDIF1703 !1704 CALL wrk_dealloc( jpi, jpj, jpk, zprt )1705 1761 CALL wrk_dealloc( jpi, jpj, zmask, zbathy, zrisfdep ) 1706 1762 CALL wrk_dealloc( jpi, jpj, zmisfdep, zmbathy ) 1707 ! 1708 IF( nn_timing == 1 ) CALL timing_stop('zgr_ zps')1709 !1710 END SUBROUTINE zgr_zps1763 1764 IF( nn_timing == 1 ) CALL timing_stop('zgr_isf') 1765 1766 END SUBROUTINE 1711 1767 1712 1768 SUBROUTINE zgr_sco -
branches/2014/dev_r4650_UKMO11_restart_functionality/NEMOGCM/NEMO/OPA_SRC/DOM/istate.F90
r5312 r5313 69 69 !! ** Purpose : Initialization of the dynamics and tracer fields. 70 70 !!---------------------------------------------------------------------- 71 ! - ML - needed for initialization of e3t_b 72 INTEGER :: ji,jj,jk ! dummy loop indices 73 REAL(wp), POINTER, DIMENSION(:,:,:,:) :: zuvd ! U & V data workspace 71 INTEGER :: ji, jj, jk ! dummy loop indices 72 REAL(wp), POINTER, DIMENSION(:,:,:,:) :: zuvd ! U & V data workspace 74 73 !!---------------------------------------------------------------------- 75 74 ! … … 84 83 IF( lk_c1d ) CALL dta_uvd_init ! Initialization of U & V input data 85 84 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 85 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 92 89 93 90 IF( ln_rstart ) THEN ! Restart from a file … … 137 134 CALL eos( tsb, rhd, rhop, gdept_0(:,:,:) ) ! before potential and in situ densities 138 135 #if ! defined key_c1d 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 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 142 143 #endif 143 144 ! -
branches/2014/dev_r4650_UKMO11_restart_functionality/NEMOGCM/NEMO/OPA_SRC/DOM/phycst.F90
r5312 r5313 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] … … 51 51 REAL(wp), PUBLIC :: rcp !: ocean specific heat [J/Kelvin] 52 52 REAL(wp), PUBLIC :: r1_rcp !: = 1. / rcp [Kelvin/J] 53 REAL(wp), PUBLIC :: rau0_rcp !: = rau0 * rcp 53 54 REAL(wp), PUBLIC :: r1_rau0_rcp !: = 1. / ( rau0 * rcp ) 54 55 … … 82 83 REAL(wp), PUBLIC :: xlic = 300.33e+6_wp !: volumetric latent heat fusion of ice [J/m3] 83 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 84 89 #endif 85 90 !!---------------------------------------------------------------------- … … 166 171 lfus = xlsn / rhosn ! latent heat of fusion of fresh ice 167 172 #endif 168 173 #if defined key_lim3 174 r1_rhoic = 1._wp / rhoic 175 r1_rhosn = 1._wp / rhosn 176 #endif 169 177 IF(lwp) THEN 170 178 WRITE(numout,*)
Note: See TracChangeset
for help on using the changeset viewer.