Changeset 3865
- Timestamp:
- 2013-04-09T18:34:38+02:00 (11 years ago)
- Location:
- branches/2013/dev_r3858_NOC_ZTC/NEMOGCM
- Files:
-
- 38 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/CONFIG/ORCA2_LIM/EXP00/namelist
r3795 r3865 638 638 ln_dynadv_cen2= .false. ! flux form - 2nd order centered scheme 639 639 ln_dynadv_ubs = .false. ! flux form - 3rd order UBS scheme 640 / 641 !----------------------------------------------------------------------- 642 &nam_vvl ! vertical coordinate options 643 !----------------------------------------------------------------------- 644 ln_vvl_zstar = .false. ! zstar vertical coordinate 645 ln_vvl_ztilde = .true. ! hybrid verticalcoordinate: only high frequency variations 646 ln_vvl_layer = .false. ! full layer vertical coordinate 647 rn_ahe3 = 0.e0 ! thickness diffusion coefficient 648 ln_vvl_dbg = .true. ! debug prints (T/F) 640 649 / 641 650 !----------------------------------------------------------------------- -
branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/LIM_SRC_2/limthd_2.F90
r3625 r3865 237 237 238 238 ! energy needed to bring ocean surface layer until its freezing 239 qcmif (ji,jj) = rau0 * rcp * fse3t_m(ji,jj ,1) * ( tfu(ji,jj) - sst_m(ji,jj) - rt0 ) * ( 1 - zinda )239 qcmif (ji,jj) = rau0 * rcp * fse3t_m(ji,jj) * ( tfu(ji,jj) - sst_m(ji,jj) - rt0 ) * ( 1 - zinda ) 240 240 241 241 ! calculate oceanic heat flux. -
branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DIA/diaar5.F90
r3862 r3865 196 196 thick0(:,:) = 0._wp 197 197 DO jk = 1, jpkm1 198 vol0 = vol0 + SUM( area (:,:) * tmask(:,:,jk) * fse3t_0(:,:,jk) )199 thick0(:,:) = thick0(:,:) + tmask_i(:,:) * tmask(:,:,jk) * fse3t_0(:,:,jk)198 vol0 = vol0 + SUM( area (:,:) * tmask(:,:,jk) * e3t_0(:,:,jk) ) 199 thick0(:,:) = thick0(:,:) + tmask_i(:,:) * tmask(:,:,jk) * e3t_0(:,:,jk) 200 200 END DO 201 201 IF( lk_mpp ) CALL mpp_sum( vol0 ) … … 212 212 ik = mbkt(ji,jj) 213 213 IF( ik > 1 ) THEN 214 zztmp = ( gdept_1d(ik) - fsdept_0(ji,jj,ik) ) / ( gdept_1d(ik) - gdept_1d(ik-1) )214 zztmp = ( gdept_1d(ik) - gdept_0(ji,jj,ik) ) / ( gdept_1d(ik) - gdept_1d(ik-1) ) 215 215 sn0(ji,jj,ik) = ( 1._wp - zztmp ) * sn0(ji,jj,ik) + zztmp * sn0(ji,jj,ik-1) 216 216 ENDIF -
branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DIA/diahth.F90
r3862 r3865 326 326 htc3(ji,jj) = htc3(ji,jj) + tsn(ji,jj,ilevel+1,jp_tem) * MIN( fse3t(ji,jj,ilevel+1), zthick(ji,jj) ) & 327 327 * tmask(ji,jj,ilevel+1) 328 htc3(ji,jj) = htc3(ji,jj) + tsn(ji,jj,ilevel+1,jp_tem) * MIN( fse3t(ji,jj,ilevel+1), zthick(ji,jj) ) &329 & * tmask(ji,jj,ilevel+1)330 328 END DO 331 329 END DO -
branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DIA/diaptr.F90
r3862 r3865 259 259 ! 260 260 #if defined key_mpp_mpi 261 ijpjjpk = jpj*jpk 261 262 ish(1) = ijpjjpk ; ish2(1) = jpj ; ish2(2) = jpk 262 263 zwork(1:ijpjjpk) = RESHAPE( p_fval, ish ) … … 314 315 END DO 315 316 #if defined key_mpp_mpi 317 ijpjjpk = jpj*jpk 316 318 ish(1) = jpj*jpk ; ish2(1) = jpj ; ish2(2) = jpk 317 319 zwork(1:ijpjjpk)= RESHAPE( p_fval, ish ) -
branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90
r3862 r3865 142 142 ENDIF 143 143 144 CALL iom_put( "toce" , tsn(:,:,:,jp_tem) ) ! temperature 145 CALL iom_put( "soce" , tsn(:,:,:,jp_sal) ) ! salinity 146 CALL iom_put( "sst" , tsn(:,:,1,jp_tem) ) ! sea surface temperature 147 CALL iom_put( "sst2" , tsn(:,:,1,jp_tem) * tsn(:,:,1,jp_tem) ) ! square of sea surface temperature 148 CALL iom_put( "sss" , tsn(:,:,1,jp_sal) ) ! sea surface salinity 149 CALL iom_put( "sss2" , tsn(:,:,1,jp_sal) * tsn(:,:,1,jp_sal) ) ! square of sea surface salinity 150 CALL iom_put( "uoce" , un ) ! i-current 151 CALL iom_put( "suoce" , un(:,:,1) ) ! surface i-current 152 CALL iom_put( "voce" , vn ) ! j-current 153 CALL iom_put( "svoce" , vn(:,:,1) ) ! surface j-current 154 155 CALL iom_put( "avt" , avt ) ! T vert. eddy diff. coef. 156 CALL iom_put( "avm" , avmu ) ! T vert. eddy visc. coef. 144 IF( lk_vvl ) THEN 145 z3d(:,:,:) = tsn(:,:,:,jp_tem) * fse3t_n(:,:,:) 146 CALL iom_put( "toce" , z3d ) ! heat content 147 CALL iom_put( "sst" , z3d(:,:,1) ) ! sea surface heat content 148 CALL iom_put( "sst2" , z3d(:,:,1) * z3d(:,:,1) ) ! sea surface content of squared temperature 149 z3d(:,:,:) = tsn(:,:,:,jp_sal) * fse3t_n(:,:,:) 150 CALL iom_put( "soce" , z3d ) ! salinity content 151 CALL iom_put( "sss" , z3d(:,:,1) ) ! sea surface salinity content 152 CALL iom_put( "sss2" , z3d(:,:,1) * z3d(:,:,1) ) ! sea surface content of squared salinity 153 ELSE 154 CALL iom_put( "toce" , tsn(:,:,:,jp_tem) ) ! temperature 155 CALL iom_put( "sst" , tsn(:,:,1,jp_tem) ) ! sea surface temperature 156 CALL iom_put( "sst2" , tsn(:,:,1,jp_tem) * tsn(:,:,1,jp_tem) ) ! square of sea surface temperature 157 CALL iom_put( "soce" , tsn(:,:,:,jp_sal) ) ! salinity 158 CALL iom_put( "sss" , tsn(:,:,1,jp_sal) ) ! sea surface salinity 159 CALL iom_put( "sss2" , tsn(:,:,1,jp_sal) * tsn(:,:,1,jp_sal) ) ! square of sea surface salinity 160 END IF 161 IF( lk_vvl .AND. (.NOT. ln_dynadv_vec) ) THEN 162 CALL iom_put( "uoce" , un(:,:,:) * fse3u_n(:,:,:) ) ! i-transport 163 CALL iom_put( "voce" , vn(:,:,:) * fse3v_n(:,:,:) ) ! j-transport 164 ELSE 165 CALL iom_put( "uoce" , un ) ! i-current 166 CALL iom_put( "voce" , vn ) ! j-current 167 END IF 168 CALL iom_put( "avt" , avt ) ! T vert. eddy diff. coef. 169 CALL iom_put( "avm" , avmu ) ! T vert. eddy visc. coef. 157 170 IF( lk_zdfddm ) THEN 158 171 CALL iom_put( "avs" , fsavs(:,:,:) ) ! S vert. eddy diff. coef. … … 250 263 ! 251 264 CALL wrk_alloc( jpi , jpj , zw2d ) 252 IF ( ln_traldf_gdia ) call wrk_alloc( jpi , jpj , jpk , zw3d )265 IF ( ln_traldf_gdia .OR. lk_vvl ) call wrk_alloc( jpi , jpj , jpk , zw3d ) 253 266 ! 254 267 ! Output the initial state and forcings … … 395 408 CALL histdef( nid_T, "vosaline", "Salinity" , "PSU" , & ! sn 396 409 & jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout ) 410 IF( lk_vvl ) THEN 411 CALL histdef( nid_T, "vovvle3t", "Level thickness" , "m" ,& ! e3t_n 412 & jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout ) 413 CALL histdef( nid_T, "vovvldep", "T point depth" , "m" ,& ! e3t_n 414 & jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout ) 415 CALL histdef( nid_T, "vovvldef", "Squared level deformation" , "%^2" ,& ! e3t_n 416 & jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout ) 417 ENDIF 397 418 ! !!! nid_T : 2D 398 419 CALL histdef( nid_T, "sosstsst", "Sea Surface temperature" , "C" , & ! sst … … 406 427 CALL histdef( nid_T, "sosfldow", "downward salt flux" , "PSU/m2/s", & ! sfx 407 428 & jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout ) 408 #if ! defined key_vvl 409 CALL histdef( nid_T, "sosst_cd", "Concentration/Dilution term on temperature"& ! emp * tsn(:,:,1,jp_tem)429 IF( .NOT. lk_vvl ) THEN 430 CALL histdef( nid_T, "sosst_cd", "Concentration/Dilution term on temperature" & ! emp * tsn(:,:,1,jp_tem) 410 431 & , "KgC/m2/s", & ! sosst_cd 411 & jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )412 CALL histdef( nid_T, "sosss_cd", "Concentration/Dilution term on salinity"& ! emp * tsn(:,:,1,jp_sal)432 & jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout ) 433 CALL histdef( nid_T, "sosss_cd", "Concentration/Dilution term on salinity" & ! emp * tsn(:,:,1,jp_sal) 413 434 & , "KgPSU/m2/s",& ! sosss_cd 414 & jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )415 #endif 435 & jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout ) 436 ENDIF 416 437 CALL histdef( nid_T, "sohefldo", "Net Downward Heat Flux" , "W/m2" , & ! qns + qsr 417 438 & jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout ) … … 585 606 ! --------------------- 586 607 587 ! ndex(1) est utilise ssi l'avant dernier argument est diff ferent de608 ! ndex(1) est utilise ssi l'avant dernier argument est different de 588 609 ! la taille du tableau en sortie. Dans ce cas , l'avant dernier argument 589 610 ! donne le nombre d'elements, et ndex la liste des indices a sortir … … 595 616 596 617 ! Write fields on T grid 597 CALL histwrite( nid_T, "votemper", it, tsn(:,:,:,jp_tem), ndim_T , ndex_T ) ! temperature 598 CALL histwrite( nid_T, "vosaline", it, tsn(:,:,:,jp_sal), ndim_T , ndex_T ) ! salinity 599 CALL histwrite( nid_T, "sosstsst", it, tsn(:,:,1,jp_tem), ndim_hT, ndex_hT ) ! sea surface temperature 600 CALL histwrite( nid_T, "sosaline", it, tsn(:,:,1,jp_sal), ndim_hT, ndex_hT ) ! sea surface salinity 618 IF( lk_vvl ) THEN 619 CALL histwrite( nid_T, "votemper", it, tsn(:,:,:,jp_tem) * fse3t_n(:,:,:) , ndim_T , ndex_T ) ! heat content 620 CALL histwrite( nid_T, "vosaline", it, tsn(:,:,:,jp_sal) * fse3t_n(:,:,:) , ndim_T , ndex_T ) ! salt content 621 CALL histwrite( nid_T, "sosstsst", it, tsn(:,:,1,jp_tem) * fse3t_n(:,:,1) , ndim_hT, ndex_hT ) ! sea surface heat content 622 CALL histwrite( nid_T, "sosaline", it, tsn(:,:,1,jp_sal) * fse3t_n(:,:,1) , ndim_hT, ndex_hT ) ! sea surface salinity content 623 ELSE 624 CALL histwrite( nid_T, "votemper", it, tsn(:,:,:,jp_tem) , ndim_T , ndex_T ) ! temperature 625 CALL histwrite( nid_T, "vosaline", it, tsn(:,:,:,jp_sal) , ndim_T , ndex_T ) ! salinity 626 CALL histwrite( nid_T, "sosstsst", it, tsn(:,:,1,jp_tem) , ndim_hT, ndex_hT ) ! sea surface temperature 627 CALL histwrite( nid_T, "sosaline", it, tsn(:,:,1,jp_sal) , ndim_hT, ndex_hT ) ! sea surface salinity 628 629 ENDIF 630 IF( lk_vvl ) THEN 631 z3d(:,:,:) = ( ( fse3t_n(:,:,:) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2 632 CALL histwrite( nid_T, "vovvle3t", it, fse3t_n (:,:,:) , ndim_T , ndex_T ) ! level thickness 633 CALL histwrite( nid_T, "vovvldep", it, fsdept_n(:,:,:) , ndim_T , ndex_T ) ! t-point depth 634 CALL histwrite( nid_T, "vovvldef", it, z3d , ndim_T , ndex_T ) ! level thickness deformation 635 ENDIF 601 636 CALL histwrite( nid_T, "sossheig", it, sshn , ndim_hT, ndex_hT ) ! sea surface height 602 637 CALL histwrite( nid_T, "sowaflup", it, ( emp-rnf ) , ndim_hT, ndex_hT ) ! upward water flux … … 604 639 ! (includes virtual salt flux beneath ice 605 640 ! in linear free surface case) 606 #if ! defined key_vvl 607 zw2d(:,:) = emp (:,:) * tsn(:,:,1,jp_tem)608 CALL histwrite( nid_T, "sosst_cd", it, zw2d, ndim_hT, ndex_hT )! c/d term on sst609 zw2d(:,:) = emp (:,:) * tsn(:,:,1,jp_sal)610 CALL histwrite( nid_T, "sosss_cd", it, zw2d, ndim_hT, ndex_hT )! c/d term on sss611 #endif 641 IF( .NOT. lk_vvl ) THEN 642 zw2d(:,:) = emp (:,:) * tsn(:,:,1,jp_tem) 643 CALL histwrite( nid_T, "sosst_cd", it, zw2d, ndim_hT, ndex_hT ) ! c/d term on sst 644 zw2d(:,:) = emp (:,:) * tsn(:,:,1,jp_sal) 645 CALL histwrite( nid_T, "sosss_cd", it, zw2d, ndim_hT, ndex_hT ) ! c/d term on sss 646 ENDIF 612 647 CALL histwrite( nid_T, "sohefldo", it, qns + qsr , ndim_hT, ndex_hT ) ! total heat flux 613 648 CALL histwrite( nid_T, "soshfldo", it, qsr , ndim_hT, ndex_hT ) ! solar heat flux … … 750 785 ! 751 786 CALL wrk_dealloc( jpi , jpj , zw2d ) 752 IF ( ln_traldf_gdia ) call wrk_dealloc( jpi , jpj , jpk , zw3d )787 IF ( ln_traldf_gdia .OR. lk_vvl ) call wrk_dealloc( jpi , jpj , jpk , zw3d ) 753 788 ! 754 789 IF( nn_timing == 1 ) CALL timing_stop('dia_wri') … … 839 874 CALL histdef( id_i, "sometauy", "Meridional Wind Stress", "N/m2" , & ! j-wind stress 840 875 & jpi, jpj, nh_i, 1 , 1, 1 , -99 , 32, clop, zsto, zout ) 876 IF( lk_vvl ) THEN 877 CALL histdef( id_i, "vovvldep", "T point depth" , "m" , & ! t-point depth 878 & jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout ) 879 END IF 841 880 842 881 #if defined key_lim2 -
branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90
r3862 r3865 134 134 !! All coordinates 135 135 !! --------------- 136 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: gdep3w !: depth of T-points (sum of e3w) (m)137 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: gdept , gdepw !: analytical depth at T-Wpoints (m)138 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e3v , e3f !: analytical vertical scale factors at V--F139 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e3t , e3u !: T--Upoints (m)140 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e3vw !: analytical vertical scale factors at VW--141 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e3w , e3uw !: W--UWpoints (m)136 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: gdep3w_0 !: depth of t-points (sum of e3w) (m) 137 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: gdept_0, gdepw_0 !: analytical (time invariant) depth at t-w points (m) 138 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e3v_0 , e3f_0 !: analytical (time invariant) vertical scale factors at v-f 139 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e3t_0 , e3u_0 !: t-u points (m) 140 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e3vw_0 !: analytical (time invariant) vertical scale factors at vw 141 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e3w_0 , e3uw_0 !: w-uw points (m) 142 142 #if defined key_vvl 143 143 LOGICAL, PUBLIC, PARAMETER :: lk_vvl = .TRUE. !: variable grid flag … … 145 145 !! All coordinates 146 146 !! --------------- 147 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: gdep3w_1 !: depth of T-points (sum of e3w) (m) 148 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: gdept_1, gdepw_1 !: analytical depth at T-W points (m) 149 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e3v_1 , e3f_1 !: analytical vertical scale factors at V--F 150 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e3t_1 , e3u_1 !: T--U points (m) 151 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e3vw_1 !: analytical vertical scale factors at VW-- 152 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e3w_1 , e3uw_1 !: W--UW points (m) 153 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e3t_b !: before - - - - T points (m) 154 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e3u_b , e3v_b !: - - - - - U--V points (m) 147 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: gdep3w_n !: now depth of T-points (sum of e3w) (m) 148 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: gdept_n, gdepw_n !: now depth at T-W points (m) 149 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e3t_n !: now vertical scale factors at t point (m) 150 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e3u_n , e3v_n !: - - - - u --v points (m) 151 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e3w_n , e3f_n !: - - - - w --f points (m) 152 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e3uw_n , e3vw_n !: - - - - uw--vw points (m) 153 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e3t_b !: before - - - - t points (m) 154 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e3u_b , e3v_b !: - - - - - u --v points (m) 155 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e3uw_b , e3vw_b !: - - - - - uw--vw points (m) 156 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e3t_a !: after - - - - t point (m) 157 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e3u_a , e3v_a !: - - - - - u --v points (m) 155 158 #else 156 159 LOGICAL, PUBLIC, PARAMETER :: lk_vvl = .FALSE. !: fixed grid flag 157 160 #endif 158 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:), TARGET :: hur , hvr !: inverse of u and v-points ocean depth (1/m) 159 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hu , hv !: depth at u- and v-points (meters) 160 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hu_0 , hv_0 !: refernce depth at u- and v-points (meters) 161 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:), TARGET :: hur , hvr !: inverse of u and v-points ocean depth (1/m) 162 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hu , hv !: depth at u- and v-points (meters) 163 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ht_0 !: reference depth at t- points (meters) 164 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hu_0 , hv_0 !: reference depth at u- and v-points (meters) 165 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: re2u_e1u !: scale factor coeffs at u points (e2u/e1u) 166 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: re1v_e2v !: scale factor coeffs at v points (e1v/e2v) 167 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: e12t , r1_e12t !: horizontal cell surface and inverse at t points 168 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: e12u , r1_e12u !: horizontal cell surface and inverse at u points 169 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: e12v , r1_e12v !: horizontal cell surface and inverse at v points 170 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: e12f , r1_e12f !: horizontal cell surface and inverse at f points 161 171 162 172 INTEGER, PUBLIC :: nla10 !: deepest W level Above ~10m (nlb10 - 1) … … 165 175 !! z-coordinate with full steps (also used in the other cases as reference z-coordinate) 166 176 !! =-----------------====------ 167 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: gdept_1d, gdepw_1d !: reference depth of t- and w-points (m)168 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: e3t_1d , e3w_1d !: reference vertical scale factors at T- and W-pts (m)169 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: e3tp , e3wp !: ocean bottom level thickness at T and W points177 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: gdept_1d, gdepw_1d !: reference depth of t- and w-points (m) 178 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: e3t_1d , e3w_1d !: reference vertical scale factors at T- and W-pts (m) 179 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: e3tp , e3wp !: ocean bottom level thickness at T and W points 170 180 171 181 !! s-coordinate and hybrid z-s-coordinate 172 182 !! =----------------======--------------- 173 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: gsigt, gsigw!: model level depth coefficient at t-, w-levels (analytic)174 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: gsi3w!: model level depth coefficient at w-level (sum of gsigw)175 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: esigt, esigw!: vertical scale factor coef. at t-, w-levels176 177 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hbatv , hbatf !: ocean depth at the vertical of V--F178 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hbatt , hbatu !: T--Upoints (m)179 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: scosrf, scobot !: ocean surface and bottom topographies180 ! ! (if deviating from coordinate surfaces in HYBRID)181 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hifv , hiff !: interface depth between stretching at V--F182 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hift , hifu !: and quasi-uniform spacing T--Upoints (m)183 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: rx1 !: Maximum grid stiffness ratio183 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: gsigt, gsigw !: model level depth coefficient at t-, w-levels (analytic) 184 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: gsi3w !: model level depth coefficient at w-level (sum of gsigw) 185 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: esigt, esigw !: vertical scale factor coef. at t-, w-levels 186 187 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hbatv , hbatf !: ocean depth at the vertical of v--f 188 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hbatt , hbatu !: t--u points (m) 189 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: scosrf, scobot !: ocean surface and bottom topographies 190 ! ! (if deviating from coordinate surfaces in HYBRID) 191 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hifv , hiff !: interface depth between stretching at v--f 192 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hift , hifu !: and quasi-uniform spacing t--u points (m) 193 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: rx1 !: Maximum grid stiffness ratio 184 194 185 195 !!---------------------------------------------------------------------- 186 196 !! masks, bathymetry 187 197 !! --------------------------------------------------------------------- 188 INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: mbathy !: number of ocean level (=0, 1, ... , jpk-1)189 INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: mbkt !: vertical index of the bottom last T- ocean level190 INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: mbku, mbkv !: vertical index of the bottom last U- and W- ocean level191 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: bathy !: ocean depth (meters)192 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: tmask_i !: interior domain T-point mask193 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: bmask !: land/ocean mask of barotropic stream function198 INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: mbathy !: number of ocean level (=0, 1, ... , jpk-1) 199 INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: mbkt !: vertical index of the bottom last T- ocean level 200 INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: mbku, mbkv !: vertical index of the bottom last U- and W- ocean level 201 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: bathy !: ocean depth (meters) 202 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: tmask_i !: interior domain T-point mask 203 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: bmask !: land/ocean mask of barotropic stream function 194 204 195 205 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: tmask, umask, vmask, fmask !: land/ocean mask at T-, U-, V- and F-pts 196 206 197 REAL(wp), PUBLIC, DIMENSION(jpiglo) :: tpol, fpol !: north fold mask (jperio= 3 or 4)207 REAL(wp), PUBLIC, DIMENSION(jpiglo) :: tpol, fpol !: north fold mask (jperio= 3 or 4) 198 208 199 209 #if defined key_noslip_accurate 200 INTEGER, PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,: ) :: npcoa!: ???201 INTEGER, PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: nicoa, njcoa!: ???210 INTEGER, PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,: ) :: npcoa !: ??? 211 INTEGER, PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: nicoa, njcoa !: ??? 202 212 #endif 203 213 … … 279 289 & glamf(jpi,jpj) , gphif(jpi,jpj) , e1f(jpi,jpj) , e2f(jpi,jpj) , ff (jpi,jpj) , STAT=ierr(3) ) 280 290 ! 281 ALLOCATE( gdep3w (jpi,jpj,jpk) , e3v(jpi,jpj,jpk) , e3f(jpi,jpj,jpk) , &282 & gdept (jpi,jpj,jpk) , e3t(jpi,jpj,jpk) , e3u(jpi,jpj,jpk) , &283 & gdepw (jpi,jpj,jpk) , e3w(jpi,jpj,jpk) , e3vw(jpi,jpj,jpk) , e3uw(jpi,jpj,jpk) , STAT=ierr(4) )291 ALLOCATE( gdep3w_0(jpi,jpj,jpk) , e3v_0(jpi,jpj,jpk) , e3f_0 (jpi,jpj,jpk) , & 292 & gdept_0 (jpi,jpj,jpk) , e3t_0(jpi,jpj,jpk) , e3u_0 (jpi,jpj,jpk) , & 293 & gdepw_0 (jpi,jpj,jpk) , e3w_0(jpi,jpj,jpk) , e3vw_0(jpi,jpj,jpk) , e3uw_0(jpi,jpj,jpk) , STAT=ierr(4) ) 284 294 ! 285 295 #if defined key_vvl 286 ALLOCATE( gdep3w_1(jpi,jpj,jpk) , e3v_1(jpi,jpj,jpk) , e3f_1 (jpi,jpj,jpk) , & 287 & gdept_1 (jpi,jpj,jpk) , e3t_1(jpi,jpj,jpk) , e3u_1 (jpi,jpj,jpk) , & 288 & gdepw_1 (jpi,jpj,jpk) , e3w_1(jpi,jpj,jpk) , e3vw_1(jpi,jpj,jpk) , e3uw_1(jpi,jpj,jpk) , & 289 & e3t_b (jpi,jpj,jpk) , e3u_b(jpi,jpj,jpk) , e3v_b (jpi,jpj,jpk) , STAT=ierr(5) ) 290 #endif 291 ! 292 ALLOCATE( hu(jpi,jpj) , hur(jpi,jpj) , hu_0(jpi,jpj) , & 293 & hv(jpi,jpj) , hvr(jpi,jpj) , hv_0(jpi,jpj) , STAT=ierr(6) ) 294 ! 295 ALLOCATE( gdept_1d(jpk) , gdepw_1d(jpk) , & 296 & e3t_1d (jpk) , e3w_1d (jpk) , e3tp (jpi,jpj), e3wp(jpi,jpj) , & 297 & gsigt (jpk) , gsigw (jpk) , gsi3w(jpk) , & 298 & esigt (jpk) , esigw (jpk) , STAT=ierr(7) ) 296 ALLOCATE( gdep3w_n(jpi,jpj,jpk) , e3t_n (jpi,jpj,jpk) , e3u_n (jpi,jpj,jpk) , & 297 & gdept_n (jpi,jpj,jpk) , e3v_n (jpi,jpj,jpk) , e3w_n (jpi,jpj,jpk) , & 298 & gdepw_n (jpi,jpj,jpk) , e3f_n (jpi,jpj,jpk) , e3vw_n(jpi,jpj,jpk) , e3uw_n(jpi,jpj,jpk) , & 299 & e3t_b (jpi,jpj,jpk) , e3u_b (jpi,jpj,jpk) , e3v_b (jpi,jpj,jpk) , & 300 & e3uw_b (jpi,jpj,jpk) , e3vw_b(jpi,jpj,jpk) , & 301 & e3t_a (jpi,jpj,jpk) , e3u_a (jpi,jpj,jpk) , e3v_a (jpi,jpj,jpk) , STAT=ierr(5) ) 302 #endif 303 ! 304 ALLOCATE( hu (jpi,jpj) , hur (jpi,jpj) , hu_0(jpi,jpj) , ht_0 (jpi,jpj) , & 305 & hv (jpi,jpj) , hvr (jpi,jpj) , hv_0(jpi,jpj) , & 306 & re2u_e1u(jpi,jpj) , re1v_e2v(jpi,jpj) , & 307 & e12t (jpi,jpj) , r1_e12t (jpi,jpj) , & 308 & e12u (jpi,jpj) , r1_e12u (jpi,jpj) , & 309 & e12v (jpi,jpj) , r1_e12v (jpi,jpj) , & 310 & e12f (jpi,jpj) , r1_e12f (jpi,jpj) , STAT=ierr(6) ) 311 ! 312 ALLOCATE( gdept_1d(jpk) , gdepw_1d(jpk) , & 313 & e3t_1d (jpk) , e3w_1d (jpk) , e3tp (jpi,jpj), e3wp(jpi,jpj) , & 314 & gsigt (jpk) , gsigw (jpk) , gsi3w(jpk) , & 315 & esigt (jpk) , esigw (jpk) , STAT=ierr(7) ) 299 316 ! 300 317 ALLOCATE( hbatv (jpi,jpj) , hbatf (jpi,jpj) , & -
branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90
r3764 r3865 86 86 CALL dom_msk ! Masks 87 87 IF( ln_sco ) CALL dom_stiff ! Maximum stiffness ratio/hydrostatic consistency 88 IF( lk_vvl ) CALL dom_vvl! Vertical variable mesh88 IF( lk_vvl ) CALL dom_vvl_init ! Vertical variable mesh 89 89 ! 90 90 IF( lk_c1d ) CALL cor_c1d ! 1D configuration: Coriolis set at T-point 91 ! 92 ! - ML - Used in dom_vvl_sf_nxt and lateral diffusion routines 93 ! but could be usefull in many other routines 94 e12t (:,:) = e1t(:,:) * e2t(:,:) 95 e12u (:,:) = e1u(:,:) * e2u(:,:) 96 e12v (:,:) = e1v(:,:) * e2v(:,:) 97 e12f (:,:) = e1f(:,:) * e2f(:,:) 98 r1_e12t (:,:) = 1._wp / e12t(:,:) 99 r1_e12u (:,:) = 1._wp / e12u(:,:) 100 r1_e12v (:,:) = 1._wp / e12v(:,:) 101 r1_e12f (:,:) = 1._wp / e12f(:,:) 102 re2u_e1u(:,:) = e2u(:,:) / e1u(:,:) 103 re1v_e2v(:,:) = e1v(:,:) / e2v(:,:) 91 104 ! 92 105 hu(:,:) = 0._wp ! Ocean depth at U- and V-points 93 106 hv(:,:) = 0._wp 94 107 DO jk = 1, jpk 95 hu(:,:) = hu(:,:) + fse3u (:,:,jk) * umask(:,:,jk)96 hv(:,:) = hv(:,:) + fse3v (:,:,jk) * vmask(:,:,jk)108 hu(:,:) = hu(:,:) + fse3u_n(:,:,jk) * umask(:,:,jk) 109 hv(:,:) = hv(:,:) + fse3v_n(:,:,jk) * vmask(:,:,jk) 97 110 END DO 98 111 ! ! Inverse of the local depth … … 345 358 DO jj = 2, jpjm1 346 359 DO jk = 1, jpkm1 347 zr1(1) = umask(ji-1,jj ,jk) *abs( (gdepw (ji ,jj ,jk )-gdepw(ji-1,jj ,jk ) &348 & +gdepw (ji ,jj ,jk+1)-gdepw(ji-1,jj ,jk+1)) &349 & /(gdepw (ji ,jj ,jk )+gdepw(ji-1,jj ,jk ) &350 & -gdepw (ji ,jj ,jk+1)-gdepw(ji-1,jj ,jk+1) + rsmall) )351 zr1(2) = umask(ji ,jj ,jk) *abs( (gdepw (ji+1,jj ,jk )-gdepw(ji ,jj ,jk ) &352 & +gdepw (ji+1,jj ,jk+1)-gdepw(ji ,jj ,jk+1)) &353 & /(gdepw (ji+1,jj ,jk )+gdepw(ji ,jj ,jk ) &354 & -gdepw (ji+1,jj ,jk+1)-gdepw(ji ,jj ,jk+1) + rsmall) )355 zr1(3) = vmask(ji ,jj ,jk) *abs( (gdepw (ji ,jj+1,jk )-gdepw(ji ,jj ,jk ) &356 & +gdepw (ji ,jj+1,jk+1)-gdepw(ji ,jj ,jk+1)) &357 & /(gdepw (ji ,jj+1,jk )+gdepw(ji ,jj ,jk ) &358 & -gdepw (ji ,jj+1,jk+1)-gdepw(ji ,jj ,jk+1) + rsmall) )359 zr1(4) = vmask(ji ,jj-1,jk) *abs( (gdepw (ji ,jj ,jk )-gdepw(ji ,jj-1,jk ) &360 & +gdepw (ji ,jj ,jk+1)-gdepw(ji ,jj-1,jk+1)) &361 & /(gdepw (ji ,jj ,jk )+gdepw(ji ,jj-1,jk ) &362 & -gdepw (ji, jj ,jk+1)-gdepw(ji ,jj-1,jk+1) + rsmall) )360 zr1(1) = umask(ji-1,jj ,jk) *abs( (gdepw_0(ji ,jj ,jk )-gdepw_0(ji-1,jj ,jk ) & 361 & +gdepw_0(ji ,jj ,jk+1)-gdepw_0(ji-1,jj ,jk+1)) & 362 & /(gdepw_0(ji ,jj ,jk )+gdepw_0(ji-1,jj ,jk ) & 363 & -gdepw_0(ji ,jj ,jk+1)-gdepw_0(ji-1,jj ,jk+1) + rsmall) ) 364 zr1(2) = umask(ji ,jj ,jk) *abs( (gdepw_0(ji+1,jj ,jk )-gdepw_0(ji ,jj ,jk ) & 365 & +gdepw_0(ji+1,jj ,jk+1)-gdepw_0(ji ,jj ,jk+1)) & 366 & /(gdepw_0(ji+1,jj ,jk )+gdepw_0(ji ,jj ,jk ) & 367 & -gdepw_0(ji+1,jj ,jk+1)-gdepw_0(ji ,jj ,jk+1) + rsmall) ) 368 zr1(3) = vmask(ji ,jj ,jk) *abs( (gdepw_0(ji ,jj+1,jk )-gdepw_0(ji ,jj ,jk ) & 369 & +gdepw_0(ji ,jj+1,jk+1)-gdepw_0(ji ,jj ,jk+1)) & 370 & /(gdepw_0(ji ,jj+1,jk )+gdepw_0(ji ,jj ,jk ) & 371 & -gdepw_0(ji ,jj+1,jk+1)-gdepw_0(ji ,jj ,jk+1) + rsmall) ) 372 zr1(4) = vmask(ji ,jj-1,jk) *abs( (gdepw_0(ji ,jj ,jk )-gdepw_0(ji ,jj-1,jk ) & 373 & +gdepw_0(ji ,jj ,jk+1)-gdepw_0(ji ,jj-1,jk+1)) & 374 & /(gdepw_0(ji ,jj ,jk )+gdepw_0(ji ,jj-1,jk ) & 375 & -gdepw_0(ji, jj ,jk+1)-gdepw_0(ji ,jj-1,jk+1) + rsmall) ) 363 376 zrxmax = MAXVAL(zr1(1:4)) 364 377 rx1(ji,jj) = MAX(rx1(ji,jj), zrxmax) -
branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90
r3294 r3865 6 6 !! History : 2.0 ! 2006-06 (B. Levier, L. Marie) original code 7 7 !! 3.1 ! 2009-02 (G. Madec, M. Leclair, R. Benshila) pure z* coordinate 8 !! ----------------------------------------------------------------------9 #if defined key_vvl 8 !! 3.3 ! 2011-10 (M. Leclair) totally rewrote domvvl: 9 !! vvl option includes z_star and z_tilde coordinates 10 10 !!---------------------------------------------------------------------- 11 11 !! 'key_vvl' variable volume 12 12 !!---------------------------------------------------------------------- 13 !! dom_vvl : defined coefficients to distribute ssh on each layers14 13 !!---------------------------------------------------------------------- 14 !! dom_vvl_init : define initial vertical scale factors, depths and column thickness 15 !! dom_vvl_sf_nxt : Compute next vertical scale factors 16 !! dom_vvl_sf_swp : Swap vertical scale factors and update the vertical grid 17 !! dom_vvl_interpol : Interpolate vertical scale factors from one grid point to another 18 !! dom_vvl_rst : read/write restart file 19 !! dom_vvl_ctl : Check the vvl options 20 !!---------------------------------------------------------------------- 21 !! * Modules used 15 22 USE oce ! ocean dynamics and tracers 16 23 USE dom_oce ! ocean space and time domain 17 USE sbc_oce ! surface boundary condition: ocean 18 USE phycst ! physical constants 24 USE sbc_oce ! ocean surface boundary condition 19 25 USE in_out_manager ! I/O manager 26 USE iom ! I/O manager library 27 USE restart ! ocean restart 20 28 USE lib_mpp ! distributed memory computing library 21 29 USE lbclnk ! ocean lateral boundary conditions (or mpp link) … … 26 34 PRIVATE 27 35 28 PUBLIC dom_vvl ! called by domain.F90 29 PUBLIC dom_vvl_2 ! called by domain.F90 30 PUBLIC dom_vvl_alloc ! called by nemogcm.F90 31 32 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: mut , muu , muv , muf !: 1/H_0 at t-,u-,v-,f-points 36 !! * Routine accessibility 37 PUBLIC dom_vvl_init ! called by domain.F90 38 PUBLIC dom_vvl_sf_nxt ! called by step.F90 39 PUBLIC dom_vvl_sf_swp ! called by step.F90 40 PUBLIC dom_vvl_interpol ! called by dynnxt.F90 41 42 !!* Namelist nam_vvl 43 LOGICAL , PUBLIC :: ln_vvl_zstar = .FALSE. ! zstar vertical coordinate 44 LOGICAL , PUBLIC :: ln_vvl_ztilde = .FALSE. ! ztilde vertical coordinate 45 LOGICAL , PUBLIC :: ln_vvl_layer = .FALSE. ! level vertical coordinate 46 LOGICAL , PUBLIC :: ln_vvl_kepe = .FALSE. ! kinetic/potential energy transfer 47 ! ! conservation: not used yet 48 REAL(wp) :: rn_ahe3 = 0.e0 ! thickness diffusion coefficient 49 LOGICAL , PUBLIC :: ln_vvl_dbg = .FALSE. ! debug control prints 50 51 !! * Module variables 52 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: un_td, vn_td ! thickness diffusion transport 53 REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: hdiv_lf ! low frequency part of hz divergence 54 REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e3t_t_b, e3t_t_n, e3t_t_a ! baroclinic scale factors 55 REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:) :: frq_rst_e3t ! retoring period for scale factors 56 REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:) :: frq_rst_hdv ! retoring period for low freq. divergence 33 57 34 58 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) :: r2dt ! vertical profile time-step, = 2 rdttra … … 39 63 # include "vectopt_loop_substitute.h90" 40 64 !!---------------------------------------------------------------------- 41 !! NEMO/OPA 4.0 , NEMO Consortium (2011)65 !! NEMO/OPA 3.3 , NEMO-Consortium (2010) 42 66 !! $Id$ 43 67 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 44 68 !!---------------------------------------------------------------------- 45 CONTAINS 69 70 CONTAINS 46 71 47 72 INTEGER FUNCTION dom_vvl_alloc() 48 73 !!---------------------------------------------------------------------- 49 !! *** ROUTINE dom_vvl_alloc *** 50 !!---------------------------------------------------------------------- 51 ! 52 ALLOCATE( mut (jpi,jpj,jpk) , muu (jpi,jpj,jpk) , muv (jpi,jpj,jpk) , muf (jpi,jpj,jpk) , & 53 & r2dt (jpk) , STAT=dom_vvl_alloc ) 54 ! 55 IF( lk_mpp ) CALL mpp_sum ( dom_vvl_alloc ) 56 IF( dom_vvl_alloc /= 0 ) CALL ctl_warn('dom_vvl_alloc: failed to allocate arrays') 57 ! 74 !! *** FUNCTION dom_vvl_alloc *** 75 !!---------------------------------------------------------------------- 76 IF( ln_vvl_zstar ) dom_vvl_alloc = 0 77 IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN 78 ALLOCATE( e3t_t_b(jpi,jpj,jpk) , e3t_t_n(jpi,jpj,jpk) , e3t_t_a(jpi,jpj,jpk) , & 79 & un_td (jpi,jpj,jpk) , vn_td (jpi,jpj,jpk) , STAT = dom_vvl_alloc ) 80 IF( lk_mpp ) CALL mpp_sum ( dom_vvl_alloc ) 81 IF( dom_vvl_alloc /= 0 ) CALL ctl_warn('dom_vvl_alloc: failed to allocate arrays') 82 ENDIF 83 IF( ln_vvl_ztilde ) THEN 84 ALLOCATE( frq_rst_e3t(jpi,jpj) , frq_rst_hdv(jpi,jpj) , hdiv_lf(jpi,jpj,jpk) , STAT= dom_vvl_alloc ) 85 IF( lk_mpp ) CALL mpp_sum ( dom_vvl_alloc ) 86 IF( dom_vvl_alloc /= 0 ) CALL ctl_warn('dom_vvl_alloc: failed to allocate arrays') 87 ENDIF 88 58 89 END FUNCTION dom_vvl_alloc 59 90 60 91 61 SUBROUTINE dom_vvl 62 !!---------------------------------------------------------------------- 63 !! *** ROUTINE dom_vvl ***92 SUBROUTINE dom_vvl_init 93 !!---------------------------------------------------------------------- 94 !! *** ROUTINE dom_vvl_init *** 64 95 !! 65 !! ** Purpose : compute mu coefficients at t-, u-, v- and f-points to 66 !! spread ssh over the whole water column (scale factors) 67 !! set the before and now ssh at u- and v-points 68 !! (also f-point in now case) 69 !!---------------------------------------------------------------------- 70 ! 71 INTEGER :: ji, jj, jk ! dummy loop indices 72 REAL(wp) :: zcoefu, zcoefv , zcoeff ! local scalars 73 REAL(wp) :: zvt , zvt_ip1, zvt_jp1, zvt_ip1jp1 ! - - 74 REAL(wp), POINTER, DIMENSION(:,:) :: zee_t, zee_u, zee_v, zee_f ! 2D workspace 75 !!---------------------------------------------------------------------- 76 ! 77 IF( nn_timing == 1 ) CALL timing_start('dom_vvl') 78 ! 79 CALL wrk_alloc( jpi, jpj, zee_t, zee_u, zee_v, zee_f ) 80 ! 81 IF(lwp) THEN 82 WRITE(numout,*) 83 WRITE(numout,*) 'dom_vvl : Variable volume initialization' 84 WRITE(numout,*) '~~~~~~~~ compute coef. used to spread ssh over each layers' 85 ENDIF 86 87 IF( dom_vvl_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'dom_vvl : unable to allocate arrays' ) 88 89 fsdept(:,:,:) = gdept (:,:,:) 90 fsdepw(:,:,:) = gdepw (:,:,:) 91 fsde3w(:,:,:) = gdep3w(:,:,:) 92 fse3t (:,:,:) = e3t (:,:,:) 93 fse3u (:,:,:) = e3u (:,:,:) 94 fse3v (:,:,:) = e3v (:,:,:) 95 fse3f (:,:,:) = e3f (:,:,:) 96 fse3w (:,:,:) = e3w (:,:,:) 97 fse3uw(:,:,:) = e3uw (:,:,:) 98 fse3vw(:,:,:) = e3vw (:,:,:) 99 100 ! !== mu computation ==! 101 zee_t(:,:) = fse3t_0(:,:,1) ! Lower bound : thickness of the first model level 102 zee_u(:,:) = fse3u_0(:,:,1) 103 zee_v(:,:) = fse3v_0(:,:,1) 104 zee_f(:,:) = fse3f_0(:,:,1) 105 DO jk = 2, jpkm1 ! Sum of the masked vertical scale factors 106 zee_t(:,:) = zee_t(:,:) + fse3t_0(:,:,jk) * tmask(:,:,jk) 107 zee_u(:,:) = zee_u(:,:) + fse3u_0(:,:,jk) * umask(:,:,jk) 108 zee_v(:,:) = zee_v(:,:) + fse3v_0(:,:,jk) * vmask(:,:,jk) 109 DO jj = 1, jpjm1 ! f-point : fmask=shlat at coasts, use the product of umask 110 zee_f(:,jj) = zee_f(:,jj) + fse3f_0(:,jj,jk) * umask(:,jj,jk) * umask(:,jj+1,jk) 111 END DO 112 END DO 113 ! ! Compute and mask the inverse of the local depth at T, U, V and F points 114 zee_t(:,:) = 1._wp / zee_t(:,:) * tmask(:,:,1) 115 zee_u(:,:) = 1._wp / zee_u(:,:) * umask(:,:,1) 116 zee_v(:,:) = 1._wp / zee_v(:,:) * vmask(:,:,1) 117 DO jj = 1, jpjm1 ! f-point case fmask cannot be used 118 zee_f(:,jj) = 1._wp / zee_f(:,jj) * umask(:,jj,1) * umask(:,jj+1,1) 96 !! ** Purpose : Initialization of all scale factors, depths 97 !! and water column heights 98 !! 99 !! ** Method : - use restart file and/or initialize 100 !! - interpolate scale factors 101 !! 102 !! ** Action : - fse3t_(n/b) and e3t_t_(n/b) 103 !! - Regrid: fse3(u/v)_n 104 !! fse3(u/v)_b 105 !! fse3w_n 106 !! fse3(u/v)w_b 107 !! fse3(u/v)w_n 108 !! fsdept_n, fsdepw_n and fsde3w_n 109 !! - h(t/u/v)_0 110 !! - frq_rst_e3t and frq_rst_hdv 111 !! 112 !! Reference : Leclair, M., and G. Madec, 2011, Ocean Modelling. 113 !!---------------------------------------------------------------------- 114 USE phycst, ONLY : rpi 115 !! * Local declarations 116 INTEGER :: jk 117 !!---------------------------------------------------------------------- 118 IF( nn_timing == 1 ) CALL timing_start('dom_vvl_init') 119 120 IF(lwp) WRITE(numout,*) 121 IF(lwp) WRITE(numout,*) 'dom_vvl_init : Variable volume activated' 122 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~' 123 124 ! choose vertical coordinate (z_star, z_tilde or layer) 125 ! ========================== 126 CALL dom_vvl_ctl 127 128 ! Allocate module arrays 129 ! ====================== 130 IF( dom_vvl_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'dom_vvl_init : unable to allocate arrays' ) 131 132 ! Read or initialize fse3t_(b/n), e3t_t_(b/n) and hdiv_lf (and e3t_a(jpk)) 133 ! ======================================================================== 134 CALL dom_vvl_rst( nit000, 'READ' ) 135 fse3t_a(:,:,jpk) = e3t_0(:,:,jpk) 136 137 ! Reconstruction of all vertical scale factors at now and before time steps 138 ! ========================================================================= 139 ! Horizontal scale factor interpolations 140 ! -------------------------------------- 141 CALL dom_vvl_interpol( fse3t_b(:,:,:), fse3u_b(:,:,:), 'U' ) 142 CALL dom_vvl_interpol( fse3t_b(:,:,:), fse3v_b(:,:,:), 'V' ) 143 CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3u_n(:,:,:), 'U' ) 144 CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3v_n(:,:,:), 'V' ) 145 CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3f_n(:,:,:), 'F' ) 146 ! Vertical scale factor interpolations 147 ! ------------------------------------ 148 CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3w_n (:,:,:), 'W' ) 149 CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3uw_n(:,:,:), 'UW' ) 150 CALL dom_vvl_interpol( fse3v_n(:,:,:), fse3vw_n(:,:,:), 'VW' ) 151 CALL dom_vvl_interpol( fse3u_b(:,:,:), fse3uw_b(:,:,:), 'UW' ) 152 CALL dom_vvl_interpol( fse3v_b(:,:,:), fse3vw_b(:,:,:), 'VW' ) 153 ! t- and w- points depth 154 ! ---------------------- 155 fsdept_n(:,:,1) = 0.5 * fse3w_n(:,:,1) 156 fsdepw_n(:,:,1) = 0.e0 157 fsde3w_n(:,:,1) = fsdept_n(:,:,1) - sshn(:,:) 158 DO jk = 2, jpk 159 fsdept_n(:,:,jk) = fsdept_n(:,:,jk-1) + fse3w_n(:,:,jk) 160 fsdepw_n(:,:,jk) = fsdepw_n(:,:,jk-1) + fse3t_n(:,:,jk-1) 161 fsde3w_n(:,:,jk) = fsdept_n(:,:,jk ) - sshn (:,:) 119 162 END DO 120 CALL lbc_lnk( zee_f, 'F', 1. ) ! lateral boundary condition on ee_f 121 ! 122 DO jk = 1, jpk ! mu coefficients 123 mut(:,:,jk) = zee_t(:,:) * tmask(:,:,jk) ! T-point at T levels 124 muu(:,:,jk) = zee_u(:,:) * umask(:,:,jk) ! U-point at T levels 125 muv(:,:,jk) = zee_v(:,:) * vmask(:,:,jk) ! V-point at T levels 126 END DO 127 DO jk = 1, jpk ! F-point : fmask=shlat at coasts, use the product of umask 128 DO jj = 1, jpjm1 129 muf(:,jj,jk) = zee_f(:,jj) * umask(:,jj,jk) * umask(:,jj+1,jk) ! at T levels 130 END DO 131 muf(:,jpj,jk) = 0._wp 132 END DO 133 CALL lbc_lnk( muf, 'F', 1. ) ! lateral boundary condition 134 135 136 hu_0(:,:) = 0.e0 ! Reference ocean depth at U- and V-points 163 ! Reference water column height at t-, u- and v- point 164 ! ---------------------------------------------------- 165 ht_0(:,:) = 0.e0 166 hu_0(:,:) = 0.e0 137 167 hv_0(:,:) = 0.e0 138 168 DO jk = 1, jpk 139 hu_0(:,:) = hu_0(:,:) + fse3u_0(:,:,jk) * umask(:,:,jk) 140 hv_0(:,:) = hv_0(:,:) + fse3v_0(:,:,jk) * vmask(:,:,jk) 169 ht_0(:,:) = ht_0(:,:) + e3t_0(:,:,jk) * tmask(:,:,jk) 170 hu_0(:,:) = hu_0(:,:) + e3u_0(:,:,jk) * umask(:,:,jk) 171 hv_0(:,:) = hv_0(:,:) + e3v_0(:,:,jk) * vmask(:,:,jk) 141 172 END DO 142 143 DO jj = 1, jpjm1 ! initialise before and now Sea Surface Height at u-, v-, f-points 144 DO ji = 1, jpim1 ! NO vector opt. 145 zcoefu = 0.50_wp / ( e1u(ji,jj) * e2u(ji,jj) ) * umask(ji,jj,1) 146 zcoefv = 0.50_wp / ( e1v(ji,jj) * e2v(ji,jj) ) * vmask(ji,jj,1) 147 zcoeff = 0.25_wp / ( e1f(ji,jj) * e2f(ji,jj) ) * umask(ji,jj,1) * umask(ji,jj+1,1) 173 174 ! Restoring frequencies for z_tilde coordinate 175 ! ============================================ 176 IF( ln_vvl_ztilde ) THEN 177 ! - ML - In the future, this should be tunable by the user (namelist) 178 frq_rst_e3t(:,:) = 2.e0_wp * rpi / ( 30.e0_wp * 86400.e0_wp ) 179 frq_rst_hdv(:,:) = 2.e0_wp * rpi / ( 5.e0_wp * 86400.e0_wp ) 180 ENDIF 181 182 IF( nn_timing == 1 ) CALL timing_stop('dom_vvl_init') 183 184 END SUBROUTINE dom_vvl_init 185 186 187 SUBROUTINE dom_vvl_sf_nxt( kt ) 188 !!---------------------------------------------------------------------- 189 !! *** ROUTINE dom_vvl_sf_nxt *** 190 !! 191 !! ** Purpose : - compute the after scale factors used in tra_zdf, dynnxt, 192 !! tranxt and dynspg routines 193 !! 194 !! ** Method : - z_star case: Repartition of ssh INCREMENT proportionnaly to the level thickness. 195 !! - z_tilde_case: after scale factor increment = 196 !! high frequency part of horizontal divergence 197 !! + retsoring towards the background grid 198 !! + thickness difusion 199 !! Then repartition of ssh INCREMENT proportionnaly 200 !! to the "baroclinic" level thickness. 201 !! 202 !! ** Action : - hdiv_lf: restoring towards full baroclinic divergence in z_tilde case 203 !! - e3t_t_a: after increment of vertical scale factor 204 !! in z_tilde case 205 !! - fse3(t/u/v)_a 206 !! 207 !! Reference : Leclair, M., and Madec, G. 2011, Ocean Modelling. 208 !!---------------------------------------------------------------------- 209 REAL(wp), POINTER, DIMENSION(:,:,:) :: ze3t 210 REAL(wp), POINTER, DIMENSION(:,: ) :: zht, z_scale, zwu, zwv, zhdiv 211 !! * Arguments 212 INTEGER, INTENT( in ) :: kt ! time step 213 !! * Local declarations 214 INTEGER :: ji, jj, jk ! dummy loop indices 215 INTEGER , DIMENSION(3) :: ijk_max, ijk_min ! temporary integers 216 REAL(wp) :: z2dt ! temporary scalars 217 REAL(wp) :: z_def_max ! temporary scalar 218 REAL(wp) :: z_tmin, z_tmax ! temporary scalars 219 !!---------------------------------------------------------------------- 220 IF( nn_timing == 1 ) CALL timing_start('dom_vvl_sf_nxt') 221 CALL wrk_alloc( jpi, jpj, zht, z_scale, zwu, zwv, zhdiv ) 222 CALL wrk_alloc( jpi, jpj, jpk, ze3t ) 223 224 IF(kt == nit000) THEN 225 IF(lwp) WRITE(numout,*) 226 IF(lwp) WRITE(numout,*) 'dom_vvl_sf_nxt : compute after scale factors' 227 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~' 228 ENDIF 229 230 ! ******************************* ! 231 ! After acale factors at t-points ! 232 ! ******************************* ! 233 234 ! ! ----------------- ! 235 IF( ln_vvl_zstar ) THEN ! z_star coordinate ! 236 ! ! ----------------- ! 237 238 z_scale(:,:) = ( ssha(:,:) - sshb(:,:) ) * tmask(:,:,1) / ( ht_0(:,:) + sshn(:,:) + 1. - tmask(:,:,1) ) 239 DO jk = 1, jpkm1 240 fse3t_a(:,:,jk) = fse3t_b(:,:,jk) + fse3t_n(:,:,jk) * z_scale(:,:) * tmask(:,:,jk) 241 END DO 242 243 ! ! --------------------------- ! 244 ELSEIF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN ! z_tilde or layer coordinate ! 245 ! ! --------------------------- ! 246 247 ! I - initialization 248 ! ================== 249 250 ! 1 - barotropic divergence 251 ! ------------------------- 252 zhdiv(:,:) = 0. 253 zht(:,:) = 0. 254 DO jk = 1, jpkm1 255 zhdiv(:,:) = zhdiv(:,:) + fse3t_n(:,:,jk) * hdivn(:,:,jk) 256 zht (:,:) = zht (:,:) + fse3t_n(:,:,jk) * tmask(:,:,jk) 257 END DO 258 zhdiv(:,:) = zhdiv(:,:) / ( zht(:,:) + 1. - tmask(:,:,1) ) 259 260 ! 2 - Low frequency baroclinic horizontal divergence (z-tilde case only) 261 ! -------------------------------------------------- 262 IF( ln_vvl_ztilde ) THEN 263 IF( kt .GT. nit000 ) THEN 264 DO jk = 1, jpkm1 265 hdiv_lf(:,:,jk) = hdiv_lf(:,:,jk) - rdt * frq_rst_hdv(:,:) & 266 & * ( hdiv_lf(:,:,jk) - fse3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) ) 267 END DO 268 ENDIF 269 END IF 270 271 ! II - after z_tilde increments of vertical scale factors 272 ! ======================================================= 273 e3t_t_a(:,:,:) = 0.e0 ! e3t_t_a used to store tendency terms 274 275 ! 1 - High frequency divergence term 276 ! ---------------------------------- 277 IF( ln_vvl_ztilde ) THEN ! z_tilde case 278 DO jk = 1, jpkm1 279 e3t_t_a(:,:,jk) = e3t_t_a(:,:,jk) - ( fse3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) - hdiv_lf(:,:,jk) ) 280 END DO 281 ELSE ! layer case 282 DO jk = 1, jpkm1 283 e3t_t_a(:,:,jk) = e3t_t_a(:,:,jk) - fse3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) 284 END DO 285 END IF 286 287 ! 2 - Restoring term (z-tilde case only) 288 ! ------------------ 289 IF( ln_vvl_ztilde ) THEN 290 DO jk = 1, jpk 291 e3t_t_a(:,:,jk) = e3t_t_a(:,:,jk) - frq_rst_e3t(:,:) * e3t_t_b(:,:,jk) 292 END DO 293 END IF 294 295 ! 3 - Thickness diffusion term 296 ! ---------------------------- 297 zwu(:,:) = 0.e0 298 zwv(:,:) = 0.e0 299 ! a - first derivative: diffusive fluxes 300 DO jk = 1, jpkm1 301 DO jj = 1, jpjm1 302 DO ji = 1, fs_jpim1 ! vector opt. 303 un_td(ji,jj,jk) = rn_ahe3 * umask(ji,jj,jk) * re2u_e1u(ji,jj) * ( e3t_t_b(ji,jj,jk) - e3t_t_b(ji+1,jj ,jk) ) 304 vn_td(ji,jj,jk) = rn_ahe3 * vmask(ji,jj,jk) * re1v_e2v(ji,jj) * ( e3t_t_b(ji,jj,jk) - e3t_t_b(ji ,jj+1,jk) ) 305 zwu(ji,jj) = zwu(ji,jj) + un_td(ji,jj,jk) 306 zwv(ji,jj) = zwv(ji,jj) + vn_td(ji,jj,jk) 307 END DO 308 END DO 309 END DO 310 ! b - correction for last oceanic u-v points 311 DO jj = 1, jpj 312 DO ji = 1, jpi 313 un_td(ji,jj,mbku(ji,jj)) = un_td(ji,jj,mbku(ji,jj)) - zwu(ji,jj) 314 vn_td(ji,jj,mbkv(ji,jj)) = vn_td(ji,jj,mbkv(ji,jj)) - zwv(ji,jj) 315 END DO 316 END DO 317 ! c - second derivative: divergence of diffusive fluxes 318 DO jk = 1, jpkm1 319 DO jj = 2, jpjm1 320 DO ji = fs_2, fs_jpim1 ! vector opt. 321 e3t_t_a(ji,jj,jk) = e3t_t_a(ji,jj,jk) + ( un_td(ji-1,jj ,jk) - un_td(ji,jj,jk) & 322 & + vn_td(ji ,jj-1,jk) - vn_td(ji,jj,jk) ) & 323 & * r1_e12t(ji,jj) 324 END DO 325 END DO 326 END DO 327 ! d - thickness diffusion transport: boundary conditions 328 ! (stored for tracer advction and continuity equation) 329 CALL lbc_lnk( un_td , 'U' , -1.) 330 CALL lbc_lnk( vn_td , 'V' , -1.) 331 332 ! 4 - Time stepping of baroclinic scale factors 333 ! --------------------------------------------- 334 ! Leapfrog time stepping 335 ! ~~~~~~~~~~~~~~~~~~~~~~ 336 IF( neuler == 0 .AND. kt == nit000 ) THEN 337 z2dt = rdt 338 ELSE 339 z2dt = 2.e0 * rdt 340 ENDIF 341 CALL lbc_lnk( e3t_t_a(:,:,:), 'T', 1. ) 342 e3t_t_a(:,:,:) = e3t_t_b(:,:,:) + z2dt * tmask(:,:,:) * e3t_t_a(:,:,:) 343 344 ! Maximum deformation control 345 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~ 346 ! - ML - Should perhaps be put in the namelist 347 z_def_max = 0.9e0 348 ze3t(:,:,jpk) = 0.e0 349 DO jk = 1, jpkm1 350 ze3t(:,:,jk) = e3t_t_a(:,:,jk) / e3t_0(:,:,jk) * tmask(:,:,jk) * tmask_i(:,:) 351 END DO 352 z_tmax = MAXVAL( ze3t(:,:,:) ) 353 z_tmin = MINVAL( ze3t(:,:,:) ) 354 ! - ML - test: for the moment, stop simulation for too large e3_t variations 355 IF( ( z_tmax .GT. z_def_max ) .OR. ( z_tmin .LT. - z_def_max ) ) THEN 356 ijk_max = MAXLOC( ze3t(:,:,:) ) 357 ijk_min = MINLOC( ze3t(:,:,:) ) 358 WRITE(numout, *) 'MAX( e3t_t_a(:,:,:) / e3t_0(:,:,:) ) =', z_tmax 359 WRITE(numout, *) 'at i, j, k=', ijk_max 360 WRITE(numout, *) 'MIN( e3t_t_a(:,:,:) / e3t_0(:,:,:) ) =', z_tmin 361 WRITE(numout, *) 'at i, j, k=', ijk_min 362 CALL ctl_stop('MAX( ABS( e3t_t_a(:,:,:) ) / e3t_0(:,:,:) ) too high') 363 ENDIF 364 ! - ML - end test 365 ! - ML - This will cause a baroclinicity error if the ctl_stop above is not used 366 e3t_t_a(:,:,:) = MIN( e3t_t_a(:,:,:), z_def_max * e3t_0(:,:,:) ) 367 e3t_t_a(:,:,:) = MAX( e3t_t_a(:,:,:), - z_def_max * e3t_0(:,:,:) ) 368 369 ! Add "tilda" part to the after scale factor 370 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 371 fse3t_a(:,:,:) = e3t_0(:,:,:) + e3t_t_a(:,:,:) 372 373 ! III - Barotropic repartition of the sea surface height over the baroclinic profile 374 ! ================================================================================== 375 ! add e3t(n-1) "star" Asselin-filtered 376 DO jk = 1, jpkm1 377 fse3t_a(:,:,jk) = fse3t_a(:,:,jk) + fse3t_b(:,:,jk) - e3t_0(:,:,jk) - e3t_t_b(:,:,jk) 378 END DO 379 ! add ( ssh increment + "baroclinicity error" ) proportionnaly to e3t(n) 380 ! - ML - baroclinicity error should be better treated in the future 381 ! i.e. locally and not spread over the water column. 382 ! (keep in mind that the idea is to reduce Eulerian velocity as much as possible) 383 zht(:,:) = 0. 384 DO jk = 1, jpkm1 385 zht(:,:) = zht(:,:) + e3t_t_a(:,:,jk) * tmask(:,:,jk) 386 END DO 387 z_scale(:,:) = ( ssha(:,:) - sshb(:,:) - zht(:,:) ) / ( ht_0(:,:) + sshn(:,:) + 1. - tmask(:,:,1) ) 388 DO jk = 1, jpkm1 389 fse3t_a(:,:,jk) = fse3t_a(:,:,jk) + fse3t_n(:,:,jk) * z_scale(:,:) * tmask(:,:,jk) 390 END DO 391 392 ENDIF 393 394 IF( ln_vvl_dbg ) THEN ! - ML - test: control prints for debuging 395 IF ( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN 396 WRITE(numout, *) 'kt =', kt 397 WRITE(numout, *) 'MAXVAL(abs(SUM(e3t_t_a))) =', & 398 & MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( zht(:,:) ) ) 399 END IF 400 zht(:,:) = 0.e0 401 DO jk = 1, jpkm1 402 zht(:,:) = zht(:,:) + fse3t_n(:,:,jk) * tmask(:,:,jk) 403 END DO 404 WRITE(numout, *) 'MAXVAL(abs(ht_0+sshn-SUM(fse3t_n))) =', & 405 & MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + sshn(:,:) - zht(:,:) ) ) 406 zht(:,:) = 0.e0 407 DO jk = 1, jpkm1 408 zht(:,:) = zht(:,:) + fse3t_a(:,:,jk) * tmask(:,:,jk) 409 END DO 410 WRITE(numout, *) 'MAXVAL(abs(ht_0+ssha-SUM(fse3t_a))) =', & 411 & MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + ssha(:,:) - zht(:,:) ) ) 412 END IF 413 414 ! *********************************** ! 415 ! After scale factors at u- v- points ! 416 ! *********************************** ! 417 418 CALL dom_vvl_interpol( fse3t_a(:,:,:), fse3u_a(:,:,:), 'U' ) 419 CALL dom_vvl_interpol( fse3t_a(:,:,:), fse3v_a(:,:,:), 'V' ) 420 421 CALL wrk_dealloc( jpi, jpj, zht, z_scale, zwu, zwv, zhdiv ) 422 CALL wrk_dealloc( jpi, jpj, jpk, ze3t ) 423 424 IF( nn_timing == 1 ) CALL timing_start('dom_vvl_sf_nxt') 425 426 END SUBROUTINE dom_vvl_sf_nxt 427 428 429 SUBROUTINE dom_vvl_sf_swp( kt ) 430 !!---------------------------------------------------------------------- 431 !! *** ROUTINE dom_vvl_sf_swp *** 432 !! 433 !! ** Purpose : compute time filter and swap of scale factors 434 !! compute all depths and related variables for next time step 435 !! write outputs and restart file 436 !! 437 !! ** Method : - swap of e3t with trick for volume/tracer conservation 438 !! - reconstruct scale factor at other grid points (interpolate) 439 !! - recompute depths and water height fields 440 !! 441 !! ** Action : - fse3t_(b/n), e3t_t_(b/n) and fse3(u/v)_n ready for next time step 442 !! - Recompute: 443 !! fse3(u/v)_b 444 !! fse3w_n 445 !! fse3(u/v)w_b 446 !! fse3(u/v)w_n 447 !! fsdept_n, fsdepw_n and fsde3w_n 448 !! h(u/v) and h(u/v)r 449 !! 450 !! Reference : Leclair, M., and G. Madec, 2009, Ocean Modelling. 451 !! Leclair, M., and G. Madec, 2011, Ocean Modelling. 452 !!---------------------------------------------------------------------- 453 !! * Arguments 454 INTEGER, INTENT( in ) :: kt ! time step 455 !! * Local declarations 456 REAL(wp), POINTER, DIMENSION(:,:,:) :: z_e3t_def 457 INTEGER :: jk ! dummy loop indices 458 !!---------------------------------------------------------------------- 459 460 IF( nn_timing == 1 ) CALL timing_start('dom_vvl_sf_swp') 461 ! 462 CALL wrk_alloc( jpi, jpj, jpk, z_e3t_def ) 463 ! 464 IF( kt == nit000 ) THEN 465 IF(lwp) WRITE(numout,*) 466 IF(lwp) WRITE(numout,*) 'dom_vvl_sf_swp : - time filter and swap of scale factors' 467 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~ - interpolate scale factors and compute depths for next time step' 468 ENDIF 469 ! 470 ! Time filter and swap of scale factors 471 ! ===================================== 472 ! - ML - fse3(t/u/v)_b are allready computed in dynnxt. 473 IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN 474 IF( neuler == 0 .AND. kt == nit000 ) THEN 475 e3t_t_b(:,:,:) = e3t_t_n(:,:,:) 476 ELSE 477 e3t_t_b(:,:,:) = e3t_t_n(:,:,:) + atfp * ( e3t_t_b(:,:,:) - 2.e0 * e3t_t_n(:,:,:) + e3t_t_a(:,:,:) ) 478 ENDIF 479 e3t_t_n(:,:,:) = e3t_t_a(:,:,:) 480 ENDIF 481 fse3t_n(:,:,:) = fse3t_a(:,:,:) 482 fse3u_n(:,:,:) = fse3u_a(:,:,:) 483 fse3v_n(:,:,:) = fse3v_a(:,:,:) 484 485 ! Compute all missing vertical scale factor and depths 486 ! ==================================================== 487 ! Horizontal scale factor interpolations 488 ! -------------------------------------- 489 ! - ML - fse3u_b and fse3v_b are allready computed in dynnxt 490 CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3f_n (:,:,:), 'F' ) 491 ! Vertical scale factor interpolations 492 ! ------------------------------------ 493 CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3w_n (:,:,:), 'W' ) 494 CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3uw_n(:,:,:), 'UW' ) 495 CALL dom_vvl_interpol( fse3v_n(:,:,:), fse3vw_n(:,:,:), 'VW' ) 496 CALL dom_vvl_interpol( fse3u_b(:,:,:), fse3uw_b(:,:,:), 'UW' ) 497 CALL dom_vvl_interpol( fse3v_b(:,:,:), fse3vw_b(:,:,:), 'VW' ) 498 ! t- and w- points depth 499 ! ---------------------- 500 fsdept_n(:,:,1) = 0.5 * fse3w_n(:,:,1) 501 fsdepw_n(:,:,1) = 0.e0 502 fsde3w_n(:,:,1) = fsdept_n(:,:,1) - sshn(:,:) 503 DO jk = 2, jpk 504 fsdept_n(:,:,jk) = fsdept_n(:,:,jk-1) + fse3w_n(:,:,jk) 505 fsdepw_n(:,:,jk) = fsdepw_n(:,:,jk-1) + fse3t_n(:,:,jk-1) 506 fsde3w_n(:,:,jk) = fsdept_n(:,:,jk ) - sshn (:,:) 507 END DO 508 ! Local depth and Inverse of the local depth of the water column at u- and v- points 509 ! ---------------------------------------------------------------------------------- 510 hu(:,:) = 0. 511 hv(:,:) = 0. 512 DO jk = 1, jpk 513 hu(:,:) = hu(:,:) + fse3u_n(:,:,jk) * umask(:,:,jk) 514 hv(:,:) = hv(:,:) + fse3v_n(:,:,jk) * vmask(:,:,jk) 515 END DO 516 ! Inverse of the local depth 517 hur(:,:) = umask(:,:,1) / ( hu(:,:) + 1. - umask(:,:,1) ) 518 hvr(:,:) = vmask(:,:,1) / ( hv(:,:) + 1. - vmask(:,:,1) ) 519 520 ! Write outputs 521 ! ============= 522 z_e3t_def(:,:,:) = ( ( fse3t_n(:,:,:) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2 523 CALL iom_put( "e3t_n" , fse3t_n (:,:,:) ) 524 CALL iom_put( "dept_n" , fsde3w_n (:,:,:) ) 525 CALL iom_put( "e3tdef" , z_e3t_def(:,:,:) ) 526 527 ! write restart file 528 ! ================== 529 IF( lrst_oce ) CALL dom_vvl_rst( kt, 'WRITE' ) 530 ! 531 CALL wrk_dealloc( jpi, jpj, jpk, z_e3t_def ) 532 ! 533 IF( nn_timing == 1 ) CALL timing_stop('dom_vvl_sf_swp') 534 535 END SUBROUTINE dom_vvl_sf_swp 536 537 538 SUBROUTINE dom_vvl_interpol( pe3_in, pe3_out, pout ) 539 !!--------------------------------------------------------------------- 540 !! *** ROUTINE dom_vvl__interpol *** 541 !! 542 !! ** Purpose : interpolate scale factors from one grid point to another 543 !! 544 !! ** Method : e3_out = e3_0 + interpolation(e3_in - e3_0) 545 !! - horizontal interpolation: grid cell surface averaging 546 !! - vertical interpolation: simple averaging 547 !!---------------------------------------------------------------------- 548 !! * Arguments 549 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( in ) :: pe3_in ! input e3 to be interpolated 550 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( inout ) :: pe3_out ! output interpolated e3 551 CHARACTER(LEN=1), INTENT( in ) :: pout ! grid point of out scale factors 552 ! ! = 'U', 'V', 'W, 'F', 'UW' or 'VW' 553 !! * Local declarations 554 INTEGER :: ji, jj, jk ! dummy loop indices 555 !!---------------------------------------------------------------------- 556 IF( nn_timing == 1 ) CALL timing_start('dom_vvl_interpol') 557 SELECT CASE ( pout ) 558 ! ! ------------------------------------- ! 559 CASE( 'U' ) ! interpolation from T-point to U-point ! 560 ! ! ------------------------------------- ! 561 ! horizontal surface weighted interpolation 562 DO jk = 1, jpk 563 DO jj = 1, jpjm1 564 DO ji = 1, fs_jpim1 ! vector opt. 565 pe3_out(ji,jj,jk) = 0.5 * umask(ji,jj,jk) * r1_e12u(ji,jj) & 566 & * ( e12t(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3t_0(ji ,jj,jk) ) & 567 & + e12t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) ) 568 END DO 569 END DO 570 END DO 571 ! boundary conditions 572 CALL lbc_lnk( pe3_out(:,:,:), 'U', 1. ) 573 pe3_out(:,:,:) = pe3_out(:,:,:) + e3u_0(:,:,:) 574 ! ! ------------------------------------- ! 575 CASE( 'V' ) ! interpolation from T-point to V-point ! 576 ! ! ------------------------------------- ! 577 ! horizontal surface weighted interpolation 578 DO jk = 1, jpk 579 DO jj = 1, jpjm1 580 DO ji = 1, fs_jpim1 ! vector opt. 581 pe3_out(ji,jj,jk) = 0.5 * vmask(ji,jj,jk) * r1_e12v(ji,jj) & 582 & * ( e12t(ji,jj ) * ( pe3_in(ji,jj ,jk) - e3t_0(ji,jj ,jk) ) & 583 & + e12t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) ) 584 END DO 585 END DO 586 END DO 587 ! boundary conditions 588 CALL lbc_lnk( pe3_out(:,:,:), 'V', 1. ) 589 pe3_out(:,:,:) = pe3_out(:,:,:) + e3v_0(:,:,:) 590 ! ! ------------------------------------- ! 591 CASE( 'F' ) ! interpolation from U-point to F-point ! 592 ! ! ------------------------------------- ! 593 ! horizontal surface weighted interpolation 594 DO jk = 1, jpk 595 DO jj = 1, jpjm1 596 DO ji = 1, fs_jpim1 ! vector opt. 597 pe3_out(ji,jj,jk) = 0.5 * umask(ji,jj,jk) * umask(ji,jj+1,jk) * r1_e12f(ji,jj) & 598 & * ( e12u(ji,jj ) * ( pe3_in(ji,jj ,jk) - e3u_0(ji,jj ,jk) ) & 599 & + e12u(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3u_0(ji,jj+1,jk) ) ) 600 END DO 601 END DO 602 END DO 603 ! boundary conditions 604 CALL lbc_lnk( pe3_out(:,:,:), 'F', 1. ) 605 pe3_out(:,:,:) = pe3_out(:,:,:) + e3f_0(:,:,:) 606 ! ! ------------------------------------- ! 607 CASE( 'W' ) ! interpolation from T-point to W-point ! 608 ! ! ------------------------------------- ! 609 ! vertical simple interpolation 610 pe3_out(:,:,1) = e3w_0(:,:,1) + pe3_in(:,:,1) - e3t_0(:,:,1) 611 ! - ML - The use of mask in this formaula enables the special treatment of the last w- point without indirect adressing 612 DO jk = 2, jpk 613 pe3_out(:,:,jk) = e3w_0(:,:,jk) + ( 1. - 0.5 * tmask(:,:,jk) ) * ( pe3_in(:,:,jk-1) - e3t_0(:,:,jk-1) ) & 614 & + 0.5 * tmask(:,:,jk) * ( pe3_in(:,:,jk ) - e3t_0(:,:,jk ) ) 615 END DO 616 ! ! -------------------------------------- ! 617 CASE( 'UW' ) ! interpolation from U-point to UW-point ! 618 ! ! -------------------------------------- ! 619 ! vertical simple interpolation 620 pe3_out(:,:,1) = e3uw_0(:,:,1) + pe3_in(:,:,1) - e3u_0(:,:,1) 621 ! - ML - The use of mask in this formaula enables the special treatment of the last w- point without indirect adressing 622 DO jk = 2, jpk 623 pe3_out(:,:,jk) = e3uw_0(:,:,jk) + ( 1. - 0.5 * umask(:,:,jk) ) * ( pe3_in(:,:,jk-1) - e3u_0(:,:,jk-1) ) & 624 & + 0.5 * umask(:,:,jk) * ( pe3_in(:,:,jk ) - e3u_0(:,:,jk ) ) 625 END DO 626 ! ! -------------------------------------- ! 627 CASE( 'VW' ) ! interpolation from V-point to VW-point ! 628 ! ! -------------------------------------- ! 629 ! vertical simple interpolation 630 pe3_out(:,:,1) = e3vw_0(:,:,1) + pe3_in(:,:,1) - e3v_0(:,:,1) 631 ! - ML - The use of mask in this formaula enables the special treatment of the last w- point without indirect adressing 632 DO jk = 2, jpk 633 pe3_out(:,:,jk) = e3vw_0(:,:,jk) + ( 1. - 0.5 * vmask(:,:,jk) ) * ( pe3_in(:,:,jk-1) - e3v_0(:,:,jk-1) ) & 634 & + 0.5 * vmask(:,:,jk) * ( pe3_in(:,:,jk ) - e3v_0(:,:,jk ) ) 635 END DO 636 END SELECT 637 638 IF( nn_timing == 1 ) CALL timing_stop('dom_vvl_interpol') 639 640 END SUBROUTINE dom_vvl_interpol 641 642 643 SUBROUTINE dom_vvl_rst( kt, cdrw ) 644 !!--------------------------------------------------------------------- 645 !! *** ROUTINE dom_vvl_rst *** 646 !! 647 !! ** Purpose : Read or write VVL file in restart file 648 !! 649 !! ** Method : use of IOM library 650 !! if the restart does not contain vertical scale factors, 651 !! they are set to the _0 values 652 !! if the restart does not contain vertical scale factors increments (z_tilde), 653 !! they are set to 0. 654 !!---------------------------------------------------------------------- 655 !! * Arguments 656 INTEGER , INTENT(in) :: kt ! ocean time-step 657 CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag 658 !! * Local declarations 659 INTEGER :: id1, id2, id3, id4, id5 ! local integers 660 !!---------------------------------------------------------------------- 661 ! 662 IF( nn_timing == 1 ) CALL timing_start('dom_vvl_rst') 663 IF( TRIM(cdrw) == 'READ' ) THEN ! Read/initialise 664 ! ! =============== 665 IF( ln_rstart ) THEN !* Read the restart file 666 id1 = iom_varid( numror, 'fse3t_b', ldstop = .FALSE. ) 667 id2 = iom_varid( numror, 'fse3t_n', ldstop = .FALSE. ) 668 id3 = iom_varid( numror, 'e3t_t_b', ldstop = .FALSE. ) 669 id4 = iom_varid( numror, 'e3t_t_n', ldstop = .FALSE. ) 670 id5 = iom_varid( numror, 'hdif_lf', ldstop = .FALSE. ) 671 ! ! --------- ! 672 ! ! all cases ! 673 ! ! --------- ! 674 IF( MIN( id1, id2 ) > 0 ) THEN ! all required arrays exist 675 CALL iom_get( numror, jpdom_autoglo, 'fse3t_b', fse3t_b(:,:,:) ) 676 CALL iom_get( numror, jpdom_autoglo, 'fse3t_n', fse3t_n(:,:,:) ) 677 IF( neuler == 0 ) THEN 678 fse3t_b(:,:,:) = fse3t_n(:,:,:) 679 ENDIF 680 ELSE ! one at least array is missing 681 CALL ctl_stop( 'dom_vvl_rst: vvl cannot restart from a non vvl run' ) 682 ENDIF 683 ! ! ----------- ! 684 IF( ln_vvl_zstar ) THEN ! z_star case ! 685 ! ! ----------- ! 686 IF( MIN( id3, id4 ) > 0 ) THEN 687 CALL ctl_stop( 'dom_vvl_rst: z_star cannot restart from a z_tilde or layer run' ) 688 ENDIF 689 ! ! ----------------------- ! 690 ELSE ! z_tilde and layer cases ! 691 ! ! ----------------------- ! 692 IF( MIN( id3, id4 ) > 0 ) THEN ! all required arrays exist 693 CALL iom_get( numror, jpdom_autoglo, 'e3t_t_b', e3t_t_b(:,:,:) ) 694 CALL iom_get( numror, jpdom_autoglo, 'e3t_t_n', e3t_t_n(:,:,:) ) 695 ELSE ! one at least array is missing 696 e3t_t_b(:,:,:) = 0.e0 697 e3t_t_n(:,:,:) = 0.e0 698 ENDIF 699 ! ! ------------ ! 700 IF( ln_vvl_ztilde ) THEN ! z_tilde case ! 701 ! ! ------------ ! 702 IF( id5 > 0 ) THEN ! required array exists 703 CALL iom_get( numror, jpdom_autoglo, 'hdiv_lf', hdiv_lf(:,:,:) ) 704 ELSE ! array is missing 705 hdiv_lf(:,:,:) = 0.e0 706 ENDIF 707 ENDIF 708 ENDIF 148 709 ! 149 zvt = e1e2t(ji ,jj ) * sshb(ji ,jj ) ! before fields 150 zvt_ip1 = e1e2t(ji+1,jj ) * sshb(ji+1,jj ) 151 zvt_jp1 = e1e2t(ji ,jj+1) * sshb(ji ,jj+1) 152 sshu_b(ji,jj) = zcoefu * ( zvt + zvt_ip1 ) 153 sshv_b(ji,jj) = zcoefv * ( zvt + zvt_jp1 ) 154 ! 155 zvt = e1e2t(ji ,jj ) * sshn(ji ,jj ) ! now fields 156 zvt_ip1 = e1e2t(ji+1,jj ) * sshn(ji+1,jj ) 157 zvt_jp1 = e1e2t(ji ,jj+1) * sshn(ji ,jj+1) 158 zvt_ip1jp1 = e1e2t(ji+1,jj+1) * sshn(ji+1,jj+1) 159 sshu_n(ji,jj) = zcoefu * ( zvt + zvt_ip1 ) 160 sshv_n(ji,jj) = zcoefv * ( zvt + zvt_jp1 ) 161 sshf_n(ji,jj) = zcoeff * ( zvt + zvt_ip1 + zvt_jp1 + zvt_ip1jp1 ) 162 END DO 163 END DO 164 CALL lbc_lnk( sshu_n, 'U', 1. ) ; CALL lbc_lnk( sshu_b, 'U', 1. ) ! lateral boundary conditions 165 CALL lbc_lnk( sshv_n, 'V', 1. ) ; CALL lbc_lnk( sshv_b, 'V', 1. ) 166 CALL lbc_lnk( sshf_n, 'F', 1. ) 167 ! 168 CALL wrk_dealloc( jpi, jpj, zee_t, zee_u, zee_v, zee_f ) 169 ! 170 IF( nn_timing == 1 ) CALL timing_stop('dom_vvl') 171 ! 172 END SUBROUTINE dom_vvl 173 174 175 SUBROUTINE dom_vvl_2( kt, pe3u_b, pe3v_b ) 176 !!---------------------------------------------------------------------- 177 !! *** ROUTINE dom_vvl_2 *** 178 !! 179 !! ** Purpose : compute the vertical scale factors at u- and v-points 180 !! in variable volume case. 181 !! 182 !! ** Method : In variable volume case (non linear sea surface) the 183 !! the vertical scale factor at velocity points is computed 184 !! as the average of the cell surface weighted e3t. 185 !! It uses the sea surface heigth so it have to be initialized 186 !! after ssh is read/set 187 !!---------------------------------------------------------------------- 188 INTEGER , INTENT(in ) :: kt ! ocean time-step index 189 REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: pe3u_b, pe3v_b ! before vertical scale factor at u- & v-pts 190 ! 191 INTEGER :: ji, jj, jk ! dummy loop indices 192 INTEGER :: iku, ikv ! local integers 193 INTEGER :: ii0, ii1, ij0, ij1 ! temporary integers 194 REAL(wp) :: zvt ! local scalars 195 !!---------------------------------------------------------------------- 196 ! 197 IF( nn_timing == 1 ) CALL timing_start('dom_vvl_2') 198 ! 199 IF( lwp .AND. kt == nit000 ) THEN 710 ELSE !* Initialize at "rest" 711 fse3t_b(:,:,:) = e3t_0(:,:,:) 712 fse3t_n(:,:,:) = e3t_0(:,:,:) 713 IF( ln_vvl_ztilde .OR. ln_vvl_layer) THEN 714 e3t_t_b(:,:,:) = 0.e0 715 e3t_t_n(:,:,:) = 0.e0 716 IF( ln_vvl_ztilde ) hdiv_lf(:,:,:) = 0.e0 717 END IF 718 ENDIF 719 720 ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN ! Create restart file 721 ! ! =================== 722 IF(lwp) WRITE(numout,*) '---- dom_vvl_rst ----' 723 ! ! --------- ! 724 ! ! all cases ! 725 ! ! --------- ! 726 CALL iom_rstput( kt, nitrst, numrow, 'fse3t_b', fse3t_b(:,:,:) ) 727 CALL iom_rstput( kt, nitrst, numrow, 'fse3t_n', fse3t_n(:,:,:) ) 728 ! ! ----------------------- ! 729 IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN ! z_tilde and layer cases ! 730 ! ! ----------------------- ! 731 CALL iom_rstput( kt, nitrst, numrow, 'e3t_t_b', e3t_t_b(:,:,:) ) 732 CALL iom_rstput( kt, nitrst, numrow, 'e3t_t_n', e3t_t_n(:,:,:) ) 733 END IF 734 ! ! -------------! 735 IF( ln_vvl_ztilde ) THEN ! z_tilde case ! 736 ! ! ------------ ! 737 CALL iom_rstput( kt, nitrst, numrow, 'hdiv_lf', hdiv_lf(:,:,:) ) 738 ENDIF 739 740 ENDIF 741 IF( nn_timing == 1 ) CALL timing_stop('dom_vvl_rst') 742 743 END SUBROUTINE dom_vvl_rst 744 745 746 SUBROUTINE dom_vvl_ctl 747 !!--------------------------------------------------------------------- 748 !! *** ROUTINE dom_vvl_ctl *** 749 !! 750 !! ** Purpose : Control the consistency between namelist options 751 !! for vertical coordinate 752 !!---------------------------------------------------------------------- 753 INTEGER :: ioptio 754 755 NAMELIST/nam_vvl/ ln_vvl_zstar, ln_vvl_ztilde, ln_vvl_layer, rn_ahe3, ln_vvl_dbg! , ln_vvl_kepe 756 !!---------------------------------------------------------------------- 757 758 REWIND ( numnam ) ! Read Namelist nam_vvl : vertical coordinate 759 READ ( numnam, nam_vvl ) 760 761 IF(lwp) THEN ! Namelist print 200 762 WRITE(numout,*) 201 WRITE(numout,*) 'dom_vvl_2 : Variable volume, fse3t_b initialization' 202 WRITE(numout,*) '~~~~~~~~~ ' 203 pe3u_b(:,:,jpk) = fse3u_0(:,:,jpk) 204 pe3v_b(:,:,jpk) = fse3u_0(:,:,jpk) 205 ENDIF 206 207 DO jk = 1, jpkm1 ! set the before scale factors at u- & v-points 208 DO jj = 2, jpjm1 209 DO ji = fs_2, fs_jpim1 210 zvt = fse3t_b(ji,jj,jk) * e1e2t(ji,jj) 211 pe3u_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji+1,jj,jk) * e1e2t(ji+1,jj) ) / ( e1u(ji,jj) * e2u(ji,jj) ) 212 pe3v_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji,jj+1,jk) * e1e2t(ji,jj+1) ) / ( e1v(ji,jj) * e2v(ji,jj) ) 213 END DO 214 END DO 215 END DO 216 217 ! Correct scale factors at locations that have been individually modified in domhgr 218 ! Such modifications break the relationship between e1e2t and e1u*e2u etc. Recompute 219 ! scale factors ignoring the modified metric. 220 ! ! ===================== 221 IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN ! ORCA R2 configuration 222 ! ! ===================== 223 IF( nn_cla == 0 ) THEN 224 ! 225 ii0 = 139 ; ii1 = 140 ! Gibraltar Strait (e2u was modified) 226 ij0 = 102 ; ij1 = 102 227 DO jk = 1, jpkm1 ! set the before scale factors at u-points 228 DO jj = mj0(ij0), mj1(ij1) 229 DO ji = mi0(ii0), mi1(ii1) 230 zvt = fse3t_b(ji,jj,jk) * e1t(ji,jj) 231 pe3u_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji+1,jj,jk) * e1t(ji+1,jj) ) / ( e1u(ji,jj) ) 232 END DO 233 END DO 234 END DO 235 ! 236 ii0 = 160 ; ii1 = 160 ! Bab el Mandeb (e2u and e1v were modified) 237 ij0 = 88 ; ij1 = 88 238 DO jk = 1, jpkm1 ! set the before scale factors at u-points 239 DO jj = mj0(ij0), mj1(ij1) 240 DO ji = mi0(ii0), mi1(ii1) 241 zvt = fse3t_b(ji,jj,jk) * e1t(ji,jj) 242 pe3u_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji+1,jj,jk) * e1t(ji+1,jj) ) / ( e1u(ji,jj) ) 243 END DO 244 END DO 245 END DO 246 DO jk = 1, jpkm1 ! set the before scale factors at v-points 247 DO jj = mj0(ij0), mj1(ij1) 248 DO ji = mi0(ii0), mi1(ii1) 249 zvt = fse3t_b(ji,jj,jk) * e2t(ji,jj) 250 pe3v_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji,jj+1,jk) * e2t(ji,jj+1) ) / ( e2v(ji,jj) ) 251 END DO 252 END DO 253 END DO 254 ENDIF 255 256 ii0 = 145 ; ii1 = 146 ! Danish Straits (e2u was modified) 257 ij0 = 116 ; ij1 = 116 258 DO jk = 1, jpkm1 ! set the before scale factors at u-points 259 DO jj = mj0(ij0), mj1(ij1) 260 DO ji = mi0(ii0), mi1(ii1) 261 zvt = fse3t_b(ji,jj,jk) * e1t(ji,jj) 262 pe3u_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji+1,jj,jk) * e1t(ji+1,jj) ) / ( e1u(ji,jj) ) 263 END DO 264 END DO 265 END DO 266 ! 267 ENDIF 268 ! ! ===================== 269 IF( cp_cfg == "orca" .AND. jp_cfg == 1 ) THEN ! ORCA R1 configuration 270 ! ! ===================== 271 272 ii0 = 281 ; ii1 = 282 ! Gibraltar Strait (e2u was modified) 273 ij0 = 200 ; ij1 = 200 274 DO jk = 1, jpkm1 ! set the before scale factors at u-points 275 DO jj = mj0(ij0), mj1(ij1) 276 DO ji = mi0(ii0), mi1(ii1) 277 zvt = fse3t_b(ji,jj,jk) * e1t(ji,jj) 278 pe3u_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji+1,jj,jk) * e1t(ji+1,jj) ) / ( e1u(ji,jj) ) 279 END DO 280 END DO 281 END DO 282 283 ii0 = 314 ; ii1 = 315 ! Bhosporus Strait (e2u was modified) 284 ij0 = 208 ; ij1 = 208 285 DO jk = 1, jpkm1 ! set the before scale factors at u-points 286 DO jj = mj0(ij0), mj1(ij1) 287 DO ji = mi0(ii0), mi1(ii1) 288 zvt = fse3t_b(ji,jj,jk) * e1t(ji,jj) 289 pe3u_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji+1,jj,jk) * e1t(ji+1,jj) ) / ( e1u(ji,jj) ) 290 END DO 291 END DO 292 END DO 293 294 ii0 = 44 ; ii1 = 44 ! Lombok Strait (e1v was modified) 295 ij0 = 124 ; ij1 = 125 296 DO jk = 1, jpkm1 ! set the before scale factors at v-points 297 DO jj = mj0(ij0), mj1(ij1) 298 DO ji = mi0(ii0), mi1(ii1) 299 zvt = fse3t_b(ji,jj,jk) * e2t(ji,jj) 300 pe3v_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji,jj+1,jk) * e2t(ji,jj+1) ) / ( e2v(ji,jj) ) 301 END DO 302 END DO 303 END DO 304 305 ii0 = 48 ; ii1 = 48 ! Sumba Strait (e1v was modified) [closed from bathy_11 on] 306 ij0 = 124 ; ij1 = 125 307 DO jk = 1, jpkm1 ! set the before scale factors at v-points 308 DO jj = mj0(ij0), mj1(ij1) 309 DO ji = mi0(ii0), mi1(ii1) 310 zvt = fse3t_b(ji,jj,jk) * e2t(ji,jj) 311 pe3v_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji,jj+1,jk) * e2t(ji,jj+1) ) / ( e2v(ji,jj) ) 312 END DO 313 END DO 314 END DO 315 316 ii0 = 53 ; ii1 = 53 ! Ombai Strait (e1v was modified) 317 ij0 = 124 ; ij1 = 125 318 DO jk = 1, jpkm1 ! set the before scale factors at v-points 319 DO jj = mj0(ij0), mj1(ij1) 320 DO ji = mi0(ii0), mi1(ii1) 321 zvt = fse3t_b(ji,jj,jk) * e2t(ji,jj) 322 pe3v_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji,jj+1,jk) * e2t(ji,jj+1) ) / ( e2v(ji,jj) ) 323 END DO 324 END DO 325 END DO 326 327 ii0 = 56 ; ii1 = 56 ! Timor Passage (e1v was modified) 328 ij0 = 124 ; ij1 = 125 329 DO jk = 1, jpkm1 ! set the before scale factors at v-points 330 DO jj = mj0(ij0), mj1(ij1) 331 DO ji = mi0(ii0), mi1(ii1) 332 zvt = fse3t_b(ji,jj,jk) * e2t(ji,jj) 333 pe3v_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji,jj+1,jk) * e2t(ji,jj+1) ) / ( e2v(ji,jj) ) 334 END DO 335 END DO 336 END DO 337 338 ii0 = 55 ; ii1 = 55 ! West Halmahera Strait (e1v was modified) 339 ij0 = 141 ; ij1 = 142 340 DO jk = 1, jpkm1 ! set the before scale factors at v-points 341 DO jj = mj0(ij0), mj1(ij1) 342 DO ji = mi0(ii0), mi1(ii1) 343 zvt = fse3t_b(ji,jj,jk) * e2t(ji,jj) 344 pe3v_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji,jj+1,jk) * e2t(ji,jj+1) ) / ( e2v(ji,jj) ) 345 END DO 346 END DO 347 END DO 348 349 ii0 = 58 ; ii1 = 58 ! East Halmahera Strait (e1v was modified) 350 ij0 = 141 ; ij1 = 142 351 DO jk = 1, jpkm1 ! set the before scale factors at v-points 352 DO jj = mj0(ij0), mj1(ij1) 353 DO ji = mi0(ii0), mi1(ii1) 354 zvt = fse3t_b(ji,jj,jk) * e2t(ji,jj) 355 pe3v_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji,jj+1,jk) * e2t(ji,jj+1) ) / ( e2v(ji,jj) ) 356 END DO 357 END DO 358 END DO 359 360 ! 361 ENDIF 362 ! ! ====================== 363 IF( cp_cfg == "orca" .AND. jp_cfg == 05 ) THEN ! ORCA R05 configuration 364 ! ! ====================== 365 ii0 = 563 ; ii1 = 564 ! Gibraltar Strait (e2u was modified) 366 ij0 = 327 ; ij1 = 327 367 DO jk = 1, jpkm1 ! set the before scale factors at u-points 368 DO jj = mj0(ij0), mj1(ij1) 369 DO ji = mi0(ii0), mi1(ii1) 370 zvt = fse3t_b(ji,jj,jk) * e1t(ji,jj) 371 pe3u_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji+1,jj,jk) * e1t(ji+1,jj) ) / ( e1u(ji,jj) ) 372 END DO 373 END DO 374 END DO 375 ! 376 ii0 = 627 ; ii1 = 628 ! Bosphore Strait (e2u was modified) 377 ij0 = 343 ; ij1 = 343 378 DO jk = 1, jpkm1 ! set the before scale factors at u-points 379 DO jj = mj0(ij0), mj1(ij1) 380 DO ji = mi0(ii0), mi1(ii1) 381 zvt = fse3t_b(ji,jj,jk) * e1t(ji,jj) 382 pe3u_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji+1,jj,jk) * e1t(ji+1,jj) ) / ( e1u(ji,jj) ) 383 END DO 384 END DO 385 END DO 386 ! 387 ii0 = 93 ; ii1 = 94 ! Sumba Strait (e2u was modified) 388 ij0 = 232 ; ij1 = 232 389 DO jk = 1, jpkm1 ! set the before scale factors at u-points 390 DO jj = mj0(ij0), mj1(ij1) 391 DO ji = mi0(ii0), mi1(ii1) 392 zvt = fse3t_b(ji,jj,jk) * e1t(ji,jj) 393 pe3u_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji+1,jj,jk) * e1t(ji+1,jj) ) / ( e1u(ji,jj) ) 394 END DO 395 END DO 396 END DO 397 ! 398 ii0 = 103 ; ii1 = 103 ! Ombai Strait (e2u was modified) 399 ij0 = 232 ; ij1 = 232 400 DO jk = 1, jpkm1 ! set the before scale factors at u-points 401 DO jj = mj0(ij0), mj1(ij1) 402 DO ji = mi0(ii0), mi1(ii1) 403 zvt = fse3t_b(ji,jj,jk) * e1t(ji,jj) 404 pe3u_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji+1,jj,jk) * e1t(ji+1,jj) ) / ( e1u(ji,jj) ) 405 END DO 406 END DO 407 END DO 408 ! 409 ii0 = 15 ; ii1 = 15 ! Palk Strait (e2u was modified) 410 ij0 = 270 ; ij1 = 270 411 DO jk = 1, jpkm1 ! set the before scale factors at u-points 412 DO jj = mj0(ij0), mj1(ij1) 413 DO ji = mi0(ii0), mi1(ii1) 414 zvt = fse3t_b(ji,jj,jk) * e1t(ji,jj) 415 pe3u_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji+1,jj,jk) * e1t(ji+1,jj) ) / ( e1u(ji,jj) ) 416 END DO 417 END DO 418 END DO 419 ! 420 ii0 = 87 ; ii1 = 87 ! Lombok Strait (e1v was modified) 421 ij0 = 232 ; ij1 = 233 422 DO jk = 1, jpkm1 ! set the before scale factors at v-points 423 DO jj = mj0(ij0), mj1(ij1) 424 DO ji = mi0(ii0), mi1(ii1) 425 zvt = fse3t_b(ji,jj,jk) * e2t(ji,jj) 426 pe3v_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji,jj+1,jk) * e2t(ji,jj+1) ) / ( e2v(ji,jj) ) 427 END DO 428 END DO 429 END DO 430 ! 431 ii0 = 662 ; ii1 = 662 ! Bab el Mandeb (e1v was modified) 432 ij0 = 276 ; ij1 = 276 433 DO jk = 1, jpkm1 ! set the before scale factors at v-points 434 DO jj = mj0(ij0), mj1(ij1) 435 DO ji = mi0(ii0), mi1(ii1) 436 zvt = fse3t_b(ji,jj,jk) * e2t(ji,jj) 437 pe3v_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji,jj+1,jk) * e2t(ji,jj+1) ) / ( e2v(ji,jj) ) 438 END DO 439 END DO 440 END DO 441 ! 442 ENDIF 443 ! End of individual corrections to scale factors 444 445 IF( ln_zps ) THEN ! minimum of the e3t at partial cell level 446 DO jj = 2, jpjm1 447 DO ji = fs_2, fs_jpim1 448 iku = mbku(ji,jj) 449 ikv = mbkv(ji,jj) 450 pe3u_b(ji,jj,iku) = MIN( fse3t_b(ji,jj,iku), fse3t_b(ji+1,jj ,iku) ) 451 pe3v_b(ji,jj,ikv) = MIN( fse3t_b(ji,jj,ikv), fse3t_b(ji ,jj+1,ikv) ) 452 END DO 453 END DO 454 ENDIF 455 456 pe3u_b(:,:,:) = pe3u_b(:,:,:) - fse3u_0(:,:,:) ! anomaly to avoid zero along closed boundary/extra halos 457 pe3v_b(:,:,:) = pe3v_b(:,:,:) - fse3v_0(:,:,:) 458 CALL lbc_lnk( pe3u_b(:,:,:), 'U', 1. ) ! lateral boundary conditions 459 CALL lbc_lnk( pe3v_b(:,:,:), 'V', 1. ) 460 pe3u_b(:,:,:) = pe3u_b(:,:,:) + fse3u_0(:,:,:) ! recover the full scale factor 461 pe3v_b(:,:,:) = pe3v_b(:,:,:) + fse3v_0(:,:,:) 462 ! 463 IF( nn_timing == 1 ) CALL timing_stop('dom_vvl_2') 464 ! 465 END SUBROUTINE dom_vvl_2 466 467 #else 468 !!---------------------------------------------------------------------- 469 !! Default option : Empty routine 470 !!---------------------------------------------------------------------- 471 CONTAINS 472 SUBROUTINE dom_vvl 473 END SUBROUTINE dom_vvl 474 SUBROUTINE dom_vvl_2(kdum, pudum, pvdum ) 475 USE par_kind 476 INTEGER , INTENT(in ) :: kdum 477 REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: pudum, pvdum 478 END SUBROUTINE dom_vvl_2 479 #endif 763 WRITE(numout,*) 'dom_vvl_ctl : choice/control of the variable vertical coordinate' 764 WRITE(numout,*) '~~~~~~~~~~~' 765 WRITE(numout,*) ' Namelist nam_vvl : chose a vertical coordinate' 766 WRITE(numout,*) ' zstar ln_vvl_zstar = ', ln_vvl_zstar 767 WRITE(numout,*) ' ztilde ln_vvl_ztilde = ', ln_vvl_ztilde 768 WRITE(numout,*) ' layer ln_vvl_layer = ', ln_vvl_layer 769 ! WRITE(numout,*) ' Namelist nam_vvl : chose kinetic-to-potential energy conservation' 770 ! WRITE(numout,*) ' ln_vvl_kepe = ', ln_vvl_kepe 771 WRITE(numout,*) ' Namelist nam_vvl : thickness diffusion coefficient' 772 WRITE(numout,*) ' rn_ahe3 = ', rn_ahe3 773 WRITE(numout,*) ' Namelist nam_vvl : debug prints' 774 WRITE(numout,*) ' ln_vvl_dbg = ', ln_vvl_dbg 775 ENDIF 776 777 ioptio = 0 ! Parameter control 778 IF( ln_vvl_zstar ) ioptio = ioptio + 1 779 IF( ln_vvl_ztilde ) ioptio = ioptio + 1 780 IF( ln_vvl_layer ) ioptio = ioptio + 1 781 782 IF( ioptio /= 1 ) CALL ctl_stop( 'Choose ONE vertical coordinate in namelist nam_vvl' ) 783 784 IF(lwp) THEN ! Print the choice 785 WRITE(numout,*) 786 IF( ln_vvl_zstar ) WRITE(numout,*) ' zstar vertical coordinate is used' 787 IF( ln_vvl_ztilde ) WRITE(numout,*) ' ztilde vertical coordinate is used' 788 IF( ln_vvl_layer ) WRITE(numout,*) ' layer vertical coordinate is used' 789 ! - ML - Option not developed yet 790 ! IF( ln_vvl_kepe ) WRITE(numout,*) ' kinetic to potential energy transfer : option used' 791 ! IF( .NOT. ln_vvl_kepe ) WRITE(numout,*) ' kinetic to potential energy transfer : option not used' 792 ENDIF 793 794 END SUBROUTINE dom_vvl_ctl 480 795 481 796 !!====================================================================== -
branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DOM/domwri.F90
r3862 r3865 183 183 CALL iom_rstput( 0, 0, inum4, 'esigw', esigw ) 184 184 ! 185 CALL iom_rstput( 0, 0, inum4, 'e3t ', e3t )! ! scale factors186 CALL iom_rstput( 0, 0, inum4, 'e3u ', e3u)187 CALL iom_rstput( 0, 0, inum4, 'e3v ', e3v)188 CALL iom_rstput( 0, 0, inum4, 'e3w ', e3w)185 CALL iom_rstput( 0, 0, inum4, 'e3t_0', e3t_0 ) ! ! scale factors 186 CALL iom_rstput( 0, 0, inum4, 'e3u_0', e3u_0 ) 187 CALL iom_rstput( 0, 0, inum4, 'e3v_0', e3v_0 ) 188 CALL iom_rstput( 0, 0, inum4, 'e3w_0', e3w_0 ) 189 189 CALL iom_rstput( 0, 0, inum4, 'rx1', rx1 ) ! ! Max. grid stiffness ratio 190 190 ! 191 CALL iom_rstput( 0, 0, inum4, 'gdept ' , gdept )! ! stretched system192 CALL iom_rstput( 0, 0, inum4, 'gdepw ' , gdepw)191 CALL iom_rstput( 0, 0, inum4, 'gdept_1d' , gdept_1d ) ! ! stretched system 192 CALL iom_rstput( 0, 0, inum4, 'gdepw_1d' , gdepw_1d ) 193 193 ENDIF 194 194 … … 196 196 ! 197 197 IF( nmsh <= 6 ) THEN ! ! 3D vertical scale factors 198 CALL iom_rstput( 0, 0, inum4, 'e3t ', e3t)199 CALL iom_rstput( 0, 0, inum4, 'e3u ', e3u)200 CALL iom_rstput( 0, 0, inum4, 'e3v ', e3v)201 CALL iom_rstput( 0, 0, inum4, 'e3w ', e3w)198 CALL iom_rstput( 0, 0, inum4, 'e3t_0', e3t_0 ) 199 CALL iom_rstput( 0, 0, inum4, 'e3u_0', e3u_0 ) 200 CALL iom_rstput( 0, 0, inum4, 'e3v_0', e3v_0 ) 201 CALL iom_rstput( 0, 0, inum4, 'e3w_0', e3w_0 ) 202 202 ELSE ! ! 2D masked bottom ocean scale factors 203 203 DO jj = 1,jpj 204 204 DO ji = 1,jpi 205 e3tp(ji,jj) = e3t (ji,jj,mbkt(ji,jj)) * tmask(ji,jj,1)206 e3wp(ji,jj) = e3w (ji,jj,mbkt(ji,jj)) * tmask(ji,jj,1)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) 207 207 END DO 208 208 END DO … … 212 212 ! 213 213 IF( nmsh <= 3 ) THEN ! ! 3D depth 214 CALL iom_rstput( 0, 0, inum4, 'gdept ', gdept, ktype = jp_r4 )214 CALL iom_rstput( 0, 0, inum4, 'gdept_0', gdept_0, ktype = jp_r4 ) 215 215 DO jk = 1,jpk 216 216 DO jj = 1, jpjm1 217 217 DO ji = 1, fs_jpim1 ! vector opt. 218 zdepu(ji,jj,jk) = MIN( gdept (ji,jj,jk) , gdept(ji+1,jj ,jk) )219 zdepv(ji,jj,jk) = MIN( gdept (ji,jj,jk) , gdept(ji ,jj+1,jk) )218 zdepu(ji,jj,jk) = MIN( gdept_0(ji,jj,jk) , gdept_0(ji+1,jj ,jk) ) 219 zdepv(ji,jj,jk) = MIN( gdept_0(ji,jj,jk) , gdept_0(ji ,jj+1,jk) ) 220 220 END DO 221 221 END DO … … 224 224 CALL iom_rstput( 0, 0, inum4, 'gdepu', zdepu, ktype = jp_r4 ) 225 225 CALL iom_rstput( 0, 0, inum4, 'gdepv', zdepv, ktype = jp_r4 ) 226 CALL iom_rstput( 0, 0, inum4, 'gdepw ', gdepw, ktype = jp_r4 )226 CALL iom_rstput( 0, 0, inum4, 'gdepw_0', gdepw_0, ktype = jp_r4 ) 227 227 ELSE ! ! 2D bottom depth 228 228 DO jj = 1,jpj 229 229 DO ji = 1,jpi 230 zprt(ji,jj) = gdept (ji,jj,mbkt(ji,jj) ) * tmask(ji,jj,1)231 zprw(ji,jj) = gdepw (ji,jj,mbkt(ji,jj)+1) * tmask(ji,jj,1)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) 232 232 END DO 233 233 END DO … … 246 246 CALL iom_rstput( 0, 0, inum4, 'gdept_1d', gdept_1d ) ! ! depth 247 247 CALL iom_rstput( 0, 0, inum4, 'gdepw_1d', gdepw_1d ) 248 CALL iom_rstput( 0, 0, inum4, 'e3t_1d' , e3t_1d ) 248 CALL iom_rstput( 0, 0, inum4, 'e3t_1d' , e3t_1d ) ! ! scale factors 249 249 CALL iom_rstput( 0, 0, inum4, 'e3w_1d' , e3w_1d ) 250 250 ENDIF -
branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90
r3862 r3865 146 146 IF( nprint == 1 .AND. lwp ) THEN 147 147 WRITE(numout,*) ' MIN val mbathy ', MINVAL( mbathy(:,:) ), ' MAX ', MAXVAL( mbathy(:,:) ) 148 WRITE(numout,*) ' MIN val depth t ', MINVAL( fsdept(:,:,:) ), &149 & ' w ', MINVAL( fsdepw(:,:,:) ), '3w ', MINVAL( fsde3w(:,:,:) )150 WRITE(numout,*) ' MIN val e3 t ', MINVAL( fse3t(:,:,:) ), ' f ', MINVAL( fse3f(:,:,:) ), &151 & ' u ', MINVAL( fse3u(:,:,:) ), ' u ', MINVAL( fse3v(:,:,:) ), &152 & ' uw', MINVAL( fse3uw(:,:,:)), ' vw', MINVAL( fse3vw(:,:,:)), &153 & ' w ', MINVAL( fse3w(:,:,:) )154 155 WRITE(numout,*) ' MAX val depth t ', MAXVAL( fsdept(:,:,:) ), &156 & ' w ', MAXVAL( fsdepw(:,:,:) ), '3w ', MAXVAL( fsde3w(:,:,:) )157 WRITE(numout,*) ' MAX val e3 t ', MAXVAL( fse3t(:,:,:) ), ' f ', MAXVAL( fse3f(:,:,:) ), &158 & ' u ', MAXVAL( fse3u(:,:,:) ), ' u ', MAXVAL( fse3v(:,:,:) ), &159 & ' uw', MAXVAL( fse3uw(:,:,:)), ' vw', MAXVAL( fse3vw(:,:,:)), &160 & ' w ', MAXVAL( fse3w(:,:,:) )148 WRITE(numout,*) ' MIN val depth t ', MINVAL( gdept_0(:,:,:) ), & 149 & ' w ', MINVAL( gdepw_0(:,:,:) ), '3w ', MINVAL( gdep3w_0(:,:,:) ) 150 WRITE(numout,*) ' MIN val e3 t ', MINVAL( e3t_0(:,:,:) ), ' f ', MINVAL( e3f_0(:,:,:) ), & 151 & ' u ', MINVAL( e3u_0(:,:,:) ), ' u ', MINVAL( e3v_0(:,:,:) ), & 152 & ' uw', MINVAL( e3uw_0(:,:,:)), ' vw', MINVAL( e3vw_0(:,:,:)), & 153 & ' w ', MINVAL( e3w_0(:,:,:) ) 154 155 WRITE(numout,*) ' MAX val depth t ', MAXVAL( gdept_0(:,:,:) ), & 156 & ' w ', MAXVAL( gdepw_0(:,:,:) ), '3w ', MAXVAL( gdep3w_0(:,:,:) ) 157 WRITE(numout,*) ' MAX val e3 t ', MAXVAL( e3t_0(:,:,:) ), ' f ', MAXVAL( e3f_0(:,:,:) ), & 158 & ' u ', MAXVAL( e3u_0(:,:,:) ), ' u ', MAXVAL( e3v_0(:,:,:) ), & 159 & ' uw', MAXVAL( e3uw_0(:,:,:)), ' vw', MAXVAL( e3vw_0(:,:,:)), & 160 & ' w ', MAXVAL( e3w_0(:,:,:) ) 161 161 ENDIF 162 162 ! … … 177 177 !! function the derivative of which gives the scale factors. 178 178 !! both depth and scale factors only depend on k (1d arrays). 179 !! w-level: gdepw_1d = fsdep(k)180 !! e3w_1d(k) = dk( fsdep)(k) = fse3(k)181 !! t-level: gdept_1d = fsdep(k+0.5)182 !! e3t_1d(k) = dk( fsdep)(k+0.5) = fse3(k+0.5)179 !! w-level: gdepw_1d = gdep(k) 180 !! e3w_1d(k) = dk(gdep)(k) = e3(k) 181 !! t-level: gdept_1d = gdep(k+0.5) 182 !! e3t_1d(k) = dk(gdep)(k+0.5) = e3(k+0.5) 183 183 !! 184 184 !! ** Action : - gdept_1d, gdepw_1d : depth of T- and W-point (m) … … 298 298 WRITE(numout,*) 299 299 WRITE(numout,*) ' Reference z-coordinate depth and scale factors:' 300 WRITE(numout, "(9x,' level gdept gdepw e3t e3w')" )300 WRITE(numout, "(9x,' level gdept_1d gdepw_1d e3t_1d e3w_1d ')" ) 301 301 WRITE(numout, "(10x, i4, 4f9.2)" ) ( jk, gdept_1d(jk), gdepw_1d(jk), e3t_1d(jk), e3w_1d(jk), jk = 1, jpk ) 302 302 ENDIF 303 303 DO jk = 1, jpk ! control positivity 304 IF( e3w_1d (jk) <= 0._wp .OR. e3t_1d (jk) <= 0._wp ) CALL ctl_stop( 'dom:zgr_z: e3w or e3t=< 0 ' )305 IF( gdepw_1d(jk) < 0._wp .OR. gdept_1d(jk) < 0._wp ) CALL ctl_stop( 'dom:zgr_z: gdepw or gdept< 0 ' )304 IF( e3w_1d (jk) <= 0._wp .OR. e3t_1d (jk) <= 0._wp ) CALL ctl_stop( 'dom:zgr_z: e3w_1d or e3t_1d =< 0 ' ) 305 IF( gdepw_1d(jk) < 0._wp .OR. gdept_1d(jk) < 0._wp ) CALL ctl_stop( 'dom:zgr_z: gdepw_1d or gdept_1d < 0 ' ) 306 306 END DO 307 307 ! … … 369 369 IF(lwp) WRITE(numout,*) ' bathymetry field: flat basin' 370 370 idta(:,:) = jpkm1 ! before last level 371 zdta(:,:) = gdepw_1d(jpk) 371 zdta(:,:) = gdepw_1d(jpk) ! last w-point depth 372 372 h_oce = gdepw_1d(jpk) 373 373 ELSE ! bump centered in the basin … … 378 378 r_bump = 50000._wp ! bump radius (meters) 379 379 h_bump = 2700._wp ! bump height (meters) 380 h_oce = gdepw_1d(jpk) 380 h_oce = gdepw_1d(jpk) ! background ocean depth (meters) 381 381 IF(lwp) WRITE(numout,*) ' bump characteristics: ' 382 382 IF(lwp) WRITE(numout,*) ' bump center (i,j) = ', ii_bump, ii_bump … … 440 440 CALL iom_close( inum ) 441 441 mbathy(:,:) = INT( bathy(:,:) ) 442 ! 442 ! ! ===================== 443 443 IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN ! ORCA R2 configuration 444 ! 444 ! ! ===================== 445 445 IF( nn_cla == 0 ) THEN 446 446 ii0 = 140 ; ii1 = 140 ! Gibraltar Strait open … … 780 780 ! 781 781 DO jk = 1, jpk 782 gdept(:,:,jk) = gdept_1d(jk)783 gdepw(:,:,jk) = gdepw_1d(jk)784 gdep3w(:,:,jk) = gdepw_1d(jk)785 e3t (:,:,jk) = e3t_1d(jk)786 e3u (:,:,jk) = e3t_1d(jk)787 e3v (:,:,jk) = e3t_1d(jk)788 e3f (:,:,jk) = e3t_1d(jk)789 e3w (:,:,jk) = e3w_1d(jk)790 e3uw(:,:,jk) = e3w_1d(jk)791 e3vw(:,:,jk) = e3w_1d(jk)782 gdept_0 (:,:,jk) = gdept_1d(jk) 783 gdepw_0 (:,:,jk) = gdepw_1d(jk) 784 gdep3w_0(:,:,jk) = gdepw_1d(jk) 785 e3t_0 (:,:,jk) = e3t_1d (jk) 786 e3u_0 (:,:,jk) = e3t_1d (jk) 787 e3v_0 (:,:,jk) = e3t_1d (jk) 788 e3f_0 (:,:,jk) = e3t_1d (jk) 789 e3w_0 (:,:,jk) = e3w_1d (jk) 790 e3uw_0 (:,:,jk) = e3w_1d (jk) 791 e3vw_0 (:,:,jk) = e3w_1d (jk) 792 792 END DO 793 793 ! … … 814 814 !! with partial steps on 3d arrays ( i, j, k ). 815 815 !! 816 !! w-level: gdepw (i,j,k) = fsdep(k)817 !! e3w (i,j,k) = dk(fsdep)(k) = fse3(i,j,k)818 !! t-level: gdept (i,j,k) = fsdep(k+0.5)819 !! e3t (i,j,k) = dk(fsdep)(k+0.5) = fse3(i,j,k+0.5)816 !! w-level: gdepw_0(i,j,k) = gdep(k) 817 !! e3w_0(i,j,k) = dk(gdep)(k) = e3(i,j,k) 818 !! t-level: gdept_0(i,j,k) = gdep(k+0.5) 819 !! e3t_0(i,j,k) = dk(gdep)(k+0.5) = e3(i,j,k+0.5) 820 820 !! 821 821 !! With the help of the bathymetric file ( bathymetry_depth_ORCA_R2.nc), … … 825 825 !! - bathy = 0 => mbathy = 0 826 826 !! - 1 < mbathy < jpkm1 827 !! - bathy > gdepw (jpk) => mbathy = jpkm1827 !! - bathy > gdepw_0(jpk) => mbathy = jpkm1 828 828 !! 829 829 !! Then, for each case, we find the new depth at t- and w- levels … … 838 838 !! 839 839 !! c a u t i o n : gdept_1d, gdepw_1d and e3._1d are positives 840 !! - - - - - - - gdept , gdepwand e3. are positives840 !! - - - - - - - gdept_0, gdepw_0 and e3. are positives 841 841 !! 842 842 !! Reference : Pacanowsky & Gnanadesikan 1997, Mon. Wea. Rev., 126, 3248-3270. … … 891 891 ! Scale factors and depth at T- and W-points 892 892 DO jk = 1, jpk ! intitialization to the reference z-coordinate 893 gdept (:,:,jk) = gdept_1d(jk)894 gdepw (:,:,jk) = gdepw_1d(jk)895 e3t (:,:,jk) = e3t_1d (jk)896 e3w (:,:,jk) = e3w_1d (jk)893 gdept_0(:,:,jk) = gdept_1d(jk) 894 gdepw_0(:,:,jk) = gdepw_1d(jk) 895 e3t_0 (:,:,jk) = e3t_1d (jk) 896 e3w_0 (:,:,jk) = e3w_1d (jk) 897 897 END DO 898 898 ! … … 906 906 ze3tp = bathy(ji,jj) - gdepw_1d(ik) 907 907 ze3wp = 0.5_wp * e3w_1d(ik) * ( 1._wp + ( ze3tp/e3t_1d(ik) ) ) 908 e3t (ji,jj,ik ) = ze3tp909 e3t (ji,jj,ik+1) = ze3tp910 e3w (ji,jj,ik ) = ze3wp911 e3w (ji,jj,ik+1) = ze3tp912 gdepw (ji,jj,ik+1) = zdepwp913 gdept (ji,jj,ik ) = gdept_1d(ik-1) + ze3wp914 gdept (ji,jj,ik+1) = gdept(ji,jj,ik) + ze3tp908 e3t_0(ji,jj,ik ) = ze3tp 909 e3t_0(ji,jj,ik+1) = ze3tp 910 e3w_0(ji,jj,ik ) = ze3wp 911 e3w_0(ji,jj,ik+1) = ze3tp 912 gdepw_0(ji,jj,ik+1) = zdepwp 913 gdept_0(ji,jj,ik ) = gdept_1d(ik-1) + ze3wp 914 gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + ze3tp 915 915 ! 916 916 ELSE ! standard case 917 IF( bathy(ji,jj) <= gdepw_1d(ik+1) ) THEN ; gdepw (ji,jj,ik+1) = bathy(ji,jj)918 ELSE ; gdepw (ji,jj,ik+1) = gdepw_1d(ik+1)917 IF( bathy(ji,jj) <= gdepw_1d(ik+1) ) THEN ; gdepw_0(ji,jj,ik+1) = bathy(ji,jj) 918 ELSE ; gdepw_0(ji,jj,ik+1) = gdepw_1d(ik+1) 919 919 ENDIF 920 920 !gm Bug? check the gdepw_1d 921 921 ! ... on ik 922 gdept (ji,jj,ik) = gdepw_1d(ik) + ( gdepw(ji,jj,ik+1) - gdepw_1d(ik) ) &923 & * ((gdept_1d( ik ) - gdepw_1d(ik) ) &924 & / ( gdepw_1d( ik+1) - gdepw_1d(ik) ))925 e3t (ji,jj,ik) = e3t_1d (ik) * ( gdepw (ji,jj,ik+1) - gdepw_1d(ik) )&926 & 927 e3w (ji,jj,ik) = 0.5_wp * ( gdepw(ji,jj,ik+1) + gdepw_1d(ik+1) - 2._wp * gdepw_1d(ik) ) &922 gdept_0(ji,jj,ik) = gdepw_1d(ik) + ( gdepw_0 (ji,jj,ik+1) - gdepw_1d(ik) ) & 923 & * ((gdept_1d( ik ) - gdepw_1d(ik) ) & 924 & / ( gdepw_1d( ik+1) - gdepw_1d(ik) )) 925 e3t_0(ji,jj,ik) = e3t_1d (ik) * ( gdepw_0 (ji,jj,ik+1) - gdepw_1d(ik) ) & 926 & / ( gdepw_1d( ik+1) - gdepw_1d(ik) ) 927 e3w_0(ji,jj,ik) = 0.5_wp * ( gdepw_0(ji,jj,ik+1) + gdepw_1d(ik+1) - 2._wp * gdepw_1d(ik) ) & 928 928 & * ( e3w_1d(ik) / ( gdepw_1d(ik+1) - gdepw_1d(ik) ) ) 929 929 ! ... on ik+1 930 e3w (ji,jj,ik+1) = e3t(ji,jj,ik)931 e3t (ji,jj,ik+1) = e3t(ji,jj,ik)932 gdept (ji,jj,ik+1) = gdept(ji,jj,ik) + e3t(ji,jj,ik)930 e3w_0 (ji,jj,ik+1) = e3t_0 (ji,jj,ik) 931 e3t_0 (ji,jj,ik+1) = e3t_0 (ji,jj,ik) 932 gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + e3t_0(ji,jj,ik) 933 933 ENDIF 934 934 ENDIF … … 941 941 ik = mbathy(ji,jj) 942 942 IF( ik > 0 ) THEN ! ocean point only 943 e3tp (ji,jj) = e3t (ji,jj,ik)944 e3wp (ji,jj) = e3w (ji,jj,ik)943 e3tp (ji,jj) = e3t_0(ji,jj,ik) 944 e3wp (ji,jj) = e3w_0(ji,jj,ik) 945 945 ! test 946 zdiff= gdepw (ji,jj,ik+1) - gdept(ji,jj,ik )946 zdiff= gdepw_0(ji,jj,ik+1) - gdept_0(ji,jj,ik ) 947 947 IF( zdiff <= 0._wp .AND. lwp ) THEN 948 948 it = it + 1 949 949 WRITE(numout,*) ' it = ', it, ' ik = ', ik, ' (i,j) = ', ji, jj 950 950 WRITE(numout,*) ' bathy = ', bathy(ji,jj) 951 WRITE(numout,*) ' gdept = ', gdept(ji,jj,ik), ' gdepw = ', gdepw(ji,jj,ik+1), ' zdiff = ', zdiff952 WRITE(numout,*) ' e3tp = ', e3t (ji,jj,ik), ' e3wp = ', e3w(ji,jj,ik )951 WRITE(numout,*) ' gdept_0 = ', gdept_0(ji,jj,ik), ' gdepw_0 = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff 952 WRITE(numout,*) ' e3tp = ', e3t_0 (ji,jj,ik), ' e3wp = ', e3w_0 (ji,jj,ik ) 953 953 ENDIF 954 954 ENDIF … … 958 958 ! Scale factors and depth at U-, V-, UW and VW-points 959 959 DO jk = 1, jpk ! initialisation to z-scale factors 960 e3u (:,:,jk) = e3t_1d(jk)961 e3v (:,:,jk) = e3t_1d(jk)962 e3uw (:,:,jk) = e3w_1d(jk)963 e3vw (:,:,jk) = e3w_1d(jk)960 e3u_0 (:,:,jk) = e3t_1d(jk) 961 e3v_0 (:,:,jk) = e3t_1d(jk) 962 e3uw_0(:,:,jk) = e3w_1d(jk) 963 e3vw_0(:,:,jk) = e3w_1d(jk) 964 964 END DO 965 965 DO jk = 1,jpk ! Computed as the minimum of neighbooring scale factors 966 966 DO jj = 1, jpjm1 967 967 DO ji = 1, fs_jpim1 ! vector opt. 968 e3u (ji,jj,jk) = MIN( e3t(ji,jj,jk), e3t(ji+1,jj,jk) )969 e3v (ji,jj,jk) = MIN( e3t(ji,jj,jk), e3t(ji,jj+1,jk) )970 e3uw (ji,jj,jk) = MIN( e3w(ji,jj,jk), e3w(ji+1,jj,jk) )971 e3vw (ji,jj,jk) = MIN( e3w(ji,jj,jk), e3w(ji,jj+1,jk) )972 END DO 973 END DO 974 END DO 975 CALL lbc_lnk( e3u , 'U', 1._wp ) ; CALL lbc_lnk( e3uw, 'U', 1._wp ) ! lateral boundary conditions976 CALL lbc_lnk( e3v , 'V', 1._wp ) ; CALL lbc_lnk( e3vw, 'V', 1._wp )968 e3u_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji+1,jj,jk) ) 969 e3v_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji,jj+1,jk) ) 970 e3uw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji+1,jj,jk) ) 971 e3vw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji,jj+1,jk) ) 972 END DO 973 END DO 974 END DO 975 CALL lbc_lnk( e3u_0 , 'U', 1._wp ) ; CALL lbc_lnk( e3uw_0, 'U', 1._wp ) ! lateral boundary conditions 976 CALL lbc_lnk( e3v_0 , 'V', 1._wp ) ; CALL lbc_lnk( e3vw_0, 'V', 1._wp ) 977 977 ! 978 978 DO jk = 1, jpk ! set to z-scale factor if zero (i.e. along closed boundaries) 979 WHERE( e3u (:,:,jk) == 0._wp ) e3u(:,:,jk) = e3t_1d(jk)980 WHERE( e3v (:,:,jk) == 0._wp ) e3v(:,:,jk) = e3t_1d(jk)981 WHERE( e3uw (:,:,jk) == 0._wp ) e3uw(:,:,jk) = e3w_1d(jk)982 WHERE( e3vw (:,:,jk) == 0._wp ) e3vw(:,:,jk) = e3w_1d(jk)979 WHERE( e3u_0 (:,:,jk) == 0._wp ) e3u_0 (:,:,jk) = e3t_1d(jk) 980 WHERE( e3v_0 (:,:,jk) == 0._wp ) e3v_0 (:,:,jk) = e3t_1d(jk) 981 WHERE( e3uw_0(:,:,jk) == 0._wp ) e3uw_0(:,:,jk) = e3w_1d(jk) 982 WHERE( e3vw_0(:,:,jk) == 0._wp ) e3vw_0(:,:,jk) = e3w_1d(jk) 983 983 END DO 984 984 985 985 ! Scale factor at F-point 986 986 DO jk = 1, jpk ! initialisation to z-scale factors 987 e3f (:,:,jk) = e3t_1d(jk)987 e3f_0(:,:,jk) = e3t_1d(jk) 988 988 END DO 989 989 DO jk = 1, jpk ! Computed as the minimum of neighbooring V-scale factors 990 990 DO jj = 1, jpjm1 991 991 DO ji = 1, fs_jpim1 ! vector opt. 992 e3f (ji,jj,jk) = MIN( e3v(ji,jj,jk), e3v(ji+1,jj,jk) )993 END DO 994 END DO 995 END DO 996 CALL lbc_lnk( e3f , 'F', 1._wp ) ! Lateral boundary conditions992 e3f_0(ji,jj,jk) = MIN( e3v_0(ji,jj,jk), e3v_0(ji+1,jj,jk) ) 993 END DO 994 END DO 995 END DO 996 CALL lbc_lnk( e3f_0, 'F', 1._wp ) ! Lateral boundary conditions 997 997 ! 998 998 DO jk = 1, jpk ! set to z-scale factor if zero (i.e. along closed boundaries) 999 WHERE( e3f (:,:,jk) == 0._wp ) e3f(:,:,jk) = e3t_1d(jk)999 WHERE( e3f_0(:,:,jk) == 0._wp ) e3f_0(:,:,jk) = e3t_1d(jk) 1000 1000 END DO 1001 1001 !!gm bug ? : must be a do loop with mj0,mj1 1002 1002 ! 1003 e3t (:,mj0(1),:) = e3t(:,mj0(2),:) ! we duplicate factor scales for jj = 1 and jj = 21004 e3w (:,mj0(1),:) = e3w(:,mj0(2),:)1005 e3u (:,mj0(1),:) = e3u(:,mj0(2),:)1006 e3v (:,mj0(1),:) = e3v(:,mj0(2),:)1007 e3f (:,mj0(1),:) = e3f(:,mj0(2),:)1003 e3t_0(:,mj0(1),:) = e3t_0(:,mj0(2),:) ! we duplicate factor scales for jj = 1 and jj = 2 1004 e3w_0(:,mj0(1),:) = e3w_0(:,mj0(2),:) 1005 e3u_0(:,mj0(1),:) = e3u_0(:,mj0(2),:) 1006 e3v_0(:,mj0(1),:) = e3v_0(:,mj0(2),:) 1007 e3f_0(:,mj0(1),:) = e3f_0(:,mj0(2),:) 1008 1008 1009 1009 ! Control of the sign 1010 IF( MINVAL( e3t (:,:,:) ) <= 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r e3t<= 0' )1011 IF( MINVAL( e3w (:,:,:) ) <= 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r e3w<= 0' )1012 IF( MINVAL( gdept (:,:,:) ) < 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r gdepw< 0' )1013 IF( MINVAL( gdepw (:,:,:) ) < 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r gdepw< 0' )1010 IF( MINVAL( e3t_0 (:,:,:) ) <= 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r e3t_0 <= 0' ) 1011 IF( MINVAL( e3w_0 (:,:,:) ) <= 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r e3w_0 <= 0' ) 1012 IF( MINVAL( gdept_0(:,:,:) ) < 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r gdept_0 < 0' ) 1013 IF( MINVAL( gdepw_0(:,:,:) ) < 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r gdepw_0 < 0' ) 1014 1014 1015 ! Compute gdep3w (vertical sum of e3w)1016 gdep3w (:,:,1) = 0.5_wp * e3w(:,:,1)1015 ! Compute gdep3w_0 (vertical sum of e3w) 1016 gdep3w_0(:,:,1) = 0.5_wp * e3w_0(:,:,1) 1017 1017 DO jk = 2, jpk 1018 gdep3w (:,:,jk) = gdep3w(:,:,jk-1) + e3w(:,:,jk)1018 gdep3w_0(:,:,jk) = gdep3w_0(:,:,jk-1) + e3w_0(:,:,jk) 1019 1019 END DO 1020 1020 … … 1025 1025 DO ji = 1, jpi 1026 1026 ik = MAX( mbathy(ji,jj), 1 ) 1027 zprt(ji,jj,1) = e3t (ji,jj,ik)1028 zprt(ji,jj,2) = e3w (ji,jj,ik)1029 zprt(ji,jj,3) = e3u (ji,jj,ik)1030 zprt(ji,jj,4) = e3v (ji,jj,ik)1031 zprt(ji,jj,5) = e3f (ji,jj,ik)1032 zprt(ji,jj,6) = gdep3w (ji,jj,ik)1027 zprt(ji,jj,1) = e3t_0 (ji,jj,ik) 1028 zprt(ji,jj,2) = e3w_0 (ji,jj,ik) 1029 zprt(ji,jj,3) = e3u_0 (ji,jj,ik) 1030 zprt(ji,jj,4) = e3v_0 (ji,jj,ik) 1031 zprt(ji,jj,5) = e3f_0 (ji,jj,ik) 1032 zprt(ji,jj,6) = gdep3w_0(ji,jj,ik) 1033 1033 END DO 1034 1034 END DO … … 1356 1356 ENDIF 1357 1357 1358 CALL lbc_lnk( e3t , 'T', 1._wp )1359 CALL lbc_lnk( e3u , 'U', 1._wp )1360 CALL lbc_lnk( e3v , 'V', 1._wp )1361 CALL lbc_lnk( e3f , 'F', 1._wp )1362 CALL lbc_lnk( e3w , 'W', 1._wp )1363 CALL lbc_lnk( e3uw , 'U', 1._wp )1364 CALL lbc_lnk( e3vw , 'V', 1._wp )1365 1366 fsdepw(:,:,:) = gdepw (:,:,:)1367 fsde3w(:,:,:) = gdep3w (:,:,:)1368 ! 1369 where (e3t (:,:,:).eq.0.0) e3t(:,:,:) = 1.01370 where (e3u (:,:,:).eq.0.0) e3u(:,:,:) = 1.01371 where (e3v (:,:,:).eq.0.0) e3v(:,:,:) = 1.01372 where (e3f (:,:,:).eq.0.0) e3f(:,:,:) = 1.01373 where (e3w (:,:,:).eq.0.0) e3w(:,:,:) = 1.01374 where (e3uw (:,:,:).eq.0.0) e3uw(:,:,:) = 1.01375 where (e3vw (:,:,:).eq.0.0) e3vw(:,:,:) = 1.01376 1377 1378 fsdept(:,:,:) = gdept (:,:,:)1379 fsdepw(:,:,:) = gdepw (:,:,:)1380 fsde3w(:,:,:) = gdep3w (:,:,:)1381 fse3t (:,:,:) = e3t (:,:,:)1382 fse3u (:,:,:) = e3u (:,:,:)1383 fse3v (:,:,:) = e3v (:,:,:)1384 fse3f (:,:,:) = e3f (:,:,:)1385 fse3w (:,:,:) = e3w (:,:,:)1386 fse3uw(:,:,:) = e3uw (:,:,:)1387 fse3vw(:,:,:) = e3vw (:,:,:)1358 CALL lbc_lnk( e3t_0 , 'T', 1._wp ) 1359 CALL lbc_lnk( e3u_0 , 'U', 1._wp ) 1360 CALL lbc_lnk( e3v_0 , 'V', 1._wp ) 1361 CALL lbc_lnk( e3f_0 , 'F', 1._wp ) 1362 CALL lbc_lnk( e3w_0 , 'W', 1._wp ) 1363 CALL lbc_lnk( e3uw_0, 'U', 1._wp ) 1364 CALL lbc_lnk( e3vw_0, 'V', 1._wp ) 1365 1366 fsdepw(:,:,:) = gdepw_0 (:,:,:) 1367 fsde3w(:,:,:) = gdep3w_0(:,:,:) 1368 ! 1369 where (e3t_0 (:,:,:).eq.0.0) e3t_0(:,:,:) = 1.0 1370 where (e3u_0 (:,:,:).eq.0.0) e3u_0(:,:,:) = 1.0 1371 where (e3v_0 (:,:,:).eq.0.0) e3v_0(:,:,:) = 1.0 1372 where (e3f_0 (:,:,:).eq.0.0) e3f_0(:,:,:) = 1.0 1373 where (e3w_0 (:,:,:).eq.0.0) e3w_0(:,:,:) = 1.0 1374 where (e3uw_0 (:,:,:).eq.0.0) e3uw_0(:,:,:) = 1.0 1375 where (e3vw_0 (:,:,:).eq.0.0) e3vw_0(:,:,:) = 1.0 1376 1377 1378 fsdept(:,:,:) = gdept_0 (:,:,:) 1379 fsdepw(:,:,:) = gdepw_0 (:,:,:) 1380 fsde3w(:,:,:) = gdep3w_0(:,:,:) 1381 fse3t (:,:,:) = e3t_0 (:,:,:) 1382 fse3u (:,:,:) = e3u_0 (:,:,:) 1383 fse3v (:,:,:) = e3v_0 (:,:,:) 1384 fse3f (:,:,:) = e3f_0 (:,:,:) 1385 fse3w (:,:,:) = e3w_0 (:,:,:) 1386 fse3uw(:,:,:) = e3uw_0 (:,:,:) 1387 fse3vw(:,:,:) = e3vw_0 (:,:,:) 1388 1388 !! 1389 1389 ! HYBRID : … … 1400 1400 1401 1401 IF( nprint == 1 .AND. lwp ) THEN ! min max values over the local domain 1402 WRITE(numout,*) ' MIN val mbathy ', MINVAL( mbathy(:,:) ), ' MAX ', MAXVAL( mbathy(:,:) )1403 WRITE(numout,*) ' MIN val depth t ', MINVAL( fsdept(:,:,:) ), &1404 & ' w ', MINVAL( fsdepw(:,:,:) ), '3w ' , MINVAL( fsde3w(:,:,:) )1405 WRITE(numout,*) ' MIN val e3 t ', MINVAL( fse3t (:,:,:) ), ' f ' , MINVAL( fse3f(:,:,:) ), &1406 & ' u ', MINVAL( fse3u (:,:,:) ), ' u ' , MINVAL( fse3v(:,:,:) ), &1407 & ' uw', MINVAL( fse3uw(:,:,:) ), ' vw' , MINVAL( fse3vw(:,:,:) ), &1408 & ' w ', MINVAL( fse3w(:,:,:) )1409 1410 WRITE(numout,*) ' MAX val depth t ', MAXVAL( fsdept(:,:,:) ), &1411 & ' w ', MAXVAL( fsdepw(:,:,:) ), '3w ' , MAXVAL( fsde3w(:,:,:) )1412 WRITE(numout,*) ' MAX val e3 t ', MAXVAL( fse3t (:,:,:) ), ' f ' , MAXVAL( fse3f(:,:,:) ), &1413 & ' u ', MAXVAL( fse3u (:,:,:) ), ' u ' , MAXVAL( fse3v(:,:,:) ), &1414 & ' uw', MAXVAL( fse3uw(:,:,:) ), ' vw' , MAXVAL( fse3vw(:,:,:) ), &1415 & ' w ', MAXVAL( fse3w(:,:,:) )1402 WRITE(numout,*) ' MIN val mbathy ', MINVAL( mbathy(:,:) ), ' MAX ', MAXVAL( mbathy(:,:) ) 1403 WRITE(numout,*) ' MIN val depth t ', MINVAL( gdept_0(:,:,:) ), & 1404 & ' w ', MINVAL( gdepw_0(:,:,:) ), '3w ' , MINVAL( gdep3w_0(:,:,:) ) 1405 WRITE(numout,*) ' MIN val e3 t ', MINVAL( e3t_0 (:,:,:) ), ' f ' , MINVAL( e3f_0 (:,:,:) ), & 1406 & ' u ', MINVAL( e3u_0 (:,:,:) ), ' u ' , MINVAL( e3v_0 (:,:,:) ), & 1407 & ' uw', MINVAL( e3uw_0 (:,:,:) ), ' vw' , MINVAL( e3vw_0 (:,:,:) ), & 1408 & ' w ', MINVAL( e3w_0 (:,:,:) ) 1409 1410 WRITE(numout,*) ' MAX val depth t ', MAXVAL( gdept_0(:,:,:) ), & 1411 & ' w ', MAXVAL( gdepw_0(:,:,:) ), '3w ' , MAXVAL( gdep3w_0(:,:,:) ) 1412 WRITE(numout,*) ' MAX val e3 t ', MAXVAL( e3t_0 (:,:,:) ), ' f ' , MAXVAL( e3f_0 (:,:,:) ), & 1413 & ' u ', MAXVAL( e3u_0 (:,:,:) ), ' u ' , MAXVAL( e3v_0 (:,:,:) ), & 1414 & ' uw', MAXVAL( e3uw_0 (:,:,:) ), ' vw' , MAXVAL( e3vw_0 (:,:,:) ), & 1415 & ' w ', MAXVAL( e3w_0 (:,:,:) ) 1416 1416 ENDIF 1417 1417 ! END DO … … 1420 1420 WRITE(numout,*) ' domzgr: vertical coordinates : point (1,1,k) bathy = ', bathy(1,1), hbatt(1,1) 1421 1421 WRITE(numout,*) ' ~~~~~~ --------------------' 1422 WRITE(numout,"(9x,' level gdept gdepw gde3w e3t e3w')")1423 WRITE(numout,"(10x,i4,4f9.2)") ( jk, fsdept(1,1,jk), fsdepw(1,1,jk), &1424 & fse3t (1,1,jk), fse3w (1,1,jk), jk=1,jpk )1422 WRITE(numout,"(9x,' level gdept_0 gdepw_0 e3t_0 e3w_0')") 1423 WRITE(numout,"(10x,i4,4f9.2)") ( jk, gdept_0(1,1,jk), gdepw_0(1,1,jk), & 1424 & e3t_0 (1,1,jk) , e3w_0 (1,1,jk) , jk=1,jpk ) 1425 1425 DO jj = mj0(20), mj1(20) 1426 1426 DO ji = mi0(20), mi1(20) … … 1428 1428 WRITE(numout,*) ' domzgr: vertical coordinates : point (20,20,k) bathy = ', bathy(ji,jj), hbatt(ji,jj) 1429 1429 WRITE(numout,*) ' ~~~~~~ --------------------' 1430 WRITE(numout,"(9x,' level gdept gdepw gde3w e3t e3w')")1431 WRITE(numout,"(10x,i4,4f9.2)") ( jk, fsdept(ji,jj,jk), fsdepw(ji,jj,jk), &1432 & fse3t (ji,jj,jk), fse3w (ji,jj,jk), jk=1,jpk )1430 WRITE(numout,"(9x,' level gdept_0 gdepw_0 e3t_0 e3w_0')") 1431 WRITE(numout,"(10x,i4,4f9.2)") ( jk, gdept_0(ji,jj,jk), gdepw_0(ji,jj,jk), & 1432 & e3t_0 (ji,jj,jk) , e3w_0 (ji,jj,jk) , jk=1,jpk ) 1433 1433 END DO 1434 1434 END DO … … 1438 1438 WRITE(numout,*) ' domzgr: vertical coordinates : point (100,74,k) bathy = ', bathy(ji,jj), hbatt(ji,jj) 1439 1439 WRITE(numout,*) ' ~~~~~~ --------------------' 1440 WRITE(numout,"(9x,' level gdept gdepw gde3w e3t e3w')")1441 WRITE(numout,"(10x,i4,4f9.2)") ( jk, fsdept(ji,jj,jk), fsdepw(ji,jj,jk), &1442 & fse3t (ji,jj,jk), fse3w (ji,jj,jk), jk=1,jpk )1440 WRITE(numout,"(9x,' level gdept_0 gdepw_0 e3t_0 e3w_0')") 1441 WRITE(numout,"(10x,i4,4f9.2)") ( jk, gdept_0(ji,jj,jk), gdepw_0(ji,jj,jk), & 1442 & e3t_0 (ji,jj,jk) , e3w_0 (ji,jj,jk) , jk=1,jpk ) 1443 1443 END DO 1444 1444 END DO … … 1558 1558 zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp) 1559 1559 zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpkm1,wp) 1560 gdept (ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsigt3(ji,jj,jk)+rn_hc*zcoeft )1561 gdepw (ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsigw3(ji,jj,jk)+rn_hc*zcoefw )1562 gdep3w (ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsi3w3(ji,jj,jk)+rn_hc*zcoeft )1560 gdept_0 (ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsigt3(ji,jj,jk)+rn_hc*zcoeft ) 1561 gdepw_0 (ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsigw3(ji,jj,jk)+rn_hc*zcoefw ) 1562 gdep3w_0(ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsi3w3(ji,jj,jk)+rn_hc*zcoeft ) 1563 1563 END DO 1564 1564 ! … … 1581 1581 & / ( hbatt(ji,jj)+hbatt(ji,jj+1) ) 1582 1582 ! 1583 e3t (ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*z_esigt3 (ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )1584 e3u (ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*z_esigtu3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )1585 e3v (ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*z_esigtv3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )1586 e3f (ji,jj,jk) = ( (hbatf(ji,jj)-rn_hc)*z_esigtf3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )1583 e3t_0(ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*z_esigt3 (ji,jj,jk) + rn_hc/REAL(jpkm1,wp) ) 1584 e3u_0(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*z_esigtu3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) ) 1585 e3v_0(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*z_esigtv3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) ) 1586 e3f_0(ji,jj,jk) = ( (hbatf(ji,jj)-rn_hc)*z_esigtf3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) ) 1587 1587 ! 1588 e3w (ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*z_esigw3 (ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )1589 e3uw (ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*z_esigwu3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )1590 e3vw (ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*z_esigwv3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )1588 e3w_0 (ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*z_esigw3 (ji,jj,jk) + rn_hc/REAL(jpkm1,wp) ) 1589 e3uw_0(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*z_esigwu3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) ) 1590 e3vw_0(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*z_esigwv3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) ) 1591 1591 END DO 1592 1592 END DO … … 1686 1686 1687 1687 DO jk = 1, jpk 1688 gdept (ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsigt3(ji,jj,jk)1689 gdepw (ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsigw3(ji,jj,jk)1690 gdep3w (ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsi3w3(ji,jj,jk)1688 gdept_0 (ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsigt3(ji,jj,jk) 1689 gdepw_0 (ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsigw3(ji,jj,jk) 1690 gdep3w_0(ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsi3w3(ji,jj,jk) 1691 1691 END DO 1692 1692 … … 1710 1710 ( hbatt(ji,jj)+hbatt(ji,jj+1) ) 1711 1711 1712 e3t (ji,jj,jk)=(scosrf(ji,jj)+hbatt(ji,jj))*z_esigt3(ji,jj,jk)1713 e3u (ji,jj,jk)=(scosrf(ji,jj)+hbatu(ji,jj))*z_esigtu3(ji,jj,jk)1714 e3v (ji,jj,jk)=(scosrf(ji,jj)+hbatv(ji,jj))*z_esigtv3(ji,jj,jk)1715 e3f (ji,jj,jk)=(scosrf(ji,jj)+hbatf(ji,jj))*z_esigtf3(ji,jj,jk)1712 e3t_0(ji,jj,jk)=(scosrf(ji,jj)+hbatt(ji,jj))*z_esigt3(ji,jj,jk) 1713 e3u_0(ji,jj,jk)=(scosrf(ji,jj)+hbatu(ji,jj))*z_esigtu3(ji,jj,jk) 1714 e3v_0(ji,jj,jk)=(scosrf(ji,jj)+hbatv(ji,jj))*z_esigtv3(ji,jj,jk) 1715 e3f_0(ji,jj,jk)=(scosrf(ji,jj)+hbatf(ji,jj))*z_esigtf3(ji,jj,jk) 1716 1716 ! 1717 e3w (ji,jj,jk)=hbatt(ji,jj)*z_esigw3(ji,jj,jk)1718 e3uw (ji,jj,jk)=hbatu(ji,jj)*z_esigwu3(ji,jj,jk)1719 e3vw (ji,jj,jk)=hbatv(ji,jj)*z_esigwv3(ji,jj,jk)1717 e3w_0(ji,jj,jk)=hbatt(ji,jj)*z_esigw3(ji,jj,jk) 1718 e3uw_0(ji,jj,jk)=hbatu(ji,jj)*z_esigwu3(ji,jj,jk) 1719 e3vw_0(ji,jj,jk)=hbatv(ji,jj)*z_esigwv3(ji,jj,jk) 1720 1720 END DO 1721 1721 … … 1723 1723 ENDDO 1724 1724 ! 1725 CALL lbc_lnk(e3t ,'T',1.) ; CALL lbc_lnk(e3u,'T',1.)1726 CALL lbc_lnk(e3v ,'T',1.) ; CALL lbc_lnk(e3f,'T',1.)1727 CALL lbc_lnk(e3w ,'T',1.)1728 CALL lbc_lnk(e3uw ,'T',1.) ; CALL lbc_lnk(e3vw,'T',1.)1725 CALL lbc_lnk(e3t_0 ,'T',1.) ; CALL lbc_lnk(e3u_0 ,'T',1.) 1726 CALL lbc_lnk(e3v_0 ,'T',1.) ; CALL lbc_lnk(e3f_0 ,'T',1.) 1727 CALL lbc_lnk(e3w_0 ,'T',1.) 1728 CALL lbc_lnk(e3uw_0,'T',1.) ; CALL lbc_lnk(e3vw_0,'T',1.) 1729 1729 ! 1730 1730 ! ! ============= … … 1784 1784 zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp) 1785 1785 zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpkm1,wp) 1786 gdept (:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*z_gsigt(jk) + hift(:,:)*zcoeft )1787 gdepw (:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*z_gsigw(jk) + hift(:,:)*zcoefw )1788 gdep3w (:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*z_gsi3w(jk) + hift(:,:)*zcoeft )1786 gdept_0 (:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*z_gsigt(jk) + hift(:,:)*zcoeft ) 1787 gdepw_0 (:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*z_gsigw(jk) + hift(:,:)*zcoefw ) 1788 gdep3w_0(:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*z_gsi3w(jk) + hift(:,:)*zcoeft ) 1789 1789 END DO 1790 1790 !!gm: e3uw, e3vw can be suppressed (modif in dynzdf, dynzdf_iso, zdfbfr) (save 2 3D arrays) … … 1792 1792 DO ji = 1, jpi 1793 1793 DO jk = 1, jpk 1794 e3t (ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*z_esigt(jk) + hift(ji,jj)/REAL(jpkm1,wp) )1795 e3u (ji,jj,jk) = ( (hbatu(ji,jj)-hifu(ji,jj))*z_esigt(jk) + hifu(ji,jj)/REAL(jpkm1,wp) )1796 e3v (ji,jj,jk) = ( (hbatv(ji,jj)-hifv(ji,jj))*z_esigt(jk) + hifv(ji,jj)/REAL(jpkm1,wp) )1797 e3f (ji,jj,jk) = ( (hbatf(ji,jj)-hiff(ji,jj))*z_esigt(jk) + hiff(ji,jj)/REAL(jpkm1,wp) )1794 e3t_0(ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*z_esigt(jk) + hift(ji,jj)/REAL(jpkm1,wp) ) 1795 e3u_0(ji,jj,jk) = ( (hbatu(ji,jj)-hifu(ji,jj))*z_esigt(jk) + hifu(ji,jj)/REAL(jpkm1,wp) ) 1796 e3v_0(ji,jj,jk) = ( (hbatv(ji,jj)-hifv(ji,jj))*z_esigt(jk) + hifv(ji,jj)/REAL(jpkm1,wp) ) 1797 e3f_0(ji,jj,jk) = ( (hbatf(ji,jj)-hiff(ji,jj))*z_esigt(jk) + hiff(ji,jj)/REAL(jpkm1,wp) ) 1798 1798 ! 1799 e3w (ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*z_esigw(jk) + hift(ji,jj)/REAL(jpkm1,wp) )1800 e3uw (ji,jj,jk) = ( (hbatu(ji,jj)-hifu(ji,jj))*z_esigw(jk) + hifu(ji,jj)/REAL(jpkm1,wp) )1801 e3vw (ji,jj,jk) = ( (hbatv(ji,jj)-hifv(ji,jj))*z_esigw(jk) + hifv(ji,jj)/REAL(jpkm1,wp) )1799 e3w_0 (ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*z_esigw(jk) + hift(ji,jj)/REAL(jpkm1,wp) ) 1800 e3uw_0(ji,jj,jk) = ( (hbatu(ji,jj)-hifu(ji,jj))*z_esigw(jk) + hifu(ji,jj)/REAL(jpkm1,wp) ) 1801 e3vw_0(ji,jj,jk) = ( (hbatv(ji,jj)-hifv(ji,jj))*z_esigw(jk) + hifv(ji,jj)/REAL(jpkm1,wp) ) 1802 1802 END DO 1803 1803 END DO -
branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr_substitute.h90
r2528 r3865 8 8 !! 3.1 ! 2009-02 (G. Madec, M. Leclair) pure z* coordinate 9 9 !!---------------------------------------------------------------------- 10 ! reference for s- or zps-coordinate (3D no time dependency) 11 # define fsdept_0(i,j,k) gdept(i,j,k) 12 # define fsdepw_0(i,j,k) gdepw(i,j,k) 13 # define fsde3w_0(i,j,k) gdep3w(i,j,k) 14 # define fse3t_0(i,j,k) e3t(i,j,k) 15 # define fse3u_0(i,j,k) e3u(i,j,k) 16 # define fse3v_0(i,j,k) e3v(i,j,k) 17 # define fse3f_0(i,j,k) e3f(i,j,k) 18 # define fse3w_0(i,j,k) e3w(i,j,k) 19 # define fse3uw_0(i,j,k) e3uw(i,j,k) 20 # define fse3vw_0(i,j,k) e3vw(i,j,k) 10 21 11 #if defined key_vvl 22 ! s* or z*-coordinate (3D + time dependency) + use of additional now arrays (..._1) 23 # define fsdept(i,j,k) gdept_1(i,j,k) 24 # define fsdepw(i,j,k) gdepw_1(i,j,k) 25 # define fsde3w(i,j,k) gdep3w_1(i,j,k) 26 # define fse3t(i,j,k) e3t_1(i,j,k) 27 # define fse3u(i,j,k) e3u_1(i,j,k) 28 # define fse3v(i,j,k) e3v_1(i,j,k) 29 # define fse3f(i,j,k) e3f_1(i,j,k) 30 # define fse3w(i,j,k) e3w_1(i,j,k) 31 # define fse3uw(i,j,k) e3uw_1(i,j,k) 32 # define fse3vw(i,j,k) e3vw_1(i,j,k) 12 ! s* or z*-coordinate (3D + time dependency) + use of additional now arrays (..._n) 33 13 34 14 # define fse3t_b(i,j,k) e3t_b(i,j,k) 35 15 # define fse3u_b(i,j,k) e3u_b(i,j,k) 36 16 # define fse3v_b(i,j,k) e3v_b(i,j,k) 37 # define fse3uw_b(i,j,k) (fse3uw_0(i,j,k)*(1.+sshu_b(i,j)*muu(i,j,k)))38 # define fse3vw_b(i,j,k) (fse3vw_0(i,j,k)*(1.+sshv_b(i,j)*muv(i,j,k)))17 # define fse3uw_b(i,j,k) e3uw_b(i,j,k) 18 # define fse3vw_b(i,j,k) e3vw_b(i,j,k) 39 19 40 # define fsdept_n(i,j,k) (fsdept_0(i,j,k)*(1.+sshn(i,j)*mut(i,j,k)))41 # define fsdepw_n(i,j,k) (fsdepw_0(i,j,k)*(1.+sshn(i,j)*mut(i,j,k)))42 # define fsde3w_n(i,j,k) (fsde3w_0(i,j,k)*(1.+sshn(i,j)*mut(i,j,k))-sshn(i,j))43 # define fse3t_n(i,j,k) (fse3t_0(i,j,k)*(1.+sshn(i,j)*mut(i,j,k)))44 # define fse3u_n(i,j,k) (fse3u_0(i,j,k)*(1.+sshu_n(i,j)*muu(i,j,k)))45 # define fse3v_n(i,j,k) (fse3v_0(i,j,k)*(1.+sshv_n(i,j)*muv(i,j,k)))46 # define fse3f_n(i,j,k) (fse3f_0(i,j,k)*(1.+sshf_n(i,j)*muf(i,j,k)))47 # define fse3w_n(i,j,k) (fse3w_0(i,j,k)*(1.+sshn(i,j)*mut(i,j,k)))48 # define fse3uw_n(i,j,k) (fse3uw_0(i,j,k)*(1.+sshu_n(i,j)*muu(i,j,k)))49 # define fse3vw_n(i,j,k) (fse3vw_0(i,j,k)*(1.+sshv_n(i,j)*muv(i,j,k)))20 # define fsdept_n(i,j,k) gdept_n(i,j,k) 21 # define fsdepw_n(i,j,k) gdepw_n(i,j,k) 22 # define fsde3w_n(i,j,k) gdep3w_n(i,j,k) 23 # define fse3t_n(i,j,k) e3t_n(i,j,k) 24 # define fse3u_n(i,j,k) e3u_n(i,j,k) 25 # define fse3v_n(i,j,k) e3v_n(i,j,k) 26 # define fse3f_n(i,j,k) e3f_n(i,j,k) 27 # define fse3w_n(i,j,k) e3w_n(i,j,k) 28 # define fse3uw_n(i,j,k) e3uw_n(i,j,k) 29 # define fse3vw_n(i,j,k) e3vw_n(i,j,k) 50 30 51 # define fse3t_m(i,j,k) (fse3t_0(i,j,k)*(1.+ssh_m(i,j)*mut(i,j,k))) 31 # define fse3t_a(i,j,k) e3t_a(i,j,k) 32 # define fse3u_a(i,j,k) e3u_a(i,j,k) 33 # define fse3v_a(i,j,k) e3v_a(i,j,k) 52 34 53 # define fse3t_a(i,j,k) (fse3t_0(i,j,k)*(1.+ssha(i,j)*mut(i,j,k))) 54 # define fse3u_a(i,j,k) (fse3u_0(i,j,k)*(1.+sshu_a(i,j)*muu(i,j,k))) 55 # define fse3v_a(i,j,k) (fse3v_0(i,j,k)*(1.+sshv_a(i,j)*muv(i,j,k))) 35 # define fse3t_m(i,j) e3t_m(i,j) 36 37 ! This part should be removed one day ... 38 ! ... In that case all occurence of the above statement functions 39 ! have to be replaced in the code by xxx_n 40 # define fsdept(i,j,k) gdept_n(i,j,k) 41 # define fsdepw(i,j,k) gdepw_n(i,j,k) 42 # define fsde3w(i,j,k) gdep3w_n(i,j,k) 43 # define fse3t(i,j,k) e3t_n(i,j,k) 44 # define fse3u(i,j,k) e3u_n(i,j,k) 45 # define fse3v(i,j,k) e3v_n(i,j,k) 46 # define fse3f(i,j,k) e3f_n(i,j,k) 47 # define fse3w(i,j,k) e3w_n(i,j,k) 48 # define fse3uw(i,j,k) e3uw_n(i,j,k) 49 # define fse3vw(i,j,k) e3vw_n(i,j,k) 56 50 57 51 #else 58 52 ! z- or s-coordinate (1D or 3D + no time dependency) use reference in all cases 59 # define fsdept(i,j,k) fsdept_0(i,j,k)60 # define fsdepw(i,j,k) fsdepw_0(i,j,k)61 # define fsde3w(i,j,k) fsde3w_0(i,j,k)62 # define fse3t(i,j,k) fse3t_0(i,j,k)63 # define fse3u(i,j,k) fse3u_0(i,j,k)64 # define fse3v(i,j,k) fse3v_0(i,j,k)65 # define fse3f(i,j,k) fse3f_0(i,j,k)66 # define fse3w(i,j,k) fse3w_0(i,j,k)67 # define fse3uw(i,j,k) fse3uw_0(i,j,k)68 # define fse3vw(i,j,k) fse3vw_0(i,j,k)69 53 70 # define fse3t_b(i,j,k) fse3t_0(i,j,k)71 # define fse3u_b(i,j,k) fse3u_0(i,j,k)72 # define fse3v_b(i,j,k) fse3v_0(i,j,k)73 # define fse3uw_b(i,j,k) fse3uw_0(i,j,k)74 # define fse3vw_b(i,j,k) fse3vw_0(i,j,k)54 # define fse3t_b(i,j,k) e3t_0(i,j,k) 55 # define fse3u_b(i,j,k) e3u_0(i,j,k) 56 # define fse3v_b(i,j,k) e3v_0(i,j,k) 57 # define fse3uw_b(i,j,k) e3uw_0(i,j,k) 58 # define fse3vw_b(i,j,k) e3vw_0(i,j,k) 75 59 76 # define fsdept_n(i,j,k) fsdept_0(i,j,k)77 # define fsdepw_n(i,j,k) fsdepw_0(i,j,k)78 # define fsde3w_n(i,j,k) fsde3w_0(i,j,k)79 # define fse3t_n(i,j,k) fse3t_0(i,j,k)80 # define fse3u_n(i,j,k) fse3u_0(i,j,k)81 # define fse3v_n(i,j,k) fse3v_0(i,j,k)82 # define fse3f_n(i,j,k) fse3f_0(i,j,k)83 # define fse3w_n(i,j,k) fse3w_0(i,j,k)84 # define fse3uw_n(i,j,k) fse3uw_0(i,j,k)85 # define fse3vw_n(i,j,k) fse3vw_0(i,j,k)60 # define fsdept_n(i,j,k) gdept_0(i,j,k) 61 # define fsdepw_n(i,j,k) gdepw_0(i,j,k) 62 # define fsde3w_n(i,j,k) gdep3w_0(i,j,k) 63 # define fse3t_n(i,j,k) e3t_0(i,j,k) 64 # define fse3u_n(i,j,k) e3u_0(i,j,k) 65 # define fse3v_n(i,j,k) e3v_0(i,j,k) 66 # define fse3f_n(i,j,k) e3f_0(i,j,k) 67 # define fse3w_n(i,j,k) e3w_0(i,j,k) 68 # define fse3uw_n(i,j,k) e3uw_0(i,j,k) 69 # define fse3vw_n(i,j,k) e3vw_0(i,j,k) 86 70 87 # define fse3t_m(i,j,k) fse3t_0(i,j,k) 71 # define fse3t_a(i,j,k) e3t_0(i,j,k) 72 # define fse3u_a(i,j,k) e3u_0(i,j,k) 73 # define fse3v_a(i,j,k) e3v_0(i,j,k) 88 74 89 # define fse3t_a(i,j,k) fse3t_0(i,j,k) 90 # define fse3u_a(i,j,k) fse3u_0(i,j,k) 91 # define fse3v_a(i,j,k) fse3v_0(i,j,k) 75 # define fse3t_m(i,j) e3t_0(i,j,1) 76 77 ! This part should be removed one day ... 78 ! ... In that case all occurence of the above statement functions 79 ! have to be replaced in the code by xxx_n 80 # define fsdept(i,j,k) gdept_0(i,j,k) 81 # define fsdepw(i,j,k) gdepw_0(i,j,k) 82 # define fsde3w(i,j,k) gdep3w_0(i,j,k) 83 # define fse3t(i,j,k) e3t_0(i,j,k) 84 # define fse3u(i,j,k) e3u_0(i,j,k) 85 # define fse3v(i,j,k) e3v_0(i,j,k) 86 # define fse3f(i,j,k) e3f_0(i,j,k) 87 # define fse3w(i,j,k) e3w_0(i,j,k) 88 # define fse3uw(i,j,k) e3uw_0(i,j,k) 89 # define fse3vw(i,j,k) e3vw_0(i,j,k) 90 92 91 #endif 93 92 !!---------------------------------------------------------------------- -
branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DOM/dtatsd.F90
r3862 r3865 222 222 DO ji = 1, jpi 223 223 DO jk = 1, jpk ! determines the intepolated T-S profiles at each (i,j) points 224 zl = fsdept_0(ji,jj,jk)224 zl = gdept_0(ji,jj,jk) 225 225 IF( zl < gdept_1d(1 ) ) THEN ! above the first level of data 226 226 ztp(jk) = ptsd(ji,jj,1 ,jp_tem) … … 260 260 ik = mbkt(ji,jj) 261 261 IF( ik > 1 ) THEN 262 zl = ( gdept_1d(ik) - fsdept_0(ji,jj,ik) ) / ( gdept_1d(ik) - gdept_1d(ik-1) )262 zl = ( gdept_1d(ik) - gdept_0(ji,jj,ik) ) / ( gdept_1d(ik) - gdept_1d(ik-1) ) 263 263 ptsd(ji,jj,ik,jp_tem) = (1.-zl) * ptsd(ji,jj,ik,jp_tem) + zl * ptsd(ji,jj,ik-1,jp_tem) 264 264 ptsd(ji,jj,ik,jp_sal) = (1.-zl) * ptsd(ji,jj,ik,jp_sal) + zl * ptsd(ji,jj,ik-1,jp_sal) -
branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DOM/istate.F90
r3862 r3865 90 90 neuler = 1 ! Set time-step indicator at nit000 (leap-frog) 91 91 CALL rst_read ! Read the restart file 92 ! ! define e3u_b, e3v_b from e3t_b read in restart file93 CALL dom_vvl_2( nit000, fse3u_b(:,:,:), fse3v_b(:,:,:) )94 92 CALL day_init ! model calendar (using both namelist and restart infos) 95 93 ELSE … … 131 129 ENDDO 132 130 ENDIF 133 ! ! define e3u_b, e3v_b from e3t_b initialized in domzgr134 CALL dom_vvl_2( nit000, fse3u_b(:,:,:), fse3v_b(:,:,:) )135 131 ! 136 132 ENDIF -
branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90
r3764 r3865 760 760 DO jj = 2, jpjm1 761 761 DO ji = 2, jpim1 762 zu(ji,jj,1) = - ( fse3u(ji,jj,1) - ssh u_n(ji,jj) * znad)763 zv(ji,jj,1) = - ( fse3v(ji,jj,1) - ssh v_n(ji,jj) * znad)762 zu(ji,jj,1) = - ( fse3u(ji,jj,1) - sshn(ji,jj) * znad) ! probable bug: changed from sshu_n for ztilde compilation 763 zv(ji,jj,1) = - ( fse3v(ji,jj,1) - sshn(ji,jj) * znad) ! probable bug: changed from sshv_n for ztilde compilation 764 764 END DO 765 765 END DO -
branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90
r3764 r3865 201 201 DO jj = 1, jpj 202 202 DO ji = 1, jpi 203 zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2.e0 * un(ji,jj,jk) + ua(ji,jj,jk) )204 zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2.e0 * vn(ji,jj,jk) + va(ji,jj,jk) )203 zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2.e0_wp * un(ji,jj,jk) + ua(ji,jj,jk) ) 204 zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2.e0_wp * vn(ji,jj,jk) + va(ji,jj,jk) ) 205 205 ! 206 206 ub(ji,jj,jk) = zuf ! ub <-- filtered velocity … … 214 214 ELSE ! Variable volume ! 215 215 ! ! ================! 216 ! Before scale factor at t-points 217 ! (used as a now filtered scale factor until the swap) 218 ! ---------------------------------------------------- 219 fse3t_b(:,:,:) = fse3t_n(:,:,:) + atfp * ( fse3t_b(:,:,:) - 2._wp * fse3t_n(:,:,:) + fse3t_a(:,:,:) ) 220 ! Add volume filter correction: compatibility with tracer advection scheme 221 ! => time filter + conservation correction (only at the first level) 222 fse3t_b(:,:,1) = fse3t_b(:,:,1) - atfp * rdt * r1_rau0 * ( emp_b(:,:) - emp(:,:) ) * tmask(:,:,1) 216 223 ! 217 DO jk = 1, jpkm1 ! Before scale factor at t-points 218 fse3t_b(:,:,jk) = fse3t_n(:,:,jk) & 219 & + atfp * ( fse3t_b(:,:,jk) + fse3t_a(:,:,jk) & 220 & - 2._wp * fse3t_n(:,:,jk) ) 221 END DO 222 zec = atfp * rdt / rau0 ! Add filter correction only at the 1st level of t-point scale factors 223 fse3t_b(:,:,1) = fse3t_b(:,:,1) - zec * ( emp_b(:,:) - emp(:,:) ) * tmask(:,:,1) 224 ! 225 IF( ln_dynadv_vec ) THEN ! vector invariant form (no thickness weighted calulation) 226 ! 227 ! ! before scale factors at u- & v-pts (computed from fse3t_b) 228 CALL dom_vvl_2( kt, fse3u_b(:,:,:), fse3v_b(:,:,:) ) 229 ! 230 DO jk = 1, jpkm1 ! Leap-Frog - Asselin filter and swap: applied on velocity 231 DO jj = 1, jpj ! -------- 224 IF( ln_dynadv_vec ) THEN 225 ! Before scale factor at (u/v)-points 226 ! ----------------------------------- 227 CALL dom_vvl_interpol( fse3t_b(:,:,:), fse3u_b(:,:,:), 'U' ) 228 CALL dom_vvl_interpol( fse3t_b(:,:,:), fse3v_b(:,:,:), 'V' ) 229 ! Leap-Frog - Asselin filter and swap: applied on velocity 230 ! ----------------------------------- 231 DO jk = 1, jpkm1 232 DO jj = 1, jpj 232 233 DO ji = 1, jpi 233 zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2. e0* un(ji,jj,jk) + ua(ji,jj,jk) )234 zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2. e0* vn(ji,jj,jk) + va(ji,jj,jk) )234 zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2._wp * un(ji,jj,jk) + ua(ji,jj,jk) ) 235 zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2._wp * vn(ji,jj,jk) + va(ji,jj,jk) ) 235 236 ! 236 237 ub(ji,jj,jk) = zuf ! ub <-- filtered velocity … … 242 243 END DO 243 244 ! 244 ELSE ! flux form (thickness weighted calulation) 245 ! 246 CALL dom_vvl_2( kt, ze3u_f, ze3v_f ) ! before scale factors at u- & v-pts (computed from fse3t_b) 247 ! 248 DO jk = 1, jpkm1 ! Leap-Frog - Asselin filter and swap: 249 DO jj = 1, jpj ! applied on thickness weighted velocity 245 ELSE 246 ! Temporary filtered scale factor at (u/v)-points (will become before scale factor) 247 !------------------------------------------------ 248 CALL dom_vvl_interpol( fse3t_b(:,:,:), ze3u_f, 'U' ) 249 CALL dom_vvl_interpol( fse3t_b(:,:,:), ze3v_f, 'V' ) 250 ! Leap-Frog - Asselin filter and swap: applied on thickness weighted velocity 251 ! ----------------------------------- =========================== 252 DO jk = 1, jpkm1 253 DO jj = 1, jpj 250 254 DO ji = 1, jpi ! --------------------------- 251 255 zue3a = ua(ji,jj,jk) * fse3u_a(ji,jj,jk) -
branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90
r3680 r3865 126 126 REAL(wp), POINTER, DIMENSION(:,:) :: zua, zva, zun, zvn, zun_e, zvn_e, zub_e, zvb_e 127 127 REAL(wp), POINTER, DIMENSION(:,:) :: zcu, zcv, zwx, zwy, zbfru, zbfrv, zu_sum, zv_sum 128 REAL(wp), POINTER, DIMENSION(:,:) :: zhu_b, zhv_b 128 129 !!---------------------------------------------------------------------- 129 130 ! … … 133 134 CALL wrk_alloc( jpi, jpj, zua, zva, zun, zvn, zun_e, zvn_e, zub_e, zvb_e ) 134 135 CALL wrk_alloc( jpi, jpj, zcu, zcv, zwx, zwy, zbfru, zbfrv, zu_sum, zv_sum ) 136 CALL wrk_alloc( jpi, jpj, zhu_b, zhv_b ) 135 137 ! 136 138 IF( kt == nit000 ) THEN !* initialisation … … 199 201 #endif 200 202 ! ! now trend 201 zua(ji,jj) = zua(ji,jj) + fse3u 202 zva(ji,jj) = zva(ji,jj) + fse3v 203 zua(ji,jj) = zua(ji,jj) + fse3u_n(ji,jj,jk) * ua(ji,jj,jk) * umask(ji,jj,jk) 204 zva(ji,jj) = zva(ji,jj) + fse3v_n(ji,jj,jk) * va(ji,jj,jk) * vmask(ji,jj,jk) 203 205 ! ! now velocity 204 zun(ji,jj) = zun(ji,jj) + fse3u 205 zvn(ji,jj) = zvn(ji,jj) + fse3v 206 zun(ji,jj) = zun(ji,jj) + fse3u_n(ji,jj,jk) * un(ji,jj,jk) 207 zvn(ji,jj) = zvn(ji,jj) + fse3v_n(ji,jj,jk) * vn(ji,jj,jk) 206 208 ! 207 209 #if defined key_vvl … … 215 217 END DO 216 218 END DO 219 220 ! before inverse water column height at u- and v- points 221 IF( lk_vvl ) THEN 222 zhu_b(:,:) = 0. 223 zhv_b(:,:) = 0. 224 DO jk = 1, jpk 225 zhu_b(:,:) = zhu_b(:,:) + fse3u_b(:,:,jk) * umask(:,:,jk) 226 zhv_b(:,:) = zhv_b(:,:) + fse3v_b(:,:,jk) * vmask(:,:,jk) 227 END DO 228 zhu_b(:,:) = umask(:,:,1) / ( zhu_b(:,:) + 1. - umask(:,:,1) ) 229 zhv_b(:,:) = vmask(:,:,1) / ( zhv_b(:,:) + 1. - vmask(:,:,1) ) 230 ELSE 231 zhu_b(:,:) = hur(:,:) 232 zhv_b(:,:) = hvr(:,:) 233 ENDIF 217 234 218 235 ! !* baroclinic momentum trend (remove the vertical mean trend) … … 355 372 vb_b(:,:) = vb_b(:,:) * hvr(:,:) 356 373 ENDIF 374 ub_b(:,:) = ub_b(:,:) * zhu_b(:,:) 375 vb_b(:,:) = vb_b(:,:) * zhv_b(:,:) 357 376 358 377 ! ----------------------------------------------------------------------- … … 683 702 CALL wrk_dealloc( jpi, jpj, zua, zva, zun, zvn, zun_e, zvn_e, zub_e, zvb_e ) 684 703 CALL wrk_dealloc( jpi, jpj, zcu, zcv, zwx, zwy, zbfru, zbfrv, zu_sum, zv_sum ) 704 CALL wrk_dealloc( jpi, jpj, zhu_b, zhv_b ) 685 705 ! 686 706 IF( nn_timing == 1 ) CALL timing_stop('dyn_spg_ts') … … 698 718 CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag 699 719 ! 720 REAL(wp), POINTER, DIMENSION(:,:) :: zzhu_b, zzhv_b 700 721 INTEGER :: ji, jk ! dummy loop indices 701 722 !!---------------------------------------------------------------------- … … 706 727 CALL iom_get( numror, jpdom_autoglo, 'vn_b' , vn_b (:,:) ) ! from barotropic loop 707 728 ELSE 729 CALL wrk_alloc( jpi, jpj, zzhu_b, zzhv_b ) 708 730 un_b (:,:) = 0._wp 709 731 vn_b (:,:) = 0._wp … … 712 734 DO jk = 1, jpkm1 713 735 DO ji = 1, jpij 714 un_b(ji,1) = un_b(ji,1) + fse3u (ji,1,jk) * un(ji,1,jk)715 vn_b(ji,1) = vn_b(ji,1) + fse3v (ji,1,jk) * vn(ji,1,jk)736 un_b(ji,1) = un_b(ji,1) + fse3u_n(ji,1,jk) * un(ji,1,jk) 737 vn_b(ji,1) = vn_b(ji,1) + fse3v_n(ji,1,jk) * vn(ji,1,jk) 716 738 END DO 717 739 END DO 718 740 ELSE ! No vector opt. 719 741 DO jk = 1, jpkm1 720 un_b(:,:) = un_b(:,:) + fse3u (:,:,jk) * un(:,:,jk)721 vn_b(:,:) = vn_b(:,:) + fse3v (:,:,jk) * vn(:,:,jk)742 un_b(:,:) = un_b(:,:) + fse3u_n(:,:,jk) * un(:,:,jk) 743 vn_b(:,:) = vn_b(:,:) + fse3v_n(:,:,jk) * vn(:,:,jk) 722 744 END DO 723 745 ENDIF … … 747 769 748 770 IF( lk_vvl ) THEN 749 ub_b (:,:) = ub_b(:,:) * umask(:,:,1) / ( hu_0(:,:) + sshu_b(:,:) + 1._wp - umask(:,:,1) ) 750 vb_b (:,:) = vb_b(:,:) * vmask(:,:,1) / ( hv_0(:,:) + sshv_b(:,:) + 1._wp - vmask(:,:,1) ) 751 ELSE 771 CALL wrk_alloc( jpi, jpj, zzhu_b, zzhv_b ) 772 ub_b (:,:) = 0. 773 vb_b (:,:) = 0. 774 zzhu_b(:,:) = 0. 775 zzhv_b(:,:) = 0. 776 ! vertical sum 777 IF( lk_vopt_loop ) THEN ! vector opt., forced unroll 778 DO jk = 1, jpkm1 779 DO ji = 1, jpij 780 ub_b (ji,1) = ub_b (ji,1) + fse3u_b(ji,1,jk) * ub (ji,1,jk) 781 vb_b (ji,1) = vb_b (ji,1) + fse3v_b(ji,1,jk) * vb (ji,1,jk) 782 zzhu_b(ji,1) = zhu_b(ji,1) + fse3u_b(ji,1,jk) * umask(ji,1,jk) 783 zzhv_b(ji,1) = zhv_b(ji,1) + fse3v_b(ji,1,jk) * vmask(ji,1,jk) 784 END DO 785 END DO 786 ELSE ! No vector opt. 787 DO jk = 1, jpkm1 788 ub_b (:,:) = ub_b (:,:) + fse3u_b(:,:,jk) * ub (:,:,jk) 789 vb_b (:,:) = vb_b (:,:) + fse3v_b(:,:,jk) * vb (:,:,jk) 790 zzhu_b(:,:) = zzhu_b(:,:) + fse3u_b(:,:,jk) * umask(:,:,jk) 791 zzhv_b(:,:) = zzhv_b(:,:) + fse3v_b(:,:,jk) * vmask(:,:,jk) 792 END DO 793 ENDIF 794 ub_b(:,:) = ub_b(:,:) / ( zzhu_b(:,:) + 1. - umask(:,:,1) ) 795 vb_b(:,:) = vb_b(:,:) / ( zzhv_b(:,:) + 1. - vmask(:,:,1) ) 796 CALL wrk_dealloc( jpi, jpj, zzhu_b, zzhv_b ) 797 ELSE 798 ub_b (:,:) = 0.e0 799 vb_b (:,:) = 0.e0 800 ! vertical sum 801 IF( lk_vopt_loop ) THEN ! vector opt., forced unroll 802 DO jk = 1, jpkm1 803 DO ji = 1, jpij 804 ub_b(ji,1) = ub_b(ji,1) + fse3u_b(ji,1,jk) * ub(ji,1,jk) 805 vb_b(ji,1) = vb_b(ji,1) + fse3v_b(ji,1,jk) * vb(ji,1,jk) 806 END DO 807 END DO 808 ELSE ! No vector opt. 809 DO jk = 1, jpkm1 810 ub_b(:,:) = ub_b(:,:) + fse3u_b(:,:,jk) * ub(:,:,jk) 811 vb_b(:,:) = vb_b(:,:) + fse3v_b(:,:,jk) * vb(:,:,jk) 812 END DO 813 ENDIF 752 814 ub_b(:,:) = ub_b(:,:) * hur(:,:) 753 815 vb_b(:,:) = vb_b(:,:) * hvr(:,:) -
branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90
r3764 r3865 8 8 !! - ! 2010-05 (K. Mogensen, A. Weaver, M. Martin, D. Lea) Assimilation interface 9 9 !! - ! 2010-09 (D.Storkey and E.O'Dea) bug fixes for BDY module 10 !!---------------------------------------------------------------------- 11 12 !!---------------------------------------------------------------------- 13 !! ssh_wzv : after ssh & now vertical velocity 14 !! ssh_nxt : filter ans swap the ssh arrays 10 !! 3.3 ! 2011-10 (M. Leclair) split former ssh_wzv routine and remove all vvl related work 11 !!---------------------------------------------------------------------- 12 13 !!---------------------------------------------------------------------- 14 !! ssh_nxt : after ssh 15 !! ssh_swp : filter ans swap the ssh arrays 16 !! wzv : compute now vertical velocity 15 17 !!---------------------------------------------------------------------- 16 18 USE oce ! ocean dynamics and tracers variables … … 20 22 USE divcur ! hor. divergence and curl (div & cur routines) 21 23 USE iom ! I/O library 24 USE restart ! only for lrst_oce 22 25 USE in_out_manager ! I/O manager 23 26 USE prtctl ! Print control … … 44 47 PRIVATE 45 48 46 PUBLIC ssh_wzv ! called by step.F9047 49 PUBLIC ssh_nxt ! called by step.F90 50 PUBLIC wzv ! called by step.F90 51 PUBLIC ssh_swp ! called by step.F90 48 52 49 53 !! * Substitutions … … 57 61 CONTAINS 58 62 59 SUBROUTINE ssh_ wzv( kt )60 !!---------------------------------------------------------------------- 61 !! *** ROUTINE ssh_ wzv***63 SUBROUTINE ssh_nxt( kt ) 64 !!---------------------------------------------------------------------- 65 !! *** ROUTINE ssh_nxt *** 62 66 !! 63 !! ** Purpose : compute the after ssh (ssha), the now vertical velocity 64 !! and update the now vertical coordinate (lk_vvl=T). 65 !! 66 !! ** Method : - Using the incompressibility hypothesis, the vertical 67 !! velocity is computed by integrating the horizontal divergence 68 !! from the bottom to the surface minus the scale factor evolution. 69 !! The boundary conditions are w=0 at the bottom (no flux) and. 67 !! ** Purpose : compute the after ssh (ssha) 68 !! 69 !! ** Method : - Using the incompressibility hypothesis, the ssh increment 70 !! is computed by integrating the horizontal divergence and multiply by 71 !! by the time step. 70 72 !! 71 73 !! ** action : ssha : after sea surface height 72 !! wn : now vertical velocity73 !! sshu_a, sshv_a, sshf_a : after sea surface height (lk_vvl=T)74 !! hu, hv, hur, hvr : ocean depth and its inverse at u-,v-points75 74 !! 76 75 !! Reference : Leclair, M., and G. Madec, 2009, Ocean Modelling. 77 76 !!---------------------------------------------------------------------- 78 INTEGER, INTENT(in) :: kt ! time step79 !80 INTEGER :: ji, jj, jk ! dummy loop indices81 REAL(wp) :: zcoefu, zcoefv, zcoeff, z2dt, z1_2dt, z1_rau0 ! local scalars82 REAL(wp), POINTER, DIMENSION(:,: ) :: z2d, zhdiv83 REAL(wp) , POINTER, DIMENSION(:,:,:) :: z3d84 !!---------------------------------------------------------------------- 85 ! 86 IF( nn_timing == 1 ) CALL timing_start('ssh_ wzv')87 ! 88 CALL wrk_alloc( jpi, jpj, z 2d, zhdiv )77 ! 78 REAL(wp), POINTER, DIMENSION(:,: ) :: zhdiv 79 INTEGER, INTENT(in) :: kt ! time step 80 ! 81 INTEGER :: jk ! dummy loop indice 82 REAL(wp) :: z2dt, z1_rau0 ! local scalars 83 !!---------------------------------------------------------------------- 84 ! 85 IF( nn_timing == 1 ) CALL timing_start('ssh_nxt') 86 ! 87 CALL wrk_alloc( jpi, jpj, zhdiv ) 89 88 ! 90 89 IF( kt == nit000 ) THEN 91 90 ! 92 91 IF(lwp) WRITE(numout,*) 93 IF(lwp) WRITE(numout,*) 'ssh_ wzv : after sea surface height and now vertical velocity'92 IF(lwp) WRITE(numout,*) 'ssh_nxt : after sea surface height' 94 93 IF(lwp) WRITE(numout,*) '~~~~~~~ ' 95 94 ! 96 wn(:,:,jpk) = 0._wp ! bottom boundary condition: w=0 (set once for all)97 !98 IF( lk_vvl ) THEN ! before and now Sea SSH at u-, v-, f-points (vvl case only)99 DO jj = 1, jpjm1100 DO ji = 1, jpim1 ! caution: use of Vector Opt. not possible101 zcoefu = 0.5 * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) )102 zcoefv = 0.5 * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) )103 zcoeff = 0.25 * umask(ji,jj,1) * umask(ji,jj+1,1)104 sshu_b(ji,jj) = zcoefu * ( e1t(ji ,jj) * e2t(ji ,jj) * sshb(ji ,jj) &105 & + e1t(ji+1,jj) * e2t(ji+1,jj) * sshb(ji+1,jj) )106 sshv_b(ji,jj) = zcoefv * ( e1t(ji,jj ) * e2t(ji,jj ) * sshb(ji,jj ) &107 & + e1t(ji,jj+1) * e2t(ji,jj+1) * sshb(ji,jj+1) )108 sshu_n(ji,jj) = zcoefu * ( e1t(ji ,jj) * e2t(ji ,jj) * sshn(ji ,jj) &109 & + e1t(ji+1,jj) * e2t(ji+1,jj) * sshn(ji+1,jj) )110 sshv_n(ji,jj) = zcoefv * ( e1t(ji,jj ) * e2t(ji,jj ) * sshn(ji,jj ) &111 & + e1t(ji,jj+1) * e2t(ji,jj+1) * sshn(ji,jj+1) )112 END DO113 END DO114 CALL lbc_lnk( sshu_b, 'U', 1. ) ; CALL lbc_lnk( sshu_n, 'U', 1. )115 CALL lbc_lnk( sshv_b, 'V', 1. ) ; CALL lbc_lnk( sshv_n, 'V', 1. )116 DO jj = 1, jpjm1117 DO ji = 1, jpim1 ! NO Vector Opt.118 sshf_n(ji,jj) = 0.5 * umask(ji,jj,1) * umask(ji,jj+1,1) &119 & / ( e1f(ji,jj ) * e2f(ji,jj ) ) &120 & * ( e1u(ji,jj ) * e2u(ji,jj ) * sshu_n(ji,jj ) &121 & + e1u(ji,jj+1) * e2u(ji,jj+1) * sshu_n(ji,jj+1) )122 END DO123 END DO124 CALL lbc_lnk( sshf_n, 'F', 1. )125 ENDIF126 !127 ENDIF128 129 ! !------------------------------------------!130 IF( lk_vvl ) THEN ! Regridding: Update Now Vertical coord. ! (only in vvl case)131 ! !------------------------------------------!132 DO jk = 1, jpkm1133 fsdept(:,:,jk) = fsdept_n(:,:,jk) ! now local depths stored in fsdep. arrays134 fsdepw(:,:,jk) = fsdepw_n(:,:,jk)135 fsde3w(:,:,jk) = fsde3w_n(:,:,jk)136 !137 fse3t (:,:,jk) = fse3t_n (:,:,jk) ! vertical scale factors stored in fse3. arrays138 fse3u (:,:,jk) = fse3u_n (:,:,jk)139 fse3v (:,:,jk) = fse3v_n (:,:,jk)140 fse3f (:,:,jk) = fse3f_n (:,:,jk)141 fse3w (:,:,jk) = fse3w_n (:,:,jk)142 fse3uw(:,:,jk) = fse3uw_n(:,:,jk)143 fse3vw(:,:,jk) = fse3vw_n(:,:,jk)144 END DO145 !146 hu(:,:) = hu_0(:,:) + sshu_n(:,:) ! now ocean depth (at u- and v-points)147 hv(:,:) = hv_0(:,:) + sshv_n(:,:)148 ! ! now masked inverse of the ocean depth (at u- and v-points)149 hur(:,:) = umask(:,:,1) / ( hu(:,:) + 1._wp - umask(:,:,1) )150 hvr(:,:) = vmask(:,:,1) / ( hv(:,:) + 1._wp - vmask(:,:,1) )151 !152 95 ENDIF 153 96 ! … … 162 105 zhdiv(:,:) = 0._wp 163 106 DO jk = 1, jpkm1 ! Horizontal divergence of barotropic transports 164 zhdiv(:,:) = zhdiv(:,:) + fse3t (:,:,jk) * hdivn(:,:,jk)107 zhdiv(:,:) = zhdiv(:,:) + fse3t_n(:,:,jk) * hdivn(:,:,jk) 165 108 END DO 166 109 ! ! Sea surface elevation time stepping … … 181 124 #if defined key_bdy 182 125 ssha(:,:) = ssha(:,:) * bdytmask(:,:) 183 CALL lbc_lnk( ssha, 'T', 1. ) ! absolutly compulsory !! (jmm)126 CALL lbc_lnk( ssha, 'T', 1. ) 184 127 #endif 185 128 #if defined key_asminc … … 190 133 ENDIF 191 134 #endif 192 ! ! Sea Surface Height at u-,v- and f-points (vvl case only) 193 IF( lk_vvl ) THEN ! (required only in key_vvl case) 194 DO jj = 1, jpjm1 195 DO ji = 1, jpim1 ! NO Vector Opt. 196 sshu_a(ji,jj) = 0.5 * umask(ji,jj,1) / ( e1u(ji ,jj) * e2u(ji ,jj) ) & 197 & * ( e1t(ji ,jj) * e2t(ji ,jj) * ssha(ji ,jj) & 198 & + e1t(ji+1,jj) * e2t(ji+1,jj) * ssha(ji+1,jj) ) 199 sshv_a(ji,jj) = 0.5 * vmask(ji,jj,1) / ( e1v(ji,jj ) * e2v(ji,jj ) ) & 200 & * ( e1t(ji,jj ) * e2t(ji,jj ) * ssha(ji,jj ) & 201 & + e1t(ji,jj+1) * e2t(ji,jj+1) * ssha(ji,jj+1) ) 135 136 ! !------------------------------! 137 ! ! outputs ! 138 ! !------------------------------! 139 CALL iom_put( "ssh" , sshn ) ! sea surface height 140 CALL iom_put( "ssh2", sshn(:,:) * sshn(:,:) ) ! square of sea surface height 141 ! 142 IF(ln_ctl) CALL prt_ctl( tab2d_1=ssha, clinfo1=' ssha - : ', mask1=tmask, ovlap=1 ) 143 ! 144 CALL wrk_dealloc( jpi, jpj, zhdiv ) 145 ! 146 IF( nn_timing == 1 ) CALL timing_stop('ssh_nxt') 147 ! 148 END SUBROUTINE ssh_nxt 149 150 151 SUBROUTINE wzv( kt ) 152 !!---------------------------------------------------------------------- 153 !! *** ROUTINE wzv *** 154 !! 155 !! ** Purpose : compute the now vertical velocity 156 !! 157 !! ** Method : - Using the incompressibility hypothesis, the vertical 158 !! velocity is computed by integrating the horizontal divergence 159 !! from the bottom to the surface minus the scale factor evolution. 160 !! The boundary conditions are w=0 at the bottom (no flux) and. 161 !! 162 !! ** action : wn : now vertical velocity 163 !! 164 !! Reference : Leclair, M., and G. Madec, 2009, Ocean Modelling. 165 !!---------------------------------------------------------------------- 166 ! 167 INTEGER, INTENT(in) :: kt ! time step 168 REAL(wp), POINTER, DIMENSION(:,: ) :: z2d 169 REAL(wp), POINTER, DIMENSION(:,:,:) :: z3d, zhdiv 170 ! 171 INTEGER :: ji, jj, jk ! dummy loop indices 172 REAL(wp) :: z1_2dt ! local scalars 173 !!---------------------------------------------------------------------- 174 175 IF( nn_timing == 1 ) CALL timing_start('wzv') 176 ! 177 CALL wrk_alloc( jpi, jpj, z2d ) 178 CALL wrk_alloc( jpi, jpj, jpk, z3d, zhdiv ) 179 ! 180 IF( kt == nit000 ) THEN 181 ! 182 IF(lwp) WRITE(numout,*) 183 IF(lwp) WRITE(numout,*) 'wzv : now vertical velocity ' 184 IF(lwp) WRITE(numout,*) '~~~ ' 185 ! 186 wn(:,:,jpk) = 0._wp ! bottom boundary condition: w=0 (set once for all) 187 ! 188 ENDIF 189 ! !------------------------------! 190 ! ! Now Vertical Velocity ! 191 ! !------------------------------! 192 z1_2dt = 1. / ( 2. * rdt ) ! set time step size (Euler/Leapfrog) 193 IF( neuler == 0 .AND. kt == nit000 ) z1_2dt = 1. / rdt 194 ! 195 IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN ! z_tilde and layer cases 196 DO jk = 1, jpkm1 197 ! horizontal divergence of thickness diffusion transport ( velocity multiplied by e3t) 198 ! - ML - note: computation allready done in dom_vvl_sf_nxt. Could be optimized (not critical and clearer this way) 199 DO jj = 2, jpjm1 200 DO ji = fs_2, fs_jpim1 ! vector opt. 201 zhdiv(ji,jj,jk) = r1_e12t(ji,jj) * ( un_td(ji,jj,jk) - un_td(ji-1,jj,jk) + vn_td(ji,jj,jk) - vn_td(ji,jj-1,jk) ) 202 END DO 202 203 END DO 203 204 END DO 204 CALL lbc_lnk( sshu_a, 'U', 1. ) ; CALL lbc_lnk( sshv_a, 'V', 1. ) ! Boundaries conditions 205 ENDIF 206 207 ! !------------------------------! 208 ! ! Now Vertical Velocity ! 209 ! !------------------------------! 210 z1_2dt = 1.e0 / z2dt 211 DO jk = jpkm1, 1, -1 ! integrate from the bottom the hor. divergence 212 ! - ML - need 3 lines here because replacement of fse3t by its expression yields too long lines otherwise 213 wn(:,:,jk) = wn(:,:,jk+1) - fse3t_n(:,:,jk) * hdivn(:,:,jk) & 214 & - ( fse3t_a(:,:,jk) - fse3t_b(:,:,jk) ) & 215 & * tmask(:,:,jk) * z1_2dt 205 CALL lbc_lnk(zhdiv, 'T', 1.) ! - ML - Perhaps not necessary: not used for horizontal "connexions" 206 ! ! Is it problematic to have a wrong vertical velocity in boundary cells? 207 ! ! Same question holds for hdivn. Perhaps just for security 208 DO jk = jpkm1, 1, -1 ! integrate from the bottom the hor. divergence 209 ! computation of w 210 wn(:,:,jk) = wn(:,:,jk+1) - ( fse3t_n(:,:,jk) * hdivn(:,:,jk) + zhdiv(:,:,jk) & 211 & + z1_2dt * ( fse3t_a(:,:,jk) - fse3t_b(:,:,jk) ) ) * tmask(:,:,jk) 212 END DO 213 ! IF( ln_vvl_layer ) wn(:,:,:) = 0.e0 214 ELSE ! z_star and linear free surface cases 215 DO jk = jpkm1, 1, -1 ! integrate from the bottom the hor. divergence 216 ! computation of w 217 wn(:,:,jk) = wn(:,:,jk+1) - ( fse3t_n(:,:,jk) * hdivn(:,:,jk) & 218 & + z1_2dt * ( fse3t_a(:,:,jk) - fse3t_b(:,:,jk) ) ) * tmask(:,:,jk) 219 END DO 220 ENDIF 221 216 222 #if defined key_bdy 217 223 wn(:,:,jk) = wn(:,:,jk) * bdytmask(:,:) 218 224 #endif 219 END DO220 225 221 226 ! !------------------------------! 222 227 ! ! outputs ! 223 228 ! !------------------------------! 224 CALL iom_put( "woce", wn ) ! vertical velocity 225 CALL iom_put( "ssh" , sshn ) ! sea surface height 226 CALL iom_put( "ssh2", sshn(:,:) * sshn(:,:) ) ! square of sea surface height 229 CALL iom_put( "woce", wn ) ! vertical velocity 227 230 IF( lk_diaar5 ) THEN ! vertical mass transport & its square value 228 231 ! Caution: in the VVL case, it only correponds to the baroclinic mass transport. 229 CALL wrk_alloc( jpi,jpj,jpk, z3d ) 230 z2d(:,:) = rau0 * e1t(:,:) * e2t(:,:) 232 z2d(:,:) = rau0 * e12t(:,:) 231 233 DO jk = 1, jpk 232 234 z3d(:,:,jk) = wn(:,:,jk) * z2d(:,:) … … 234 236 CALL iom_put( "w_masstr" , z3d ) 235 237 CALL iom_put( "w_masstr2", z3d(:,:,:) * z3d(:,:,:) ) 236 CALL wrk_dealloc( jpi,jpj,jpk, z3d ) 237 ENDIF 238 ! 239 IF(ln_ctl) CALL prt_ctl( tab2d_1=ssha, clinfo1=' ssha - : ', mask1=tmask, ovlap=1 ) 240 ! 241 CALL wrk_dealloc( jpi, jpj, z2d, zhdiv ) 242 ! 243 IF( nn_timing == 1 ) CALL timing_stop('ssh_wzv') 244 ! 245 END SUBROUTINE ssh_wzv 246 247 248 SUBROUTINE ssh_nxt( kt ) 238 ENDIF 239 ! 240 CALL wrk_dealloc( jpi, jpj, z2d ) 241 CALL wrk_dealloc( jpi, jpj, jpk, z3d, zhdiv ) 242 ! 243 IF( nn_timing == 1 ) CALL timing_stop('wzv') 244 245 246 END SUBROUTINE wzv 247 248 249 SUBROUTINE ssh_swp( kt ) 249 250 !!---------------------------------------------------------------------- 250 251 !! *** ROUTINE ssh_nxt *** … … 252 253 !! ** Purpose : achieve the sea surface height time stepping by 253 254 !! applying Asselin time filter and swapping the arrays 254 !! ssha already computed in ssh_ wzv255 !! ssha already computed in ssh_nxt 255 256 !! 256 257 !! ** Method : - apply Asselin time fiter to now ssh (excluding the forcing … … 266 267 !!---------------------------------------------------------------------- 267 268 INTEGER, INTENT(in) :: kt ! ocean time-step index 268 !! 269 INTEGER :: ji, jj ! dummy loop indices 270 REAL(wp) :: zec ! temporary scalar 271 !!---------------------------------------------------------------------- 272 ! 273 IF( nn_timing == 1 ) CALL timing_start('ssh_nxt') 269 !!---------------------------------------------------------------------- 270 ! 271 IF( nn_timing == 1 ) CALL timing_start('ssh_swp') 274 272 ! 275 273 IF( kt == nit000 ) THEN 276 274 IF(lwp) WRITE(numout,*) 277 IF(lwp) WRITE(numout,*) 'ssh_ nxt : next sea surface height (Asselin time filter + swap)'275 IF(lwp) WRITE(numout,*) 'ssh_swp : Asselin time filter and swap of sea surface height' 278 276 IF(lwp) WRITE(numout,*) '~~~~~~~ ' 279 277 ENDIF 280 278 281 ! !--------------------------! 282 IF( lk_vvl ) THEN ! Variable volume levels ! (ssh at t-, u-, v, f-points) 283 ! !--------------------------! 284 ! 285 IF( neuler == 0 .AND. kt == nit000 ) THEN !** Euler time-stepping at first time-step : no filter 286 sshn (:,:) = ssha (:,:) ! now <-- after (before already = now) 287 sshu_n(:,:) = sshu_a(:,:) 288 sshv_n(:,:) = sshv_a(:,:) 289 DO jj = 1, jpjm1 ! ssh now at f-point 290 DO ji = 1, jpim1 ! NO Vector Opt. 291 sshf_n(ji,jj) = 0.5 * umask(ji,jj,1) * umask(ji,jj+1,1) & 292 & / ( e1f(ji,jj ) * e2f(ji,jj ) ) & 293 & * ( e1u(ji,jj ) * e2u(ji,jj ) * sshu_n(ji,jj ) & 294 & + e1u(ji,jj+1) * e2u(ji,jj+1) * sshu_n(ji,jj+1) ) 295 END DO 296 END DO 297 CALL lbc_lnk( sshf_n, 'F', 1. ) ! Boundaries conditions 298 ! 299 ELSE !** Leap-Frog time-stepping: Asselin filter + swap 300 zec = atfp * rdt / rau0 301 DO jj = 1, jpj 302 DO ji = 1, jpi ! before <-- now filtered 303 sshb (ji,jj) = sshn (ji,jj) + atfp * ( sshb(ji,jj) - 2 * sshn(ji,jj) + ssha(ji,jj) ) & 304 & - zec * ( emp_b(ji,jj) - emp(ji,jj) ) * tmask(ji,jj,1) 305 sshn (ji,jj) = ssha (ji,jj) ! now <-- after 306 sshu_n(ji,jj) = sshu_a(ji,jj) 307 sshv_n(ji,jj) = sshv_a(ji,jj) 308 END DO 309 END DO 310 DO jj = 1, jpjm1 ! ssh now at f-point 311 DO ji = 1, jpim1 ! NO Vector Opt. 312 sshf_n(ji,jj) = 0.5 * umask(ji,jj,1) * umask(ji,jj+1,1) & 313 & / ( e1f(ji,jj ) * e2f(ji,jj ) ) & 314 & * ( e1u(ji,jj ) * e2u(ji,jj ) * sshu_n(ji,jj ) & 315 & + e1u(ji,jj+1) * e2u(ji,jj+1) * sshu_n(ji,jj+1) ) 316 END DO 317 END DO 318 CALL lbc_lnk( sshf_n, 'F', 1. ) ! Boundaries conditions 319 ! 320 DO jj = 1, jpjm1 ! ssh before at u- & v-points 321 DO ji = 1, jpim1 ! NO Vector Opt. 322 sshu_b(ji,jj) = 0.5 * umask(ji,jj,1) / ( e1u(ji ,jj) * e2u(ji ,jj) ) & 323 & * ( e1t(ji ,jj) * e2t(ji ,jj) * sshb(ji ,jj) & 324 & + e1t(ji+1,jj) * e2t(ji+1,jj) * sshb(ji+1,jj) ) 325 sshv_b(ji,jj) = 0.5 * vmask(ji,jj,1) / ( e1v(ji,jj ) * e2v(ji,jj ) ) & 326 & * ( e1t(ji,jj ) * e2t(ji,jj ) * sshb(ji,jj ) & 327 & + e1t(ji,jj+1) * e2t(ji,jj+1) * sshb(ji,jj+1) ) 328 END DO 329 END DO 330 CALL lbc_lnk( sshu_b, 'U', 1. ) 331 CALL lbc_lnk( sshv_b, 'V', 1. ) ! Boundaries conditions 332 ! 333 ENDIF 334 ! !--------------------------! 335 ELSE ! fixed levels ! (ssh at t-point only) 336 ! !--------------------------! 337 ! 338 IF( neuler == 0 .AND. kt == nit000 ) THEN !** Euler time-stepping at first time-step : no filter 339 sshn(:,:) = ssha(:,:) ! now <-- after (before already = now) 340 ! 341 ELSE ! Leap-Frog time-stepping: Asselin filter + swap 342 DO jj = 1, jpj 343 DO ji = 1, jpi ! before <-- now filtered 344 sshb(ji,jj) = sshn(ji,jj) + atfp * ( sshb(ji,jj) - 2 * sshn(ji,jj) + ssha(ji,jj) ) 345 sshn(ji,jj) = ssha(ji,jj) ! now <-- after 346 END DO 347 END DO 348 ENDIF 349 ! 279 IF( neuler == 0 .AND. kt == nit000 ) THEN !** Euler time-stepping at first time-step : no filter 280 sshn(:,:) = ssha(:,:) ! now <-- after (before already = now) 281 ELSE !** Leap-Frog time-stepping: Asselin filter + swap 282 sshb(:,:) = sshn(:,:) + atfp * ( sshb(:,:) - 2 * sshn(:,:) + ssha(:,:) ) ! before <-- now filtered 283 IF( lk_vvl ) sshb(:,:) = sshb(:,:) - atfp * rdt / rau0 * ( emp_b(:,:) - emp(:,:) ) * tmask(:,:,1) 284 sshn(:,:) = ssha(:,:) ! now <-- after 350 285 ENDIF 351 286 ! … … 357 292 IF(ln_ctl) CALL prt_ctl( tab2d_1=sshb, clinfo1=' sshb - : ', mask1=tmask, ovlap=1 ) 358 293 ! 359 IF( nn_timing == 1 ) CALL timing_stop('ssh_ nxt')360 ! 361 END SUBROUTINE ssh_ nxt294 IF( nn_timing == 1 ) CALL timing_stop('ssh_swp') 295 ! 296 END SUBROUTINE ssh_swp 362 297 363 298 !!====================================================================== -
branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90
r3680 r3865 91 91 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: sss_m !: mean (nn_fsbc time-step) surface sea salinity [psu] 92 92 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ssh_m !: mean (nn_fsbc time-step) sea surface height [m] 93 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: e3t_m !: mean (nn_fsbc time-step) sea surface height [m] 93 94 94 95 !! * Substitutions … … 105 106 !! *** FUNCTION sbc_oce_alloc *** 106 107 !!--------------------------------------------------------------------- 107 INTEGER :: ierr( 4)108 INTEGER :: ierr(5) 108 109 !!--------------------------------------------------------------------- 109 110 ierr(:) = 0 … … 126 127 & ssu_m (jpi,jpj) , sst_m(jpi,jpj) , & 127 128 & ssv_m (jpi,jpj) , sss_m (jpi,jpj), ssh_m(jpi,jpj) , STAT=ierr(4) ) 129 ! 130 #if defined key_vvl 131 ALLOCATE( e3t_m(jpi,jpj) , STAT=ierr(5) ) 132 #endif 128 133 ! 129 134 sbc_oce_alloc = MAXVAL( ierr ) -
branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90
r3625 r3865 412 412 ! Freezing/melting potential 413 413 ! Calculated over NEMO leapfrog timestep (hence 2*dt) 414 nfrzmlt(:,:)=rau0*rcp*fse3t_m(:,: ,1)*(Tocnfrz-sst_m(:,:))/(2.0*dt)414 nfrzmlt(:,:)=rau0*rcp*fse3t_m(:,:)*(Tocnfrz-sst_m(:,:))/(2.0*dt) 415 415 416 416 ztmp(:,:) = nfrzmlt(:,:) -
branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/SBC/sbcssm.F90
r3680 r3865 24 24 PRIVATE 25 25 26 PUBLIC sbc_ssm ! routine called by step.F90 27 PUBLIC sbc_ssm_init ! routine called by sbcmod.F90 28 26 PUBLIC sbc_ssm ! routine called by step.F90 27 29 28 LOGICAL, SAVE :: l_ssm_mean = .FALSE. ! keep track of whether means have been read 30 29 ! from restart file 31 30 32 31 !! * Substitutions 33 32 # include "domzgr_substitute.h90" … … 67 66 ELSE ; ssh_m(:,:) = sshn(:,:) 68 67 ENDIF 69 68 ! 69 IF( lk_vvl ) fse3t_m(:,:) = fse3t_n(:,:,1) 70 70 ! 71 71 ELSE … … 84 84 ELSE ; ssh_m(:,:) = zcoef * sshn(:,:) 85 85 ENDIF 86 IF( lk_vvl ) fse3t_m(:,:) = zcoef * fse3t_n(:,:,1) 86 87 ! ! ---------------------------------------- ! 87 88 ELSEIF( MOD( kt - 2 , nn_fsbc ) == 0 ) THEN ! Initialisation: New mean computation ! … … 92 93 sss_m(:,:) = 0.e0 93 94 ssh_m(:,:) = 0.e0 95 IF( lk_vvl ) fse3t_m(:,:) = 0.e0 94 96 ENDIF 95 97 ! ! ---------------------------------------- ! … … 104 106 ELSE ; ssh_m(:,:) = ssh_m(:,:) + sshn(:,:) 105 107 ENDIF 108 IF( lk_vvl ) fse3t_m(:,:) = fse3t_m(:,:) + fse3t_n(:,:,1) 106 109 107 110 ! ! ---------------------------------------- ! … … 114 117 ssv_m(:,:) = ssv_m(:,:) * zcoef ! 115 118 ssh_m(:,:) = ssh_m(:,:) * zcoef ! mean SSH [m] 119 IF( lk_vvl ) fse3t_m(:,:) = fse3t_m(:,:) * zcoef ! mean vertical scale factor [m] 116 120 ! 117 121 ENDIF … … 130 134 CALL iom_rstput( kt, nitrst, numrow, 'sss_m' , sss_m ) 131 135 CALL iom_rstput( kt, nitrst, numrow, 'ssh_m' , ssh_m ) 136 IF( lk_vvl ) THEN 137 CALL iom_rstput( kt, nitrst, numrow, 'fse3t_m' , fse3t_m(:,:) ) 138 END IF 132 139 ! 133 140 ENDIF … … 168 175 CALL iom_get( numror, jpdom_autoglo, 'sss_m' , sss_m ) ! " " salinity (T-point) 169 176 CALL iom_get( numror, jpdom_autoglo, 'ssh_m' , ssh_m ) ! " " height (T-point) 177 IF( lk_vvl ) CALL iom_get( numror, jpdom_autoglo, 'fse3t_m', fse3t_m(:,:) ) 170 178 ! 171 179 IF( zf_sbc /= REAL( nn_fsbc, wp ) ) THEN ! nn_fsbc has changed between 2 runs … … 178 186 sss_m(:,:) = zcoef * sss_m(:,:) 179 187 ssh_m(:,:) = zcoef * ssh_m(:,:) 188 IF( lk_vvl ) fse3t_m(:,:) = zcoef * fse3t_m(:,:) 180 189 ELSE 181 190 IF(lwp) WRITE(numout,*) '~~~~~~~ mean fields read in the ocean restart file' -
branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/SBC/sbcwave.F90
r3680 r3865 155 155 DO jj = 1, jpj-1 156 156 DO ji = 1, jpi-1 157 usd3d(ji,jj,jk) = usd2d(ji,jj)*exp(2.0*uwavenum(ji,jj)*(-MIN( gdept (ji,jj,jk) , gdept(ji+1,jj ,jk))))158 vsd3d(ji,jj,jk) = vsd2d(ji,jj)*exp(2.0*vwavenum(ji,jj)*(-MIN( gdept (ji,jj,jk) , gdept(ji ,jj+1,jk))))157 usd3d(ji,jj,jk) = usd2d(ji,jj)*exp(2.0*uwavenum(ji,jj)*(-MIN( gdept_0(ji,jj,jk) , gdept_0(ji+1,jj ,jk)))) 158 vsd3d(ji,jj,jk) = vsd2d(ji,jj)*exp(2.0*vwavenum(ji,jj)*(-MIN( gdept_0(ji,jj,jk) , gdept_0(ji ,jj+1,jk)))) 159 159 END DO 160 160 END DO 161 usd3d(jpi,:,jk) = usd2d(jpi,:)*exp( 2.0*uwavenum(jpi,:)*(-gdept (jpi,:,jk)) )162 vsd3d(:,jpj,jk) = vsd2d(:,jpj)*exp( 2.0*vwavenum(:,jpj)*(-gdept (:,jpj,jk)) )161 usd3d(jpi,:,jk) = usd2d(jpi,:)*exp( 2.0*uwavenum(jpi,:)*(-gdept_0(jpi,:,jk)) ) 162 vsd3d(:,jpj,jk) = vsd2d(:,jpj)*exp( 2.0*vwavenum(:,jpj)*(-gdept_0(:,jpj,jk)) ) 163 163 END DO 164 164 -
branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/TRA/eosbn2.F90
r3625 r3865 74 74 CONTAINS 75 75 76 SUBROUTINE eos_insitu( pts, prd )76 SUBROUTINE eos_insitu( pts, prd, pdep ) 77 77 !!---------------------------------------------------------------------- 78 78 !! *** ROUTINE eos_insitu *** … … 114 114 ! ! 2 : salinity [psu] 115 115 REAL(wp), DIMENSION(:,:,:) , INTENT( out) :: prd ! in situ density [-] 116 REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: pdep ! depth [m] 116 117 !! 117 118 INTEGER :: ji, jj, jk ! dummy loop indices … … 140 141 zt = pts (ji,jj,jk,jp_tem) 141 142 zs = pts (ji,jj,jk,jp_sal) 142 zh = fsdept(ji,jj,jk) ! depth143 zh = pdep(ji,jj,jk) ! depth 143 144 zsr= zws (ji,jj,jk) ! square root salinity 144 145 ! … … 198 199 199 200 200 SUBROUTINE eos_insitu_pot( pts, prd, prhop )201 SUBROUTINE eos_insitu_pot( pts, prd, prhop, pdep ) 201 202 !!---------------------------------------------------------------------- 202 203 !! *** ROUTINE eos_insitu_pot *** … … 249 250 REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT( out) :: prd ! in situ density [-] 250 251 REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT( out) :: prhop ! potential density (surface referenced) 252 REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pdep ! depth [m] 251 253 ! 252 254 INTEGER :: ji, jj, jk ! dummy loop indices … … 271 273 zt = pts (ji,jj,jk,jp_tem) 272 274 zs = pts (ji,jj,jk,jp_sal) 273 zh = fsdept(ji,jj,jk) ! depth275 zh = pdep(ji,jj,jk) ! depth 274 276 zsr= zws (ji,jj,jk) ! square root salinity 275 277 ! -
branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/TRA/traadv.F90
r3718 r3865 14 14 USE oce ! ocean dynamics and active tracers 15 15 USE dom_oce ! ocean space and time domain 16 USE domvvl ! variable vertical scale factors 16 17 USE traadv_cen2 ! 2nd order centered scheme (tra_adv_cen2 routine) 17 18 USE traadv_tvd ! TVD scheme (tra_adv_tvd routine) … … 92 93 zwn(:,:,jk) = e1t(:,:) * e2t(:,:) * wn(:,:,jk) 93 94 END DO 95 ! 96 IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN 97 zun(:,:,:) = zun(:,:,:) + un_td(:,:,:) 98 zvn(:,:,:) = zvn(:,:,:) + vn_td(:,:,:) 99 ENDIF 100 ! 94 101 zun(:,:,jpk) = 0._wp ! no transport trough the bottom 95 102 zvn(:,:,jpk) = 0._wp ! no transport trough the bottom -
branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/TRA/trabbl.F90
r3764 r3865 66 66 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ahu_bbl_0, ahv_bbl_0 ! diffusive bbl flux coefficients at u and v-points 67 67 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC :: e3u_bbl_0, e3v_bbl_0 ! thichness of the bbl (e3) at u and v-points (PUBLIC for TAM) 68 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC :: r1_e1e2t ! inverse of the cell surface at t-point [1/m2] (PUBLIC for TAM)69 68 70 69 !! * Substitutions … … 85 84 & vtr_bbl (jpi,jpj) , ahv_bbl (jpi,jpj) , mbkv_d (jpi,jpj) , mgrhv(jpi,jpj) , & 86 85 & ahu_bbl_0(jpi,jpj) , ahv_bbl_0(jpi,jpj) , & 87 & e3u_bbl_0(jpi,jpj) , e3v_bbl_0(jpi,jpj) , r1_e1e2t(jpi,jpj) , STAT= tra_bbl_alloc)86 & e3u_bbl_0(jpi,jpj) , e3v_bbl_0(jpi,jpj) , STAT= tra_bbl_alloc ) 88 87 ! 89 88 IF( lk_mpp ) CALL mpp_sum ( tra_bbl_alloc ) … … 217 216 # endif 218 217 ik = mbkt(ji,jj) ! bottom T-level index 219 zbtr = r1_e1 e2t(ji,jj) / fse3t(ji,jj,ik)218 zbtr = r1_e12t(ji,jj) / fse3t(ji,jj,ik) 220 219 pta(ji,jj,ik,jn) = pta(ji,jj,ik,jn) & 221 220 & + ( ahu_bbl(ji ,jj ) * ( zptb(ji+1,jj ) - zptb(ji ,jj ) ) & … … 279 278 ! 280 279 ! ! up -slope T-point (shelf bottom point) 281 zbtr = r1_e1 e2t(iis,jj) / fse3t(iis,jj,ikus)280 zbtr = r1_e12t(iis,jj) / fse3t(iis,jj,ikus) 282 281 ztra = zu_bbl * ( ptb(iid,jj,ikus,jn) - ptb(iis,jj,ikus,jn) ) * zbtr 283 282 pta(iis,jj,ikus,jn) = pta(iis,jj,ikus,jn) + ztra 284 283 ! 285 284 DO jk = ikus, ikud-1 ! down-slope upper to down T-point (deep column) 286 zbtr = r1_e1 e2t(iid,jj) / fse3t(iid,jj,jk)285 zbtr = r1_e12t(iid,jj) / fse3t(iid,jj,jk) 287 286 ztra = zu_bbl * ( ptb(iid,jj,jk+1,jn) - ptb(iid,jj,jk,jn) ) * zbtr 288 287 pta(iid,jj,jk,jn) = pta(iid,jj,jk,jn) + ztra 289 288 END DO 290 289 ! 291 zbtr = r1_e1 e2t(iid,jj) / fse3t(iid,jj,ikud)290 zbtr = r1_e12t(iid,jj) / fse3t(iid,jj,ikud) 292 291 ztra = zu_bbl * ( ptb(iis,jj,ikus,jn) - ptb(iid,jj,ikud,jn) ) * zbtr 293 292 pta(iid,jj,ikud,jn) = pta(iid,jj,ikud,jn) + ztra … … 301 300 ! 302 301 ! up -slope T-point (shelf bottom point) 303 zbtr = r1_e1 e2t(ji,ijs) / fse3t(ji,ijs,ikvs)302 zbtr = r1_e12t(ji,ijs) / fse3t(ji,ijs,ikvs) 304 303 ztra = zv_bbl * ( ptb(ji,ijd,ikvs,jn) - ptb(ji,ijs,ikvs,jn) ) * zbtr 305 304 pta(ji,ijs,ikvs,jn) = pta(ji,ijs,ikvs,jn) + ztra 306 305 ! 307 306 DO jk = ikvs, ikvd-1 ! down-slope upper to down T-point (deep column) 308 zbtr = r1_e1 e2t(ji,ijd) / fse3t(ji,ijd,jk)307 zbtr = r1_e12t(ji,ijd) / fse3t(ji,ijd,jk) 309 308 ztra = zv_bbl * ( ptb(ji,ijd,jk+1,jn) - ptb(ji,ijd,jk,jn) ) * zbtr 310 309 pta(ji,ijd,jk,jn) = pta(ji,ijd,jk,jn) + ztra 311 310 END DO 312 311 ! ! down-slope T-point (deep bottom point) 313 zbtr = r1_e1 e2t(ji,ijd) / fse3t(ji,ijd,ikvd)312 zbtr = r1_e12t(ji,ijd) / fse3t(ji,ijd,ikvd) 314 313 ztra = zv_bbl * ( ptb(ji,ijs,ikvs,jn) - ptb(ji,ijd,ikvd,jn) ) * zbtr 315 314 pta(ji,ijd,ikvd,jn) = pta(ji,ijd,ikvd,jn) + ztra … … 423 422 ztb (ji,jj) = tsb(ji,jj,ik,jp_tem) * tmask(ji,jj,1) ! bottom before T and S 424 423 zsb (ji,jj) = tsb(ji,jj,ik,jp_sal) * tmask(ji,jj,1) 425 zdep(ji,jj) = fsdept_0(ji,jj,ik)! bottom T-level reference depth424 zdep(ji,jj) = gdept_0(ji,jj,ik) ! bottom T-level reference depth 426 425 ! 427 426 zub(ji,jj) = un(ji,jj,mbku(ji,jj)) ! bottom velocity … … 594 593 IF( nn_eos /= 0 ) CALL ctl_stop ( ' bbl parameterisation requires eos = 0. We stop.' ) 595 594 596 597 ! !* inverse of surface of T-cells598 r1_e1e2t(:,:) = 1._wp / ( e1t(:,:) * e2t(:,:) )599 600 595 ! !* vertical index of "deep" bottom u- and v-points 601 596 DO jj = 1, jpjm1 ! (the "shelf" bottom k-indices are mbku and mbkv) … … 605 600 END DO 606 601 END DO 607 ! convert einto REAL to use lbc_lnk ; impose a min value of 1 as a zero can be set in lbclnk602 ! convert into REAL to use lbc_lnk ; impose a min value of 1 as a zero can be set in lbclnk 608 603 zmbk(:,:) = REAL( mbku_d(:,:), wp ) ; CALL lbc_lnk(zmbk,'U',1.) ; mbku_d(:,:) = MAX( INT( zmbk(:,:) ), 1 ) 609 604 zmbk(:,:) = REAL( mbkv_d(:,:), wp ) ; CALL lbc_lnk(zmbk,'V',1.) ; mbkv_d(:,:) = MAX( INT( zmbk(:,:) ), 1 ) 610 605 611 606 !* sign of grad(H) at u- and v-points 612 607 mgrhu(jpi,:) = 0. ; mgrhu(:,jpj) = 0. ; mgrhv(jpi,:) = 0. ; mgrhv(:,jpj) = 0. 613 608 DO jj = 1, jpjm1 614 609 DO ji = 1, jpim1 615 mgrhu(ji,jj) = INT( SIGN( 1.e0, fsdept_0(ji+1,jj,mbkt(ji+1,jj)) - fsdept_0(ji,jj,mbkt(ji,jj)) ) )616 mgrhv(ji,jj) = INT( SIGN( 1.e0, fsdept_0(ji,jj+1,mbkt(ji,jj+1)) - fsdept_0(ji,jj,mbkt(ji,jj)) ) )610 mgrhu(ji,jj) = INT( SIGN( 1.e0, gdept_0(ji+1,jj,mbkt(ji+1,jj)) - gdept_0(ji,jj,mbkt(ji,jj)) ) ) 611 mgrhv(ji,jj) = INT( SIGN( 1.e0, gdept_0(ji,jj+1,mbkt(ji,jj+1)) - gdept_0(ji,jj,mbkt(ji,jj)) ) ) 617 612 END DO 618 613 END DO 619 614 620 615 DO jj = 1, jpjm1 !* bbl thickness at u- (v-) point 621 DO ji = 1, jpim1 622 e3u_bbl_0(ji,jj) = MIN( fse3u_0(ji,jj,mbkt(ji+1,jj )), fse3u_0(ji,jj,mbkt(ji,jj)) )623 e3v_bbl_0(ji,jj) = MIN( fse3v_0(ji,jj,mbkt(ji ,jj+1)), fse3v_0(ji,jj,mbkt(ji,jj)) )616 DO ji = 1, jpim1 ! minimum of top & bottom e3u_0 (e3v_0) 617 e3u_bbl_0(ji,jj) = MIN( e3u_0(ji,jj,mbkt(ji+1,jj )), e3u_0(ji,jj,mbkt(ji,jj)) ) 618 e3v_bbl_0(ji,jj) = MIN( e3v_0(ji,jj,mbkt(ji ,jj+1)), e3v_0(ji,jj,mbkt(ji,jj)) ) 624 619 END DO 625 620 END DO -
branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_bilap.F90
r3294 r3865 110 110 DO jj = 1, jpjm1 111 111 DO ji = 1, fs_jpim1 ! vector opt. 112 zeeu(ji,jj) = e2u(ji,jj) * fse3u(ji,jj,jk) / e1u(ji,jj) * umask(ji,jj,jk)113 zeev(ji,jj) = e1v(ji,jj) * fse3v(ji,jj,jk) / e2v(ji,jj) * vmask(ji,jj,jk)112 zeeu(ji,jj) = re2u_e1u(ji,jj) * fse3u_n(ji,jj,jk) * umask(ji,jj,jk) 113 zeev(ji,jj) = re1v_e2v(ji,jj) * fse3v_n(ji,jj,jk) * vmask(ji,jj,jk) 114 114 END DO 115 115 END DO … … 133 133 DO jj = 2, jpjm1 ! Second derivative (divergence) time the eddy diffusivity coefficient 134 134 DO ji = fs_2, fs_jpim1 ! vector opt. 135 zbtr = 1.0 / ( e1 t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )135 zbtr = 1.0 / ( e12t(ji,jj) * fse3t_n(ji,jj,jk) ) 136 136 zlt(ji,jj) = fsahtt(ji,jj,jk) * zbtr * ( ztu(ji,jj,jk) - ztu(ji-1,jj,jk) & 137 137 & + ztv(ji,jj,jk) - ztv(ji,jj-1,jk) ) … … 151 151 DO ji = fs_2, fs_jpim1 ! vector opt. 152 152 ! horizontal diffusive trends 153 zbtr = 1.0 / ( e1 t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )153 zbtr = 1.0 / ( e12t(ji,jj) * fse3t_n(ji,jj,jk) ) 154 154 ztra = zbtr * ( ztu(ji,jj,jk) - ztu(ji-1,jj,jk) + ztv(ji,jj,jk) - ztv(ji,jj-1,jk) ) 155 155 ! add it to the general tracer trends -
branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_bilapg.F90
r3805 r3865 210 210 DO jj = 1, jpjm1 211 211 DO ji = 1, jpim1 212 zabe1 = e2u(ji,jj) * fse3u(ji,jj,jk) / e1u(ji,jj)213 zabe2 = e1v(ji,jj) * fse3v(ji,jj,jk) / e2v(ji,jj)212 zabe1 = re2u_e1u(ji,jj) * fse3u_n(ji,jj,jk) 213 zabe2 = re1v_e2v(ji,jj) * fse3v_n(ji,jj,jk) 214 214 215 215 zmku = 1./MAX( tmask(ji+1,jj,jk )+tmask(ji,jj,jk+1) & … … 279 279 DO jk = 2, jpkm1 280 280 DO ji = 2, jpim1 281 zcof0 = e1 t(ji,jj) * e2t(ji,jj) / fse3w(ji,jj,jk) &281 zcof0 = e12t(ji,jj) / fse3w_n(ji,jj,jk) & 282 282 & * ( wslpi(ji,jj,jk) * wslpi(ji,jj,jk) & 283 283 & + wslpj(ji,jj,jk) * wslpj(ji,jj,jk) ) … … 310 310 DO ji = 2, jpim1 311 311 ! eddy coef. divided by the volume element 312 zbtr = 1.0 / ( e1 t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )312 zbtr = 1.0 / ( e12t(ji,jj) * fse3t_n(ji,jj,jk) ) 313 313 ! vertical divergence 314 314 ztav = fsahtt(ji,jj,jk) * ( zftw(ji,jk) - zftw(ji,jk+1) ) … … 322 322 DO ji = 2, jpim1 323 323 ! inverse of the volume element 324 zbtr = 1.0 / ( e1 t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )324 zbtr = 1.0 / ( e12t(ji,jj) * fse3t_n(ji,jj,jk) ) 325 325 ! vertical divergence 326 326 ztav = zftw(ji,jk) - zftw(ji,jk+1) -
branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_iso.F90
r3805 r3865 176 176 DO jj = 1 , jpjm1 177 177 DO ji = 1, fs_jpim1 ! vector opt. 178 zabe1 = ( fsahtu(ji,jj,jk) + pahtb0 ) * e2u(ji,jj) * fse3u(ji,jj,jk) / e1u(ji,jj)179 zabe2 = ( fsahtv(ji,jj,jk) + pahtb0 ) * e1v(ji,jj) * fse3v(ji,jj,jk) / e2v(ji,jj)178 zabe1 = ( fsahtu(ji,jj,jk) + pahtb0 ) * re2u_e1u(ji,jj) * fse3u_n(ji,jj,jk) 179 zabe2 = ( fsahtv(ji,jj,jk) + pahtb0 ) * re1v_e2v(ji,jj) * fse3v_n(ji,jj,jk) 180 180 ! 181 181 zmsku = 1. / MAX( tmask(ji+1,jj,jk ) + tmask(ji,jj,jk+1) & … … 201 201 DO jj = 2 , jpjm1 202 202 DO ji = fs_2, fs_jpim1 ! vector opt. 203 zbtr = 1.0 / ( e1 t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )203 zbtr = 1.0 / ( e12t(ji,jj) * fse3t_n(ji,jj,jk) ) 204 204 ztra = zbtr * ( zftu(ji,jj,jk) - zftu(ji-1,jj,jk) + zftv(ji,jj,jk) - zftv(ji,jj-1,jk) ) 205 205 pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra … … 288 288 DO jj = 2, jpjm1 289 289 DO ji = fs_2, fs_jpim1 ! vector opt. 290 zbtr = 1.0 / ( e1 t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )290 zbtr = 1.0 / ( e12t(ji,jj) * fse3t_n(ji,jj,jk) ) 291 291 ztra = ( ztfw(ji,jj,jk) - ztfw(ji,jj,jk+1) ) * zbtr 292 292 pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra -
branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_lap.F90
r3294 r3865 31 31 32 32 PUBLIC tra_ldf_lap ! routine called by step.F90 33 34 REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:,:) :: e1ur, e2vr ! scale factor coefficients35 33 36 34 !! * Substitutions … … 85 83 IF(lwp) WRITE(numout,*) 'tra_ldf_lap : iso-level laplacian diffusion on ', cdtype 86 84 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ ' 87 !88 IF( .NOT. ALLOCATED( e1ur ) ) THEN89 ! This routine may be called for both active and passive tracers.90 ! Allocate and set saved arrays on first call only.91 ALLOCATE( e1ur(jpi,jpj), e2vr(jpi,jpj), STAT=ierr )92 IF( lk_mpp ) CALL mpp_sum( ierr )93 IF( ierr /= 0 ) CALL ctl_stop( 'STOP', 'tra_ldf_lap : unable to allocate arrays' )94 !95 e1ur(:,:) = e2u(:,:) / e1u(:,:)96 e2vr(:,:) = e1v(:,:) / e2v(:,:)97 ENDIF98 85 ENDIF 99 86 … … 107 94 DO jj = 1, jpjm1 108 95 DO ji = 1, fs_jpim1 ! vector opt. 109 zabe1 = fsahtu(ji,jj,jk) * umask(ji,jj,jk) * e1ur(ji,jj) * fse3u(ji,jj,jk)110 zabe2 = fsahtv(ji,jj,jk) * vmask(ji,jj,jk) * e2vr(ji,jj) * fse3v(ji,jj,jk)96 zabe1 = fsahtu(ji,jj,jk) * umask(ji,jj,jk) * re2u_e1u(ji,jj) * fse3u_n(ji,jj,jk) 97 zabe2 = fsahtv(ji,jj,jk) * vmask(ji,jj,jk) * re1v_e2v(ji,jj) * fse3v_n(ji,jj,jk) 111 98 ztu(ji,jj,jk) = zabe1 * ( ptb(ji+1,jj ,jk,jn) - ptb(ji,jj,jk,jn) ) 112 99 ztv(ji,jj,jk) = zabe2 * ( ptb(ji ,jj+1,jk,jn) - ptb(ji,jj,jk,jn) ) … … 120 107 ikv = mbkv(ji,jj) 121 108 IF( iku == jk ) THEN 122 zabe1 = fsahtu(ji,jj,iku) * umask(ji,jj,iku) * e1ur(ji,jj) * fse3u(ji,jj,iku)109 zabe1 = fsahtu(ji,jj,iku) * umask(ji,jj,iku) * re2u_e1u(ji,jj) * fse3u_n(ji,jj,iku) 123 110 ztu(ji,jj,jk) = zabe1 * pgu(ji,jj,jn) 124 111 ENDIF 125 112 IF( ikv == jk ) THEN 126 zabe2 = fsahtv(ji,jj,ikv) * vmask(ji,jj,ikv) * e2vr(ji,jj) * fse3v(ji,jj,ikv)113 zabe2 = fsahtv(ji,jj,ikv) * vmask(ji,jj,ikv) * re1v_e2v(ji,jj) * fse3v_n(ji,jj,ikv) 127 114 ztv(ji,jj,jk) = zabe2 * pgv(ji,jj,jn) 128 115 ENDIF … … 136 123 DO jj = 2, jpjm1 137 124 DO ji = fs_2, fs_jpim1 ! vector opt. 138 zbtr = 1._wp / ( e1 t(ji,jj) *e2t(ji,jj) * fse3t(ji,jj,jk) )125 zbtr = 1._wp / ( e12t(ji,jj) * fse3t_n(ji,jj,jk) ) 139 126 ! horizontal diffusive trends added to the general tracer trends 140 127 pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + zbtr * ( ztu(ji,jj,jk) - ztu(ji-1,jj,jk) & -
branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/TRA/traqsr.F90
r3862 r3865 448 448 !CDIR NOVERRCHK 449 449 DO ji = 1, jpi 450 zc0 = ze0(ji,jj,jk-1) * EXP( - fse3t_0(ji,jj,jk-1) * xsi0r )451 zc1 = ze1(ji,jj,jk-1) * EXP( - fse3t_0(ji,jj,jk-1) * zekb(ji,jj) )452 zc2 = ze2(ji,jj,jk-1) * EXP( - fse3t_0(ji,jj,jk-1) * zekg(ji,jj) )453 zc3 = ze3(ji,jj,jk-1) * EXP( - fse3t_0(ji,jj,jk-1) * zekr(ji,jj) )450 zc0 = ze0(ji,jj,jk-1) * EXP( - e3t_0(ji,jj,jk-1) * xsi0r ) 451 zc1 = ze1(ji,jj,jk-1) * EXP( - e3t_0(ji,jj,jk-1) * zekb(ji,jj) ) 452 zc2 = ze2(ji,jj,jk-1) * EXP( - e3t_0(ji,jj,jk-1) * zekg(ji,jj) ) 453 zc3 = ze3(ji,jj,jk-1) * EXP( - e3t_0(ji,jj,jk-1) * zekr(ji,jj) ) 454 454 ze0(ji,jj,jk) = zc0 455 455 ze1(ji,jj,jk) = zc1 -
branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/oce.F90
r3625 r3865 36 36 !! ------------ ! fields ! fields ! trends ! 37 37 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:), TARGET :: sshb , sshn , ssha !: sea surface height at t-point [m] 38 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: sshu_b , sshu_n , sshu_a !: sea surface height at u-point [m]39 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: sshv_b , sshv_n , sshv_a !: sea surface height at u-point [m]40 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: sshf_n !: sea surface height at f-point [m]41 38 ! 42 39 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: spgu, spgv !: horizontal surface pressure gradient … … 77 74 & rn2b (jpi,jpj,jpk) , rn2 (jpi,jpj,jpk) , STAT=ierr(1) ) 78 75 ! 79 ALLOCATE( rhd (jpi,jpj,jpk) , & 80 & rhop(jpi,jpj,jpk) , & 81 & sshb (jpi,jpj) , sshn (jpi,jpj) , ssha (jpi,jpj) , & 82 & sshu_b(jpi,jpj) , sshu_n(jpi,jpj) , sshu_a(jpi,jpj) , & 83 & sshv_b(jpi,jpj) , sshv_n(jpi,jpj) , sshv_a(jpi,jpj) , & 84 & sshf_n(jpi,jpj) , & 85 & spgu (jpi,jpj) , spgv(jpi,jpj) , & 86 & gtsu(jpi,jpj,jpts), gtsv(jpi,jpj,jpts), & 87 & gru(jpi,jpj) , grv(jpi,jpj) , STAT=ierr(2) ) 76 ALLOCATE(rhd (jpi,jpj,jpk) , & 77 & rhop(jpi,jpj,jpk) , & 78 & sshb (jpi,jpj) , sshn (jpi,jpj) , ssha (jpi,jpj) , & 79 & spgu (jpi,jpj) , spgv(jpi,jpj) , & 80 & gtsu(jpi,jpj,jpts), gtsv(jpi,jpj,jpts), & 81 & gru(jpi,jpj) , grv(jpi,jpj) , STAT=ierr(2) ) 88 82 ! 89 83 ALLOCATE( snwice_mass(jpi,jpj) , snwice_mass_b(jpi,jpj), & -
branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/step.F90
r3769 r3865 97 97 98 98 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 99 ! Ocean dynamics : ssh, wn, hdiv, rot ! 100 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 101 CALL ssh_wzv( kstp ) ! after ssh & vertical velocity 99 ! Ocean dynamics : hdiv, rot, ssh, e3, wn 100 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 101 CALL ssh_nxt ( kstp ) ! after ssh 102 IF( lk_vvl ) CALL dom_vvl_sf_nxt( kstp ) ! after vertical scale factors 103 CALL wzv ( kstp ) ! now cross-level velocity 102 104 103 105 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> … … 139 141 ! 140 142 IF( lk_ldfslp ) THEN ! slope of lateral mixing 141 CALL eos( tsb, rhd )! before in situ density143 CALL eos( tsb, rhd, gdept_0(:,:,:) ) ! before in situ density 142 144 IF( ln_zps ) CALL zps_hde( kstp, jpts, tsb, gtsu, gtsv, & ! Partial steps: before horizontal gradient 143 145 & rhd, gru , grv ) ! of t, s, rd at the last ocean level … … 202 204 IF( ln_zdfnpc ) CALL tra_npc( kstp ) ! update after fields by non-penetrative convection 203 205 CALL tra_nxt( kstp ) ! tracer fields at next time step 204 CALL eos ( tsa, rhd, rhop )! Time-filtered in situ density for hpg computation206 CALL eos ( tsa, rhd, rhop, fsdept_n(:,:,:) ) ! Time-filtered in situ density for hpg computation 205 207 IF( ln_zps ) CALL zps_hde( kstp, jpts, tsa, gtsu, gtsv, & ! zps: time filtered hor. derivative 206 208 & rhd, gru , grv ) ! of t, s, rd at the last ocean level 207 209 208 210 ELSE ! centered hpg (eos then time stepping) 209 CALL eos ( tsn, rhd, rhop )! now in situ density for hpg computation211 CALL eos ( tsn, rhd, rhop, fsdept_n(:,:,:) ) ! now in situ density for hpg computation 210 212 IF( ln_zps ) CALL zps_hde( kstp, jpts, tsn, gtsu, gtsv, & ! zps: now hor. derivative 211 213 & rhd, gru , grv ) ! of t, s, rd at the last ocean level … … 238 240 CALL dyn_nxt( kstp ) ! lateral velocity at next time step 239 241 240 CALL ssh_nxt( kstp ) ! sea surface height at next time step 242 CALL ssh_swp( kstp ) ! swap of sea surface height 243 IF( lk_vvl ) CALL dom_vvl_sf_swp( kstp ) ! swap of vertical scale factors 241 244 242 245 IF( ln_diahsb ) CALL dia_hsb( kstp ) ! - ML - global conservation diagnostics -
branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/step_oce.F90
r3769 r3865 64 64 USE bdydyn3d ! bdy cond. for baroclinic vel. (bdy_dyn3d routine) 65 65 66 USE sshwzv ! vertical velocity and ssh (ssh_wzv routine) 66 USE sshwzv ! vertical velocity and ssh (ssh_nxt routine) 67 ! (ssh_swp routine) 68 ! (wzv routine) 69 USE domvvl ! variable vertical scale factors (dom_vvl_sf_nxt routine) 70 ! (dom_vvl_sf_swp routine) 67 71 68 72 USE ldfslp ! iso-neutral slopes (ldf_slp routine) -
branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/TOP_SRC/C14b/trcwri_c14b.F90
r3680 r3865 36 36 DO jn = jp_c14b0, jp_c14b1 37 37 cltra = TRIM( ctrcnm(jn) ) ! short title for tracer 38 CALL iom_put( cltra, trn(:,:,:,jn) ) 38 IF( lk_vvl ) THEN 39 CALL iom_put( cltra, trn(:,:,:,jn) * fse3t_n(:,:,:) ) 40 ELSE 41 CALL iom_put( cltra, trn(:,:,:,jn) ) 42 ENDIF 39 43 END DO 40 44 ! -
branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/TOP_SRC/CFC/trcwri_cfc.F90
r3680 r3865 36 36 DO jn = jp_cfc0, jp_cfc1 37 37 cltra = TRIM( ctrcnm(jn) ) ! short title for tracer 38 CALL iom_put( cltra, trn(:,:,:,jn) ) 38 IF( lk_vvl ) THEN 39 CALL iom_put( cltra, trn(:,:,:,jn) * fse3t_n(:,:,:) ) 40 ELSE 41 CALL iom_put( cltra, trn(:,:,:,jn) ) 42 ENDIF 39 43 END DO 40 44 ! -
branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/TOP_SRC/MY_TRC/trcwri_my_trc.F90
r3680 r3865 36 36 DO jn = jp_myt0, jp_myt1 37 37 cltra = TRIM( ctrcnm(jn) ) ! short title for tracer 38 CALL iom_put( cltra, trn(:,:,:,jn) ) 38 IF( lk_vvl ) THEN 39 CALL iom_put( cltra, trn(:,:,:,jn) * fse3t_n(:,:,:) ) 40 ELSE 41 CALL iom_put( cltra, trn(:,:,:,jn) ) 42 ENDIF 39 43 END DO 40 44 ! -
branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/TOP_SRC/PISCES/trcwri_pisces.F90
r3680 r3865 39 39 DO jn = jp_pcs0, jp_pcs1 40 40 cltra = TRIM( ctrcnm(jn) ) ! short title for tracer 41 IF( lk_vvl ) THEN 42 CALL iom_put( cltra, trn(:,:,:,jn) * fse3t_n(:,:,:) ) 43 ELSE 44 CALL iom_put( cltra, trn(:,:,:,jn) ) 45 ENDIF 41 46 CALL iom_put( cltra, trn(:,:,:,jn) * zrfact ) 42 47 END DO … … 47 52 IF( jn == jppo4 ) zrfact = po4r * 1.0e+6 48 53 cltra = TRIM( ctrcnm(jn) ) ! short title for tracer 49 CALL iom_put( cltra, trn(:,:,:,jn) * zrfact ) 54 IF( lk_vvl ) THEN 55 CALL iom_put( cltra, trn(:,:,:,jn) * fse3t_n(:,:,:) * zrfact ) 56 ELSE 57 CALL iom_put( cltra, trn(:,:,:,jn) * zrfact ) 58 ENDIF 50 59 END DO 51 60 #endif -
branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/TOP_SRC/oce_trc.F90
r3680 r3865 77 77 USE oce , ONLY : sshb => sshb !: sea surface height at t-point [m] 78 78 USE oce , ONLY : ssha => ssha !: sea surface height at t-point [m] 79 USE oce , ONLY : sshu_n => sshu_n !: sea surface height at u-point [m]80 USE oce , ONLY : sshu_b => sshu_b !: sea surface height at u-point [m]81 USE oce , ONLY : sshu_a => sshu_a !: sea surface height at u-point [m]82 USE oce , ONLY : sshv_n => sshv_n !: sea surface height at v-point [m]83 USE oce , ONLY : sshv_b => sshv_b !: sea surface height at v-point [m]84 USE oce , ONLY : sshv_a => sshv_a !: sea surface height at v-point [m]85 USE oce , ONLY : sshf_n => sshf_n !: sea surface height at v-point [m]86 79 USE oce , ONLY : l_traldf_rot => l_traldf_rot !: rotated laplacian operator for lateral diffusion 87 80 #if defined key_offline
Note: See TracChangeset
for help on using the changeset viewer.