Changeset 8093 for branches/2017
- Timestamp:
- 2017-05-30T10:13:14+02:00 (7 years ago)
- Location:
- branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM
- Files:
-
- 1 added
- 22 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/SHARED/namelist_ref
r7990 r8093 601 601 602 602 !!====================================================================== 603 !! *** Bottom boundary condition *** 604 !!====================================================================== 605 !! nambfr bottom friction 606 !! nambbc bottom temperature boundary condition 603 !! *** top/Bottom boundary condition *** !! 604 !!====================================================================== 605 !! nambfr bottom friction (default: NONE) 606 !! namtfr top friction (default: NONE) 607 !! nambbc bottom temperature boundary condition (default: NO) 607 608 !! nambbl bottom boundary layer scheme (default: NO) 608 609 !!====================================================================== … … 611 612 &nambfr ! bottom friction (default: linear) 612 613 !----------------------------------------------------------------------- 613 nn_bfr = 1 ! type of bottom friction : = 0 : free slip, = 1 : linear friction 614 ! = 2 : nonlinear friction 614 nn_bfr = 1 ! type of top/bottom drag: free slip (=0), linear drag (=1), 615 ! ! nonlinear drag (=2), nonlinear with logarithmic formulation (=3) 616 ln_bfrimp = .true. ! implicit bottom friction (requires ln_zdfexp = .false. if true) 617 ln_loglayer = .false. ! logarithmic formulation (non linear case only) 618 ! 615 619 rn_bfri1 = 4.e-4 ! bottom drag coefficient (linear case) 616 620 rn_bfri2 = 1.e-3 ! bottom drag coefficient (non linear case). Minimum coeft if ln_loglayer=T … … 619 623 rn_bfrz0 = 3.e-3 ! bottom roughness [m] if ln_loglayer=T 620 624 ln_bfr2d = .false. ! horizontal variation of the bottom friction coef (read a 2D mask file ) 621 rn_bfrien = 50. ! local multiplying factor of bfr (ln_bfr2d=T) 625 rn_bfrien = 50. ! local boost factor 626 ! 622 627 rn_tfri1 = 4.e-4 ! top drag coefficient (linear case) 623 628 rn_tfri2 = 2.5e-3 ! top drag coefficient (non linear case). Minimum coeft if ln_loglayer=T … … 626 631 rn_tfrz0 = 3.e-3 ! top roughness [m] if ln_loglayer=T 627 632 ln_tfr2d = .false. ! horizontal variation of the top friction coef (read a 2D mask file ) 628 rn_tfrien = 50. ! local multiplying factor of tfr (ln_tfr2d=T) 629 630 ln_bfrimp = .true. ! implicit bottom friction (requires ln_zdfexp = .false. if true) 631 ln_loglayer = .false. ! logarithmic formulation (non linear case) 632 / 633 rn_tfrien = 50. ! local boost factor 634 / 635 636 !----------------------------------------------------------------------- 637 &namdrg ! top/bottom drag coeeficient (default: NO selection) 638 !----------------------------------------------------------------------- 639 ln_NONE = .false. ! free-slip : Cd = 0 (F => fill namdrg_bot 640 ln_lin = .false. ! linear drag: Cd = Cd0 Uc0 & namdrg_top) 641 ln_non_lin = .false. ! non-linear drag: Cd = Cd0 |U| 642 ln_loglayer= .false. ! logarithmic drag: Cd = vkarmn/log(z/z0) |U| 643 ! 644 ln_drgimp = .true. ! implicit top/bottom friction flag 645 / 646 647 !----------------------------------------------------------------------- 648 &namdrg_top ! TOP friction (ln_isfcav=T) 649 !----------------------------------------------------------------------- 650 rn_Cd0 = 1.e-3 ! drag coefficient [-] 651 rn_Uc0 = 0.4 ! ref. velocity [m/s] (linear drag=Cd0*Uc0) 652 rn_Cdmax = 0.1 ! drag value maximum [-] (logarithmic drag) 653 rn_ke0 = 2.5e-3 ! background kinetic energy [m2/s2] (non-linear cases) 654 rn_z0 = 3.0e-3 ! roughness [m] (ln_loglayer=T) 655 ln_boost = .false. ! =T regional boost of Cd0 ; =F constant 656 rn_boost= 50. ! local boost factor [-] 657 / 658 659 !----------------------------------------------------------------------- 660 &namdrg_bot ! BOTTOM friction 661 !----------------------------------------------------------------------- 662 rn_Cd0 = 1.e-3 ! drag coefficient [-] 663 rn_Uc0 = 0.4 ! ref. velocity [m/s] (linear drag=Cd0*Uc0) 664 rn_Cdmax = 0.1 ! drag value maximum [-] (logarithmic drag) 665 rn_ke0 = 2.5e-3 ! background kinetic energy [m2/s2] (non-linear cases) 666 rn_z0 = 3.e-3 ! roughness [m] (ln_loglayer=T) 667 ln_boost = .false. ! =T regional boost of Cd0 ; =F constant 668 rn_boost= 50. ! local boost factor [-] 669 / 670 633 671 !----------------------------------------------------------------------- 634 672 &nambbc ! bottom temperature boundary condition (default: NO) … … 935 973 rn_emin0 = 1.e-4 ! surface minimum value of tke [m2/s2] 936 974 rn_bshear = 1.e-20 ! background shear (>0) currently a numerical threshold (do not change it) 975 nn_pdl = 1 ! Prandtl number function of richarson number (=1, avt=pdl(Ri)*avm) or not (=0, avt=avm) 937 976 nn_mxl = 2 ! mixing length: = 0 bounded by the distance to surface and bottom 938 977 ! = 1 bounded by the local vertical scale factor 939 978 ! = 2 first vertical derivative of mixing length bounded by 1 940 979 ! = 3 as =2 with distinct disspipative an mixing length scale 941 nn_pdl = 1 ! Prandtl number function of richarson number (=1, avt=pdl(Ri)*avm) or not (=0, avt=avm)942 980 ln_mxl0 = .true. ! surface mixing length scale = F(wind stress) (T) or not (F) 943 981 rn_mxl0 = 0.04 ! surface buoyancy lenght scale minimum value 982 !!gm ln_drg = .false. ! top/bottom friction added as boundary condition of TKE 944 983 ln_lc = .true. ! Langmuir cell parameterisation (Axell 2002) 945 984 rn_lc = 0.15 ! coef. associated to Langmuir cells -
branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/CRS/crs.F90
r7953 r8093 140 140 ! Physical and dynamical ocean fields for output or passing to TOP, time-mean fields 141 141 REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: tsn_crs 142 REAL(wp), DIMENSION(:,:,:) , ALLOCATABLE :: un_crs, vn_crs, wn_crs , rke_crs142 REAL(wp), DIMENSION(:,:,:) , ALLOCATABLE :: un_crs, vn_crs, wn_crs 143 143 REAL(wp), DIMENSION(:,:,:) , ALLOCATABLE :: hdivn_crs 144 144 REAL(wp), DIMENSION(:,:) , ALLOCATABLE :: sshn_crs … … 227 227 228 228 229 ALLOCATE( un_crs(jpi_crs,jpj_crs,jpk) , vn_crs(jpi_crs,jpj_crs,jpk) , & 230 & wn_crs(jpi_crs,jpj_crs,jpk) , hdivn_crs(jpi_crs,jpj_crs,jpk),& 231 & rke_crs(jpi_crs,jpj_crs,jpk), STAT=ierr(11)) 229 ALLOCATE( un_crs(jpi_crs,jpj_crs,jpk) , vn_crs(jpi_crs,jpj_crs,jpk) , & 230 & wn_crs(jpi_crs,jpj_crs,jpk) , hdivn_crs(jpi_crs,jpj_crs,jpk), STAT=ierr(11)) 232 231 233 232 ALLOCATE( sshn_crs(jpi_crs,jpj_crs), emp_crs (jpi_crs,jpj_crs), emp_b_crs(jpi_crs,jpj_crs), & -
branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/CRS/crsfld.F90
r7953 r8093 58 58 INTEGER :: ji, jj, jk ! dummy loop indices 59 59 REAL(wp) :: z2dcrsu, z2dcrsv ! local scalars 60 REAL(wp) :: zztmp ! - - 60 61 ! 61 62 REAL(wp), POINTER, DIMENSION(:,:,:) :: ze3t, ze3u, ze3v, ze3w ! 3D workspace for e3 62 REAL(wp), POINTER, DIMENSION(:,:,:) :: zt, zt_crs 63 REAL(wp), POINTER, DIMENSION(:,:,:) :: zt, zt_crs, z3d 63 64 REAL(wp), POINTER, DIMENSION(:,:,:) :: zs, zs_crs 64 65 !!---------------------------------------------------------------------- … … 69 70 CALL wrk_alloc( jpi,jpj,jpk, ze3t, ze3w ) 70 71 CALL wrk_alloc( jpi,jpj,jpk, ze3u, ze3v ) 71 CALL wrk_alloc( jpi,jpj,jpk, zt , zs )72 CALL wrk_alloc( jpi,jpj,jpk, zt , zs ,z3d ) 72 73 ! 73 74 CALL wrk_alloc( jpi_crs,jpj_crs,jpk, zt_crs, zs_crs ) … … 86 87 avs_crs (:,:,: ) = 0._wp ! avt 87 88 hdivn_crs(:,:,: ) = 0._wp ! hdiv 88 rke_crs (:,:,: ) = 0._wp ! rke89 89 sshn_crs (:,: ) = 0._wp ! ssh 90 90 utau_crs (:,: ) = 0._wp ! taux … … 158 158 CALL iom_put( "voces" , zs_crs ) ! vS 159 159 160 161 ! Kinetic energy 162 CALL crs_dom_ope( rke, 'VOL', 'T', tmask, rke_crs, p_e12=e1e2t, p_e3=ze3t, psgn=1.0 ) 163 CALL iom_put( "eken", rke_crs ) 164 160 IF( iom_use( "eken") ) THEN ! kinetic energy 161 z3d(:,:,jk) = 0._wp 162 DO jk = 1, jpkm1 163 DO jj = 2, jpjm1 164 DO ji = fs_2, fs_jpim1 ! vector opt. 165 zztmp = r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) 166 z3d(ji,jj,jk) = 0.25_wp * zztmp * ( & 167 & un(ji-1,jj,jk)**2 * e2u(ji-1,jj) * e3u_n(ji-1,jj,jk) & 168 & + un(ji ,jj,jk)**2 * e2u(ji ,jj) * e3u_n(ji ,jj,jk) & 169 & + vn(ji,jj-1,jk)**2 * e1v(ji,jj-1) * e3v_n(ji,jj-1,jk) & 170 & + vn(ji,jj ,jk)**2 * e1v(ji,jj ) * e3v_n(ji,jj ,jk) ) 171 END DO 172 END DO 173 END DO 174 CALL lbc_lnk( z3d, 'T', 1. ) 175 ! 176 CALL crs_dom_ope( z3d, 'VOL', 'T', tmask, zt_crs, p_e12=e1e2t, p_e3=ze3t, psgn=1.0 ) 177 CALL iom_put( "eken", zt_crs ) 178 ENDIF 165 179 ! Horizontal divergence ( following OPA_SRC/DYN/divhor.F90 ) 166 180 DO jk = 1, jpkm1 … … 175 189 hdivn_crs(ji,jj,jk) = ( z2dcrsu + z2dcrsv ) / crs_volt_wgt(ji,jj,jk) 176 190 ENDIF 177 END DO178 END DO179 END DO191 END DO 192 END DO 193 END DO 180 194 CALL crs_lbc_lnk( hdivn_crs, 'T', 1.0 ) 181 195 ! -
branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90
r7953 r8093 29 29 USE dynadv, ONLY: ln_dynadv_vec 30 30 USE zdf_oce ! ocean vertical physics 31 USE zdfdrg ! ocean vertical physics: top/bottom friction 31 32 USE ldftra ! lateral physics: eddy diffusivity coef. 32 33 USE ldfdyn ! lateral physics: eddy viscosity coef. … … 119 120 !! ** Method : use iom_put 120 121 !!---------------------------------------------------------------------- 121 !!122 122 INTEGER, INTENT( in ) :: kt ! ocean time-step index 123 123 !! 124 INTEGER :: ji, jj, jk! dummy loop indices125 INTEGER :: jkbot !126 REAL(wp) :: zztmp, zztmpx, zztmpy !127 !!128 REAL(wp), POINTER, DIMENSION(:,:) :: z2d! 2D workspace129 REAL(wp), POINTER, DIMENSION(:,:,:) :: z3d! 3D workspace124 INTEGER :: ji, jj, jk ! dummy loop indices 125 INTEGER :: ikbot ! local integer 126 REAL(wp):: zztmp , zztmpx ! local scalar 127 REAL(wp):: zztmp2, zztmpy ! - - 128 REAL(wp), DIMENSION(jpi,jpj) :: z2d ! 2D workspace 129 REAL(wp), DIMENSION(jpi,jpj,jpk) :: z3d ! 3D workspace 130 130 !!---------------------------------------------------------------------- 131 131 ! 132 132 IF( nn_timing == 1 ) CALL timing_start('dia_wri') 133 133 ! 134 CALL wrk_alloc( jpi , jpj , z2d )135 CALL wrk_alloc( jpi , jpj, jpk , z3d )136 !137 134 ! Output the initial state and forcings 138 135 IF( ninist == 1 ) THEN … … 162 159 DO jj = 1, jpj 163 160 DO ji = 1, jpi 164 jkbot = mbkt(ji,jj)165 z2d(ji,jj) = tsn(ji,jj, jkbot,jp_tem)161 ikbot = mbkt(ji,jj) 162 z2d(ji,jj) = tsn(ji,jj,ikbot,jp_tem) 166 163 END DO 167 164 END DO … … 174 171 DO jj = 1, jpj 175 172 DO ji = 1, jpi 176 jkbot = mbkt(ji,jj)177 z2d(ji,jj) = tsn(ji,jj, jkbot,jp_sal)173 ikbot = mbkt(ji,jj) 174 z2d(ji,jj) = tsn(ji,jj,ikbot,jp_sal) 178 175 END DO 179 176 END DO … … 182 179 183 180 IF ( iom_use("taubot") ) THEN ! bottom stress 181 zztmp = rau0 * 0.25 184 182 z2d(:,:) = 0._wp 185 183 DO jj = 2, jpjm1 186 184 DO ji = fs_2, fs_jpim1 ! vector opt. 185 !!gm old 186 !!gm BUG missing x 0.5 187 187 zztmpx = ( bfrua(ji ,jj) * un(ji ,jj,mbku(ji ,jj)) & 188 188 & + bfrua(ji-1,jj) * un(ji-1,jj,mbku(ji-1,jj)) ) … … 190 190 & + bfrva(ji,jj-1) * vn(ji,jj-1,mbkv(ji,jj-1)) ) 191 191 z2d(ji,jj) = rau0 * SQRT( zztmpx * zztmpx + zztmpy * zztmpy ) * tmask(ji,jj,1) 192 !!gm 193 zztmp2 = ( ( rCdU_bot(ji+1,jj)+rCdU_bot(ji ,jj) ) * un(ji ,jj,mbku(ji ,jj)) )**2 & 194 & + ( ( rCdU_bot(ji ,jj)+rCdU_bot(ji-1,jj) ) * un(ji-1,jj,mbku(ji-1,jj)) )**2 & 195 & + ( ( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj ) ) * vn(ji,jj ,mbkv(ji,jj )) )**2 & 196 & + ( ( rCdU_bot(ji,jj )+rCdU_bot(ji,jj-1) ) * vn(ji,jj-1,mbkv(ji,jj-1)) )**2 197 z2d(ji,jj) = zztmp * SQRT( zztmp2 ) * tmask(ji,jj,1) 198 !!gm new end 192 199 ! 193 200 ENDDO … … 202 209 DO jj = 1, jpj 203 210 DO ji = 1, jpi 204 jkbot = mbku(ji,jj)205 z2d(ji,jj) = un(ji,jj, jkbot)211 ikbot = mbku(ji,jj) 212 z2d(ji,jj) = un(ji,jj,ikbot) 206 213 END DO 207 214 END DO … … 214 221 DO jj = 1, jpj 215 222 DO ji = 1, jpi 216 jkbot = mbkv(ji,jj)217 z2d(ji,jj) = vn(ji,jj, jkbot)223 ikbot = mbkv(ji,jj) 224 z2d(ji,jj) = vn(ji,jj,ikbot) 218 225 END DO 219 226 END DO … … 281 288 ! 282 289 IF ( iom_use("eken") ) THEN 283 rke(:,:,jk) = 0._wp ! kinetic energy290 z3d(:,:,jk) = 0._wp ! kinetic energy 284 291 DO jk = 1, jpkm1 285 292 DO jj = 2, jpjm1 286 293 DO ji = fs_2, fs_jpim1 ! vector opt. 287 zztmp = 1._wp / ( e1e2t(ji,jj) * e3t_n(ji,jj,jk) ) 288 zztmpx = 0.5 * ( un(ji-1,jj,jk) * un(ji-1,jj,jk) * e2u(ji-1,jj) * e3u_n(ji-1,jj,jk) & 289 & + un(ji ,jj,jk) * un(ji ,jj,jk) * e2u(ji ,jj) * e3u_n(ji ,jj,jk) ) & 290 & * zztmp 291 ! 292 zztmpy = 0.5 * ( vn(ji,jj-1,jk) * vn(ji,jj-1,jk) * e1v(ji,jj-1) * e3v_n(ji,jj-1,jk) & 293 & + vn(ji,jj ,jk) * vn(ji,jj ,jk) * e1v(ji,jj ) * e3v_n(ji,jj ,jk) ) & 294 & * zztmp 295 ! 296 rke(ji,jj,jk) = 0.5_wp * ( zztmpx + zztmpy ) 297 ! 298 ENDDO 299 ENDDO 300 ENDDO 301 CALL lbc_lnk( rke, 'T', 1. ) 302 CALL iom_put( "eken", rke ) 294 zztmp = r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) 295 z3d(ji,jj,jk) = 0.25_wp * zztmp * ( & 296 & un(ji-1,jj,jk)**2 * e2u(ji-1,jj) * e3u_n(ji-1,jj,jk) & 297 & + un(ji ,jj,jk)**2 * e2u(ji ,jj) * e3u_n(ji ,jj,jk) & 298 & + vn(ji,jj-1,jk)**2 * e1v(ji,jj-1) * e3v_n(ji,jj-1,jk) & 299 & + vn(ji,jj ,jk)**2 * e1v(ji,jj ) * e3v_n(ji,jj ,jk) ) 300 END DO 301 END DO 302 END DO 303 CALL lbc_lnk( z3d, 'T', 1. ) 304 CALL iom_put( "eken", z3d ) 303 305 ENDIF 304 306 ! … … 407 409 CALL iom_put( "bn2", rn2 ) !Brunt-Vaisala buoyancy frequency (N^2) 408 410 ! 409 CALL wrk_dealloc( jpi , jpj , z2d ) 410 CALL wrk_dealloc( jpi , jpj, jpk , z3d ) 411 ! 412 ! If we want tmb values 413 414 IF (ln_diatmb) THEN 411 412 IF (ln_diatmb) THEN ! If we want tmb values 415 413 CALL dia_tmb 416 414 ENDIF -
branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/DYN/dynbfr.F90
r7753 r8093 15 15 USE dom_oce ! ocean space and time domain variables 16 16 USE zdf_oce ! ocean vertical physics variables 17 !!gm new 18 USE zdfdrg ! vertical physics: top/bottom drag coef. 19 !!gm old 17 20 USE zdfbfr ! ocean bottom friction variables 21 !!gm 18 22 USE trd_oce ! trends: ocean variables 19 23 USE trddyn ! trend manager: dynamics 24 ! 20 25 USE in_out_manager ! I/O manager 21 26 USE prtctl ! Print control … … 50 55 INTEGER :: ikbu, ikbv ! local integers 51 56 REAL(wp) :: zm1_2dt ! local scalar 52 REAL(wp), POINTER, DIMENSION(:,:,:) :: ztrdu, ztrdv 57 REAL(wp) :: zCdu, zCdv ! - - 58 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ztrdu, ztrdv 53 59 !!--------------------------------------------------------------------- 54 60 ! … … 61 67 62 68 !!gm bug : time step is only rdt (not 2 rdt if euler start !) 63 zm1_2dt = - 1._wp / ( 2._wp * rdt )69 zm1_2dt = - 1._wp / ( 2._wp * rdt ) 64 70 65 IF( l_trddyn ) THEN ! trends: store the input trends66 CALL wrk_alloc( jpi,jpj,jpk, ztrdu, ztrdv)67 ztrdu(:,:,:) = ua(:,:,:)68 ztrdv(:,:,:) = va(:,:,:)69 ENDIF71 IF( l_trddyn ) THEN ! trends: store the input trends 72 ALLOCATE( ztrdu(jpi,jpj,jpk) , ztrdv(jpi,jpj,jpk) ) 73 ztrdu(:,:,:) = ua(:,:,:) 74 ztrdv(:,:,:) = va(:,:,:) 75 ENDIF 70 76 71 77 72 DO jj = 2, jpjm1 73 DO ji = 2, jpim1 74 ikbu = mbku(ji,jj) ! deepest ocean u- & v-levels 75 ikbv = mbkv(ji,jj) 76 ! 77 ! Apply stability criteria on absolute value : abs(bfr/e3) < 1/(2dt) => bfr/e3 > -1/(2dt) 78 ua(ji,jj,ikbu) = ua(ji,jj,ikbu) + MAX( bfrua(ji,jj) / e3u_n(ji,jj,ikbu) , zm1_2dt ) * ub(ji,jj,ikbu) 79 va(ji,jj,ikbv) = va(ji,jj,ikbv) + MAX( bfrva(ji,jj) / e3v_n(ji,jj,ikbv) , zm1_2dt ) * vb(ji,jj,ikbv) 80 END DO 81 END DO 82 ! 83 IF( ln_isfcav ) THEN ! ocean cavities 84 DO jj = 2, jpjm1 85 DO ji = 2, jpim1 86 ! (ISF) stability criteria for top friction 87 ikbu = miku(ji,jj) ! first wet ocean u- & v-levels 88 ikbv = mikv(ji,jj) 89 ! 90 ! Apply stability criteria on absolute value : abs(bfr/e3) < 1/(2dt) => bfr/e3 > -1/(2dt) 78 DO jj = 2, jpjm1 79 DO ji = 2, jpim1 80 ikbu = mbku(ji,jj) ! deepest wet ocean u- & v-levels 81 ikbv = mbkv(ji,jj) 82 ! 83 ! Apply stability criteria on absolute value : abs(bfr/e3) < 1/(2dt) => bfr/e3 > -1/(2dt) 84 !!gm old 85 ua(ji,jj,ikbu) = ua(ji,jj,ikbu) + MAX( bfrua(ji,jj) / e3u_n(ji,jj,ikbu) , zm1_2dt ) * ub(ji,jj,ikbu) 86 va(ji,jj,ikbv) = va(ji,jj,ikbv) + MAX( bfrva(ji,jj) / e3v_n(ji,jj,ikbv) , zm1_2dt ) * vb(ji,jj,ikbv) 87 !!gm new 88 ! zCdu = 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) / e3u_n(ji,jj,ikbu) 89 ! zCdv = 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) / e3v_n(ji,jj,ikbv) 90 ! ! 91 ! ua(ji,jj,ikbu) = ua(ji,jj,ikbu) + MAX( zCdu , zm1_2dt ) * ub(ji,jj,ikbu) 92 ! va(ji,jj,ikbv) = va(ji,jj,ikbv) + MAX( zCdv , zm1_2dt ) * vb(ji,jj,ikbv) 93 !!gm 94 END DO 95 END DO 96 ! 97 IF( ln_isfcav ) THEN ! ocean cavities 98 DO jj = 2, jpjm1 99 DO ji = 2, jpim1 100 ikbu = miku(ji,jj) ! first wet ocean u- & v-levels 101 ikbv = mikv(ji,jj) 102 ! 103 ! Apply stability criteria on absolute value : abs(bfr/e3) < 1/(2dt) => bfr/e3 > -1/(2dt) 104 !!gm old 91 105 ua(ji,jj,ikbu) = ua(ji,jj,ikbu) + MAX( tfrua(ji,jj) / e3u_n(ji,jj,ikbu) , zm1_2dt ) * ub(ji,jj,ikbu) & 92 106 & * (1.-umask(ji,jj,1)) 93 107 va(ji,jj,ikbv) = va(ji,jj,ikbv) + MAX( tfrva(ji,jj) / e3v_n(ji,jj,ikbv) , zm1_2dt ) * vb(ji,jj,ikbv) & 94 108 & * (1.-vmask(ji,jj,1)) 95 ! (ISF) 109 !!gm new 110 ! zCdu = 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) / e3u_n(ji,jj,ikbu) ! NB: Cdtop masked 111 ! zCdv = 0.5*( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) / e3v_n(ji,jj,ikbv) 112 ! ! 113 ! ua(ji,jj,ikbu) = ua(ji,jj,ikbu) + MAX( zCdu , zm1_2dt ) * ub(ji,jj,ikbu) 114 ! va(ji,jj,ikbv) = va(ji,jj,ikbv) + MAX( zCdv , zm1_2dt ) * vb(ji,jj,ikbv) 115 !!gm 96 116 END DO 97 117 END DO 98 END IF118 END IF 99 119 ! 100 120 IF( l_trddyn ) THEN ! trends: send trends to trddyn for further diagnostics … … 102 122 ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 103 123 CALL trd_dyn( ztrdu(:,:,:), ztrdv(:,:,:), jpdyn_bfr, kt ) 104 CALL wrk_dealloc( jpi,jpj,jpk,ztrdu, ztrdv )124 DEALLOCATE( ztrdu, ztrdv ) 105 125 ENDIF 106 126 ! ! print mean trends (used for debugging) -
branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90
r7761 r8093 1457 1457 !! *** ROUTINE interp1 *** 1458 1458 !! 1459 !! ** Purpose : Calculate the first order of deri avtive of1459 !! ** Purpose : Calculate the first order of derivative of 1460 1460 !! a cubic spline function y=a+b*x+c*x^2+d*x^3 1461 1461 !! -
branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90
r7831 r8093 28 28 USE sbc_oce ! surface boundary condition: ocean 29 29 USE zdf_oce ! Bottom friction coefts 30 !!gm new 31 USE zdfdrg ! vertical physics: top/bottom drag coef. 32 !!gm 30 33 USE sbcisf ! ice shelf variable (fwfisf) 31 34 USE sbcapr ! surface boundary condition: atmospheric pressure … … 146 149 REAL(wp) :: zraur, z1_2dt_b, z2dt_bf ! local scalars 147 150 REAL(wp) :: zx1, zy1, zx2, zy2 ! - - 148 REAL(wp) :: z1_12, z1_8, z1_4, z1_2 ! - -151 REAL(wp) :: z1_12, z1_8, z1_4, z1_2 ! - - 149 152 REAL(wp) :: zu_spg, zv_spg ! - - 150 REAL(wp) :: zhura, zhvra ! - - 151 REAL(wp) :: za0, za1, za2, za3 ! - - 152 ! 153 REAL(wp), POINTER, DIMENSION(:,:) :: zsshp2_e 154 REAL(wp), POINTER, DIMENSION(:,:) :: zu_trd, zv_trd, zu_frc, zv_frc, zssh_frc 155 REAL(wp), POINTER, DIMENSION(:,:) :: zwx, zwy, zhdiv 156 REAL(wp), POINTER, DIMENSION(:,:) :: zhup2_e, zhvp2_e, zhust_e, zhvst_e 157 REAL(wp), POINTER, DIMENSION(:,:) :: zsshu_a, zsshv_a 158 REAL(wp), POINTER, DIMENSION(:,:) :: zhf 159 REAL(wp), POINTER, DIMENSION(:,:) :: zcpx, zcpy ! Wetting/Dying gravity filter coef. 153 REAL(wp) :: zhura, zhvra ! - - 154 REAL(wp) :: za0, za1, za2, za3 ! - - 155 REAL(wp) :: zztmp ! - - 156 REAL(wp), DIMENSION(jpi,jpj) :: zsshp2_e 157 REAL(wp), DIMENSION(jpi,jpj) :: zu_trd, zv_trd, zu_frc, zv_frc, zssh_frc 158 REAL(wp), DIMENSION(jpi,jpj) :: zwx, zwy, zhdiv 159 REAL(wp), DIMENSION(jpi,jpj) :: zhup2_e, zhvp2_e, zhust_e, zhvst_e 160 REAL(wp), DIMENSION(jpi,jpj) :: zsshu_a, zsshv_a 161 REAL(wp), DIMENSION(jpi,jpj) :: zhf 162 REAL(wp), DIMENSION(jpi,jpj) :: zCdU_u, zCdU_v ! top/bottom stress at u- & v-points 163 ! 164 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zcpx, zcpy ! Wetting/Dying gravity filter coef. 160 165 !!---------------------------------------------------------------------- 161 166 ! 162 IF( nn_timing == 1 ) CALL timing_start('dyn_spg_ts') 163 ! 164 ! !* Allocate temporary arrays 165 CALL wrk_alloc( jpi,jpj, zsshp2_e, zhdiv ) 166 CALL wrk_alloc( jpi,jpj, zu_trd, zv_trd) 167 CALL wrk_alloc( jpi,jpj, zwx, zwy, zssh_frc, zu_frc, zv_frc) 168 CALL wrk_alloc( jpi,jpj, zhup2_e, zhvp2_e, zhust_e, zhvst_e) 169 CALL wrk_alloc( jpi,jpj, zsshu_a, zsshv_a ) 170 CALL wrk_alloc( jpi,jpj, zhf ) 171 IF( ln_wd ) CALL wrk_alloc( jpi, jpj, zcpx, zcpy ) 167 IF( nn_timing == 1 ) CALL timing_start('dyn_spg_ts') 168 ! 169 IF( ln_wd ) ALLOCATE( zcpx(jpi,jpj), zcpy(jpi,jpj) ) 172 170 ! 173 171 zmdi=1.e+20 ! missing data indicator for masking … … 211 209 ! 212 210 ENDIF 211 ! 212 !!gm old/new 213 IF( ln_isfcav ) THEN ! top+bottom friction (ocean cavities) 214 zCdU_u(:,:) = bfrua(:,:) + tfrua(:,:) 215 zCdU_v(:,:) = bfrva(:,:) + tfrva(:,:) 216 ELSE ! bottom friction only 217 zCdU_u(:,:) = bfrua(:,:) 218 zCdU_v(:,:) = bfrva(:,:) 219 ENDIF 220 !!gm new 221 ! IF( ln_isfcav ) THEN ! top+bottom friction (ocean cavities) 222 ! DO jj = 2, jpjm1 223 ! DO ji = fs_2, fs_jpim1 ! vector opt. 224 ! zCdU_u(ji,jj) = 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) + rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) 225 ! zCdU_v(ji,jj) = 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) + rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) 226 ! END DO 227 ! END DO 228 ! ELSE ! bottom friction only 229 ! DO jj = 2, jpjm1 230 ! DO ji = fs_2, fs_jpim1 ! vector opt. 231 ! zCdU_u(ji,jj) = 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) 232 ! zCdU_v(ji,jj) = 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) 233 ! END DO 234 ! END DO 235 ! ENDIF 236 !!gm 213 237 ! 214 238 ! Set arrays to remove/compute coriolis trend. … … 415 439 ! no worries about sshn(ji+1,jj) - sshn(ji ,jj) = 0, it won't happen ! here 416 440 zcpx(ji,jj) = ABS( (sshn(ji+1,jj) + ht_wd(ji+1,jj) - sshn(ji,jj) - ht_wd(ji,jj)) & 417 &/ (sshn(ji+1,jj) - sshn(ji ,jj)) )441 & / (sshn(ji+1,jj) - sshn(ji ,jj)) ) 418 442 ELSE 419 443 zcpx(ji,jj) = 0._wp … … 491 515 ! Note that the "unclipped" bottom friction parameter is used even with explicit drag 492 516 IF( ln_wd ) THEN 493 zu_frc(:,:) = zu_frc(:,:) + MAX(r1_hu_n(:,:) * bfrua(:,:),-1._wp / rdtbt) * zwx(:,:) 494 zv_frc(:,:) = zv_frc(:,:) + MAX(r1_hv_n(:,:) * bfrva(:,:),-1._wp / rdtbt) * zwy(:,:) 517 zztmp = - 1._wp / rdtbt 518 DO jj = 2, jpjm1 519 DO ji = fs_2, fs_jpim1 ! vector opt. 520 !!gm old 521 zu_frc(ji,jj) = zu_frc(ji,jj) + MAX( r1_hu_n(ji,jj) * bfrua(ji,jj) , zztmp ) * zwx(ji,jj) 522 zv_frc(ji,jj) = zv_frc(ji,jj) + MAX( r1_hv_n(ji,jj) * bfrva(ji,jj) , zztmp ) * zwy(ji,jj) 523 !!gm new 524 ! zu_frc(ji,jj) = zu_frc(ji,jj) + MAX( r1_hu_n(ji,jj) * 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) , zztmp ) * zwx(ji,jj) 525 ! zv_frc(ji,jj) = zv_frc(ji,jj) + MAX( r1_hv_n(ji,jj) * 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) , zztmp ) * zwy(ji,jj) 526 !!gm 527 END DO 528 END DO 495 529 ELSE 496 zu_frc(:,:) = zu_frc(:,:) + r1_hu_n(:,:) * bfrua(:,:) * zwx(:,:) 497 zv_frc(:,:) = zv_frc(:,:) + r1_hv_n(:,:) * bfrva(:,:) * zwy(:,:) 530 DO jj = 2, jpjm1 531 DO ji = fs_2, fs_jpim1 ! vector opt. 532 !!gm old 533 zu_frc(ji,jj) = zu_frc(ji,jj) + r1_hu_n(ji,jj) * bfrua(ji,jj) * zwx(ji,jj) 534 zv_frc(ji,jj) = zv_frc(ji,jj) + r1_hv_n(ji,jj) * bfrva(ji,jj) * zwy(ji,jj) 535 !!gm new 536 ! zu_frc(ji,jj) = zu_frc(ji,jj) + r1_hu_n(ji,jj) * 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) * zwx(ji,jj) 537 ! zv_frc(ji,jj) = zv_frc(ji,jj) + r1_hv_n(ji,jj) * 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) * zwy(ji,jj) 538 !!gm 539 END DO 540 END DO 498 541 END IF 499 542 ! … … 520 563 ! 521 564 ! Note that the "unclipped" top friction parameter is used even with explicit drag 522 zu_frc(:,:) = zu_frc(:,:) + r1_hu_n(:,:) * tfrua(:,:) * zwx(:,:) 523 zv_frc(:,:) = zv_frc(:,:) + r1_hv_n(:,:) * tfrva(:,:) * zwy(:,:) 565 DO jj = 2, jpjm1 566 DO ji = fs_2, fs_jpim1 ! vector opt. 567 zu_frc(ji,jj) = zu_frc(ji,jj) + r1_hu_n(ji,jj) * tfrua(ji,jj) * zwx(ji,jj) 568 zv_frc(ji,jj) = zv_frc(ji,jj) + r1_hv_n(ji,jj) * tfrva(ji,jj) * zwy(ji,jj) 569 END DO 570 END DO 524 571 ! 525 572 IF (ln_bt_fw) THEN ! Add wind forcing 526 zu_frc(:,:) = zu_frc(:,:) + zraur * utau(:,:) * r1_hu_n(:,:) 527 zv_frc(:,:) = zv_frc(:,:) + zraur * vtau(:,:) * r1_hv_n(:,:) 573 DO jj = 2, jpjm1 574 DO ji = fs_2, fs_jpim1 ! vector opt. 575 zu_frc(ji,jj) = zu_frc(ji,jj) + zraur * utau(ji,jj) * r1_hu_n(ji,jj) 576 zv_frc(ji,jj) = zv_frc(ji,jj) + zraur * vtau(ji,jj) * r1_hv_n(ji,jj) 577 END DO 578 END DO 528 579 ELSE 529 zu_frc(:,:) = zu_frc(:,:) + zraur * z1_2 * ( utau_b(:,:) + utau(:,:) ) * r1_hu_n(:,:) 530 zv_frc(:,:) = zv_frc(:,:) + zraur * z1_2 * ( vtau_b(:,:) + vtau(:,:) ) * r1_hv_n(:,:) 580 DO jj = 2, jpjm1 581 DO ji = fs_2, fs_jpim1 ! vector opt. 582 zu_frc(ji,jj) = zu_frc(ji,jj) + zraur * z1_2 * ( utau_b(ji,jj) + utau(ji,jj) ) * r1_hu_n(ji,jj) 583 zv_frc(ji,jj) = zv_frc(ji,jj) + zraur * z1_2 * ( vtau_b(ji,jj) + vtau(ji,jj) ) * r1_hv_n(ji,jj) 584 END DO 585 END DO 531 586 ENDIF 532 587 ! … … 885 940 ENDIF 886 941 ! 887 ! Add bottom stresses: 888 zu_trd(:,:) = zu_trd(:,:) + bfrua(:,:) * un_e(:,:) * hur_e(:,:) 889 zv_trd(:,:) = zv_trd(:,:) + bfrva(:,:) * vn_e(:,:) * hvr_e(:,:) 890 ! 891 ! Add top stresses: 892 zu_trd(:,:) = zu_trd(:,:) + tfrua(:,:) * un_e(:,:) * hur_e(:,:) 893 zv_trd(:,:) = zv_trd(:,:) + tfrva(:,:) * vn_e(:,:) * hvr_e(:,:) 942 DO jj = 2, jpjm1 943 DO ji = fs_2, fs_jpim1 ! vector opt. 944 ! Add top/bottom stresses: 945 !!gm old/new 946 zu_trd(ji,jj) = zu_trd(ji,jj) + zCdU_u(ji,jj) * un_e(ji,jj) * hur_e(ji,jj) 947 zv_trd(ji,jj) = zv_trd(ji,jj) + zCdU_v(ji,jj) * vn_e(ji,jj) * hvr_e(ji,jj) 948 !!gm 949 END DO 950 END DO 894 951 ! 895 952 ! Surface pressure trend: … … 1091 1148 IF( lrst_oce .AND.ln_bt_fw ) CALL ts_rst( kt, 'WRITE' ) 1092 1149 ! 1093 CALL wrk_dealloc( jpi,jpj, zsshp2_e, zhdiv ) 1094 CALL wrk_dealloc( jpi,jpj, zu_trd, zv_trd ) 1095 CALL wrk_dealloc( jpi,jpj, zwx, zwy, zssh_frc, zu_frc, zv_frc ) 1096 CALL wrk_dealloc( jpi,jpj, zhup2_e, zhvp2_e, zhust_e, zhvst_e ) 1097 CALL wrk_dealloc( jpi,jpj, zsshu_a, zsshv_a ) 1098 CALL wrk_dealloc( jpi,jpj, zhf ) 1099 IF( ln_wd ) CALL wrk_dealloc( jpi, jpj, zcpx, zcpy ) 1150 IF( ln_wd ) DEALLOCATE( zcpx, zcpy ) 1100 1151 ! 1101 1152 IF ( ln_diatmb ) THEN -
branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_imp.F90
r7990 r8093 22 22 USE dynadv , ONLY: ln_dynadv_vec ! Momentum advection form 23 23 USE zdf_oce ! ocean vertical physics 24 !!gm new 25 USE zdfdrg ! vertical physics: top/bottom drag coef. 26 !!gm old 24 27 USE zdfbfr ! Bottom friction setup 28 !!gm 25 29 ! 26 30 USE in_out_manager ! I/O manager … … 116 120 ikbu = mbku(ji,jj) ! ocean bottom level at u- and v-points 117 121 ikbv = mbkv(ji,jj) ! (deepest ocean u- and v-points) 122 !!gm old 118 123 avmu(ji,jj,ikbu+1) = -bfrua(ji,jj) * e3uw_n(ji,jj,ikbu+1) 119 124 avmv(ji,jj,ikbv+1) = -bfrva(ji,jj) * e3vw_n(ji,jj,ikbv+1) 125 !!gm new 126 ! avmu(ji,jj,ikbu+1) = -0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) * e3uw_n(ji,jj,ikbu+1) 127 ! avmv(ji,jj,ikbv+1) = -0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) * e3vw_n(ji,jj,ikbv+1) 128 !!gm 120 129 END DO 121 130 END DO … … 125 134 ikbu = miku(ji,jj) ! ocean top level at u- and v-points 126 135 ikbv = mikv(ji,jj) ! (first wet ocean u- and v-points) 136 !!gm old 127 137 IF( ikbu >= 2 ) avmu(ji,jj,ikbu) = -tfrua(ji,jj) * e3uw_n(ji,jj,ikbu) 128 138 IF( ikbv >= 2 ) avmv(ji,jj,ikbv) = -tfrva(ji,jj) * e3vw_n(ji,jj,ikbv) 139 !!gm new 140 ! top Cd is masked (=0 outside cavities) no need of test on mik>=2 141 ! avmu(ji,jj,ikbu) = -0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * e3uw_n(ji,jj,ikbu) 142 ! avmv(ji,jj,ikbv) = -0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * e3vw_n(ji,jj,ikbv) 143 !!gm 129 144 END DO 130 145 END DO … … 148 163 ze3ua = ( 1._wp - r_vvl ) * e3u_n(ji,jj,ikbu) + r_vvl * e3u_a(ji,jj,ikbu) 149 164 ze3va = ( 1._wp - r_vvl ) * e3v_n(ji,jj,ikbv) + r_vvl * e3v_a(ji,jj,ikbv) 165 !!gm old 150 166 ua(ji,jj,ikbu) = ua(ji,jj,ikbu) + p2dt * bfrua(ji,jj) * ua_b(ji,jj) / ze3ua 151 167 va(ji,jj,ikbv) = va(ji,jj,ikbv) + p2dt * bfrva(ji,jj) * va_b(ji,jj) / ze3va 168 !!gm new 169 ! ua(ji,jj,ikbu) = ua(ji,jj,ikbu) + p2dt * 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) * ua_b(ji,jj) / ze3ua 170 ! va(ji,jj,ikbv) = va(ji,jj,ikbv) + p2dt * 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) * va_b(ji,jj) / ze3va 171 !!gm 152 172 END DO 153 173 END DO … … 159 179 ze3ua = ( 1._wp - r_vvl ) * e3u_n(ji,jj,ikbu) + r_vvl * e3u_a(ji,jj,ikbu) 160 180 ze3va = ( 1._wp - r_vvl ) * e3v_n(ji,jj,ikbv) + r_vvl * e3v_a(ji,jj,ikbv) 181 !!gm old 161 182 ua(ji,jj,ikbu) = ua(ji,jj,ikbu) + p2dt * tfrua(ji,jj) * ua_b(ji,jj) / ze3ua 162 183 va(ji,jj,ikbv) = va(ji,jj,ikbv) + p2dt * tfrva(ji,jj) * va_b(ji,jj) / ze3va 184 !!gm new 185 ua(ji,jj,ikbu) = ua(ji,jj,ikbu) + p2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * ua_b(ji,jj) / ze3ua 186 va(ji,jj,ikbv) = va(ji,jj,ikbv) + p2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * va_b(ji,jj) / ze3va 187 !!gm 163 188 END DO 164 189 END DO -
branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/TRA/eosbn2.F90
r7753 r8093 928 928 pn2(ji,jj,jk) = grav * ( zaw * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) ) & 929 929 & - zbw * ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) ) ) & 930 & / e3w_n(ji,jj,jk) * tmask(ji,jj,jk)930 & / e3w_n(ji,jj,jk) * wmask(ji,jj,jk) 931 931 END DO 932 932 END DO -
branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdf_oce.F90
r7990 r8093 46 46 47 47 48 REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:) :: av t , avs !: vertical T & S diffusivity coef at w-pt[m2/s]49 REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:) :: avm !: vertical momentum viscosity coef at w-pt [m2/s] 48 REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:) :: avm, avt, avs !: vertical mixing coefficient (w-point) [m2/s] 49 !!gm 50 50 REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:) :: avmu , avmv !: vertical viscosity coef at uw- & vw-pts [m2/s] 51 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: avt_k , avm_k !: avt, avs computed by turbulent closure alone 51 !!gm 52 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: avm_k, avt_k !: avm, avt computed by turbulent closure alone 52 53 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: en !: now turbulent kinetic energy [m2/s2] 54 REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:) :: avtb_2d !: horizontal shape of background Kz profile 53 55 REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:) :: avmb , avtb !: background profile of avm and avt 54 REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:) :: avtb_2d !: horizontal shape of background Kz profile 56 !!gm 55 57 REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:) :: bfrua, bfrva !: bottom friction coefficients 56 58 REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:) :: tfrua, tfrva !: top friction coefficients 59 !!gm 57 60 58 61 !!---------------------------------------------------------------------- 59 !! NEMO/OPA 4.0 , NEMO Consortium (201 1)62 !! NEMO/OPA 4.0 , NEMO Consortium (2017) 60 63 !! $Id$ 61 64 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) -
branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfddm.F90
r7990 r8093 39 39 CONTAINS 40 40 41 SUBROUTINE zdf_ddm( kt )41 SUBROUTINE zdf_ddm( kt, p_avm, p_avt, p_avs ) 42 42 !!---------------------------------------------------------------------- 43 43 !! *** ROUTINE zdf_ddm *** … … 69 69 !! References : Merryfield et al., JPO, 29, 1124-1142, 1999. 70 70 !!---------------------------------------------------------------------- 71 INTEGER, INTENT(in) :: kt ! ocean time-step indexocean time step 71 INTEGER, INTENT(in ) :: kt ! ocean time-step indexocean time step 72 REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: p_avm ! Kz on momentum (w-points) 73 REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: p_avt ! Kz on temperature (w-points) 74 REAL(wp), DIMENSION(:,:,:), INTENT( out) :: p_avs ! Kz on salinity (w-points) 72 75 ! 73 76 INTEGER :: ji, jj , jk ! dummy loop indices … … 91 94 !!gm and many acces in memory 92 95 93 DO jj = 1, jpj ! R=zrau = (alpha / beta) (dk[t] / dk[s])96 DO jj = 1, jpj !== R=zrau = (alpha / beta) (dk[t] / dk[s]) ==! 94 97 DO ji = 1, jpi 95 98 zrw = ( gdepw_n(ji,jj,jk ) - gdept_n(ji,jj,jk) ) & … … 109 112 END DO 110 113 111 DO jj = 1, jpj ! indicators:114 DO jj = 1, jpj !== indicators ==! 112 115 DO ji = 1, jpi 113 116 ! stability indicator: msks=1 if rn2>0; 0 elsewhere … … 154 157 & + 0.15 * zrau(ji,jj) * zmskd2(ji,jj) ) 155 158 ! add to the eddy viscosity coef. previously computed 156 avs(ji,jj,jk) =avt(ji,jj,jk) + zavfs + zavds157 avt(ji,jj,jk) =avt(ji,jj,jk) + zavft + zavdt158 avm(ji,jj,jk) =avm(ji,jj,jk) + MAX( zavft + zavdt, zavfs + zavds )159 p_avs(ji,jj,jk) = p_avt(ji,jj,jk) + zavfs + zavds 160 p_avt(ji,jj,jk) = p_avt(ji,jj,jk) + zavft + zavdt 161 p_avm(ji,jj,jk) = p_avm(ji,jj,jk) + MAX( zavft + zavdt, zavfs + zavds ) 159 162 END DO 160 163 END DO -
branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfevd.F90
r7990 r8093 38 38 CONTAINS 39 39 40 SUBROUTINE zdf_evd( kt )40 SUBROUTINE zdf_evd( kt, p_avm, p_avt ) 41 41 !!---------------------------------------------------------------------- 42 42 !! *** ROUTINE zdf_evd *** … … 55 55 !! ** Action : avt, avm enhanced where static instability occurs 56 56 !!---------------------------------------------------------------------- 57 INTEGER, INTENT(in) :: kt ! ocean time-step indexocean time step 57 INTEGER , INTENT(in ) :: kt ! ocean time-step indexocean time step 58 REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: p_avm, p_avt ! momentum and tracer Kz (w-points) 58 59 ! 59 60 INTEGER :: ji, jj, jk ! dummy loop indices … … 71 72 ! 72 73 ! 73 zavt_evd(:,:,:) = avt(:,:,:)! set avt prior to evd application74 zavt_evd(:,:,:) = p_avt(:,:,:) ! set avt prior to evd application 74 75 ! 75 76 SELECT CASE ( nn_evdm ) … … 77 78 CASE ( 1 ) !== enhance tracer & momentum Kz ==! (if rn2<-1.e-12) 78 79 ! 79 zavm_evd(:,:,:) = avm(:,:,:)! set avm prior to evd application80 zavm_evd(:,:,:) = p_avm(:,:,:) ! set avm prior to evd application 80 81 ! 81 82 !! change last digits results 82 83 ! WHERE( MAX( rn2(2:jpi,2:jpj,2:jpkm1), rn2b(2:jpi,2:jpj,2:jpkm1) ) <= -1.e-12 ) THEN 83 ! avt(2:jpi,2:jpj,2:jpkm1) = rn_evd * wmask(2:jpi,2:jpj,2:jpkm1)84 ! avm(2:jpi,2:jpj,2:jpkm1) = rn_evd * wmask(2:jpi,2:jpj,2:jpkm1)84 ! p_avt(2:jpi,2:jpj,2:jpkm1) = rn_evd * wmask(2:jpi,2:jpj,2:jpkm1) 85 ! p_avm(2:jpi,2:jpj,2:jpkm1) = rn_evd * wmask(2:jpi,2:jpj,2:jpkm1) 85 86 ! END WHERE 86 87 ! … … 89 90 DO ji = 2, jpim1 90 91 IF( MIN( rn2(ji,jj,jk), rn2b(ji,jj,jk) ) <= -1.e-12 ) THEN 91 avt(ji,jj,jk) = rn_evd * wmask(ji,jj,jk)92 avm(ji,jj,jk) = rn_evd * wmask(ji,jj,jk)92 p_avt(ji,jj,jk) = rn_evd * wmask(ji,jj,jk) 93 p_avm(ji,jj,jk) = rn_evd * wmask(ji,jj,jk) 93 94 ENDIF 94 95 END DO … … 96 97 END DO 97 98 ! 98 zavm_evd(:,:,:) = avm(:,:,:) - zavm_evd(:,:,:) ! change in avm due to evd99 CALL iom_put( "avm_evd", zavm_evd ) ! output this change99 zavm_evd(:,:,:) = p_avm(:,:,:) - zavm_evd(:,:,:) ! change in avm due to evd 100 CALL iom_put( "avm_evd", zavm_evd ) ! output this change 100 101 ! 101 102 CASE DEFAULT !== enhance tracer Kz ==! (if rn2<-1.e-12) 102 103 !! change last digits results 103 104 ! WHERE( MAX( rn2(2:jpi,2:jpj,2:jpkm1), rn2b(2:jpi,2:jpj,2:jpkm1) ) <= -1.e-12 ) 104 ! avt(2:jpi,2:jpj,2:jpkm1) = rn_evd * wmask(2:jpi,2:jpj,2:jpkm1)105 ! p_avt(2:jpi,2:jpj,2:jpkm1) = rn_evd * wmask(2:jpi,2:jpj,2:jpkm1) 105 106 ! END WHERE 106 107 … … 109 110 DO ji = 2, jpim1 110 111 IF( MIN( rn2(ji,jj,jk), rn2b(ji,jj,jk) ) <= -1.e-12 ) & 111 avt(ji,jj,jk) = rn_evd * wmask(ji,jj,jk)112 p_avt(ji,jj,jk) = rn_evd * wmask(ji,jj,jk) 112 113 END DO 113 114 END DO … … 116 117 END SELECT 117 118 ! 118 zavt_evd(:,:,:) = avt(:,:,:) - zavt_evd(:,:,:) ! change in avt due to evd119 zavt_evd(:,:,:) = p_avt(:,:,:) - zavt_evd(:,:,:) ! change in avt due to evd 119 120 CALL iom_put( "avt_evd", zavt_evd ) ! output this change 120 121 IF( l_trdtra ) CALL trd_tra( kt, 'TRA', jp_tem, jptra_evd, zavt_evd ) -
branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfgls.F90
r7991 r8093 19 19 USE domvvl ! ocean space and time domain : variable volume layer 20 20 USE zdf_oce ! ocean vertical physics 21 USE zdfsh2 ! vertical physics: shear production term of TKE 22 USE zdfbfr ! bottom friction (only for rn_bfrz0) 21 !!gm old 22 USE zdfbfr , ONLY : rn_tfrz0, rn_bfrz0 ! top/bottom roughness 23 !!gm new 24 USE zdfdrg , ONLY : r_z0_top , r_z0_bot ! top/bottom roughness 25 USE zdfdrg , ONLY : rCdU_top , rCdU_bot ! top/bottom friction 26 !!gm 23 27 USE sbc_oce ! surface boundary condition: ocean 24 28 USE phycst ! physical constants 25 29 USE zdfmxl ! mixed layer 26 USE sbcwave , ONLY: hsw ! significant wave height30 USE sbcwave , ONLY : hsw ! significant wave height 27 31 ! 28 32 USE lbclnk ! ocean lateral boundary conditions (or mpp link) … … 45 49 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: hmxl_n !: now mixing length 46 50 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: zwall !: wall function 47 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ustars2 !: Squared surface velocity scale at T-points 48 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ustarb2 !: Squared bottom velocity scale at T-points 51 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ustar2_surf !: Squared surface velocity scale at T-points 52 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ustar2_top !: Squared top velocity scale at T-points 53 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ustar2_bot !: Squared bottom velocity scale at T-points 49 54 50 55 ! !! ** Namelist namzdf_gls ** … … 101 106 REAL(wp) :: rsc_tke, rsc_psi, rpsi1, rpsi2, rpsi3, rsc_psi0 ! - - - - 102 107 REAL(wp) :: rpsi3m, rpsi3p, rpp, rmm, rnn ! - - - - 108 ! 109 REAL(wp) :: r2_3 = 2._wp/3._wp ! constant=2/3 103 110 104 111 !! * Substitutions … … 115 122 !! *** FUNCTION zdf_gls_alloc *** 116 123 !!---------------------------------------------------------------------- 117 ALLOCATE( hmxl_n(jpi,jpj,jpk) , ustar s2(jpi,jpj) ,&118 & zwall (jpi,jpj,jpk) , ustar b2(jpi,jpj) ,STAT= zdf_gls_alloc )124 ALLOCATE( hmxl_n(jpi,jpj,jpk) , ustar2_surf(jpi,jpj) , & 125 & zwall (jpi,jpj,jpk) , ustar2_top (jpi,jpj) , ustar2_bot(jpi,jpj) , STAT= zdf_gls_alloc ) 119 126 ! 120 127 IF( lk_mpp ) CALL mpp_sum ( zdf_gls_alloc ) … … 123 130 124 131 125 SUBROUTINE zdf_gls( kt )132 SUBROUTINE zdf_gls( kt, p_sh2, p_avm, p_avt ) 126 133 !!---------------------------------------------------------------------- 127 134 !! *** ROUTINE zdf_gls *** … … 130 137 !! coefficients using the GLS turbulent closure scheme. 131 138 !!---------------------------------------------------------------------- 132 INTEGER, INTENT(in) :: kt ! ocean time step 133 INTEGER :: ji, jj, jk, ibot, ibotm1, dir ! dummy loop arguments 134 REAL(wp) :: zesh2, zsigpsi, zcoef, zex1, zex2 ! local scalars 135 REAL(wp) :: ztx2, zty2, zup, zdown, zcof ! - - 136 REAL(wp) :: zratio, zrn2, zflxb, sh ! - - 139 INTEGER , INTENT(in ) :: kt ! ocean time step 140 REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: p_sh2 ! shear production term 141 REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: p_avm, p_avt ! momentum and tracer Kz (w-points) 142 ! 143 INTEGER :: ji, jj, jk ! dummy loop arguments 144 INTEGER :: ibot, ibotm1 ! local integers 145 INTEGER :: itop, itopp1 ! - - 146 REAL(wp) :: zesh2, zsigpsi, zcoef, zex1 , zex2 ! local scalars 147 REAL(wp) :: ztx2, zty2, zup, zdown, zcof, zdir ! - - 148 REAL(wp) :: zratio, zrn2, zflxb, sh , z_en ! - - 137 149 REAL(wp) :: prod, buoy, diss, zdiss, sm ! - - 138 REAL(wp) :: gh, gm, shr, dif, zsqen, zav ! - - 139 REAL(wp), POINTER, DIMENSION(:,: ) :: zdep 140 REAL(wp), POINTER, DIMENSION(:,: ) :: zkar 141 REAL(wp), POINTER, DIMENSION(:,: ) :: zflxs ! Turbulence fluxed induced by internal waves 142 REAL(wp), POINTER, DIMENSION(:,: ) :: zhsro ! Surface roughness (surface waves) 143 REAL(wp), POINTER, DIMENSION(:,:,:) :: eb ! tke at time before 144 REAL(wp), POINTER, DIMENSION(:,:,:) :: hmxl_b ! mixing length at time before 145 REAL(wp), POINTER, DIMENSION(:,:,:) :: shear ! vertical shear 146 REAL(wp), POINTER, DIMENSION(:,:,:) :: eps ! dissipation rate 147 REAL(wp), POINTER, DIMENSION(:,:,:) :: zwall_psi ! Wall function use in the wb case (ln_sigpsi) 148 REAL(wp), POINTER, DIMENSION(:,:,:) :: psi ! psi at time now 149 REAL(wp), POINTER, DIMENSION(:,:,:) :: z_elem_a ! element of the first matrix diagonal 150 REAL(wp), POINTER, DIMENSION(:,:,:) :: z_elem_b ! element of the second matrix diagonal 151 REAL(wp), POINTER, DIMENSION(:,:,:) :: z_elem_c ! element of the third matrix diagonal 152 REAL(wp), POINTER, DIMENSION(:,:,:) :: zstt, zstm ! stability function on tracer and momentum 150 REAL(wp) :: gh, gm, shr, dif, zsqen, zavt, zavm ! - - 151 REAL(wp), DIMENSION(jpi,jpj) :: zdep 152 REAL(wp), DIMENSION(jpi,jpj) :: zkar 153 REAL(wp), DIMENSION(jpi,jpj) :: zflxs ! Turbulence fluxed induced by internal waves 154 REAL(wp), DIMENSION(jpi,jpj) :: zhsro ! Surface roughness (surface waves) 155 REAL(wp), DIMENSION(jpi,jpj,jpk) :: eb ! tke at time before 156 REAL(wp), DIMENSION(jpi,jpj,jpk) :: hmxl_b ! mixing length at time before 157 REAL(wp), DIMENSION(jpi,jpj,jpk) :: eps ! dissipation rate 158 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zwall_psi ! Wall function use in the wb case (ln_sigpsi) 159 REAL(wp), DIMENSION(jpi,jpj,jpk) :: psi ! psi at time now 160 REAL(wp), DIMENSION(jpi,jpj,jpk) :: z_elem_a ! element of the first matrix diagonal 161 REAL(wp), DIMENSION(jpi,jpj,jpk) :: z_elem_b ! element of the second matrix diagonal 162 REAL(wp), DIMENSION(jpi,jpj,jpk) :: z_elem_c ! element of the third matrix diagonal 163 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zstt, zstm ! stability function on tracer and momentum 153 164 !!-------------------------------------------------------------------- 154 165 ! 155 166 IF( nn_timing == 1 ) CALL timing_start('zdf_gls') 156 167 ! 157 CALL wrk_alloc( jpi,jpj, zdep, zkar, zflxs, zhsro )158 CALL wrk_alloc( jpi,jpj,jpk, eb, hmxl_b, shear, eps, zwall_psi, z_elem_a, z_elem_b, z_elem_c, psi )159 CALL wrk_alloc( jpi,jpj,jpk, zstt, zstm )160 161 168 ! Preliminary computing 162 169 163 ustar s2 = 0._wp ; ustarb2 = 0._wp ; psi = 0._wp ; zwall_psi = 0._wp164 165 ! restore before values computed GLS alone166 avt(:,:,:) = avt_k(:,:,:)167 avm(:,:,:) = avm_k(:,:,:) 168 169 ! Compute surface and bottom friction at T-points170 ustar2_surf(:,:) = 0._wp ; psi(:,:,:) = 0._wp 171 ustar2_top (:,:) = 0._wp ; zwall_psi(:,:,:) = 0._wp 172 ustar2_bot (:,:) = 0._wp 173 174 175 176 ! Compute surface, top and bottom friction at T-points 170 177 DO jj = 2, jpjm1 171 178 DO ji = fs_2, fs_jpim1 ! vector opt. 172 179 ! 173 180 ! surface friction 174 ustar s2(ji,jj) = r1_rau0 * taum(ji,jj) * tmask(ji,jj,1)181 ustar2_surf(ji,jj) = r1_rau0 * taum(ji,jj) * tmask(ji,jj,1) 175 182 ! 183 !!gm old 176 184 ! bottom friction (explicit before friction) 177 185 ! Note that we chose here not to bound the friction as in dynbfr) … … 180 188 zty2 = ( bfrva(ji,jj) * vb(ji,jj,mbkv(ji,jj)) + bfrva(ji,jj-1) * vb(ji,jj-1,mbkv(ji,jj-1)) ) & 181 189 & * ( 1._wp - 0.5_wp * vmask(ji,jj,1) * vmask(ji,jj-1,1) ) 182 ustar b2(ji,jj) = SQRT( ztx2 * ztx2 + zty2 * zty2 ) * tmask(ji,jj,1)190 ustar2_bot(ji,jj) = SQRT( ztx2 * ztx2 + zty2 * zty2 ) * tmask(ji,jj,1) 183 191 END DO 184 192 END DO 185 193 !!gm new 194 !!gm Rq we may add here r_ke0(_top/_bot) ? ==>> think about that... 195 ! ! bottom friction (explicit before friction) 196 ! zmsku = ( 2._wp - umask(ji-1,jj,mbkt(ji,jj)) * umask(ji,jj,mbkt(ji,jj)) ) 197 ! zmskv = ( 2._wp - vmask(ji,jj-1,mbkt(ji,jj)) * vmask(ji,jj,mbkt(ji,jj)) ) ! (CAUTION: CdU<0) 198 ! ustar2_bot(ji,jj) = - rCdU_bot(ji,jj) * SQRT( ( zmsku*( ub(ji,jj,mbkt(ji,jj))+ub(ji-1,jj,mbkt(ji,jj)) ) )**2 199 ! & + ( zmskv*( vb(ji,jj,mbkt(ji,jj))+vb(ji,jj-1,mbkt(ji,jj)) ) )**2 ) 200 ! END DO 201 ! END DO 202 ! IF( ln_isfcav ) THEN !top friction 203 ! DO jj = 2, jpjm1 204 ! DO ji = fs_2, fs_jpim1 ! vector opt. 205 ! zmsku = ( 2. - umask(ji-1,jj,mikt(ji,jj)) * umask(ji,jj,mikt(ji,jj)) ) 206 ! zmskv = ( 2. - vmask(ji,jj-1,mikt(ji,jj)) * vmask(ji,jj,mikt(ji,jj)) ) ! (CAUTION: CdU<0) 207 ! ustar2_top(ji,jj) = - rCdU_top(ji,jj) * SQRT( ( zmsku*( ub(ji,jj,mikt(ji,jj))+ub(ji-1,jj,mikt(ji,jj)) ) )**2 208 ! & + ( zmskv*( vb(ji,jj,mikt(ji,jj))+vb(ji,jj-1,mikt(ji,jj)) ) )**2 ) 209 ! END DO 210 ! END DO 211 ! ENDIF 212 !!gm 186 213 SELECT CASE ( nn_z0_met ) !== Set surface roughness length ==! 187 214 CASE ( 0 ) ! Constant roughness 188 215 zhsro(:,:) = rn_hsro 189 216 CASE ( 1 ) ! Standard Charnock formula 190 zhsro(:,:) = MAX( rsbc_zs1 * ustars2(:,:), rn_hsro)217 zhsro(:,:) = MAX( rsbc_zs1 * ustar2_surf(:,:) , rn_hsro ) 191 218 CASE ( 2 ) ! Roughness formulae according to Rascle et al., Ocean Modelling (2008) 192 219 !!gm zcof = 2._wp * 0.6_wp / 28._wp 193 !!gm zdep(:,:) = 30._wp * TANH( zcof/ SQRT( MAX(ustar s2(:,:),rsmall) ) ) ! Wave age (eq. 10)194 zdep (:,:) = 30.*TANH(2.*0.3/(28.*SQRT(MAX(ustars2(:,:),rsmall))))! Wave age (eq. 10)195 zhsro(:,:) = MAX(rsbc_zs2 * ustar s2(:,:) * zdep(:,:)**1.5, rn_hsro) ! zhsro = rn_frac_hs * Hsw (eq. 11)220 !!gm zdep(:,:) = 30._wp * TANH( zcof/ SQRT( MAX(ustar2_surf(:,:),rsmall) ) ) ! Wave age (eq. 10) 221 zdep (:,:) = 30.*TANH( 2.*0.3/(28.*SQRT(MAX(ustar2_surf(:,:),rsmall))) ) ! Wave age (eq. 10) 222 zhsro(:,:) = MAX(rsbc_zs2 * ustar2_surf(:,:) * zdep(:,:)**1.5, rn_hsro) ! zhsro = rn_frac_hs * Hsw (eq. 11) 196 223 CASE ( 3 ) ! Roughness given by the wave model (coupled or read in file) 197 224 !!gm BUG missing a multiplicative coefficient.... 198 225 zhsro(:,:) = hsw(:,:) 199 226 END SELECT 200 201 ! !== Compute shear and dissipation rate ==! 202 CALL zdf_sh2( shear ) 203 204 227 ! 205 228 DO jk = 2, jpkm1 !== Compute dissipation rate ==! 206 229 DO jj = 1, jpjm1 … … 245 268 DO ji = 2, jpim1 246 269 ! 247 buoy = - avt(ji,jj,jk) * rn2(ji,jj,jk)! stratif. destruction270 buoy = - p_avt(ji,jj,jk) * rn2(ji,jj,jk) ! stratif. destruction 248 271 ! 249 272 diss = eps(ji,jj,jk) ! dissipation 250 273 ! 251 dir = 0.5_wp + SIGN( 0.5_wp, shear(ji,jj,jk) + buoy ) ! dir =1(=0) if shear(ji,jj,jk)+buoy >0(<0) 252 ! 253 zesh2 = dir*(shear(ji,jj,jk)+buoy)+(1._wp-dir)*shear(ji,jj,jk) ! production term 254 zdiss = dir*(diss/en(ji,jj,jk)) +(1._wp-dir)*(diss-buoy)/en(ji,jj,jk) ! dissipation term 274 zdir = 0.5_wp + SIGN( 0.5_wp, p_sh2(ji,jj,jk) + buoy ) ! zdir =1(=0) if shear(ji,jj,jk)+buoy >0(<0) 275 ! 276 zesh2 = zdir*(p_sh2(ji,jj,jk)+buoy)+(1._wp-zdir)*p_sh2(ji,jj,jk) ! production term 277 zdiss = zdir*(diss/en(ji,jj,jk)) +(1._wp-zdir)*(diss-buoy)/en(ji,jj,jk) ! dissipation term 278 !!gm better coding, identical results 279 ! zesh2 = p_sh2(ji,jj,jk) + zdir*buoy ! production term 280 ! zdiss = ( diss - (1._wp-zdir)*buoy ) / en(ji,jj,jk) ! dissipation term 281 !!gm 255 282 ! 256 283 ! Compute a wall function from 1. to rsc_psi*zwall/rsc_psi0 … … 269 296 zcof = rfact_tke * tmask(ji,jj,jk) 270 297 ! ! lower diagonal 271 z_elem_a(ji,jj,jk) = zcof * ( avm(ji,jj,jk ) +avm(ji,jj,jk-1) ) / ( e3t_n(ji,jj,jk-1) * e3w_n(ji,jj,jk) )298 z_elem_a(ji,jj,jk) = zcof * ( p_avm(ji,jj,jk ) + p_avm(ji,jj,jk-1) ) / ( e3t_n(ji,jj,jk-1) * e3w_n(ji,jj,jk) ) 272 299 ! ! upper diagonal 273 z_elem_c(ji,jj,jk) = zcof * ( avm(ji,jj,jk+1) +avm(ji,jj,jk ) ) / ( e3t_n(ji,jj,jk ) * e3w_n(ji,jj,jk) )300 z_elem_c(ji,jj,jk) = zcof * ( p_avm(ji,jj,jk+1) + p_avm(ji,jj,jk ) ) / ( e3t_n(ji,jj,jk ) * e3w_n(ji,jj,jk) ) 274 301 ! ! diagonal 275 302 z_elem_b(ji,jj,jk) = 1._wp - z_elem_a(ji,jj,jk) - z_elem_c(ji,jj,jk) + rdt * zdiss * wmask(ji,jj,jk) … … 293 320 CASE ( 0 ) ! Dirichlet case 294 321 ! First level 295 en(:,:,1) = rc02r * ustar s2(:,:) * (1._wp + rsbc_tke1)**(2._wp/3._wp)322 en(:,:,1) = rc02r * ustar2_surf(:,:) * (1._wp + rsbc_tke1)**r2_3 296 323 en(:,:,1) = MAX(en(:,:,1), rn_emin) 297 324 z_elem_a(:,:,1) = en(:,:,1) … … 300 327 ! 301 328 ! One level below 302 en(:,:,2) = rc02r * ustar s2(:,:) * (1._wp + rsbc_tke1 * ((zhsro(:,:)+gdepw_n(:,:,2)) &329 en(:,:,2) = rc02r * ustar2_surf(:,:) * (1._wp + rsbc_tke1 * ((zhsro(:,:)+gdepw_n(:,:,2)) & 303 330 & / zhsro(:,:) )**(1.5_wp*ra_sf))**(2._wp/3._wp) 304 331 en(:,:,2) = MAX(en(:,:,2), rn_emin ) … … 311 338 ! 312 339 ! Dirichlet conditions at k=1 313 en(:,:,1) = rc02r * ustar s2(:,:) * (1._wp + rsbc_tke1)**(2._wp/3._wp)340 en(:,:,1) = rc02r * ustar2_surf(:,:) * (1._wp + rsbc_tke1)**r2_3 314 341 en(:,:,1) = MAX(en(:,:,1), rn_emin) 315 342 z_elem_a(:,:,1) = en(:,:,1) … … 322 349 z_elem_a(:,:,2) = 0._wp 323 350 zkar(:,:) = (rl_sf + (vkarmn-rl_sf)*(1.-exp(-rtrans*gdept_n(:,:,1)/zhsro(:,:)) )) 324 zflxs(:,:) = rsbc_tke2 * ustar s2(:,:)**1.5_wp * zkar(:,:) &351 zflxs(:,:) = rsbc_tke2 * ustar2_surf(:,:)**1.5_wp * zkar(:,:) & 325 352 & * ((zhsro(:,:)+gdept_n(:,:,1)) / zhsro(:,:) )**(1.5_wp*ra_sf) 326 353 … … 340 367 DO jj = 2, jpjm1 341 368 DO ji = fs_2, fs_jpim1 ! vector opt. 369 !!gm This means that bottom and ocean w-level above have a specified "en" value. Sure ???? 370 !! With thick deep ocean level thickness, this may be quite large, no ??? 371 !! in particular in ocean cavities where top stratification can be large... 342 372 ibot = mbkt(ji,jj) + 1 ! k bottom level of w-point 343 373 ibotm1 = mbkt(ji,jj) ! k-1 bottom level of w-point but >=1 344 374 ! 345 ! Bottom level Dirichlet condition: 346 z_elem_a(ji,jj,ibot ) = 0._wp 347 z_elem_c(ji,jj,ibot ) = 0._wp 348 z_elem_b(ji,jj,ibot ) = 1._wp 349 en(ji,jj,ibot ) = MAX( rc02r * ustarb2(ji,jj), rn_emin ) 350 ! 351 ! Just above last level, Dirichlet condition again 352 z_elem_a(ji,jj,ibotm1) = 0._wp 353 z_elem_c(ji,jj,ibotm1) = 0._wp 354 z_elem_b(ji,jj,ibotm1) = 1._wp 355 en(ji,jj,ibotm1) = MAX( rc02r * ustarb2(ji,jj), rn_emin ) 356 END DO 357 END DO 375 z_en = MAX( rc02r * ustar2_bot(ji,jj), rn_emin ) 376 ! 377 ! Dirichlet condition applied at: 378 ! Bottom level (ibot) & Just above it (ibotm1) 379 z_elem_a(ji,jj,ibot) = 0._wp ; z_elem_a(ji,jj,ibotm1) = 0._wp 380 z_elem_c(ji,jj,ibot) = 0._wp ; z_elem_c(ji,jj,ibotm1) = 0._wp 381 z_elem_b(ji,jj,ibot) = 1._wp ; z_elem_b(ji,jj,ibotm1) = 1._wp 382 en (ji,jj,ibot) = z_en ; en (ji,jj,ibotm1) = z_en 383 END DO 384 END DO 385 !!gm new 386 ! IF( ln_isfcav) THEN ! top boundary (ocean cavity) 387 ! DO jj = 2, jpjm1 388 ! DO ji = fs_2, fs_jpim1 ! vector opt. 389 ! itop = mikt(ji,jj) ! k top w-point 390 ! itopp1 = mikt(ji,jj) + 1 ! k+1 1st w-point below the top one 391 ! ! ! mask at the ocean surface points 392 ! z_en = MAX( rc02r * ustar2_top(ji,jj), rn_emin ) * ( 1._wp - tmask(ji,jj,1) ) 393 ! ! 394 ! ! Dirichlet condition applied at: 395 ! ! top level (itop) & Just below it (itopp1) 396 ! z_elem_a(ji,jj,itop) = 0._wp ; z_elem_a(ji,jj,ipotm1) = 0._wp 397 ! z_elem_c(ji,jj,itop) = 0._wp ; z_elem_c(ji,jj,itopp1) = 0._wp 398 ! z_elem_b(ji,jj,itop) = 1._wp ; z_elem_b(ji,jj,itopp1) = 1._wp 399 ! en (ji,jj,itop) = zen ; en (ji,jj,itopp1) = z_en 400 ! END DO 401 ! END DO 402 ! ENDIF 403 !!gm 358 404 ! 359 405 CASE ( 1 ) ! Neumman boundary condition … … 364 410 ibotm1 = mbkt(ji,jj) ! k-1 bottom level of w-point but >=1 365 411 ! 412 z_en = MAX( rc02r * ustar2_bot(ji,jj), rn_emin ) 413 ! 366 414 ! Bottom level Dirichlet condition: 367 z_elem_a(ji,jj,ibot) = 0._wp 368 z_elem_c(ji,jj,ibot) = 0._wp 369 z_elem_b(ji,jj,ibot) = 1._wp 370 en(ji,jj,ibot) = MAX( rc02r * ustarb2(ji,jj), rn_emin ) 371 ! 372 ! Just above last level: Neumann condition 373 z_elem_b(ji,jj,ibotm1) = z_elem_b(ji,jj,ibotm1) + z_elem_c(ji,jj,ibotm1) ! Remove z_elem_c from z_elem_b 374 z_elem_c(ji,jj,ibotm1) = 0._wp 375 END DO 376 END DO 415 ! Bottom level (ibot) & Just above it (ibotm1) 416 ! Dirichlet ! Neumann 417 z_elem_a(ji,jj,ibot) = 0._wp ! ! Remove z_elem_c from z_elem_b 418 z_elem_b(ji,jj,ibot) = 1._wp ; z_elem_b(ji,jj,ibotm1) = z_elem_b(ji,jj,ibotm1) + z_elem_c(ji,jj,ibotm1) 419 z_elem_c(ji,jj,ibot) = 0._wp ; z_elem_c(ji,jj,ibotm1) = 0._wp 420 END DO 421 END DO 422 !!gm new 423 ! IF( ln_isfcav) THEN ! top boundary (ocean cavity) 424 ! DO jj = 2, jpjm1 425 ! DO ji = fs_2, fs_jpim1 ! vector opt. 426 ! itop = mikt(ji,jj) ! k top w-point 427 ! itopp1 = mikt(ji,jj) + 1 ! k+1 1st w-point below the top one 428 ! ! ! mask at the ocean surface points 429 ! z_en = MAX( rc02r * ustar2_top(ji,jj), rn_emin ) * ( 1._wp - tmask(ji,jj,1) ) 430 ! ! 431 ! ! Bottom level Dirichlet condition: 432 ! ! Bottom level (ibot) & Just above it (ibotm1) 433 ! ! Dirichlet ! Neumann 434 ! z_elem_a(ji,jj,itop) = 0._wp ! ! Remove z_elem_c from z_elem_b 435 ! z_elem_b(ji,jj,itop) = 1._wp ; z_elem_b(ji,jj,itopp1) = z_elem_b(ji,jj,itopp1) + z_elem_c(ji,jj,itopp1) 436 ! z_elem_c(ji,jj,itop) = 0._wp ; z_elem_c(ji,jj,itopp1) = 0._wp 437 ! END DO 438 ! END DO 439 ! ENDIF 440 !!gm 377 441 ! 378 442 END SELECT … … 465 529 zratio = psi(ji,jj,jk) / eb(ji,jj,jk) 466 530 ! 467 ! psi3+ : stable : B=-KhN²<0 => N²>0 if rn2>0 dir = 1 (stable) otherwisedir = 0 (unstable)468 dir = 0.5_wp + SIGN( 0.5_wp, rn2(ji,jj,jk) )469 ! 470 rpsi3 = dir * rpsi3m + ( 1._wp -dir ) * rpsi3p531 ! psi3+ : stable : B=-KhN²<0 => N²>0 if rn2>0 zdir = 1 (stable) otherwise zdir = 0 (unstable) 532 zdir = 0.5_wp + SIGN( 0.5_wp, rn2(ji,jj,jk) ) 533 ! 534 rpsi3 = zdir * rpsi3m + ( 1._wp - zdir ) * rpsi3p 471 535 ! 472 536 ! shear prod. - stratif. destruction 473 prod = rpsi1 * zratio * shear(ji,jj,jk)537 prod = rpsi1 * zratio * p_sh2(ji,jj,jk) 474 538 ! 475 539 ! stratif. destruction 476 buoy = rpsi3 * zratio * (- avt(ji,jj,jk) * rn2(ji,jj,jk) )540 buoy = rpsi3 * zratio * (- p_avt(ji,jj,jk) * rn2(ji,jj,jk) ) 477 541 ! 478 542 ! shear prod. - stratif. destruction 479 543 diss = rpsi2 * zratio * zwall(ji,jj,jk) * eps(ji,jj,jk) 480 544 ! 481 dir = 0.5_wp + SIGN( 0.5_wp, prod + buoy ) !dir =1(=0) if shear(ji,jj,jk)+buoy >0(<0)482 ! 483 zesh2 = dir * ( prod + buoy ) + (1._wp -dir ) * prod ! production term484 zdiss = dir * ( diss / psi(ji,jj,jk) ) + (1._wp -dir ) * (diss-buoy) / psi(ji,jj,jk) ! dissipation term545 zdir = 0.5_wp + SIGN( 0.5_wp, prod + buoy ) ! zdir =1(=0) if shear(ji,jj,jk)+buoy >0(<0) 546 ! 547 zesh2 = zdir * ( prod + buoy ) + (1._wp - zdir ) * prod ! production term 548 zdiss = zdir * ( diss / psi(ji,jj,jk) ) + (1._wp - zdir ) * (diss-buoy) / psi(ji,jj,jk) ! dissipation term 485 549 ! 486 550 ! building the matrix 487 551 zcof = rfact_psi * zwall_psi(ji,jj,jk) * tmask(ji,jj,jk) 488 552 ! ! lower diagonal 489 z_elem_a(ji,jj,jk) = zcof * ( avm(ji,jj,jk ) +avm(ji,jj,jk-1) ) / ( e3t_n(ji,jj,jk-1) * e3w_n(ji,jj,jk) )553 z_elem_a(ji,jj,jk) = zcof * ( p_avm(ji,jj,jk ) + p_avm(ji,jj,jk-1) ) / ( e3t_n(ji,jj,jk-1) * e3w_n(ji,jj,jk) ) 490 554 ! ! upper diagonal 491 z_elem_c(ji,jj,jk) = zcof * ( avm(ji,jj,jk+1) +avm(ji,jj,jk ) ) / ( e3t_n(ji,jj,jk ) * e3w_n(ji,jj,jk) )555 z_elem_c(ji,jj,jk) = zcof * ( p_avm(ji,jj,jk+1) + p_avm(ji,jj,jk ) ) / ( e3t_n(ji,jj,jk ) * e3w_n(ji,jj,jk) ) 492 556 ! ! diagonal 493 557 z_elem_b(ji,jj,jk) = 1._wp - z_elem_a(ji,jj,jk) - z_elem_c(ji,jj,jk) + rdt * zdiss * wmask(ji,jj,jk) … … 506 570 ! 507 571 CASE ( 0 ) ! Dirichlet boundary conditions 508 ! 509 ! Surface value 510 zdep(:,:) = zhsro(:,:) * rl_sf ! Cosmetic 511 psi (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 512 z_elem_a(:,:,1) = psi(:,:,1) 513 z_elem_c(:,:,1) = 0._wp 514 z_elem_b(:,:,1) = 1._wp 515 ! 516 ! One level below 517 zkar(:,:) = (rl_sf + (vkarmn-rl_sf)*(1._wp-exp(-rtrans*gdepw_n(:,:,2)/zhsro(:,:) ))) 518 zdep(:,:) = (zhsro(:,:) + gdepw_n(:,:,2)) * zkar(:,:) 519 psi (:,:,2) = rc0**rpp * en(:,:,2)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 520 z_elem_a(:,:,2) = 0._wp 521 z_elem_c(:,:,2) = 0._wp 522 z_elem_b(:,:,2) = 1._wp 523 ! 524 ! 572 ! 573 ! Surface value 574 zdep (:,:) = zhsro(:,:) * rl_sf ! Cosmetic 575 psi (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 576 z_elem_a(:,:,1) = psi(:,:,1) 577 z_elem_c(:,:,1) = 0._wp 578 z_elem_b(:,:,1) = 1._wp 579 ! 580 ! One level below 581 zkar (:,:) = (rl_sf + (vkarmn-rl_sf)*(1._wp-exp(-rtrans*gdepw_n(:,:,2)/zhsro(:,:) ))) 582 zdep (:,:) = (zhsro(:,:) + gdepw_n(:,:,2)) * zkar(:,:) 583 psi (:,:,2) = rc0**rpp * en(:,:,2)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 584 z_elem_a(:,:,2) = 0._wp 585 z_elem_c(:,:,2) = 0._wp 586 z_elem_b(:,:,2) = 1._wp 587 ! 525 588 CASE ( 1 ) ! Neumann boundary condition on d(psi)/dz 526 ! 527 ! Surface value: Dirichlet 528 zdep(:,:) = zhsro(:,:) * rl_sf 529 psi (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 530 z_elem_a(:,:,1) = psi(:,:,1) 531 z_elem_c(:,:,1) = 0._wp 532 z_elem_b(:,:,1) = 1._wp 533 ! 534 ! Neumann condition at k=2 535 z_elem_b(:,:,2) = z_elem_b(:,:,2) + z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b 536 z_elem_a(:,:,2) = 0._wp 537 ! 538 ! Set psi vertical flux at the surface: 539 zkar(:,:) = rl_sf + (vkarmn-rl_sf)*(1._wp-exp(-rtrans*gdept_n(:,:,1)/zhsro(:,:) )) ! Lengh scale slope 540 zdep(:,:) = ((zhsro(:,:) + gdept_n(:,:,1)) / zhsro(:,:))**(rmm*ra_sf) 541 zflxs(:,:) = (rnn + rsbc_tke1 * (rnn + rmm*ra_sf) * zdep(:,:))*(1._wp + rsbc_tke1*zdep(:,:))**(2._wp*rmm/3._wp-1_wp) 542 zdep(:,:) = rsbc_psi1 * (zwall_psi(:,:,1)*avm(:,:,1)+zwall_psi(:,:,2)*avm(:,:,2)) * & 543 & ustars2(:,:)**rmm * zkar(:,:)**rnn * (zhsro(:,:) + gdept_n(:,:,1))**(rnn-1.) 544 zflxs(:,:) = zdep(:,:) * zflxs(:,:) 545 psi(:,:,2) = psi(:,:,2) + zflxs(:,:) / e3w_n(:,:,2) 546 ! 547 ! 589 ! 590 ! Surface value: Dirichlet 591 zdep (:,:) = zhsro(:,:) * rl_sf 592 psi (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 593 z_elem_a(:,:,1) = psi(:,:,1) 594 z_elem_c(:,:,1) = 0._wp 595 z_elem_b(:,:,1) = 1._wp 596 ! 597 ! Neumann condition at k=2 598 z_elem_b(:,:,2) = z_elem_b(:,:,2) + z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b 599 z_elem_a(:,:,2) = 0._wp 600 ! 601 ! Set psi vertical flux at the surface: 602 zkar (:,:) = rl_sf + (vkarmn-rl_sf)*(1._wp-exp(-rtrans*gdept_n(:,:,1)/zhsro(:,:) )) ! Lengh scale slope 603 zdep (:,:) = ((zhsro(:,:) + gdept_n(:,:,1)) / zhsro(:,:))**(rmm*ra_sf) 604 zflxs(:,:) = (rnn + rsbc_tke1 * (rnn + rmm*ra_sf) * zdep(:,:))*(1._wp + rsbc_tke1*zdep(:,:))**(2._wp*rmm/3._wp-1_wp) 605 zdep (:,:) = rsbc_psi1 * (zwall_psi(:,:,1)*avm(:,:,1)+zwall_psi(:,:,2)*avm(:,:,2)) * & 606 & ustar2_surf(:,:)**rmm * zkar(:,:)**rnn * (zhsro(:,:) + gdept_n(:,:,1))**(rnn-1.) 607 zflxs(:,:) = zdep(:,:) * zflxs(:,:) 608 psi (:,:,2) = psi(:,:,2) + zflxs(:,:) / e3w_n(:,:,2) 609 ! 548 610 END SELECT 549 611 … … 552 614 ! 553 615 SELECT CASE ( nn_bc_bot ) 554 !555 616 ! 556 617 CASE ( 0 ) ! Dirichlet … … 597 658 ! Set psi vertical flux at the bottom: 598 659 zdep(ji,jj) = rn_bfrz0 + 0.5_wp*e3t_n(ji,jj,ibotm1) 599 zflxb = rsbc_psi2 * ( avm(ji,jj,ibot) +avm(ji,jj,ibotm1) ) &660 zflxb = rsbc_psi2 * ( p_avm(ji,jj,ibot) + p_avm(ji,jj,ibotm1) ) & 600 661 & * (0.5_wp*(en(ji,jj,ibot)+en(ji,jj,ibotm1)))**rmm * zdep(ji,jj)**(rnn-1._wp) 601 662 psi(ji,jj,ibotm1) = psi(ji,jj,ibotm1) + zflxb / e3w_n(ji,jj,ibotm1) … … 730 791 gh = gh * rf6 731 792 ! Gm = M²l²/q² Shear number 732 shr = shear(ji,jj,jk) / MAX(avm(ji,jj,jk), rsmall )793 shr = p_sh2(ji,jj,jk) / MAX( p_avm(ji,jj,jk), rsmall ) 733 794 gm = MAX( shr * zcof , 1.e-10 ) 734 795 gm = gm * rf6 … … 765 826 DO jj = 2, jpjm1 766 827 DO ji = fs_2, fs_jpim1 ! vector opt. 767 zsqen 768 zav 769 avt(ji,jj,jk) = MAX( zav, avtb(jk) ) * wmask(ji,jj,jk) ! apply mask for zdfmxl routine770 zav = zsqen * zstm(ji,jj,jk)771 avm(ji,jj,jk) = MAX( zav, avmb(jk) ) ! Note that avm is not masked at the surface and the bottom828 zsqen = SQRT( 2._wp * en(ji,jj,jk) ) * hmxl_n(ji,jj,jk) 829 zavt = zsqen * zstt(ji,jj,jk) 830 zavm = zsqen * zstm(ji,jj,jk) 831 p_avt(ji,jj,jk) = MAX( zavt, avtb(jk) ) * wmask(ji,jj,jk) ! apply mask for zdfmxl routine 832 p_avm(ji,jj,jk) = MAX( zavm, avmb(jk) ) ! Note that avm is not masked at the surface and the bottom 772 833 END DO 773 834 END DO 774 835 END DO 775 836 avt(:,:,1) = 0._wp 776 !!gm I'm sure this is needed to compute the shear term at next time-step777 CALL lbc_lnk( avm, 'W', 1. ) ; CALL lbc_lnk( avt, 'W', 1. )778 837 ! 779 838 IF(ln_ctl) THEN … … 781 840 CALL prt_ctl( tab3d_1=avm, clinfo1=' gls - m: ', ovlap=1, kdim=jpk ) 782 841 ENDIF 783 !784 avt_k(:,:,:) = avt(:,:,:) !== store avt, avm values computed by GLS only ==!785 avm_k(:,:,:) = avm(:,:,:)786 !787 IF( lrst_oce ) CALL gls_rst( kt, 'WRITE' ) ! write the GLS info in the restart file788 !789 CALL wrk_dealloc( jpi,jpj, zdep, zkar, zflxs, zhsro )790 CALL wrk_dealloc( jpi,jpj,jpk, eb, hmxl_b, shear, eps, zwall_psi, z_elem_a, z_elem_b, z_elem_c, psi )791 CALL wrk_dealloc( jpi,jpj,jpk, zstt, zstm )792 842 ! 793 843 IF( nn_timing == 1 ) CALL timing_stop('zdf_gls') … … 837 887 IF(lwp) THEN !* Control print 838 888 WRITE(numout,*) 839 WRITE(numout,*) 'zdf_gls_init : glsturbulent closure scheme'889 WRITE(numout,*) 'zdf_gls_init : GLS turbulent closure scheme' 840 890 WRITE(numout,*) '~~~~~~~~~~~~' 841 891 WRITE(numout,*) ' Namelist namzdf_gls : set gls mixing parameters' … … 854 904 WRITE(numout,*) ' Type of closure nn_clos = ', nn_clos 855 905 WRITE(numout,*) ' Surface roughness (m) rn_hsro = ', rn_hsro 906 !!gm old 907 WRITE(numout,*) ' top roughness (m) (nambfr namelist) rn_tfrz0 = ', rn_tfrz0 856 908 WRITE(numout,*) ' Bottom roughness (m) (nambfr namelist) rn_bfrz0 = ', rn_bfrz0 909 !!gm new 910 ! WRITE(numout,*) ' Namelist namdrg_top/_bot used values:' 911 ! WRITE(numout,*) ' top ocean cavity roughness (m) rn_z0(_top) = ', r_z0_top 912 ! WRITE(numout,*) ' Bottom seafloor roughness (m) rn_z0(_bot) = ', r_z0_bot 913 !!gm 914 WRITE(numout,*) 857 915 ENDIF 858 916 … … 861 919 862 920 ! !* Check of some namelist values 863 IF( nn_bc_surf < 0 .OR. nn_bc_surf > 1 )CALL ctl_stop( 'zdf_gls_init: bad flag: nn_bc_surf is 0 or 1' )864 IF( nn_bc_surf < 0 .OR. nn_bc_surf > 1 )CALL ctl_stop( 'zdf_gls_init: bad flag: nn_bc_surf is 0 or 1' )865 IF( nn_z0_met < 0 .OR. nn_z0_met > 3 )CALL ctl_stop( 'zdf_gls_init: bad flag: nn_z0_met is 0, 1, 2 or 3' )866 IF( nn_z0_met == 3 .AND. .NOT.ln_wave )CALL ctl_stop( 'zdf_gls_init: nn_z0_met=3 requires ln_wave=T' )867 IF( nn_stab_func < 0 .OR. nn_stab_func > 3 )CALL ctl_stop( 'zdf_gls_init: bad flag: nn_stab_func is 0, 1, 2 and 3' )868 IF( nn_clos < 0 .OR. nn_clos > 3 )CALL ctl_stop( 'zdf_gls_init: bad flag: nn_clos is 0, 1, 2 or 3' )921 IF( nn_bc_surf < 0 .OR. nn_bc_surf > 1 ) CALL ctl_stop( 'zdf_gls_init: bad flag: nn_bc_surf is 0 or 1' ) 922 IF( nn_bc_surf < 0 .OR. nn_bc_surf > 1 ) CALL ctl_stop( 'zdf_gls_init: bad flag: nn_bc_surf is 0 or 1' ) 923 IF( nn_z0_met < 0 .OR. nn_z0_met > 3 ) CALL ctl_stop( 'zdf_gls_init: bad flag: nn_z0_met is 0, 1, 2 or 3' ) 924 IF( nn_z0_met == 3 .AND. .NOT.ln_wave ) CALL ctl_stop( 'zdf_gls_init: nn_z0_met=3 requires ln_wave=T' ) 925 IF( nn_stab_func < 0 .OR. nn_stab_func > 3 ) CALL ctl_stop( 'zdf_gls_init: bad flag: nn_stab_func is 0, 1, 2 and 3' ) 926 IF( nn_clos < 0 .OR. nn_clos > 3 ) CALL ctl_stop( 'zdf_gls_init: bad flag: nn_clos is 0, 1, 2 or 3' ) 869 927 870 928 SELECT CASE ( nn_clos ) !* set the parameters for the chosen closure … … 872 930 CASE( 0 ) ! k-kl (Mellor-Yamada) 873 931 ! 874 IF(lwp) WRITE(numout,*) 'The choosen closure is k-kl closed to the classical Mellor-Yamada' 932 IF(lwp) WRITE(numout,*) ' ==>> k-kl closure chosen (i.e. closed to the classical Mellor-Yamada)' 933 IF(lwp) WRITE(numout,*) 875 934 rpp = 0._wp 876 935 rmm = 1._wp … … 890 949 CASE( 1 ) ! k-eps 891 950 ! 892 IF(lwp) WRITE(numout,*) 'The choosen closure is k-eps' 951 IF(lwp) WRITE(numout,*) ' ==>> k-eps closure chosen' 952 IF(lwp) WRITE(numout,*) 893 953 rpp = 3._wp 894 954 rmm = 1.5_wp … … 908 968 CASE( 2 ) ! k-omega 909 969 ! 910 IF(lwp) WRITE(numout,*) 'The choosen closure is k-omega' 970 IF(lwp) WRITE(numout,*) ' ==>> k-omega closure chosen' 971 IF(lwp) WRITE(numout,*) 911 972 rpp = -1._wp 912 973 rmm = 0.5_wp … … 926 987 CASE( 3 ) ! generic 927 988 ! 928 IF(lwp) WRITE(numout,*) 'The choosen closure is generic' 989 IF(lwp) WRITE(numout,*) ' ==>> generic closure chosen' 990 IF(lwp) WRITE(numout,*) 929 991 rpp = 2._wp 930 992 rmm = 1._wp … … 949 1011 CASE ( 0 ) ! Galperin stability functions 950 1012 ! 951 IF(lwp) WRITE(numout,*) ' Stability functions from Galperin'1013 IF(lwp) WRITE(numout,*) ' ==>> Stability functions from Galperin' 952 1014 rc2 = 0._wp 953 1015 rc3 = 0._wp … … 961 1023 CASE ( 1 ) ! Kantha-Clayson stability functions 962 1024 ! 963 IF(lwp) WRITE(numout,*) ' Stability functions from Kantha-Clayson'1025 IF(lwp) WRITE(numout,*) ' ==>> Stability functions from Kantha-Clayson' 964 1026 rc2 = 0.7_wp 965 1027 rc3 = 0.2_wp … … 973 1035 CASE ( 2 ) ! Canuto A stability functions 974 1036 ! 975 IF(lwp) WRITE(numout,*) ' Stability functions from Canuto A'1037 IF(lwp) WRITE(numout,*) ' ==>> Stability functions from Canuto A' 976 1038 rs0 = 1.5_wp * rl1 * rl5*rl5 977 1039 rs1 = -rl4*(rl6+rl7) + 2._wp*rl4*rl5*(rl1-(1._wp/3._wp)*rl2-rl3) + 1.5_wp*rl1*rl5*rl8 … … 997 1059 CASE ( 3 ) ! Canuto B stability functions 998 1060 ! 999 IF(lwp) WRITE(numout,*) ' Stability functions from Canuto B'1061 IF(lwp) WRITE(numout,*) ' ==>> Stability functions from Canuto B' 1000 1062 rs0 = 1.5_wp * rm1 * rm5*rm5 1001 1063 rs1 = -rm4 * (rm6+rm7) + 2._wp * rm4*rm5*(rm1-(1._wp/3._wp)*rm2-rm3) + 1.5_wp * rm1*rm5*rm8 … … 1052 1114 IF(lwp) THEN !* Control print 1053 1115 WRITE(numout,*) 1054 WRITE(numout,*) 'Limit values' 1055 WRITE(numout,*) '~~~~~~~~~~~~' 1056 WRITE(numout,*) 'Parameter m = ',rmm 1057 WRITE(numout,*) 'Parameter n = ',rnn 1058 WRITE(numout,*) 'Parameter p = ',rpp 1059 WRITE(numout,*) 'rpsi1 = ',rpsi1 1060 WRITE(numout,*) 'rpsi2 = ',rpsi2 1061 WRITE(numout,*) 'rpsi3m = ',rpsi3m 1062 WRITE(numout,*) 'rpsi3p = ',rpsi3p 1063 WRITE(numout,*) 'rsc_tke = ',rsc_tke 1064 WRITE(numout,*) 'rsc_psi = ',rsc_psi 1065 WRITE(numout,*) 'rsc_psi0 = ',rsc_psi0 1066 WRITE(numout,*) 'rc0 = ',rc0 1116 WRITE(numout,*) ' Limit values :' 1117 WRITE(numout,*) ' Parameter m = ', rmm 1118 WRITE(numout,*) ' Parameter n = ', rnn 1119 WRITE(numout,*) ' Parameter p = ', rpp 1120 WRITE(numout,*) ' rpsi1 = ', rpsi1 1121 WRITE(numout,*) ' rpsi2 = ', rpsi2 1122 WRITE(numout,*) ' rpsi3m = ', rpsi3m 1123 WRITE(numout,*) ' rpsi3p = ', rpsi3p 1124 WRITE(numout,*) ' rsc_tke = ', rsc_tke 1125 WRITE(numout,*) ' rsc_psi = ', rsc_psi 1126 WRITE(numout,*) ' rsc_psi0 = ', rsc_psi0 1127 WRITE(numout,*) ' rc0 = ', rc0 1067 1128 WRITE(numout,*) 1068 WRITE(numout,*) 'Shear free turbulence parameters:' 1069 WRITE(numout,*) 'rcm_sf = ',rcm_sf 1070 WRITE(numout,*) 'ra_sf = ',ra_sf 1071 WRITE(numout,*) 'rl_sf = ',rl_sf 1072 WRITE(numout,*) 1129 WRITE(numout,*) ' Shear free turbulence parameters:' 1130 WRITE(numout,*) ' rcm_sf = ', rcm_sf 1131 WRITE(numout,*) ' ra_sf = ', ra_sf 1132 WRITE(numout,*) ' rl_sf = ', rl_sf 1073 1133 ENDIF 1074 1134 … … 1090 1150 ! 1091 1151 ! !* Wall proximity function 1092 zwall (:,:,:) = 1._wp * tmask(:,:,:) 1152 !!gm tmask or wmask ???? 1153 zwall(:,:,:) = 1._wp * tmask(:,:,:) 1093 1154 1094 1155 ! !* read or initialize all required files … … 1133 1194 CALL iom_get( numror, jpdom_autoglo, 'hmxl_n', hmxl_n ) 1134 1195 ELSE 1135 IF(lwp) WRITE(numout,*) ' ===>>>> : previous run without gls scheme, en and hmxl_n computed by iterative loop' 1196 IF(lwp) WRITE(numout,*) 1197 IF(lwp) WRITE(numout,*) ' ==>> previous run without GLS scheme, set en and hmxl_n to background values' 1136 1198 en (:,:,:) = rn_emin 1137 1199 hmxl_n(:,:,:) = 0.05_wp 1138 avt_k (:,:,:) = avt(:,:,:) 1139 avm_k (:,:,:) = avm(:,:,:) 1140 DO jit = nit000 + 1, nit000 + 10 ; CALL zdf_gls( jit ) ; END DO 1200 ! avt_k, avm_k already set to the background value in zdf_phy_init 1141 1201 ENDIF 1142 1202 ELSE !* Start from rest 1143 IF(lwp) WRITE(numout,*) ' ===>>>> : Initialisation of en and hmxl_n by background values' 1203 IF(lwp) WRITE(numout,*) 1204 IF(lwp) WRITE(numout,*) ' ==>> start from rest, set en and hmxl_n by background values' 1144 1205 en (:,:,:) = rn_emin 1145 1206 hmxl_n(:,:,:) = 0.05_wp 1146 DO jk = 1, jpk 1147 avt_k(:,:,jk) = avtb(jk) * wmask(:,:,jk) 1148 avm_k(:,:,jk) = avmb(jk) * wmask(:,:,jk) 1149 END DO 1207 ! avt_k, avm_k already set to the background value in zdf_phy_init 1150 1208 ENDIF 1151 1209 ! -
branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfiwm.F90
r7990 r8093 35 35 PUBLIC zdf_iwm ! called in step module 36 36 PUBLIC zdf_iwm_init ! called in nemogcm module 37 PUBLIC zdf_iwm_alloc ! called in nemogcm module38 37 39 38 ! !!* Namelist namzdf_iwm : internal wave-driven mixing * … … 78 77 79 78 80 SUBROUTINE zdf_iwm( kt )79 SUBROUTINE zdf_iwm( kt, p_avm, p_avt, p_avs ) 81 80 !!---------------------------------------------------------------------- 82 81 !! *** ROUTINE zdf_iwm *** … … 127 126 !! References : de Lavergne et al. 2015, JPO; 2016, in prep. 128 127 !!---------------------------------------------------------------------- 129 INTEGER, INTENT(in) :: kt ! ocean time-step 128 INTEGER , INTENT(in ) :: kt ! ocean time step 129 REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: p_avm ! momentum Kz (w-points) 130 REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: p_avt, p_avs ! tracer Kz (w-points) 130 131 ! 131 132 INTEGER :: ji, jj, jk ! dummy loop indices 132 REAL(wp) :: z tpc! scalar workspace133 REAL(wp), DIMENSION( :,:) , POINTER:: zfact ! Used for vertical structure134 REAL(wp), DIMENSION( :,:) , POINTER:: zhdep ! Ocean depth135 REAL(wp), DIMENSION( :,:,:), POINTER:: zwkb ! WKB-stretched height above bottom136 REAL(wp), DIMENSION( :,:,:), POINTER:: zweight ! Weight for high mode vertical distribution137 REAL(wp), DIMENSION( :,:,:), POINTER:: znu_t ! Molecular kinematic viscosity (T grid)138 REAL(wp), DIMENSION( :,:,:), POINTER:: znu_w ! Molecular kinematic viscosity (W grid)139 REAL(wp), DIMENSION( :,:,:), POINTER:: zReb ! Turbulence intensity parameter133 REAL(wp) :: zztmp ! scalar workspace 134 REAL(wp), DIMENSION(jpi,jpj) :: zfact ! Used for vertical structure 135 REAL(wp), DIMENSION(jpi,jpj) :: zhdep ! Ocean depth 136 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zwkb ! WKB-stretched height above bottom 137 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zweight ! Weight for high mode vertical distribution 138 REAL(wp), DIMENSION(jpi,jpj,jpk) :: znu_t ! Molecular kinematic viscosity (T grid) 139 REAL(wp), DIMENSION(jpi,jpj,jpk) :: znu_w ! Molecular kinematic viscosity (W grid) 140 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zReb ! Turbulence intensity parameter 140 141 !!---------------------------------------------------------------------- 141 142 ! 142 143 IF( nn_timing == 1 ) CALL timing_start('zdf_iwm') 143 144 ! 144 CALL wrk_alloc( jpi,jpj, zfact, zhdep ) 145 CALL wrk_alloc( jpi,jpj,jpk, zwkb, zweight, znu_t, znu_w, zReb ) 146 147 ! ! ----------------------------- ! 148 ! ! Internal wave-driven mixing ! (compute zav_wave) 149 ! ! ----------------------------- ! 145 ! ! ----------------------------- ! 146 ! ! Internal wave-driven mixing ! (compute zav_wave) 147 ! ! ----------------------------- ! 150 148 ! 151 149 ! !* Critical slope mixing: distribute energy over the time-varying ocean depth, … … 162 160 emix_iwm(:,:,jk) = zfact(:,:) * ( EXP( ( gde3w_n(:,:,jk ) - zhdep(:,:) ) / hcri_iwm(:,:) ) & 163 161 & - EXP( ( gde3w_n(:,:,jk-1) - zhdep(:,:) ) / hcri_iwm(:,:) ) ) * wmask(:,:,jk) & 162 & / ( gde3w_n(:,:,jk) - gde3w_n(:,:,jk-1) ) 163 164 164 !!gm delta(gde3w_n) = e3t_n !! Please verify the grid-point position w versus t-point 165 & / ( gde3w_n(:,:,jk) - gde3w_n(:,:,jk-1) ) 165 !!gm it seems to me that only 1/hcri_iwm is used ==> compute it one for all 166 166 167 END DO 167 168 168 169 ! !* Pycnocline-intensified mixing: distribute energy over the time-varying 169 170 ! !* ocean depth as proportional to sqrt(rn2)^nn_zpyc 170 ! 171 ! ! (NB: N2 is masked, so no use of wmask here) 171 172 SELECT CASE ( nn_zpyc ) 172 173 ! … … 210 211 ! !* ocean depth as proportional to rn2 * exp(-z_wkb/rn_hbot) 211 212 ! 212 zwkb (:,:,:) = 0._wp213 zfact(:,:) = 0._wp213 zwkb (:,:,:) = 0._wp 214 zfact(:,:) = 0._wp 214 215 DO jk = 2, jpkm1 215 216 zfact(:,:) = zfact(:,:) + e3w_n(:,:,jk) * SQRT( MAX( 0._wp, rn2(:,:,jk) ) ) * wmask(:,:,jk) 216 217 zwkb(:,:,jk) = zfact(:,:) 217 218 END DO 219 !!gm even better: 220 ! DO jk = 2, jpkm1 221 ! zwkb(:,:) = zwkb(:,:) + e3w_n(:,:,jk) * SQRT( MAX( 0._wp, rn2(:,:,jk) ) ) 222 ! END DO 223 ! zfact(:,:) = zwkb(:,:,jpkm1) 224 !!gm or just use zwkb(k=jpk-1) instead of zfact... 225 !!gm 218 226 ! 219 227 DO jk = 2, jpkm1 … … 221 229 DO ji = 1, jpi 222 230 IF( zfact(ji,jj) /= 0 ) zwkb(ji,jj,jk) = zhdep(ji,jj) * ( zfact(ji,jj) - zwkb(ji,jj,jk) ) & 223 & * tmask(ji,jj,jk) / zfact(ji,jj)224 END DO 225 END DO 226 END DO 227 zwkb(:,:,1) = zhdep(:,:) * tmask(:,:,1)231 & * wmask(ji,jj,jk) / zfact(ji,jj) 232 END DO 233 END DO 234 END DO 235 zwkb(:,:,1) = zhdep(:,:) * wmask(:,:,1) 228 236 ! 229 237 zweight(:,:,:) = 0._wp … … 247 255 emix_iwm(:,:,jk) = emix_iwm(:,:,jk) + zweight(:,:,jk) * zfact(:,:) * wmask(:,:,jk) & 248 256 & / ( gde3w_n(:,:,jk) - gde3w_n(:,:,jk-1) ) 249 END DO 250 ! 257 !!gm use of e3t_n just above? 258 END DO 259 ! 260 !!gm this is to be replaced by just a constant value znu=1.e-6 m2/s 251 261 ! Calculate molecular kinematic viscosity 252 262 znu_t(:,:,:) = 1.e-4_wp * ( 17.91_wp - 0.53810_wp * tsn(:,:,:,jp_tem) + 0.00694_wp * tsn(:,:,:,jp_tem) * tsn(:,:,:,jp_tem) & … … 255 265 znu_w(:,:,jk) = 0.5_wp * ( znu_t(:,:,jk-1) + znu_t(:,:,jk) ) * wmask(:,:,jk) 256 266 END DO 267 !!gm end 257 268 ! 258 269 ! Calculate turbulence intensity parameter Reb … … 285 296 ! 286 297 IF( kt == nit000 ) THEN !* Control print at first time-step: diagnose the energy consumed by zav_wave 287 z tpc= 0._wp298 zztmp = 0._wp 288 299 !!gm used of glosum 3D.... 289 300 DO jk = 2, jpkm1 290 301 DO jj = 1, jpj 291 302 DO ji = 1, jpi 292 z tpc = ztpc+ e3w_n(ji,jj,jk) * e1e2t(ji,jj) &293 & * MAX( 0._wp, rn2(ji,jj,jk) ) * zav_wave(ji,jj,jk) * wmask(ji,jj,jk) * tmask_i(ji,jj)303 zztmp = zztmp + e3w_n(ji,jj,jk) * e1e2t(ji,jj) & 304 & * MAX( 0._wp, rn2(ji,jj,jk) ) * zav_wave(ji,jj,jk) * wmask(ji,jj,jk) * tmask_i(ji,jj) 294 305 END DO 295 306 END DO 296 307 END DO 297 IF( lk_mpp ) CALL mpp_sum( z tpc)298 z tpc = rau0 * ztpc! Global integral of rauo * Kz * N^2 = power contributing to mixing308 IF( lk_mpp ) CALL mpp_sum( zztmp ) 309 zztmp = rau0 * zztmp ! Global integral of rauo * Kz * N^2 = power contributing to mixing 299 310 ! 300 311 IF(lwp) THEN … … 303 314 WRITE(numout,*) '~~~~~~~ ' 304 315 WRITE(numout,*) 305 WRITE(numout,*) ' Total power consumption by av_wave : ztpc = ', ztpc* 1.e-12_wp, 'TW'316 WRITE(numout,*) ' Total power consumption by av_wave = ', zztmp * 1.e-12_wp, 'TW' 306 317 ENDIF 307 318 ENDIF … … 323 334 CALL iom_put( "av_ratio", zav_ratio ) 324 335 DO jk = 2, jpkm1 !* update momentum & tracer diffusivity with wave-driven mixing 325 avs(:,:,jk) =avs(:,:,jk) + zav_wave(:,:,jk) * zav_ratio(:,:,jk)326 avt(:,:,jk) =avt(:,:,jk) + zav_wave(:,:,jk)327 avm(:,:,jk) =avm(:,:,jk) + zav_wave(:,:,jk)336 p_avs(:,:,jk) = p_avs(:,:,jk) + zav_wave(:,:,jk) * zav_ratio(:,:,jk) 337 p_avt(:,:,jk) = p_avt(:,:,jk) + zav_wave(:,:,jk) 338 p_avm(:,:,jk) = p_avm(:,:,jk) + zav_wave(:,:,jk) 328 339 END DO 329 340 ! 330 341 ELSE !* update momentum & tracer diffusivity with wave-driven mixing 331 342 DO jk = 2, jpkm1 332 avs(:,:,jk) =avs(:,:,jk) + zav_wave(:,:,jk)333 avt(:,:,jk) =avt(:,:,jk) + zav_wave(:,:,jk)334 avm(:,:,jk) =avm(:,:,jk) + zav_wave(:,:,jk)343 p_avs(:,:,jk) = p_avs(:,:,jk) + zav_wave(:,:,jk) 344 p_avt(:,:,jk) = p_avt(:,:,jk) + zav_wave(:,:,jk) 345 p_avm(:,:,jk) = p_avm(:,:,jk) + zav_wave(:,:,jk) 335 346 END DO 336 347 ENDIF … … 353 364 CALL iom_put( "emix_iwm", emix_iwm ) 354 365 355 CALL wrk_dealloc( jpi,jpj, zfact, zhdep )356 CALL wrk_dealloc( jpi,jpj,jpk, zwkb, zweight, znu_t, znu_w, zReb )357 358 366 IF(ln_ctl) CALL prt_ctl(tab3d_1=zav_wave , clinfo1=' iwm - av_wave: ', tab3d_2=avt, clinfo2=' avt: ', ovlap=1, kdim=jpk) 359 367 ! -
branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfphy.F90
r7990 r8093 11 11 !! zdf_phy : upadate at each time-step the vertical mixing coeff. 12 12 !!---------------------------------------------------------------------- 13 USE par_oce ! ocean parameters13 USE oce ! ocean dynamics and tracers variables 14 14 USE zdf_oce ! vertical physics: shared variables 15 USE zdfdrg ! vertical physics: top/bottom drag coef. 16 !!gm old 15 17 USE zdfbfr ! vertical physics: bottom friction 16 USE zdfric ! vertical physics: Richardson vertical mixing 18 !!gm 19 USE zdfsh2 ! vertical physics: shear production term of TKE 20 USE zdfric ! vertical physics: RIChardson dependent vertical mixing 17 21 USE zdftke ! vertical physics: TKE vertical mixing 18 22 USE zdfgls ! vertical physics: GLS vertical mixing … … 44 48 INTEGER, PARAMETER :: np_TKE = 3 ! Turbulente Kinetic Eenergy closure scheme for Kz 45 49 INTEGER, PARAMETER :: np_GLS = 4 ! Generic Length Scale closure scheme for Kz 46 50 51 LOGICAL :: l_zdfsh2 ! shear production term flag (=F for CST, =T otherwise (i.e. TKE, GLS, RIC)) 52 47 53 !!---------------------------------------------------------------------- 48 54 !! NEMO/OPA 4.0 , NEMO Consortium (2017) … … 61 67 !! set horizontal shape and vertical profile of background mixing coef. 62 68 !!---------------------------------------------------------------------- 69 INTEGER :: jk ! dummy loop indices 63 70 INTEGER :: ioptio, ios ! local integers 64 71 !! … … 90 97 WRITE(numout,*) ' vertical closure scheme' 91 98 WRITE(numout,*) ' constant vertical mixing coefficient ln_zdfcst = ', ln_zdfcst 92 WRITE(numout,*) ' constant vertical mixing coefficientln_zdfric = ', ln_zdfric93 WRITE(numout,*) ' constant vertical mixing coefficientln_zdftke = ', ln_zdftke94 WRITE(numout,*) ' constant vertical mixing coefficientln_zdfgls = ', ln_zdfgls99 WRITE(numout,*) ' Richardson number dependent closure ln_zdfric = ', ln_zdfric 100 WRITE(numout,*) ' Turbulent Kinetic Energy closure (TKE) ln_zdftke = ', ln_zdftke 101 WRITE(numout,*) ' Generic Length Scale closure (GLS) ln_zdfgls = ', ln_zdfgls 95 102 WRITE(numout,*) ' convection: ' 96 103 WRITE(numout,*) ' enhanced vertical diffusion ln_zdfevd = ', ln_zdfevd … … 137 144 ENDIF 138 145 ! 139 ! ! set 1st/last level Av to zero one for all 140 avt(:,:,1) = 0._wp ; avs(:,:,1) = 0._wp ; avm(:,:,1) = 0._wp 146 DO jk = 1, jpk ! set turbulent closure Kz to the background value (avt_k, avm_k) 147 avt_k(:,:,jk) = avtb_2d(:,:) * avtb(jk) * wmask (:,:,jk) 148 avm_k(:,:,jk) = avmb(jk) * wmask (:,:,jk) 149 END DO 150 !!gm to be tested only the 1st & last levels 151 ! avt (:,:, 1 ) = 0._wp ; avs(:,:, 1 ) = 0._wp ; avm (:,:, 1 ) = 0._wp 152 ! avt (:,:,jpk) = 0._wp ; avs(:,:,jpk) = 0._wp ; avm (:,:,jpk) = 0._wp 153 !!gm 154 avt (:,:,:) = 0._wp ; avs(:,:,:) = 0._wp ; avm (:,:,:) = 0._wp 141 155 142 156 ! !== Convection ==! … … 170 184 IF( ln_zdfric .OR. ln_zdfgls ) CALL ctl_stop( 'zdf_phy_init: zdfric and zdfgls never tested with ice shelves cavities ' ) 171 185 ENDIF 186 ! ! shear production term flag 187 IF( ln_zdfcst ) THEN ; l_zdfsh2 = .FALSE. 188 ELSE ; l_zdfsh2 = .TRUE. 189 ENDIF 172 190 173 191 ! !== gravity wave-driven mixing ==! … … 175 193 IF( ln_zdfswm ) CALL zdf_swm_init ! surface wave-driven mixing 176 194 177 ! !== bottom friction ==! 195 ! !== top/bottom friction ==! 196 CALL zdf_drg_init 197 !!gm old 178 198 CALL zdf_bfr_init 199 !!gm 179 200 ! 180 201 ! !== time-stepping ==! … … 200 221 ! 201 222 INTEGER :: ji, jj, jk ! dummy loop indice 223 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zsh2 ! shear production 202 224 !! --------------------------------------------------------------------- 203 225 ! 226 !!gm old 204 227 CALL zdf_bfr( kt ) !* bottom friction (if quadratic) 205 ! 228 !!gm 229 ! 230 ! IF( l_zdfdrg ) THEN !== update top/bottom drag ==! (non-linear cases) 231 ! ! 232 ! ! !* bottom drag 233 ! CALL zdf_drg( kt, mbkt , r_Cdmin_bot, r_Cdmax_bot, & ! <<== in 234 ! & r_z0_bot, r_ke0_bot, rCd0_bot, & 235 ! & rCdU_bot ) ! ==>> out : bottom drag [m/s] 236 ! IF( ln_isfcav ) THEN !* top drag (ocean cavities) 237 ! CALL zdf_drg( kt, mikt , r_Cdmin_top, r_Cdmax_top, & ! <<== in 238 ! & r_z0_top, r_ke0_top, rCd0_top, & 239 ! & rCdU_top ) ! ==>> out : bottom drag [m/s] 240 ! ENDIF 241 ! ENDIF 242 ! 243 ! !== Kz from chosen turbulent closure ==! (avm_k, avt_k) 244 ! 245 IF( l_zdfsh2 ) & !* shear production at w-points (energy conserving form) 246 CALL zdf_sh2( ub, vb, un, vn, avm_k, & ! <<== in 247 & zsh2 ) ! ==>> out : shear production 248 ! 206 249 SELECT CASE ( nzdf_phy ) !* Vertical eddy viscosity and diffusivity coefficients at w-points 207 CASE( np_RIC ) ; CALL zdf_ric( kt ) ! Richardson number dependent Kz 208 CASE( np_TKE ) ; CALL zdf_tke( kt ) ! TKE closure scheme for Kz 209 CASE( np_GLS ) ; CALL zdf_gls( kt ) ! GLS closure scheme for Kz 210 CASE( np_CST ) ! Constant Kz (reset avt, avm to the background value) 211 avt(2:jpim1,2:jpjm1,1:jpkm1) = rn_avt0 * wmask(2:jpim1,2:jpjm1,1:jpkm1) 212 avm(2:jpim1,2:jpjm1,1:jpkm1) = rn_avm0 * wmask(2:jpim1,2:jpjm1,1:jpkm1) 250 CASE( np_RIC ) ; CALL zdf_ric( kt, gdept_n, zsh2, avm_k, avt_k ) ! Richardson number dependent Kz 251 CASE( np_TKE ) ; CALL zdf_tke( kt , zsh2, avm_k, avt_k ) ! TKE closure scheme for Kz 252 CASE( np_GLS ) ; CALL zdf_gls( kt , zsh2, avm_k, avt_k ) ! GLS closure scheme for Kz 253 ! CASE( np_CST ) ! Constant Kz (reset avt, avm to the background value) 254 ! ! avt_k and avm_k set one for all at initialisation phase 255 !!gm avt(2:jpim1,2:jpjm1,1:jpkm1) = rn_avt0 * wmask(2:jpim1,2:jpjm1,1:jpkm1) 256 !!gm avm(2:jpim1,2:jpjm1,1:jpkm1) = rn_avm0 * wmask(2:jpim1,2:jpjm1,1:jpkm1) 213 257 END SELECT 258 ! 259 ! !== ocean Kz ==! (avt, avs, avm) 260 ! 261 ! !* start from turbulent closure values 262 avt(:,:,2:jpkm1) = avt_k(:,:,2:jpkm1) 263 avm(:,:,2:jpkm1) = avm_k(:,:,2:jpkm1) 214 264 ! 215 265 IF( ln_rnf_mouth ) THEN !* increase diffusivity at rivers mouths … … 219 269 ENDIF 220 270 ! 221 IF( ln_zdfevd ) CALL zdf_evd( kt )!* convection: enhanced vertical eddy diffusivity271 IF( ln_zdfevd ) CALL zdf_evd( kt, avm, avt ) !* convection: enhanced vertical eddy diffusivity 222 272 ! 223 273 ! !* double diffusive mixing 224 274 IF( ln_zdfddm ) THEN ! update avt and compute avs 225 CALL zdf_ddm( kt )275 CALL zdf_ddm( kt, avm, avt, avs ) 226 276 ELSE ! same mixing on all tracers 227 277 avs(2:jpim1,2:jpjm1,1:jpkm1) = avt(2:jpim1,2:jpjm1,1:jpkm1) … … 229 279 ! 230 280 ! !* wave-induced mixing 231 IF( ln_zdfswm ) CALL zdf_swm( kt )! surface wave (Qiao et al. 2004)232 IF( ln_zdfiwm ) CALL zdf_iwm( kt )! internal wave (de Lavergne et al 2017)281 IF( ln_zdfswm ) CALL zdf_swm( kt, avm, avt, avs ) ! surface wave (Qiao et al. 2004) 282 IF( ln_zdfiwm ) CALL zdf_iwm( kt, avm, avt, avs ) ! internal wave (de Lavergne et al 2017) 233 283 234 284 235 285 ! !* Lateral boundary conditions (sign unchanged) 236 CALL lbc_lnk( avm, 'W', 1. ) 237 CALL lbc_lnk( avt, 'W', 1. ) 238 CALL lbc_lnk( avs, 'W', 1. ) 239 ! 240 !!gm ===>>> TO BE REMOVED 241 DO jk = 1, jpkm1 !* vertical eddy viscosity at wu- and wv-points 242 DO jj = 2, jpjm1 243 DO ji = 2, jpim1 244 avmu(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji+1,jj ,jk) ) * wumask(ji,jj,jk) 245 avmv(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji ,jj+1,jk) ) * wvmask(ji,jj,jk) 246 END DO 247 END DO 248 END DO 249 CALL lbc_lnk( avmu, 'U', 1. ) ; CALL lbc_lnk( avmv, 'V', 1. ) ! Lateral boundary conditions 250 !!gm end 286 CALL lbc_lnk( avm_k, 'W', 1. ) ! needed to compute the shear production term 287 CALL lbc_lnk( avt_k, 'W', 1. ) !!gm a priori useless ==>> to be tested 288 CALL lbc_lnk( avm , 'W', 1. ) ! needed to compute avm at u- and v-points 289 CALL lbc_lnk( avt , 'W', 1. ) !!gm a priori only avm_k and avm are required 290 CALL lbc_lnk( avs , 'W', 1. ) !!gm To be tested 291 ! 251 292 252 293 253 294 CALL zdf_mxl( kt ) !* mixed layer depth, and level 254 295 255 !!gm !==>> to be moved in zdftke & zdfgls 256 ! write TKE or GLS information in the restart file 257 IF( lrst_oce .AND. ln_zdftke ) CALL tke_rst( kt, 'WRITE' ) 258 IF( lrst_oce .AND. ln_zdfgls ) CALL gls_rst( kt, 'WRITE' ) 259 ! 296 IF( lrst_oce ) THEN !* write TKE, GLS or RIC fields in the restart file 297 IF( ln_zdftke ) CALL tke_rst( kt, 'WRITE' ) 298 IF( ln_zdfgls ) CALL gls_rst( kt, 'WRITE' ) 299 IF( ln_zdfric ) CALL ric_rst( kt, 'WRITE' ) 300 ENDIF 301 ! 260 302 END SUBROUTINE zdf_phy 261 303 -
branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfric.F90
r7991 r8093 16 16 17 17 !!---------------------------------------------------------------------- 18 !! zdf_ric_init : initialization, namelist read, & parameters control 18 19 !! zdf_ric : update momentum and tracer Kz from the Richardson number 19 !! zdf_ric_init : initialization, namelist read, & parameters control20 !! ric_rst : read/write RIC restart in ocean restart file 20 21 !!---------------------------------------------------------------------- 21 22 USE oce ! ocean dynamics and tracers variables 22 23 USE dom_oce ! ocean space and time domain variables 23 24 USE zdf_oce ! vertical physics: variables 24 USE zdfsh2 ! vertical physics: shear production term of TKE25 25 USE phycst ! physical constants 26 26 USE sbc_oce, ONLY : taum 27 USE eosbn2 , ONLY : neos28 27 ! 29 28 USE in_out_manager ! I/O manager 30 USE lbclnk ! ocean lateral boundary condition (or mpp link) 31 USE lib_mpp ! MPP library 29 USE iom ! I/O manager library 32 30 USE timing ! Timing 33 31 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) … … 37 35 PRIVATE 38 36 39 PUBLIC zdf_ric ! called by step.F90 40 PUBLIC zdf_ric_init ! called by opa.F90 37 PUBLIC zdf_ric ! called by zdfphy.F90 38 PUBLIC ric_rst ! called by zdfphy.F90 39 PUBLIC zdf_ric_init ! called by nemogcm.F90 41 40 42 41 ! !!* Namelist namzdf_ric : Richardson number dependent Kz * … … 54 53 # include "vectopt_loop_substitute.h90" 55 54 !!---------------------------------------------------------------------- 56 !! NEMO/OPA 4.0 , NEMO Consortium (201 1)55 !! NEMO/OPA 4.0 , NEMO Consortium (2017) 57 56 !! $Id$ 58 57 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 59 58 !!---------------------------------------------------------------------- 60 59 CONTAINS 61 62 SUBROUTINE zdf_ric( kt )63 !!----------------------------------------------------------------------64 !! *** ROUTINE zdfric ***65 !!66 !! ** Purpose : Compute the before eddy viscosity and diffusivity as67 !! a function of the local richardson number.68 !!69 !! ** Method : Local richardson number dependent formulation of the70 !! vertical eddy viscosity and diffusivity coefficients.71 !! The eddy coefficients are given by:72 !! avm = avm0 + avmb73 !! avt = avm0 / (1 + rn_alp*ri)74 !! with ri = N^2 / dz(u)**275 !! = e3w**2 * rn2/[ mi( dk(ub) )+mj( dk(vb) ) ]76 !! avm0= rn_avmri / (1 + rn_alp*ri)**nn_ric77 !! Where ri is the before local Richardson number,78 !! rn_avmri is the maximum value reaches by avm and avt79 !! avmb and avtb are the background (or minimum) values80 !! and rn_alp, nn_ric are adjustable parameters.81 !! Typical values used are : avm0=1.e-2 m2/s, avmb=1.e-6 m2/s82 !! avtb=1.e-7 m2/s, rn_alp=5. and nn_ric=2.83 !! a numerical threshold is impose on the vertical shear (1.e-20)84 !! As second step compute Ekman depth from wind stress forcing85 !! and apply namelist provided vertical coeff within this depth.86 !! The Ekman depth is:87 !! Ustar = SQRT(Taum/rho0)88 !! ekd= rn_ekmfc * Ustar / f089 !! Large et al. (1994, eq.29) suggest rn_ekmfc=0.7; however, the derivation90 !! of the above equation indicates the value is somewhat arbitrary; therefore91 !! we allow the freedom to increase or decrease this value, if the92 !! Ekman depth estimate appears too shallow or too deep, respectively.93 !! Ekd is then limited by rn_mldmin and rn_mldmax provided in the94 !! namelist95 !! N.B. the mask are required for implicit scheme, and surface96 !! and bottom value already set in zdfphy.F9097 !!98 !! ** Action : avm, avt mixing coeff (inner domain values only)99 !!100 !! References : Pacanowski & Philander 1981, JPO, 1441-1451.101 !! PFJ Lermusiaux 2001.102 !!----------------------------------------------------------------------103 INTEGER, INTENT(in) :: kt ! ocean time-step104 !!105 INTEGER :: ji, jj, jk ! dummy loop indices106 LOGICAL :: ll_av_weight = .TRUE. ! local logical107 REAL(wp) :: zcfRi, zav, zustar ! local scalars108 REAL(wp), DIMENSION(jpi,jpj) :: zh_ekm ! 2D workspace109 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zsh2 ! 3D -110 !!----------------------------------------------------------------------111 !112 IF( nn_timing == 1 ) CALL timing_start('zdf_ric')113 !114 ! !== avm and avt = F(Richardson number) ==!115 !116 ! !* Shear production at uw- and vw-points (energy conserving form)117 CALL zdf_sh2( zsh2 )118 !119 DO jk = 2, jpkm1 !* Vertical eddy viscosity and diffusivity coefficients120 DO jj = 1, jpjm1121 DO ji = 1, jpim1122 ! ! coefficient = F(richardson number) (avm-weighted Ri)123 zcfRi = 1._wp / ( 1._wp + rn_alp * MAX( 0._wp , avm(ji,jj,jk) * rn2(ji,jj,jk) / ( zsh2(ji,jj,jk) + 1.e-20 ) ) )124 zav = rn_avmri * zcfRi**nn_ric125 !126 ! ! avm and avt coefficients127 avm(ji,jj,jk) = MAX( zav , avmb(jk) ) * wmask(ji,jj,jk)128 avt(ji,jj,jk) = MAX( zav * zcfRi , avtb(jk) ) * wmask(ji,jj,jk)129 END DO130 END DO131 END DO132 !133 !!gm BUG <<<<==== This param can't work at low latitude134 !!gm it provides there much to thick mixed layer ( summer 150m in GYRE configuration !!! )135 !136 IF( ln_mldw ) THEN !== set a minimum value in the Ekman layer ==!137 !138 DO jj = 2, jpjm1 !* Ekman depth139 DO ji = 2, jpim1140 zustar = SQRT( taum(ji,jj) * r1_rau0 )141 zh_ekm(ji,jj) = rn_ekmfc * zustar / ( ABS( ff_t(ji,jj) ) + rsmall ) ! Ekman depth142 zh_ekm(ji,jj) = MAX( rn_mldmin , MIN( zh_ekm(ji,jj) , rn_mldmax ) ) ! set allowed rang143 END DO144 END DO145 DO jk = 2, jpkm1 !* minimum mixing coeff. within the Ekman layer146 DO jj = 2, jpjm1147 DO ji = 2, jpim1148 IF( gdept_n(ji,jj,jk) < zh_ekm(ji,jj) ) THEN149 avm(ji,jj,jk) = MAX( avm(ji,jj,jk), rn_wvmix ) * wmask(ji,jj,jk)150 avt(ji,jj,jk) = MAX( avt(ji,jj,jk), rn_wtmix ) * wmask(ji,jj,jk)151 ENDIF152 END DO153 END DO154 END DO155 ENDIF156 !157 IF( nn_timing == 1 ) CALL timing_stop('zdf_ric')158 !159 END SUBROUTINE zdf_ric160 161 60 162 61 SUBROUTINE zdf_ric_init … … 205 104 ENDIF 206 105 ! 207 DO jk = 1, jpk ! Initialization of vertical eddy coef. to the background value 208 avt(:,:,jk) = avtb(jk) * tmask(:,:,jk) 209 avm(:,:,jk) = avmb(jk) * umask(:,:,jk) 106 CALL ric_rst( nit000, 'READ' ) !* read or initialize all required files 107 ! 108 END SUBROUTINE zdf_ric_init 109 110 111 SUBROUTINE zdf_ric( kt, pdept, p_sh2, p_avm, p_avt ) 112 !!---------------------------------------------------------------------- 113 !! *** ROUTINE zdfric *** 114 !! 115 !! ** Purpose : Compute the before eddy viscosity and diffusivity as 116 !! a function of the local richardson number. 117 !! 118 !! ** Method : Local richardson number dependent formulation of the 119 !! vertical eddy viscosity and diffusivity coefficients. 120 !! The eddy coefficients are given by: 121 !! avm = avm0 + avmb 122 !! avt = avm0 / (1 + rn_alp*ri) 123 !! with ri = N^2 / dz(u)**2 124 !! = e3w**2 * rn2/[ mi( dk(ub) )+mj( dk(vb) ) ] 125 !! avm0= rn_avmri / (1 + rn_alp*Ri)**nn_ric 126 !! where ri is the before local Richardson number, 127 !! rn_avmri is the maximum value reaches by avm and avt 128 !! and rn_alp, nn_ric are adjustable parameters. 129 !! Typical values : rn_alp=5. and nn_ric=2. 130 !! 131 !! As second step compute Ekman depth from wind stress forcing 132 !! and apply namelist provided vertical coeff within this depth. 133 !! The Ekman depth is: 134 !! Ustar = SQRT(Taum/rho0) 135 !! ekd= rn_ekmfc * Ustar / f0 136 !! Large et al. (1994, eq.24) suggest rn_ekmfc=0.7; however, the derivation 137 !! of the above equation indicates the value is somewhat arbitrary; therefore 138 !! we allow the freedom to increase or decrease this value, if the 139 !! Ekman depth estimate appears too shallow or too deep, respectively. 140 !! Ekd is then limited by rn_mldmin and rn_mldmax provided in the 141 !! namelist 142 !! N.B. the mask are required for implicit scheme, and surface 143 !! and bottom value already set in zdfphy.F90 144 !! 145 !! ** Action : avm, avt mixing coeff (inner domain values only) 146 !! 147 !! References : Pacanowski & Philander 1981, JPO, 1441-1451. 148 !! PFJ Lermusiaux 2001. 149 !!---------------------------------------------------------------------- 150 INTEGER , INTENT(in ) :: kt ! ocean time-step 151 REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: pdept ! depth of t-point [m] 152 REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: p_sh2 ! shear production term 153 REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: p_avm, p_avt ! momentum and tracer Kz (w-points) 154 !! 155 INTEGER :: ji, jj, jk ! dummy loop indices 156 REAL(wp) :: zcfRi, zav, zustar, zhek ! local scalars 157 REAL(wp), DIMENSION(jpi,jpj) :: zh_ekm ! 2D workspace 158 !!---------------------------------------------------------------------- 159 ! 160 IF( nn_timing == 1 ) CALL timing_start('zdf_ric') 161 ! 162 ! !== avm and avt = F(Richardson number) ==! 163 DO jk = 2, jpkm1 164 DO jj = 1, jpjm1 165 DO ji = 1, jpim1 ! coefficient = F(richardson number) (avm-weighted Ri) 166 zcfRi = 1._wp / ( 1._wp + rn_alp * MAX( 0._wp , avm(ji,jj,jk) * rn2(ji,jj,jk) / ( p_sh2(ji,jj,jk) + 1.e-20 ) ) ) 167 zav = rn_avmri * zcfRi**nn_ric 168 ! ! avm and avt coefficients 169 p_avm(ji,jj,jk) = MAX( zav , avmb(jk) ) * wmask(ji,jj,jk) 170 p_avt(ji,jj,jk) = MAX( zav * zcfRi , avtb(jk) ) * wmask(ji,jj,jk) 171 END DO 172 END DO 210 173 END DO 211 174 ! 212 END SUBROUTINE zdf_ric_init 175 !!gm BUG <<<<==== This param can't work at low latitude 176 !!gm it provides there much to thick mixed layer ( summer 150m in GYRE configuration !!! ) 177 ! 178 IF( ln_mldw ) THEN !== set a minimum value in the Ekman layer ==! 179 ! 180 DO jj = 2, jpjm1 !* Ekman depth 181 DO ji = 2, jpim1 182 zustar = SQRT( taum(ji,jj) * r1_rau0 ) 183 zhek = rn_ekmfc * zustar / ( ABS( ff_t(ji,jj) ) + rsmall ) ! Ekman depth 184 zh_ekm(ji,jj) = MAX( rn_mldmin , MIN( zhek , rn_mldmax ) ) ! set allowed range 185 END DO 186 END DO 187 DO jk = 2, jpkm1 !* minimum mixing coeff. within the Ekman layer 188 DO jj = 2, jpjm1 189 DO ji = 2, jpim1 190 IF( pdept(ji,jj,jk) < zh_ekm(ji,jj) ) THEN 191 p_avm(ji,jj,jk) = MAX( p_avm(ji,jj,jk), rn_wvmix ) * wmask(ji,jj,jk) 192 p_avt(ji,jj,jk) = MAX( p_avt(ji,jj,jk), rn_wtmix ) * wmask(ji,jj,jk) 193 ENDIF 194 END DO 195 END DO 196 END DO 197 ENDIF 198 ! 199 IF( nn_timing == 1 ) CALL timing_stop('zdf_ric') 200 ! 201 END SUBROUTINE zdf_ric 202 203 204 SUBROUTINE ric_rst( kt, cdrw ) 205 !!--------------------------------------------------------------------- 206 !! *** ROUTINE ric_rst *** 207 !! 208 !! ** Purpose : Read or write TKE file (en) in restart file 209 !! 210 !! ** Method : use of IOM library 211 !! if the restart does not contain TKE, en is either 212 !! set to rn_emin or recomputed 213 !!---------------------------------------------------------------------- 214 INTEGER , INTENT(in) :: kt ! ocean time-step 215 CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag 216 ! 217 INTEGER :: jit, jk ! dummy loop indices 218 INTEGER :: id1, id2 ! local integers 219 !!---------------------------------------------------------------------- 220 ! 221 IF( TRIM(cdrw) == 'READ' ) THEN ! Read/initialise 222 ! ! --------------- 223 ! !* Read the restart file 224 IF( ln_rstart ) THEN 225 id1 = iom_varid( numror, 'avt_k', ldstop = .FALSE. ) 226 id2 = iom_varid( numror, 'avm_k', ldstop = .FALSE. ) 227 ! 228 IF( MIN( id1, id2 ) > 0 ) THEN ! restart exists => read it 229 CALL iom_get( numror, jpdom_autoglo, 'avt_k', avt_k ) 230 CALL iom_get( numror, jpdom_autoglo, 'avm_k', avm_k ) 231 ENDIF 232 ENDIF 233 ! !* otherwise Kz already set to the background value in zdf_phy_init 234 ! 235 ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN ! Create restart file 236 ! ! ------------------- 237 IF(lwp) WRITE(numout,*) '---- ric-rst ----' 238 CALL iom_rstput( kt, nitrst, numrow, 'avt_k', avt_k ) 239 CALL iom_rstput( kt, nitrst, numrow, 'avm_k', avm_k ) 240 ! 241 ENDIF 242 ! 243 END SUBROUTINE ric_rst 213 244 214 245 !!====================================================================== -
branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfsh2.F90
r7997 r8093 11 11 !! zdf_sh2 : compute mixing the shear production term of TKE 12 12 !!---------------------------------------------------------------------- 13 USE oce ! ocean: shared variables14 13 USE dom_oce ! domain: ocean 15 USE zdf_oce ! vertical physics: variables16 14 ! 17 15 USE in_out_manager ! I/O manager … … 31 29 CONTAINS 32 30 33 SUBROUTINE zdf_sh2( psh2 )31 SUBROUTINE zdf_sh2( pub, pvb, pun, pvn, p_avm, p_sh2 ) 34 32 !!---------------------------------------------------------------------- 35 33 !! *** ROUTINE zdf_sh2 *** … … 46 44 !! NB: wet-point only horizontal averaging of shear 47 45 !! 48 !! ** Action : - p sh2 shear prod. term at w-point (interior oceandomain only)49 !! 46 !! ** Action : - p_sh2 shear prod. term at w-point (inner domain only) 47 !! ***** 50 48 !! References : Bruchard, OM 2002 51 49 !! --------------------------------------------------------------------- 52 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( out) :: psh2 ! shear production of TKE (w-points) 50 REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: pub, pvb, pun, pvn ! before, now horizontal velocities 51 REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: p_avm ! vertical eddy viscosity (w-points) 52 REAL(wp), DIMENSION(:,:,:) , INTENT( out) :: p_sh2 ! shear production of TKE (w-points) 53 53 ! 54 54 INTEGER :: ji, jj, jk ! dummy loop arguments 55 55 REAL(wp), DIMENSION(jpi,jpj) :: zsh2u, zsh2v ! 2D workspace 56 56 !!-------------------------------------------------------------------- 57 ! 57 58 IF( nn_timing == 1 ) CALL timing_start('zdf_sh2') 58 59 ! … … 60 61 DO jj = 1, jpjm1 !* 2 x shear production at uw- and vw-points (energy conserving form) 61 62 DO ji = 1, jpim1 62 zsh2u(ji,jj) = ( avm(ji+1,jj,jk) +avm(ji,jj,jk) ) &63 & * ( un(ji,jj,jk-1) -un(ji,jj,jk) ) &64 & * ( ub(ji,jj,jk-1) - ub(ji,jj,jk) ) / ( e3uw_n(ji,jj,jk) * e3uw_b(ji,jj,jk) ) * wumask(ji,jj,jk)65 zsh2v(ji,jj) = ( avm(ji,jj+1,jk) +avm(ji,jj,jk) ) &66 & * ( vn(ji,jj,jk-1) -vn(ji,jj,jk) ) &67 & * ( vb(ji,jj,jk-1) -vb(ji,jj,jk) ) / ( e3vw_n(ji,jj,jk) * e3vw_b(ji,jj,jk) ) * wvmask(ji,jj,jk)63 zsh2u(ji,jj) = ( p_avm(ji+1,jj,jk) + p_avm(ji,jj,jk) ) & 64 & * ( pun(ji,jj,jk-1) - pun(ji,jj,jk) ) & 65 & * ( pub(ji,jj,jk-1) - pub(ji,jj,jk) ) / ( e3uw_n(ji,jj,jk) * e3uw_b(ji,jj,jk) ) * wumask(ji,jj,jk) 66 zsh2v(ji,jj) = ( p_avm(ji,jj+1,jk) + p_avm(ji,jj,jk) ) & 67 & * ( pvn(ji,jj,jk-1) - pvn(ji,jj,jk) ) & 68 & * ( pvb(ji,jj,jk-1) - pvb(ji,jj,jk) ) / ( e3vw_n(ji,jj,jk) * e3vw_b(ji,jj,jk) ) * wvmask(ji,jj,jk) 68 69 END DO 69 70 END DO 70 71 DO jj = 2, jpjm1 !* shear production at w-point 71 72 DO ji = 2, jpim1 ! coast mask: =2 at the coast ; =1 otherwise (NB: wmask useless as zsh2 are masked) 72 p sh2(ji,jj,jk) = 0.25 * ( ( zsh2u(ji-1,jj) + zsh2u(ji,jj) ) * ( 2. - umask(ji-1,jj,jk) * umask(ji,jj,jk) ) &73 & + ( zsh2v(ji,jj-1) + zsh2v(ji,jj) ) * ( 2. - vmask(ji,jj-1,jk) * vmask(ji,jj,jk) ) )73 p_sh2(ji,jj,jk) = 0.25 * ( ( zsh2u(ji-1,jj) + zsh2u(ji,jj) ) * ( 2. - umask(ji-1,jj,jk) * umask(ji,jj,jk) ) & 74 & + ( zsh2v(ji,jj-1) + zsh2v(ji,jj) ) * ( 2. - vmask(ji,jj-1,jk) * vmask(ji,jj,jk) ) ) 74 75 END DO 75 76 END DO -
branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfswm.F90
r7990 r8093 34 34 CONTAINS 35 35 36 SUBROUTINE zdf_swm( kt )36 SUBROUTINE zdf_swm( kt, p_avm, p_avt, p_avs ) 37 37 !!--------------------------------------------------------------------- 38 38 !! *** ROUTINE zdf_swm *** … … 51 51 !! reference : Qiao et al. GRL, 2004 52 52 !!--------------------------------------------------------------------- 53 INTEGER, INTENT(in) :: kt ! ocean time step 53 INTEGER , INTENT(in ) :: kt ! ocean time step 54 REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: p_avm ! momentum Kz (w-points) 55 REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: p_avt, p_avs ! tracer Kz (w-points) 54 56 ! 55 57 INTEGER :: ji, jj, jk ! dummy loop indices … … 63 65 zqb = zcoef * hsw(ji,jj) * tsd2d(ji,jj) * EXP( -3. * wnum(ji,jj) * gdepw_n(ji,jj,jk) ) * wmask(ji,jj,jk) 64 66 ! 65 avt(ji,jj,jk) =avt(ji,jj,jk) + zqb66 avs(ji,jj,jk) =avs(ji,jj,jk) + zqb67 avm(ji,jj,jk) =avm(ji,jj,jk) + zqb67 p_avt(ji,jj,jk) = p_avt(ji,jj,jk) + zqb 68 p_avs(ji,jj,jk) = p_avs(ji,jj,jk) + zqb 69 p_avm(ji,jj,jk) = p_avm(ji,jj,jk) + zqb 68 70 END DO 69 71 END DO -
branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90
r7991 r8093 43 43 USE sbc_oce ! surface boundary condition: ocean 44 44 USE zdf_oce ! vertical physics: ocean variables 45 USE zdfsh2 ! vertical physics: shear production term of TKE 45 !!gm new 46 USE zdfdrg ! vertical physics: top/bottom drag coef. 47 !!gm 46 48 USE zdfmxl ! vertical physics: mixed layer 47 USE lbclnk ! ocean lateral boundary conditions (or mpp link)48 USE prtctl ! Print control49 USE in_out_manager ! I/O manager50 USE iom ! I/O manager library51 USE lib_mpp ! MPP library52 USE wrk_nemo ! work arrays53 USE timing ! Timing54 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)55 49 #if defined key_agrif 56 50 USE agrif_opa_interp 57 51 USE agrif_opa_update 58 52 #endif 53 ! 54 USE in_out_manager ! I/O manager 55 USE iom ! I/O manager library 56 USE lib_mpp ! MPP library 57 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 58 USE prtctl ! Print control 59 USE wrk_nemo ! work arrays 60 USE timing ! Timing 61 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 59 62 60 63 IMPLICIT NONE … … 87 90 REAL(wp) :: rhftau_scl = 1.0_wp ! scale factor applied to HF part of taum (nn_etau=3) 88 91 89 REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:) :: htau! depth of tke penetration (nn_htau)90 REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: dissl! now mixing lenght of dissipation91 REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: apdlr! now mixing lenght of dissipation92 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: htau ! depth of tke penetration (nn_htau) 93 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: dissl ! now mixing lenght of dissipation 94 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: apdlr ! now mixing lenght of dissipation 92 95 93 96 !! * Substitutions … … 112 115 113 116 114 SUBROUTINE zdf_tke( kt )117 SUBROUTINE zdf_tke( kt, p_sh2, p_avm, p_avt ) 115 118 !!---------------------------------------------------------------------- 116 119 !! *** ROUTINE zdf_tke *** … … 157 160 !! Bruchard OM 2002 158 161 !!---------------------------------------------------------------------- 159 INTEGER, INTENT(in) :: kt ! ocean time step 162 INTEGER , INTENT(in ) :: kt ! ocean time step 163 REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: p_sh2 ! shear production term 164 REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: p_avm, p_avt ! momentum and tracer Kz (w-points) 160 165 !!---------------------------------------------------------------------- 161 166 ! … … 165 170 #endif 166 171 ! 167 avt(:,:,:) = avt_k(:,:,:) ! restore before Kz computed with TKE alone 168 avm(:,:,:) = avm_k(:,:,:) 169 ! 170 CALL tke_tke ! now tke (en) 171 ! 172 CALL tke_avn ! now avt, avm, dissl 173 ! 174 avt_k(:,:,:) = avt(:,:,:) ! store Kz computed with TKE alone 175 avm_k(:,:,:) = avm(:,:,:) 172 CALL tke_tke( gdepw_n, e3t_n, e3w_n, p_sh2, p_avm, p_avt ) ! now tke (en) 173 ! 174 CALL tke_avn( gdepw_n, e3t_n, e3w_n, p_avm, p_avt ) ! now avt, avm, dissl 176 175 ! 177 176 #if defined key_agrif … … 179 178 IF( .NOT.Agrif_Root() ) CALL Agrif_Update_Tke( kt ) ! children only 180 179 #endif 181 !182 IF( lrst_oce ) CALL tke_rst( kt, 'WRITE' )183 180 ! 184 181 END SUBROUTINE zdf_tke 185 182 186 183 187 SUBROUTINE tke_tke 184 SUBROUTINE tke_tke( pdepw, p_e3t, p_e3w, p_sh2 & 185 & , p_avm, p_avt ) 188 186 !!---------------------------------------------------------------------- 189 187 !! *** ROUTINE tke_tke *** … … 200 198 !! ** Action : - en : now turbulent kinetic energy) 201 199 !! --------------------------------------------------------------------- 202 INTEGER :: ji, jj, jk ! dummy loop arguments 203 !!bfr INTEGER :: ikbu, ikbv, ikbum1, ikbvm1 ! local integers 204 !!bfr INTEGER :: ikbt, ikbumm1, ikbvmm1 ! - - 205 !!bfr REAL(wp) :: zebot ! local scalars 200 REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: pdepw ! depth of w-points 201 REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: p_e3t, p_e3w ! level thickness (t- & w-points) 202 REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: p_sh2 ! shear production term 203 REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: p_avm, p_avt ! vertical eddy viscosity & diffusivity (w-points) 204 ! 205 INTEGER :: ji, jj, jk ! dummy loop arguments 206 !!bfr REAL(wp) :: zebot, zmshu, zmskv ! local scalars 206 207 REAL(wp) :: zrhoa = 1.22 ! Air density kg/m3 207 208 REAL(wp) :: zcdrag = 1.5e-3 ! drag coefficient … … 212 213 REAL(wp) :: zus , zwlc , zind ! - - 213 214 REAL(wp) :: zzd_up, zzd_lw ! - - 214 INTEGER , DIMENSION(jpi,jpj) :: imlc 215 REAL(wp), POINTER, DIMENSION(:,: ) :: zhlc 216 REAL(wp), POINTER, DIMENSION(:,:,:) :: zpelc, zdiag, zd_up, zd_lw 217 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zsh2 215 INTEGER , DIMENSION(jpi,jpj) :: imlc 216 REAL(wp), DIMENSION(jpi,jpj) :: zhlc 217 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zpelc, zdiag, zd_up, zd_lw 218 218 !!-------------------------------------------------------------------- 219 219 ! 220 220 IF( nn_timing == 1 ) CALL timing_start('tke_tke') 221 !222 CALL wrk_alloc( jpi,jpj, zhlc )223 CALL wrk_alloc( jpi,jpj,jpk, zpelc, zdiag, zd_up, zd_lw )224 221 ! 225 222 zbbrau = rn_ebb / rau0 ! Local constant initialisation … … 256 253 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 257 254 ! en(bot) = (rn_ebb0/rau0)*0.5*sqrt(u_botfr^2+v_botfr^2) (min value rn_emin) 255 !!gm old 258 256 !! DO jj = 2, jpjm1 259 257 !! DO ji = fs_2, fs_jpim1 ! vector opt. … … 263 261 !! bfrva(ji,jj-1) * vb(ji,jj-1,mbkv(ji,jj-1) ) 264 262 !! zebot = 0.001875_wp * SQRT( ztx2 * ztx2 + zty2 * zty2 ) ! where 0.001875 = (rn_ebb0/rau0) * 0.5 = 3.75*0.5/1000. 265 !! en (ji,jj,mbkt(ji,jj)+1) = MAX( zebot, rn_emin ) * tmask(ji,jj,1)263 !! en(ji,jj,mbkt(ji,jj)+1) = MAX( zebot, rn_emin ) * ssmask(ji,jj) 266 264 !! END DO 267 265 !! END DO 266 !!gm new 267 !! 268 !! ====>>>> add below an wet-only calculation of u and v at t-point like in zdfsh2: 269 !! zmsku = ( 2. - umask(ji-1,jj,mbkt(ji,jj)) * umask(ji,jj,mbkt(ji,jj)) ) 270 !! zmskv = ( 2. - vmask(ji,jj-1,mbkt(ji,jj)) * vmask(ji,jj,mbkt(ji,jj)) ) 271 !! 272 !! 273 !! DO jj = 2, jpjm1 274 !! DO ji = fs_2, fs_jpim1 ! vector opt. 275 !! zmsku = ( 2. - umask(ji-1,jj,mbkt(ji,jj)) * umask(ji,jj,mbkt(ji,jj)) ) 276 !! zmskv = ( 2. - vmask(ji,jj-1,mbkt(ji,jj)) * vmask(ji,jj,mbkt(ji,jj)) ) 277 !! ! ! where 0.001875 = (rn_ebb0/rau0) * 0.5 = 3.75*0.5/1000. (CAUTION CdU<0) 278 !! zebot = - 0.001875_wp * rCdU_bot(ji,jj) * SQRT( ( zmsku*( ub(ji,jj,mbkt(ji,jj))+ub(ji-1,jj,mbkt(ji,jj)) ) )**2 279 !! & + ( zmskv*( vb(ji,jj,mbkt(ji,jj))+vb(ji,jj-1,mbkt(ji,jj)) ) )**2 ) 280 !! en(ji,jj,mbkt(ji,jj)+1) = MAX( zebot, rn_emin ) * ssmask(ji,jj) 281 !! END DO 282 !! END DO 283 !! IF( ln_isfcav ) THEN !top friction 284 !! DO jj = 2, jpjm1 285 !! DO ji = fs_2, fs_jpim1 ! vector opt. 286 !! zmsku = ( 2. - umask(ji-1,jj,mikt(ji,jj)) * umask(ji,jj,mikt(ji,jj)) ) 287 !! zmskv = ( 2. - vmask(ji,jj-1,mikt(ji,jj)) * vmask(ji,jj,mikt(ji,jj)) ) 288 !! ! ! where 0.001875 = (rn_ebb0/rau0) * 0.5 = 3.75*0.5/1000. (CAUTION CdU<0) 289 !! zebot = - 0.001875_wp * rCdU_top(ji,jj) * SQRT( ( zmsku*( ub(ji,jj,mikt(ji,jj))+ub(ji-1,jj,mikt(ji,jj)) ) )**2 290 !! & + ( zmskv*( vb(ji,jj,mikt(ji,jj))+vb(ji,jj-1,mikt(ji,jj)) ) )**2 ) 291 !! en(ji,jj,mikt(ji,jj)+1) = MAX( zebot, rn_emin ) * (1._wp - tmask(ji,jj,1)) ! masked at ocean surface 292 !! END DO 293 !! END DO 294 !! ENDIF 295 !! 268 296 !!bfr - end commented area 269 297 ! 270 298 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 271 IF( ln_lc ) THEN ! Langmuir circulation source term added to tke 299 IF( ln_lc ) THEN ! Langmuir circulation source term added to tke ! (Axell JGR 2002) 272 300 ! !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 273 301 ! 274 302 ! !* total energy produce by LC : cumulative sum over jk 275 zpelc(:,:,1) = MAX( rn2b(:,:,1), 0._wp ) * gdepw_n(:,:,1) * e3w_n(:,:,1)303 zpelc(:,:,1) = MAX( rn2b(:,:,1), 0._wp ) * pdepw(:,:,1) * p_e3w(:,:,1) 276 304 DO jk = 2, jpk 277 zpelc(:,:,jk) = zpelc(:,:,jk-1) + MAX( rn2b(:,:,jk), 0._wp ) * gdepw_n(:,:,jk) * e3w_n(:,:,jk)305 zpelc(:,:,jk) = zpelc(:,:,jk-1) + MAX( rn2b(:,:,jk), 0._wp ) * pdepw(:,:,jk) * p_e3w(:,:,jk) 278 306 END DO 279 307 ! !* finite Langmuir Circulation depth … … 291 319 DO jj = 1, jpj 292 320 DO ji = 1, jpi 293 zhlc(ji,jj) = gdepw_n(ji,jj,imlc(ji,jj))321 zhlc(ji,jj) = pdepw(ji,jj,imlc(ji,jj)) 294 322 END DO 295 323 END DO … … 300 328 zus = zcof * SQRT( taum(ji,jj) ) ! Stokes drift 301 329 ! ! vertical velocity due to LC 302 zind = 0.5 - SIGN( 0.5, gdepw_n(ji,jj,jk) - zhlc(ji,jj) )303 zwlc = zind * rn_lc * zus * SIN( rpi * gdepw_n(ji,jj,jk) / zhlc(ji,jj) )330 zind = 0.5 - SIGN( 0.5, pdepw(ji,jj,jk) - zhlc(ji,jj) ) 331 zwlc = zind * rn_lc * zus * SIN( rpi * pdepw(ji,jj,jk) / zhlc(ji,jj) ) 304 332 ! ! TKE Langmuir circulation source term 305 333 en(ji,jj,jk) = en(ji,jj,jk) + rdt * MAX(0.,1._wp - 2.*fr_i(ji,jj) ) * ( zwlc * zwlc * zwlc ) & … … 318 346 ! ! zdiag : diagonal zd_up : upper diagonal zd_lw : lower diagonal 319 347 ! 320 ! !* Shear production at uw- and vw-points (energy conserving form)321 CALL zdf_sh2( zsh2 )322 !323 348 IF( nn_pdl == 1 ) THEN !* Prandtl number = F( Ri ) 324 349 DO jk = 2, jpkm1 … … 326 351 DO ji = 2, jpim1 327 352 ! ! local Richardson number 328 zri = MAX( rn2b(ji,jj,jk), 0._wp ) * avm(ji,jj,jk) / ( zsh2(ji,jj,jk) + rn_bshear )353 zri = MAX( rn2b(ji,jj,jk), 0._wp ) * p_avm(ji,jj,jk) / ( p_sh2(ji,jj,jk) + rn_bshear ) 329 354 ! ! inverse of Prandtl number 330 355 apdlr(ji,jj,jk) = MAX( 0.1_wp, ri_cri / MAX( ri_cri , zri ) ) … … 340 365 ! ! A minimum of 2.e-5 m2/s is imposed on TKE vertical 341 366 ! ! eddy coefficient (ensure numerical stability) 342 zzd_up = zcof * MAX( avm(ji,jj,jk+1) +avm(ji,jj,jk ) , 2.e-5_wp ) & ! upper diagonal343 & / ( e3t_n(ji,jj,jk ) * e3w_n(ji,jj,jk ) )344 zzd_lw = zcof * MAX( avm(ji,jj,jk ) +avm(ji,jj,jk-1) , 2.e-5_wp ) & ! lower diagonal345 & / ( e3t_n(ji,jj,jk-1) * e3w_n(ji,jj,jk ) )367 zzd_up = zcof * MAX( p_avm(ji,jj,jk+1) + p_avm(ji,jj,jk ) , 2.e-5_wp ) & ! upper diagonal 368 & / ( p_e3t(ji,jj,jk ) * p_e3w(ji,jj,jk ) ) 369 zzd_lw = zcof * MAX( p_avm(ji,jj,jk ) + p_avm(ji,jj,jk-1) , 2.e-5_wp ) & ! lower diagonal 370 & / ( p_e3t(ji,jj,jk-1) * p_e3w(ji,jj,jk ) ) 346 371 ! 347 372 zd_up(ji,jj,jk) = zzd_up ! Matrix (zdiag, zd_up, zd_lw) … … 350 375 ! 351 376 ! ! right hand side in en 352 en(ji,jj,jk) = en(ji,jj,jk) + rdt * ( zsh2(ji,jj,jk)& ! shear353 & - avt(ji,jj,jk) * rn2(ji,jj,jk)& ! stratification377 en(ji,jj,jk) = en(ji,jj,jk) + rdt * ( p_sh2(ji,jj,jk) & ! shear 378 & - p_avt(ji,jj,jk) * rn2(ji,jj,jk) & ! stratification 354 379 & + zfact3 * dissl(ji,jj,jk) * en(ji,jj,jk) & ! dissipation 355 380 & ) * wmask(ji,jj,jk) … … 407 432 DO jj = 2, jpjm1 408 433 DO ji = fs_2, fs_jpim1 ! vector opt. 409 en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( - gdepw_n(ji,jj,jk) / htau(ji,jj) ) &434 en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -pdepw(ji,jj,jk) / htau(ji,jj) ) & 410 435 & * MAX(0.,1._wp - 2.*fr_i(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1) 411 436 END DO … … 416 441 DO ji = fs_2, fs_jpim1 ! vector opt. 417 442 jk = nmln(ji,jj) 418 en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( - gdepw_n(ji,jj,jk) / htau(ji,jj) ) &443 en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -pdepw(ji,jj,jk) / htau(ji,jj) ) & 419 444 & * MAX(0.,1._wp - 2.*fr_i(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1) 420 445 END DO … … 429 454 zdif = taum(ji,jj) - ztau ! mean of modulus - modulus of the mean 430 455 zdif = rhftau_scl * MAX( 0._wp, zdif + rhftau_add ) ! apply some modifications... 431 en(ji,jj,jk) = en(ji,jj,jk) + zbbrau * zdif * EXP( - gdepw_n(ji,jj,jk) / htau(ji,jj) ) &456 en(ji,jj,jk) = en(ji,jj,jk) + zbbrau * zdif * EXP( -pdepw(ji,jj,jk) / htau(ji,jj) ) & 432 457 & * MAX(0.,1._wp - 2.*fr_i(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1) 433 458 END DO … … 435 460 END DO 436 461 ENDIF 437 !!gm not sure we need this lbc ....438 CALL lbc_lnk( en, 'W', 1. ) ! Lateral boundary conditions (sign unchanged)439 !440 CALL wrk_dealloc( jpi,jpj, zhlc )441 CALL wrk_dealloc( jpi,jpj,jpk, zpelc, zdiag, zd_up, zd_lw )442 462 ! 443 463 IF( nn_timing == 1 ) CALL timing_stop('tke_tke') … … 446 466 447 467 448 SUBROUTINE tke_avn 468 SUBROUTINE tke_avn( pdepw, p_e3t, p_e3w, p_avm, p_avt ) 449 469 !!---------------------------------------------------------------------- 450 470 !! *** ROUTINE tke_avn *** … … 480 500 !! ** Action : - avt, avm : now vertical eddy diffusivity and viscosity (w-point) 481 501 !!---------------------------------------------------------------------- 502 REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: pdepw ! depth (w-points) 503 REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: p_e3t, p_e3w ! level thickness (t- & w-points) 504 REAL(wp), DIMENSION(:,:,:), INTENT( out) :: p_avm, p_avt ! vertical eddy viscosity & diffusivity (w-points) 505 ! 482 506 INTEGER :: ji, jj, jk ! dummy loop indices 483 507 REAL(wp) :: zrn2, zraug, zcoef, zav ! local scalars 484 508 REAL(wp) :: zdku, zdkv, zsqen ! - - 485 509 REAL(wp) :: zemxl, zemlm, zemlp ! - - 486 REAL(wp), POINTER, DIMENSION(:,:,:) :: zmpdl,zmxlm, zmxld510 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zmxlm, zmxld 487 511 !!-------------------------------------------------------------------- 488 512 ! 489 513 IF( nn_timing == 1 ) CALL timing_start('tke_avn') 490 491 CALL wrk_alloc( jpi,jpj,jpk, zmpdl, zmxlm, zmxld )492 514 493 515 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< … … 502 524 ! 503 525 IF( ln_mxl0 ) THEN ! surface mixing length = F(stress) : l=vkarmn*2.e5*taum/(rau0*g) 526 zraug = vkarmn * 2.e5_wp / ( rau0 * grav ) 504 527 DO jj = 2, jpjm1 505 528 DO ji = fs_2, fs_jpim1 506 zraug = vkarmn * 2.e5_wp / ( rau0 * grav )507 529 zmxlm(ji,jj,1) = MAX( rn_mxl0, zraug * taum(ji,jj) * tmask(ji,jj,1) ) 508 530 END DO … … 529 551 ! 530 552 !!gm Not sure of that coding for ISF.... 531 ! where wmask = 0 set zmxlm == e3w_n553 ! where wmask = 0 set zmxlm == p_e3w 532 554 CASE ( 0 ) ! bounded by the distance to surface and bottom 533 555 DO jk = 2, jpkm1 534 556 DO jj = 2, jpjm1 535 557 DO ji = fs_2, fs_jpim1 ! vector opt. 536 zemxl = MIN( gdepw_n(ji,jj,jk) - gdepw_n(ji,jj,mikt(ji,jj)), zmxlm(ji,jj,jk), &537 & gdepw_n(ji,jj,mbkt(ji,jj)+1) - gdepw_n(ji,jj,jk) )558 zemxl = MIN( pdepw(ji,jj,jk) - pdepw(ji,jj,mikt(ji,jj)), zmxlm(ji,jj,jk), & 559 & pdepw(ji,jj,mbkt(ji,jj)+1) - pdepw(ji,jj,jk) ) 538 560 ! wmask prevent zmxlm = 0 if jk = mikt(ji,jj) 539 zmxlm(ji,jj,jk) = zemxl * wmask(ji,jj,jk) + MIN( zmxlm(ji,jj,jk),e3w_n(ji,jj,jk)) * (1 - wmask(ji,jj,jk))540 zmxld(ji,jj,jk) = zemxl * wmask(ji,jj,jk) + MIN( zmxlm(ji,jj,jk),e3w_n(ji,jj,jk)) * (1 - wmask(ji,jj,jk))561 zmxlm(ji,jj,jk) = zemxl * wmask(ji,jj,jk) + MIN( zmxlm(ji,jj,jk) , p_e3w(ji,jj,jk) ) * (1 - wmask(ji,jj,jk)) 562 zmxld(ji,jj,jk) = zemxl * wmask(ji,jj,jk) + MIN( zmxlm(ji,jj,jk) , p_e3w(ji,jj,jk) ) * (1 - wmask(ji,jj,jk)) 541 563 END DO 542 564 END DO … … 547 569 DO jj = 2, jpjm1 548 570 DO ji = fs_2, fs_jpim1 ! vector opt. 549 zemxl = MIN( e3w_n(ji,jj,jk), zmxlm(ji,jj,jk) )571 zemxl = MIN( p_e3w(ji,jj,jk), zmxlm(ji,jj,jk) ) 550 572 zmxlm(ji,jj,jk) = zemxl 551 573 zmxld(ji,jj,jk) = zemxl … … 558 580 DO jj = 2, jpjm1 559 581 DO ji = fs_2, fs_jpim1 ! vector opt. 560 zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk-1) + e3t_n(ji,jj,jk-1), zmxlm(ji,jj,jk) )582 zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk-1) + p_e3t(ji,jj,jk-1), zmxlm(ji,jj,jk) ) 561 583 END DO 562 584 END DO … … 565 587 DO jj = 2, jpjm1 566 588 DO ji = fs_2, fs_jpim1 ! vector opt. 567 zemxl = MIN( zmxlm(ji,jj,jk+1) + e3t_n(ji,jj,jk+1), zmxlm(ji,jj,jk) )589 zemxl = MIN( zmxlm(ji,jj,jk+1) + p_e3t(ji,jj,jk+1), zmxlm(ji,jj,jk) ) 568 590 zmxlm(ji,jj,jk) = zemxl 569 591 zmxld(ji,jj,jk) = zemxl … … 576 598 DO jj = 2, jpjm1 577 599 DO ji = fs_2, fs_jpim1 ! vector opt. 578 zmxld(ji,jj,jk) = MIN( zmxld(ji,jj,jk-1) + e3t_n(ji,jj,jk-1), zmxlm(ji,jj,jk) )600 zmxld(ji,jj,jk) = MIN( zmxld(ji,jj,jk-1) + p_e3t(ji,jj,jk-1), zmxlm(ji,jj,jk) ) 579 601 END DO 580 602 END DO … … 583 605 DO jj = 2, jpjm1 584 606 DO ji = fs_2, fs_jpim1 ! vector opt. 585 zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk+1) + e3t_n(ji,jj,jk+1), zmxlm(ji,jj,jk) )607 zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk+1) + p_e3t(ji,jj,jk+1), zmxlm(ji,jj,jk) ) 586 608 END DO 587 609 END DO … … 609 631 zsqen = SQRT( en(ji,jj,jk) ) 610 632 zav = rn_ediff * zmxlm(ji,jj,jk) * zsqen 611 avm(ji,jj,jk) = MAX( zav, avmb(jk) ) * wmask(ji,jj,jk)612 avt(ji,jj,jk) = MAX( zav, avtb_2d(ji,jj) * avtb(jk) ) * wmask(ji,jj,jk)633 p_avm(ji,jj,jk) = MAX( zav, avmb(jk) ) * wmask(ji,jj,jk) 634 p_avt(ji,jj,jk) = MAX( zav, avtb_2d(ji,jj) * avtb(jk) ) * wmask(ji,jj,jk) 613 635 dissl(ji,jj,jk) = zsqen / zmxld(ji,jj,jk) 614 636 END DO … … 621 643 DO jj = 2, jpjm1 622 644 DO ji = fs_2, fs_jpim1 ! vector opt. 623 avt(ji,jj,jk) = MAX( apdlr(ji,jj,jk) *avt(ji,jj,jk), avtb_2d(ji,jj) * avtb(jk) ) * tmask(ji,jj,jk)645 p_avt(ji,jj,jk) = MAX( apdlr(ji,jj,jk) * p_avt(ji,jj,jk), avtb_2d(ji,jj) * avtb(jk) ) * tmask(ji,jj,jk) 624 646 END DO 625 647 END DO 626 648 END DO 627 649 ENDIF 628 !!gm I'm sure this is needed to compute the shear term at next time-step629 CALL lbc_lnk( avm, 'W', 1. ) ! Lateral boundary conditions (sign unchanged)630 CALL lbc_lnk( avt, 'W', 1. ) ! Lateral boundary conditions on avt (sign unchanged)631 650 632 651 IF(ln_ctl) THEN … … 634 653 CALL prt_ctl( tab3d_1=avm, clinfo1=' tke - m: ', ovlap=1, kdim=jpk ) 635 654 ENDIF 636 !637 CALL wrk_dealloc( jpi,jpj,jpk, zmpdl, zmxlm, zmxld )638 655 ! 639 656 IF( nn_timing == 1 ) CALL timing_stop('tke_avn') … … 737 754 ! !* set vertical eddy coef. to the background value 738 755 DO jk = 1, jpk 739 avt (:,:,jk) = avtb(jk) * wmask(:,:,jk)740 avm (:,:,jk) = avmb(jk) * wmask(:,:,jk)756 avt(:,:,jk) = avtb(jk) * wmask(:,:,jk) 757 avm(:,:,jk) = avmb(jk) * wmask(:,:,jk) 741 758 END DO 742 759 dissl(:,:,:) = 1.e-12_wp … … 748 765 749 766 SUBROUTINE tke_rst( kt, cdrw ) 750 !!--------------------------------------------------------------------- 751 !! *** ROUTINE tke_rst *** 752 !! 753 !! ** Purpose : Read or write TKE file (en) in restart file 754 !! 755 !! ** Method : use of IOM library 756 !! if the restart does not contain TKE, en is either 757 !! set to rn_emin or recomputed 758 !!---------------------------------------------------------------------- 759 INTEGER , INTENT(in) :: kt ! ocean time-step 760 CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag 761 ! 762 INTEGER :: jit, jk ! dummy loop indices 763 INTEGER :: id1, id2, id3, id4 ! local integers 764 !!---------------------------------------------------------------------- 765 ! 766 IF( TRIM(cdrw) == 'READ' ) THEN ! Read/initialise 767 ! ! --------------- 768 IF( ln_rstart ) THEN !* Read the restart file 769 id1 = iom_varid( numror, 'en' , ldstop = .FALSE. ) 770 id2 = iom_varid( numror, 'avt_k', ldstop = .FALSE. ) 771 id3 = iom_varid( numror, 'avm_k', ldstop = .FALSE. ) 772 id4 = iom_varid( numror, 'dissl', ldstop = .FALSE. ) 773 ! 774 IF( id1 > 0 ) THEN ! 'en' exists 775 CALL iom_get( numror, jpdom_autoglo, 'en', en ) 776 IF( MIN( id2, id3, id4 ) > 0 ) THEN ! all required arrays exist 777 CALL iom_get( numror, jpdom_autoglo, 'avt_k', avt_k ) 778 CALL iom_get( numror, jpdom_autoglo, 'avm_k', avm_k ) 779 CALL iom_get( numror, jpdom_autoglo, 'dissl', dissl ) 780 ELSE ! one at least array is missing 781 CALL tke_avn ! compute avt, avm, and dissl (approximation) 782 ENDIF 783 ELSE ! No TKE array found: initialisation 784 IF(lwp) WRITE(numout,*) ' ===>>>> : previous run without tke scheme, en computed by iterative loop' 785 en (:,:,:) = rn_emin * wmask(:,:,:) 786 CALL tke_avn ! recompute avt, avm, and dissl (approximation) 787 avt_k(:,:,:) = avt(:,:,:) 788 avm_k(:,:,:) = avm(:,:,:) 789 ! 790 DO jit = nit000 + 1, nit000 + 10 ; CALL zdf_tke( jit ) ; END DO 791 avt_k(:,:,:) = avt(:,:,:) 792 avm_k(:,:,:) = avm(:,:,:) 793 ENDIF 794 ELSE !* Start from rest 795 en(:,:,:) = rn_emin * tmask(:,:,:) 796 DO jk = 1, jpk ! set the Kz to the background value 797 avt_k(:,:,jk) = avtb(jk) * wmask (:,:,jk) 798 avm_k(:,:,jk) = avmb(jk) * wmask (:,:,jk) 799 END DO 800 ENDIF 801 ! 802 ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN ! Create restart file 803 ! ! ------------------- 804 IF(lwp) WRITE(numout,*) '---- tke-rst ----' 805 CALL iom_rstput( kt, nitrst, numrow, 'en' , en ) 806 CALL iom_rstput( kt, nitrst, numrow, 'avt_k', avt_k ) 807 CALL iom_rstput( kt, nitrst, numrow, 'avm_k', avm_k ) 808 CALL iom_rstput( kt, nitrst, numrow, 'dissl', dissl ) 809 ! 810 ENDIF 811 ! 767 !!--------------------------------------------------------------------- 768 !! *** ROUTINE tke_rst *** 769 !! 770 !! ** Purpose : Read or write TKE file (en) in restart file 771 !! 772 !! ** Method : use of IOM library 773 !! if the restart does not contain TKE, en is either 774 !! set to rn_emin or recomputed 775 !!---------------------------------------------------------------------- 776 INTEGER , INTENT(in) :: kt ! ocean time-step 777 CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag 778 ! 779 INTEGER :: jit, jk ! dummy loop indices 780 INTEGER :: id1, id2, id3, id4 ! local integers 781 !!---------------------------------------------------------------------- 782 ! 783 IF( TRIM(cdrw) == 'READ' ) THEN ! Read/initialise 784 ! ! --------------- 785 IF( ln_rstart ) THEN !* Read the restart file 786 id1 = iom_varid( numror, 'en' , ldstop = .FALSE. ) 787 id2 = iom_varid( numror, 'avt_k', ldstop = .FALSE. ) 788 id3 = iom_varid( numror, 'avm_k', ldstop = .FALSE. ) 789 id4 = iom_varid( numror, 'dissl', ldstop = .FALSE. ) 790 ! 791 IF( MIN( id1, id2, id3, id4 ) > 0 ) THEN ! fields exist 792 CALL iom_get( numror, jpdom_autoglo, 'en', en ) 793 CALL iom_get( numror, jpdom_autoglo, 'avt_k', avt_k ) 794 CALL iom_get( numror, jpdom_autoglo, 'avm_k', avm_k ) 795 CALL iom_get( numror, jpdom_autoglo, 'dissl', dissl ) 796 ELSE ! start TKE from rest 797 IF(lwp) WRITE(numout,*) ' ==>> previous run without TKE scheme, set en to background values' 798 en(:,:,:) = rn_emin * wmask(:,:,:) 799 ! avt_k, avm_k already set to the background value in zdf_phy_init 800 ENDIF 801 ELSE !* Start from rest 802 IF(lwp) WRITE(numout,*) ' ==>> start from rest: set en to the background value' 803 en(:,:,:) = rn_emin * wmask(:,:,:) 804 ! avt_k, avm_k already set to the background value in zdf_phy_init 805 ENDIF 806 ! 807 ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN ! Create restart file 808 ! ! ------------------- 809 IF(lwp) WRITE(numout,*) '---- tke-rst ----' 810 CALL iom_rstput( kt, nitrst, numrow, 'en' , en ) 811 CALL iom_rstput( kt, nitrst, numrow, 'avt_k', avt_k ) 812 CALL iom_rstput( kt, nitrst, numrow, 'avm_k', avm_k ) 813 CALL iom_rstput( kt, nitrst, numrow, 'dissl', dissl ) 814 ! 815 ENDIF 816 ! 812 817 END SUBROUTINE tke_rst 813 818 -
branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/oce.F90
r7646 r8093 63 63 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: riceload 64 64 65 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: rke !: kinetic energy66 67 65 !! arrays relating to embedding ice in the ocean. These arrays need to be declared 68 66 !! even if no ice model is required. In the no ice model or traditional levitating … … 99 97 & rhd (jpi,jpj,jpk) , rhop (jpi,jpj,jpk) , STAT=ierr(1) ) 100 98 ! 101 ALLOCATE(rke(jpi,jpj,jpk) , & 102 & sshb(jpi,jpj) , sshn(jpi,jpj) , ssha(jpi,jpj) , & 103 & ub_b(jpi,jpj) , un_b(jpi,jpj) , ua_b(jpi,jpj) , & 104 & vb_b(jpi,jpj) , vn_b(jpi,jpj) , va_b(jpi,jpj) , & 105 & spgu (jpi,jpj) , spgv(jpi,jpj) , & 106 & gtsu(jpi,jpj,jpts), gtsv(jpi,jpj,jpts), & 107 & gru(jpi,jpj) , grv(jpi,jpj) , & 108 & gtui(jpi,jpj,jpts), gtvi(jpi,jpj,jpts), & 109 & grui(jpi,jpj) , grvi(jpi,jpj) , & 110 & riceload(jpi,jpj), STAT=ierr(2) ) 99 ALLOCATE( sshb(jpi,jpj) , sshn(jpi,jpj) , ssha(jpi,jpj) , & 100 & ub_b(jpi,jpj) , un_b(jpi,jpj) , ua_b(jpi,jpj) , & 101 & vb_b(jpi,jpj) , vn_b(jpi,jpj) , va_b(jpi,jpj) , & 102 & spgu (jpi,jpj) , spgv(jpi,jpj) , & 103 & gtsu(jpi,jpj,jpts), gtsv(jpi,jpj,jpts) , & 104 & gru(jpi,jpj) , grv(jpi,jpj) , & 105 & gtui(jpi,jpj,jpts), gtvi(jpi,jpj,jpts) , & 106 & grui(jpi,jpj) , grvi(jpi,jpj) , & 107 & riceload(jpi,jpj) , STAT=ierr(2) ) 111 108 ! 112 109 ALLOCATE( snwice_mass(jpi,jpj) , snwice_mass_b(jpi,jpj), snwice_fmass(jpi,jpj) , STAT=ierr(3) ) -
branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/step.F90
r7954 r8093 36 36 !!---------------------------------------------------------------------- 37 37 USE step_oce ! time stepping definition modules 38 !!gm to be removed when removing avmu, avmv 39 USE zdf_oce ! ocean vertical physics variables 40 !!gm 38 41 ! 39 42 USE iom ! xIOs server … … 128 131 CALL zdf_phy( kstp ) ! vertical physics update (bfr, avt, avs, avm + MLD) 129 132 133 134 !!gm ===>>> TO BE REMOVED (require changes in zdf_oce, dynzdf(_imp,_exp), dynldf_iso, diawri) 135 DO jk = 1, jpkm1 !* vertical eddy viscosity at wu- and wv-points 136 DO jj = 2, jpjm1 137 DO ji = 2, jpim1 138 avmu(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji+1,jj ,jk) ) * wumask(ji,jj,jk) 139 avmv(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji ,jj+1,jk) ) * wvmask(ji,jj,jk) 140 END DO 141 END DO 142 END DO 143 CALL lbc_lnk( avmu, 'U', 1. ) ; CALL lbc_lnk( avmv, 'V', 1. ) ! Lateral boundary conditions 144 !!gm end 145 146 130 147 ! LATERAL PHYSICS 131 148 ! -
branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/SETTE/sette.sh
r7931 r8093 1098 1098 export TEST_NAME="SHORT" 1099 1099 cd ${CONFIG_DIR0} 1100 . ./makenemo -m ${CMP_NAM} -n ORCA2AGUL_1_2 -r ORCA2_LIM3_PISCES -d "OPA_SRC LIM_SRC_3 NST_SRC" -j 8 add_key "key_agrif" del_key "key_ zdftmx key_top"1100 . ./makenemo -m ${CMP_NAM} -n ORCA2AGUL_1_2 -r ORCA2_LIM3_PISCES -d "OPA_SRC LIM_SRC_3 NST_SRC" -j 8 add_key "key_agrif" del_key "key_top" 1101 1101 cd ${SETTE_DIR} 1102 1102 . ./param.cfg … … 1145 1145 export TEST_NAME="SHORT_NOZOOM" 1146 1146 cd ${CONFIG_DIR0} 1147 . ./makenemo -m ${CMP_NAM} -n ORCA2AGUL_2_2 -r ORCA2_LIM3_PISCES -d "OPA_SRC LIM_SRC_3 NST_SRC" -j 8 add_key "key_agrif" del_key "key_ zdftmx key_top"1147 . ./makenemo -m ${CMP_NAM} -n ORCA2AGUL_2_2 -r ORCA2_LIM3_PISCES -d "OPA_SRC LIM_SRC_3 NST_SRC" -j 8 add_key "key_agrif" del_key "key_top" 1148 1148 cd ${SETTE_DIR} 1149 1149 . ./param.cfg … … 1182 1182 export TEST_NAME="SHORT_NOAGRIF" 1183 1183 cd ${CONFIG_DIR0} 1184 . ./makenemo -m ${CMP_NAM} -n ORCA2AGUL_2_2_NAG -r ORCA2_LIM3_PISCES -d "OPA_SRC LIM_SRC_3 NST_SRC" -j 8 del_key "key_ zdftmx key_top"1184 . ./makenemo -m ${CMP_NAM} -n ORCA2AGUL_2_2_NAG -r ORCA2_LIM3_PISCES -d "OPA_SRC LIM_SRC_3 NST_SRC" -j 8 del_key "key_top" 1185 1185 cd ${SETTE_DIR} 1186 1186 . ./param.cfg … … 1220 1220 export TEST_NAME="LONG" 1221 1221 cd ${CONFIG_DIR0} 1222 . ./makenemo -m ${CMP_NAM} -n ORCA2AGUL_LONG -r ORCA2_LIM3_PISCES -d "OPA_SRC LIM_SRC_3 NST_SRC" -j 8 add_key "key_agrif" del_key "key_ zdftmx key_top"1222 . ./makenemo -m ${CMP_NAM} -n ORCA2AGUL_LONG -r ORCA2_LIM3_PISCES -d "OPA_SRC LIM_SRC_3 NST_SRC" -j 8 add_key "key_agrif" del_key "key_top" 1223 1223 cd ${SETTE_DIR} 1224 1224 . ./param.cfg … … 1319 1319 export TEST_NAME="REPRO_4_4" 1320 1320 cd ${CONFIG_DIR0} 1321 . ./makenemo -m ${CMP_NAM} -n ORCA2AGUL_16 -r ORCA2_LIM3_PISCES -d "OPA_SRC LIM_SRC_3 NST_SRC" -j 8 add_key "key_agrif" del_key "key_ zdftmx key_top"1321 . ./makenemo -m ${CMP_NAM} -n ORCA2AGUL_16 -r ORCA2_LIM3_PISCES -d "OPA_SRC LIM_SRC_3 NST_SRC" -j 8 add_key "key_agrif" del_key "key_top" 1322 1322 cd ${SETTE_DIR} 1323 1323 . ./param.cfg
Note: See TracChangeset
for help on using the changeset viewer.