- Timestamp:
- 2017-12-01T18:44:09+01:00 (6 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfgls.F90
r7646 r8882 5 5 !! turbulent closure parameterization 6 6 !!====================================================================== 7 !! History : 3.0 ! 2009-09 (G. Reffray) Original code 8 !! 3.3 ! 2010-10 (C. Bricaud) Add in the reference 7 !! History : 3.0 ! 2009-09 (G. Reffray) Original code 8 !! 3.3 ! 2010-10 (C. Bricaud) Add in the reference 9 !! 4.0 ! 2017-04 (G. Madec) remove CPP keys & avm at t-point only 10 !! - ! 2017-05 (G. Madec) add top friction as boundary condition 9 11 !!---------------------------------------------------------------------- 10 #if defined key_zdfgls 11 !!---------------------------------------------------------------------- 12 !! 'key_zdfgls' Generic Length Scale vertical physics 12 13 13 !!---------------------------------------------------------------------- 14 14 !! zdf_gls : update momentum and tracer Kz from a gls scheme … … 20 20 USE domvvl ! ocean space and time domain : variable volume layer 21 21 USE zdf_oce ! ocean vertical physics 22 USE zdfbfr ! bottom friction (only for rn_bfrz0) 22 USE zdfdrg , ONLY : r_z0_top , r_z0_bot ! top/bottom roughness 23 USE zdfdrg , ONLY : rCdU_top , rCdU_bot ! top/bottom friction 23 24 USE sbc_oce ! surface boundary condition: ocean 24 25 USE phycst ! physical constants 25 26 USE zdfmxl ! mixed layer 26 USE sbcwave , ONLY: hsw ! significant wave height27 USE sbcwave , ONLY : hsw ! significant wave height 27 28 ! 28 29 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 29 30 USE lib_mpp ! MPP manager 30 USE wrk_nemo ! work arrays31 31 USE prtctl ! Print control 32 32 USE in_out_manager ! I/O manager … … 38 38 PRIVATE 39 39 40 PUBLIC zdf_gls ! routine called in step module 41 PUBLIC zdf_gls_init ! routine called in opa module 42 PUBLIC gls_rst ! routine called in step module 43 44 LOGICAL , PUBLIC, PARAMETER :: lk_zdfgls = .TRUE. !: TKE vertical mixing flag 40 PUBLIC zdf_gls ! called in zdfphy 41 PUBLIC zdf_gls_init ! called in zdfphy 42 PUBLIC gls_rst ! called in zdfphy 43 45 44 ! 46 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: mxln !: now mixing length45 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: hmxl_n !: now mixing length 47 46 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: zwall !: wall function 48 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ustars2 !: Squared surface velocity scale at T-points 49 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ustarb2 !: Squared bottom velocity scale at T-points 47 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ustar2_surf !: Squared surface velocity scale at T-points 48 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ustar2_top !: Squared top velocity scale at T-points 49 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ustar2_bot !: Squared bottom velocity scale at T-points 50 50 51 51 ! !! ** Namelist namzdf_gls ** … … 102 102 REAL(wp) :: rsc_tke, rsc_psi, rpsi1, rpsi2, rpsi3, rsc_psi0 ! - - - - 103 103 REAL(wp) :: rpsi3m, rpsi3p, rpp, rmm, rnn ! - - - - 104 ! 105 REAL(wp) :: r2_3 = 2._wp/3._wp ! constant=2/3 104 106 105 107 !! * Substitutions … … 116 118 !! *** FUNCTION zdf_gls_alloc *** 117 119 !!---------------------------------------------------------------------- 118 ALLOCATE( mxln(jpi,jpj,jpk), zwall(jpi,jpj,jpk) ,&119 & ustars2(jpi,jpj) , ustarb2(jpi,jpj), STAT= zdf_gls_alloc )120 ALLOCATE( hmxl_n(jpi,jpj,jpk) , ustar2_surf(jpi,jpj) , & 121 & zwall (jpi,jpj,jpk) , ustar2_top (jpi,jpj) , ustar2_bot(jpi,jpj) , STAT= zdf_gls_alloc ) 120 122 ! 121 123 IF( lk_mpp ) CALL mpp_sum ( zdf_gls_alloc ) … … 124 126 125 127 126 SUBROUTINE zdf_gls( kt )128 SUBROUTINE zdf_gls( kt, p_sh2, p_avm, p_avt ) 127 129 !!---------------------------------------------------------------------- 128 130 !! *** ROUTINE zdf_gls *** … … 131 133 !! coefficients using the GLS turbulent closure scheme. 132 134 !!---------------------------------------------------------------------- 133 INTEGER, INTENT(in) :: kt ! ocean time step 134 INTEGER :: ji, jj, jk, ibot, ibotm1, dir ! dummy loop arguments 135 REAL(wp) :: zesh2, zsigpsi, zcoef, zex1, zex2 ! local scalars 136 REAL(wp) :: ztx2, zty2, zup, zdown, zcof ! - - 137 REAL(wp) :: zratio, zrn2, zflxb, sh ! - - 135 INTEGER , INTENT(in ) :: kt ! ocean time step 136 REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: p_sh2 ! shear production term 137 REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: p_avm, p_avt ! momentum and tracer Kz (w-points) 138 ! 139 INTEGER :: ji, jj, jk ! dummy loop arguments 140 INTEGER :: ibot, ibotm1 ! local integers 141 INTEGER :: itop, itopp1 ! - - 142 REAL(wp) :: zesh2, zsigpsi, zcoef, zex1 , zex2 ! local scalars 143 REAL(wp) :: ztx2, zty2, zup, zdown, zcof, zdir ! - - 144 REAL(wp) :: zratio, zrn2, zflxb, sh , z_en ! - - 138 145 REAL(wp) :: prod, buoy, diss, zdiss, sm ! - - 139 REAL(wp) :: gh, gm, shr, dif, zsqen, zav ! - - 140 REAL(wp), POINTER, DIMENSION(:,: ) :: zdep 141 REAL(wp), POINTER, DIMENSION(:,: ) :: zkar 142 REAL(wp), POINTER, DIMENSION(:,: ) :: zflxs ! Turbulence fluxed induced by internal waves 143 REAL(wp), POINTER, DIMENSION(:,: ) :: zhsro ! Surface roughness (surface waves) 144 REAL(wp), POINTER, DIMENSION(:,:,:) :: eb ! tke at time before 145 REAL(wp), POINTER, DIMENSION(:,:,:) :: mxlb ! mixing length at time before 146 REAL(wp), POINTER, DIMENSION(:,:,:) :: shear ! vertical shear 147 REAL(wp), POINTER, DIMENSION(:,:,:) :: eps ! dissipation rate 148 REAL(wp), POINTER, DIMENSION(:,:,:) :: zwall_psi ! Wall function use in the wb case (ln_sigpsi) 149 REAL(wp), POINTER, DIMENSION(:,:,:) :: psi ! psi at time now 150 REAL(wp), POINTER, DIMENSION(:,:,:) :: z_elem_a ! element of the first matrix diagonal 151 REAL(wp), POINTER, DIMENSION(:,:,:) :: z_elem_b ! element of the second matrix diagonal 152 REAL(wp), POINTER, DIMENSION(:,:,:) :: z_elem_c ! element of the third matrix diagonal 146 REAL(wp) :: gh, gm, shr, dif, zsqen, zavt, zavm ! - - 147 REAL(wp) :: zmsku, zmskv ! - - 148 REAL(wp), DIMENSION(jpi,jpj) :: zdep 149 REAL(wp), DIMENSION(jpi,jpj) :: zkar 150 REAL(wp), DIMENSION(jpi,jpj) :: zflxs ! Turbulence fluxed induced by internal waves 151 REAL(wp), DIMENSION(jpi,jpj) :: zhsro ! Surface roughness (surface waves) 152 REAL(wp), DIMENSION(jpi,jpj,jpk) :: eb ! tke at time before 153 REAL(wp), DIMENSION(jpi,jpj,jpk) :: hmxl_b ! mixing length at time before 154 REAL(wp), DIMENSION(jpi,jpj,jpk) :: eps ! dissipation rate 155 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zwall_psi ! Wall function use in the wb case (ln_sigpsi) 156 REAL(wp), DIMENSION(jpi,jpj,jpk) :: psi ! psi at time now 157 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zd_lw, zd_up, zdiag ! lower, upper and diagonal of the matrix 158 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zstt, zstm ! stability function on tracer and momentum 153 159 !!-------------------------------------------------------------------- 154 160 ! 155 IF( nn_timing == 1 ) CALL timing_start('zdf_gls') 156 ! 157 CALL wrk_alloc( jpi,jpj, zdep, zkar, zflxs, zhsro ) 158 CALL wrk_alloc( jpi,jpj,jpk, eb, mxlb, shear, eps, zwall_psi, z_elem_a, z_elem_b, z_elem_c, psi ) 159 161 IF( ln_timing ) CALL timing_start('zdf_gls') 162 ! 160 163 ! Preliminary computing 161 164 162 ustars2 = 0._wp ; ustarb2 = 0._wp ; psi = 0._wp ; zwall_psi = 0._wp 163 164 IF( kt /= nit000 ) THEN ! restore before value to compute tke 165 avt (:,:,:) = avt_k (:,:,:) 166 avm (:,:,:) = avm_k (:,:,:) 167 avmu(:,:,:) = avmu_k(:,:,:) 168 avmv(:,:,:) = avmv_k(:,:,:) 169 ENDIF 170 171 ! Compute surface and bottom friction at T-points 165 ustar2_surf(:,:) = 0._wp ; psi(:,:,:) = 0._wp 166 ustar2_top (:,:) = 0._wp ; zwall_psi(:,:,:) = 0._wp 167 ustar2_bot (:,:) = 0._wp 168 169 ! Compute surface, top and bottom friction at T-points 172 170 DO jj = 2, jpjm1 173 171 DO ji = fs_2, fs_jpim1 ! vector opt. 174 172 ! 175 173 ! surface friction 176 ustar s2(ji,jj) = r1_rau0 * taum(ji,jj) * tmask(ji,jj,1)174 ustar2_surf(ji,jj) = r1_rau0 * taum(ji,jj) * tmask(ji,jj,1) 177 175 ! 178 ! bottom friction (explicit before friction) 179 ! Note that we chose here not to bound the friction as in dynbfr) 180 ztx2 = ( bfrua(ji,jj) * ub(ji,jj,mbku(ji,jj)) + bfrua(ji-1,jj) * ub(ji-1,jj,mbku(ji-1,jj)) ) & 181 & * ( 1._wp - 0.5_wp * umask(ji,jj,1) * umask(ji-1,jj,1) ) 182 zty2 = ( bfrva(ji,jj) * vb(ji,jj,mbkv(ji,jj)) + bfrva(ji,jj-1) * vb(ji,jj-1,mbkv(ji,jj-1)) ) & 183 & * ( 1._wp - 0.5_wp * vmask(ji,jj,1) * vmask(ji,jj-1,1) ) 184 ustarb2(ji,jj) = SQRT( ztx2 * ztx2 + zty2 * zty2 ) * tmask(ji,jj,1) 185 END DO 186 END DO 187 188 ! Set surface roughness length 189 SELECT CASE ( nn_z0_met ) 190 ! 191 CASE ( 0 ) ! Constant roughness 176 !!gm Rq we may add here r_ke0(_top/_bot) ? ==>> think about that... 177 ! bottom friction (explicit before friction) 178 zmsku = ( 2._wp - umask(ji-1,jj,mbkt(ji,jj)) * umask(ji,jj,mbkt(ji,jj)) ) 179 zmskv = ( 2._wp - vmask(ji,jj-1,mbkt(ji,jj)) * vmask(ji,jj,mbkt(ji,jj)) ) ! (CAUTION: CdU<0) 180 ustar2_bot(ji,jj) = - rCdU_bot(ji,jj) * SQRT( ( zmsku*( ub(ji,jj,mbkt(ji,jj))+ub(ji-1,jj,mbkt(ji,jj)) ) )**2 & 181 & + ( zmskv*( vb(ji,jj,mbkt(ji,jj))+vb(ji,jj-1,mbkt(ji,jj)) ) )**2 ) 182 END DO 183 END DO 184 IF( ln_isfcav ) THEN !top friction 185 DO jj = 2, jpjm1 186 DO ji = fs_2, fs_jpim1 ! vector opt. 187 zmsku = ( 2. - umask(ji-1,jj,mikt(ji,jj)) * umask(ji,jj,mikt(ji,jj)) ) 188 zmskv = ( 2. - vmask(ji,jj-1,mikt(ji,jj)) * vmask(ji,jj,mikt(ji,jj)) ) ! (CAUTION: CdU<0) 189 ustar2_top(ji,jj) = - rCdU_top(ji,jj) * SQRT( ( zmsku*( ub(ji,jj,mikt(ji,jj))+ub(ji-1,jj,mikt(ji,jj)) ) )**2 & 190 & + ( zmskv*( vb(ji,jj,mikt(ji,jj))+vb(ji,jj-1,mikt(ji,jj)) ) )**2 ) 191 END DO 192 END DO 193 ENDIF 194 195 SELECT CASE ( nn_z0_met ) !== Set surface roughness length ==! 196 CASE ( 0 ) ! Constant roughness 192 197 zhsro(:,:) = rn_hsro 193 198 CASE ( 1 ) ! Standard Charnock formula 194 zhsro(:,:) = MAX( rsbc_zs1 * ustars2(:,:), rn_hsro)199 zhsro(:,:) = MAX( rsbc_zs1 * ustar2_surf(:,:) , rn_hsro ) 195 200 CASE ( 2 ) ! Roughness formulae according to Rascle et al., Ocean Modelling (2008) 196 zdep(:,:) = 30.*TANH(2.*0.3/(28.*SQRT(MAX(ustars2(:,:),rsmall)))) ! Wave age (eq. 10) 197 zhsro(:,:) = MAX(rsbc_zs2 * ustars2(:,:) * zdep(:,:)**1.5, rn_hsro) ! zhsro = rn_frac_hs * Hsw (eq. 11) 201 !!gm faster coding : the 2 comment lines should be used 202 !!gm zcof = 2._wp * 0.6_wp / 28._wp 203 !!gm zdep(:,:) = 30._wp * TANH( zcof/ SQRT( MAX(ustar2_surf(:,:),rsmall) ) ) ! Wave age (eq. 10) 204 zdep (:,:) = 30.*TANH( 2.*0.3/(28.*SQRT(MAX(ustar2_surf(:,:),rsmall))) ) ! Wave age (eq. 10) 205 zhsro(:,:) = MAX(rsbc_zs2 * ustar2_surf(:,:) * zdep(:,:)**1.5, rn_hsro) ! zhsro = rn_frac_hs * Hsw (eq. 11) 198 206 CASE ( 3 ) ! Roughness given by the wave model (coupled or read in file) 199 zhsro(:,:) = hsw(:,:)207 zhsro(:,:) = rn_frac_hs * hsw(:,:) ! (rn_frac_hs=1.6 see Eq. (5) of Rascle et al. 2008 ) 200 208 END SELECT 201 202 ! Compute shear and dissipation rate 203 DO jk = 2, jpkm1 204 DO jj = 2, jpjm1 205 DO ji = fs_2, fs_jpim1 ! vector opt. 206 avmu(ji,jj,jk) = avmu(ji,jj,jk) * ( un(ji,jj,jk-1) - un(ji,jj,jk) ) & 207 & * ( ub(ji,jj,jk-1) - ub(ji,jj,jk) ) & 208 & / ( e3uw_n(ji,jj,jk) * e3uw_b(ji,jj,jk) ) 209 avmv(ji,jj,jk) = avmv(ji,jj,jk) * ( vn(ji,jj,jk-1) - vn(ji,jj,jk) ) & 210 & * ( vb(ji,jj,jk-1) - vb(ji,jj,jk) ) & 211 & / ( e3vw_n(ji,jj,jk) * e3vw_b(ji,jj,jk) ) 212 eps(ji,jj,jk) = rc03 * en(ji,jj,jk) * SQRT(en(ji,jj,jk)) / mxln(ji,jj,jk) 209 ! 210 DO jk = 2, jpkm1 !== Compute dissipation rate ==! 211 DO jj = 1, jpjm1 212 DO ji = 1, jpim1 213 eps(ji,jj,jk) = rc03 * en(ji,jj,jk) * SQRT( en(ji,jj,jk) ) / hmxl_n(ji,jj,jk) 213 214 END DO 214 215 END DO 215 216 END DO 216 !217 ! Lateral boundary conditions (avmu,avmv) (sign unchanged)218 CALL lbc_lnk( avmu, 'U', 1. ) ; CALL lbc_lnk( avmv, 'V', 1. )219 217 220 218 ! Save tke at before time step 221 eb (:,:,:) = en(:,:,:)222 mxlb(:,:,:) = mxln(:,:,:)219 eb (:,:,:) = en (:,:,:) 220 hmxl_b(:,:,:) = hmxl_n(:,:,:) 223 221 224 222 IF( nn_clos == 0 ) THEN ! Mellor-Yamada … … 226 224 DO jj = 2, jpjm1 227 225 DO ji = fs_2, fs_jpim1 ! vector opt. 228 zup = mxln(ji,jj,jk) * gdepw_n(ji,jj,mbkt(ji,jj)+1)226 zup = hmxl_n(ji,jj,jk) * gdepw_n(ji,jj,mbkt(ji,jj)+1) 229 227 zdown = vkarmn * gdepw_n(ji,jj,jk) * ( -gdepw_n(ji,jj,jk) + gdepw_n(ji,jj,mbkt(ji,jj)+1) ) 230 228 zcoef = ( zup / MAX( zdown, rsmall ) ) … … 245 243 ! The surface boundary condition are set after 246 244 ! The bottom boundary condition are also set after. In standard e(bottom)=0. 247 ! z _elem_b : diagonal z_elem_c : upper diagonal z_elem_a: lower diagonal245 ! zdiag : diagonal zd_up : upper diagonal zd_lw : lower diagonal 248 246 ! Warning : after this step, en : right hand side of the matrix 249 247 250 248 DO jk = 2, jpkm1 251 249 DO jj = 2, jpjm1 252 DO ji = fs_2, fs_jpim1 ! vector opt. 253 ! 254 ! shear prod. at w-point weightened by mask 255 shear(ji,jj,jk) = ( avmu(ji-1,jj,jk) + avmu(ji,jj,jk) ) / MAX( 1.e0 , umask(ji-1,jj,jk) + umask(ji,jj,jk) ) & 256 & + ( avmv(ji,jj-1,jk) + avmv(ji,jj,jk) ) / MAX( 1.e0 , vmask(ji,jj-1,jk) + vmask(ji,jj,jk) ) 257 ! 258 ! stratif. destruction 259 buoy = - avt(ji,jj,jk) * rn2(ji,jj,jk) 260 ! 261 ! shear prod. - stratif. destruction 262 diss = eps(ji,jj,jk) 263 ! 264 dir = 0.5_wp + SIGN( 0.5_wp, shear(ji,jj,jk) + buoy ) ! dir =1(=0) if shear(ji,jj,jk)+buoy >0(<0) 265 ! 266 zesh2 = dir*(shear(ji,jj,jk)+buoy)+(1._wp-dir)*shear(ji,jj,jk) ! production term 267 zdiss = dir*(diss/en(ji,jj,jk)) +(1._wp-dir)*(diss-buoy)/en(ji,jj,jk) ! dissipation term 250 DO ji = 2, jpim1 251 ! 252 buoy = - p_avt(ji,jj,jk) * rn2(ji,jj,jk) ! stratif. destruction 253 ! 254 diss = eps(ji,jj,jk) ! dissipation 255 ! 256 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) 257 ! 258 zesh2 = zdir*(p_sh2(ji,jj,jk)+buoy)+(1._wp-zdir)*p_sh2(ji,jj,jk) ! production term 259 zdiss = zdir*(diss/en(ji,jj,jk)) +(1._wp-zdir)*(diss-buoy)/en(ji,jj,jk) ! dissipation term 260 !!gm better coding, identical results 261 ! zesh2 = p_sh2(ji,jj,jk) + zdir*buoy ! production term 262 ! zdiss = ( diss - (1._wp-zdir)*buoy ) / en(ji,jj,jk) ! dissipation term 263 !!gm 268 264 ! 269 265 ! Compute a wall function from 1. to rsc_psi*zwall/rsc_psi0 … … 281 277 ! building the matrix 282 278 zcof = rfact_tke * tmask(ji,jj,jk) 283 ! 284 ! lower diagonal 285 z_elem_a(ji,jj,jk) = zcof * ( avm (ji,jj,jk ) + avm (ji,jj,jk-1) ) & 286 & / ( e3t_n(ji,jj,jk-1) * e3w_n(ji,jj,jk ) ) 287 ! 288 ! upper diagonal 289 z_elem_c(ji,jj,jk) = zcof * ( avm (ji,jj,jk+1) + avm (ji,jj,jk ) ) & 290 & / ( e3t_n(ji,jj,jk ) * e3w_n(ji,jj,jk) ) 291 ! 292 ! diagonal 293 z_elem_b(ji,jj,jk) = 1._wp - z_elem_a(ji,jj,jk) - z_elem_c(ji,jj,jk) & 294 & + rdt * zdiss * tmask(ji,jj,jk) 295 ! 296 ! right hand side in en 297 en(ji,jj,jk) = en(ji,jj,jk) + rdt * zesh2 * tmask(ji,jj,jk) 279 ! ! lower diagonal 280 zd_lw(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) ) 281 ! ! upper diagonal 282 zd_up(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) ) 283 ! ! diagonal 284 zdiag(ji,jj,jk) = 1._wp - zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) + rdt * zdiss * wmask(ji,jj,jk) 285 ! ! right hand side in en 286 en(ji,jj,jk) = en(ji,jj,jk) + rdt * zesh2 * wmask(ji,jj,jk) 298 287 END DO 299 288 END DO 300 289 END DO 301 290 ! 302 z _elem_b(:,:,jpk) = 1._wp291 zdiag(:,:,jpk) = 1._wp 303 292 ! 304 293 ! Set surface condition on zwall_psi (1 at the bottom) 305 zwall_psi(:,:, 1) = zwall_psi(:,:,2)306 zwall_psi(:,:,jpk) = 1. 294 zwall_psi(:,:, 1 ) = zwall_psi(:,:,2) 295 zwall_psi(:,:,jpk) = 1._wp 307 296 ! 308 297 ! Surface boundary condition on tke … … 311 300 SELECT CASE ( nn_bc_surf ) 312 301 ! 313 CASE ( 0 ) ! Dirichlet case302 CASE ( 0 ) ! Dirichlet boundary condition (set e at k=1 & 2) 314 303 ! First level 315 en(:,:,1) = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1)**(2._wp/3._wp) 316 en(:,:,1) = MAX(en(:,:,1), rn_emin) 317 z_elem_a(:,:,1) = en(:,:,1) 318 z_elem_c(:,:,1) = 0._wp 319 z_elem_b(:,:,1) = 1._wp 304 en (:,:,1) = MAX( rn_emin , rc02r * ustar2_surf(:,:) * (1._wp + rsbc_tke1)**r2_3 ) 305 zd_lw(:,:,1) = en(:,:,1) 306 zd_up(:,:,1) = 0._wp 307 zdiag(:,:,1) = 1._wp 320 308 ! 321 309 ! One level below 322 en(:,:,2) = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1 * ((zhsro(:,:)+gdepw_n(:,:,2)) & 323 & / zhsro(:,:) )**(1.5_wp*ra_sf))**(2._wp/3._wp) 324 en(:,:,2) = MAX(en(:,:,2), rn_emin ) 325 z_elem_a(:,:,2) = 0._wp 326 z_elem_c(:,:,2) = 0._wp 327 z_elem_b(:,:,2) = 1._wp 328 ! 329 ! 330 CASE ( 1 ) ! Neumann boundary condition on d(e)/dz 310 en (:,:,2) = MAX( rc02r * ustar2_surf(:,:) * ( 1._wp + rsbc_tke1 * ((zhsro(:,:)+gdepw_n(:,:,2)) & 311 & / zhsro(:,:) )**(1.5_wp*ra_sf) )**(2._wp/3._wp) , rn_emin ) 312 zd_lw(:,:,2) = 0._wp 313 zd_up(:,:,2) = 0._wp 314 zdiag(:,:,2) = 1._wp 315 ! 316 ! 317 CASE ( 1 ) ! Neumann boundary condition (set d(e)/dz) 331 318 ! 332 319 ! Dirichlet conditions at k=1 333 en(:,:,1) = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1)**(2._wp/3._wp) 334 en(:,:,1) = MAX(en(:,:,1), rn_emin) 335 z_elem_a(:,:,1) = en(:,:,1) 336 z_elem_c(:,:,1) = 0._wp 337 z_elem_b(:,:,1) = 1._wp 320 en (:,:,1) = MAX( rc02r * ustar2_surf(:,:) * (1._wp + rsbc_tke1)**r2_3 , rn_emin ) 321 zd_lw(:,:,1) = en(:,:,1) 322 zd_up(:,:,1) = 0._wp 323 zdiag(:,:,1) = 1._wp 338 324 ! 339 325 ! at k=2, set de/dz=Fw 340 326 !cbr 341 z _elem_b(:,:,2) = z_elem_b(:,:,2) + z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b342 z _elem_a(:,:,2) = 0._wp343 zkar (:,:) = (rl_sf + (vkarmn-rl_sf)*(1.-exp(-rtrans*gdept_n(:,:,1)/zhsro(:,:)) ))344 zflxs(:,:) = rsbc_tke2 * ustars2(:,:)**1.5_wp * zkar(:,:) &345 & * ((zhsro(:,:)+gdept_n(:,:,1)) / zhsro(:,:))**(1.5_wp*ra_sf)346 347 en(:,:,2) = en(:,:,2) + zflxs(:,:) /e3w_n(:,:,2)327 zdiag(:,:,2) = zdiag(:,:,2) + zd_lw(:,:,2) ! Remove zd_lw from zdiag 328 zd_lw(:,:,2) = 0._wp 329 zkar (:,:) = (rl_sf + (vkarmn-rl_sf)*(1.-EXP(-rtrans*gdept_n(:,:,1)/zhsro(:,:)) )) 330 zflxs(:,:) = rsbc_tke2 * ustar2_surf(:,:)**1.5_wp * zkar(:,:) & 331 & * ( ( zhsro(:,:)+gdept_n(:,:,1) ) / zhsro(:,:) )**(1.5_wp*ra_sf) 332 !!gm why not : * ( 1._wp + gdept_n(:,:,1) / zhsro(:,:) )**(1.5_wp*ra_sf) 333 en(:,:,2) = en(:,:,2) + zflxs(:,:) / e3w_n(:,:,2) 348 334 ! 349 335 ! … … 356 342 ! 357 343 CASE ( 0 ) ! Dirichlet 358 ! ! en(ibot) = u*^2 / Co2 and mxln(ibot) = rn_lmin344 ! ! en(ibot) = u*^2 / Co2 and hmxl_n(ibot) = rn_lmin 359 345 ! ! Balance between the production and the dissipation terms 346 DO jj = 2, jpjm1 347 DO ji = fs_2, fs_jpim1 ! vector opt. 348 !!gm This means that bottom and ocean w-level above have a specified "en" value. Sure ???? 349 !! With thick deep ocean level thickness, this may be quite large, no ??? 350 !! in particular in ocean cavities where top stratification can be large... 351 ibot = mbkt(ji,jj) + 1 ! k bottom level of w-point 352 ibotm1 = mbkt(ji,jj) ! k-1 bottom level of w-point but >=1 353 ! 354 z_en = MAX( rc02r * ustar2_bot(ji,jj), rn_emin ) 355 ! 356 ! Dirichlet condition applied at: 357 ! Bottom level (ibot) & Just above it (ibotm1) 358 zd_lw(ji,jj,ibot) = 0._wp ; zd_lw(ji,jj,ibotm1) = 0._wp 359 zd_up(ji,jj,ibot) = 0._wp ; zd_up(ji,jj,ibotm1) = 0._wp 360 zdiag(ji,jj,ibot) = 1._wp ; zdiag(ji,jj,ibotm1) = 1._wp 361 en (ji,jj,ibot) = z_en ; en (ji,jj,ibotm1) = z_en 362 END DO 363 END DO 364 ! 365 IF( ln_isfcav) THEN ! top boundary (ocean cavity) 366 DO jj = 2, jpjm1 367 DO ji = fs_2, fs_jpim1 ! vector opt. 368 itop = mikt(ji,jj) ! k top w-point 369 itopp1 = mikt(ji,jj) + 1 ! k+1 1st w-point below the top one 370 ! ! mask at the ocean surface points 371 z_en = MAX( rc02r * ustar2_top(ji,jj), rn_emin ) * ( 1._wp - tmask(ji,jj,1) ) 372 ! 373 !!gm TO BE VERIFIED !!! 374 ! Dirichlet condition applied at: 375 ! top level (itop) & Just below it (itopp1) 376 zd_lw(ji,jj,itop) = 0._wp ; zd_lw(ji,jj,itopp1) = 0._wp 377 zd_up(ji,jj,itop) = 0._wp ; zd_up(ji,jj,itopp1) = 0._wp 378 zdiag(ji,jj,itop) = 1._wp ; zdiag(ji,jj,itopp1) = 1._wp 379 en (ji,jj,itop) = z_en ; en (ji,jj,itopp1) = z_en 380 END DO 381 END DO 382 ENDIF 383 ! 384 CASE ( 1 ) ! Neumman boundary condition 385 ! 360 386 DO jj = 2, jpjm1 361 387 DO ji = fs_2, fs_jpim1 ! vector opt. … … 363 389 ibotm1 = mbkt(ji,jj) ! k-1 bottom level of w-point but >=1 364 390 ! 391 z_en = MAX( rc02r * ustar2_bot(ji,jj), rn_emin ) 392 ! 365 393 ! Bottom level Dirichlet condition: 366 z_elem_a(ji,jj,ibot ) = 0._wp 367 z_elem_c(ji,jj,ibot ) = 0._wp 368 z_elem_b(ji,jj,ibot ) = 1._wp 369 en(ji,jj,ibot ) = MAX( rc02r * ustarb2(ji,jj), rn_emin ) 370 ! 371 ! Just above last level, Dirichlet condition again 372 z_elem_a(ji,jj,ibotm1) = 0._wp 373 z_elem_c(ji,jj,ibotm1) = 0._wp 374 z_elem_b(ji,jj,ibotm1) = 1._wp 375 en(ji,jj,ibotm1) = MAX( rc02r * ustarb2(ji,jj), rn_emin ) 376 END DO 377 END DO 378 ! 379 CASE ( 1 ) ! Neumman boundary condition 380 ! 381 DO jj = 2, jpjm1 382 DO ji = fs_2, fs_jpim1 ! vector opt. 383 ibot = mbkt(ji,jj) + 1 ! k bottom level of w-point 384 ibotm1 = mbkt(ji,jj) ! k-1 bottom level of w-point but >=1 385 ! 386 ! Bottom level Dirichlet condition: 387 z_elem_a(ji,jj,ibot) = 0._wp 388 z_elem_c(ji,jj,ibot) = 0._wp 389 z_elem_b(ji,jj,ibot) = 1._wp 390 en(ji,jj,ibot) = MAX( rc02r * ustarb2(ji,jj), rn_emin ) 391 ! 392 ! Just above last level: Neumann condition 393 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 394 z_elem_c(ji,jj,ibotm1) = 0._wp 395 END DO 396 END DO 394 ! Bottom level (ibot) & Just above it (ibotm1) 395 ! Dirichlet ! Neumann 396 zd_lw(ji,jj,ibot) = 0._wp ! ! Remove zd_up from zdiag 397 zdiag(ji,jj,ibot) = 1._wp ; zdiag(ji,jj,ibotm1) = zdiag(ji,jj,ibotm1) + zd_up(ji,jj,ibotm1) 398 zd_up(ji,jj,ibot) = 0._wp ; zd_up(ji,jj,ibotm1) = 0._wp 399 END DO 400 END DO 401 IF( ln_isfcav) THEN ! top boundary (ocean cavity) 402 DO jj = 2, jpjm1 403 DO ji = fs_2, fs_jpim1 ! vector opt. 404 itop = mikt(ji,jj) ! k top w-point 405 itopp1 = mikt(ji,jj) + 1 ! k+1 1st w-point below the top one 406 ! ! mask at the ocean surface points 407 z_en = MAX( rc02r * ustar2_top(ji,jj), rn_emin ) * ( 1._wp - tmask(ji,jj,1) ) 408 ! 409 ! Bottom level Dirichlet condition: 410 ! Bottom level (ibot) & Just above it (ibotm1) 411 ! Dirichlet ! Neumann 412 zd_lw(ji,jj,itop) = 0._wp ! ! Remove zd_up from zdiag 413 zdiag(ji,jj,itop) = 1._wp ; zdiag(ji,jj,itopp1) = zdiag(ji,jj,itopp1) + zd_up(ji,jj,itopp1) 414 zd_up(ji,jj,itop) = 0._wp ; zd_up(ji,jj,itopp1) = 0._wp 415 END DO 416 END DO 417 ENDIF 397 418 ! 398 419 END SELECT … … 404 425 DO jj = 2, jpjm1 405 426 DO ji = fs_2, fs_jpim1 ! vector opt. 406 z _elem_b(ji,jj,jk) = z_elem_b(ji,jj,jk) - z_elem_a(ji,jj,jk) * z_elem_c(ji,jj,jk-1) / z_elem_b(ji,jj,jk-1)427 zdiag(ji,jj,jk) = zdiag(ji,jj,jk) - zd_lw(ji,jj,jk) * zd_up(ji,jj,jk-1) / zdiag(ji,jj,jk-1) 407 428 END DO 408 429 END DO … … 411 432 DO jj = 2, jpjm1 412 433 DO ji = fs_2, fs_jpim1 ! vector opt. 413 z _elem_a(ji,jj,jk) = en(ji,jj,jk) - z_elem_a(ji,jj,jk) / z_elem_b(ji,jj,jk-1) * z_elem_a(ji,jj,jk-1)434 zd_lw(ji,jj,jk) = en(ji,jj,jk) - zd_lw(ji,jj,jk) / zdiag(ji,jj,jk-1) * zd_lw(ji,jj,jk-1) 414 435 END DO 415 436 END DO … … 418 439 DO jj = 2, jpjm1 419 440 DO ji = fs_2, fs_jpim1 ! vector opt. 420 en(ji,jj,jk) = ( z _elem_a(ji,jj,jk) - z_elem_c(ji,jj,jk) * en(ji,jj,jk+1) ) / z_elem_b(ji,jj,jk)441 en(ji,jj,jk) = ( zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) * en(ji,jj,jk+1) ) / zdiag(ji,jj,jk) 421 442 END DO 422 443 END DO … … 437 458 DO jj = 2, jpjm1 438 459 DO ji = fs_2, fs_jpim1 ! vector opt. 439 psi(ji,jj,jk) = eb(ji,jj,jk) * mxlb(ji,jj,jk)460 psi(ji,jj,jk) = eb(ji,jj,jk) * hmxl_b(ji,jj,jk) 440 461 END DO 441 462 END DO … … 455 476 DO jj = 2, jpjm1 456 477 DO ji = fs_2, fs_jpim1 ! vector opt. 457 psi(ji,jj,jk) = SQRT( eb(ji,jj,jk) ) / ( rc0 * mxlb(ji,jj,jk) )478 psi(ji,jj,jk) = SQRT( eb(ji,jj,jk) ) / ( rc0 * hmxl_b(ji,jj,jk) ) 458 479 END DO 459 480 END DO … … 464 485 DO jj = 2, jpjm1 465 486 DO ji = fs_2, fs_jpim1 ! vector opt. 466 psi(ji,jj,jk) = rc02 * eb(ji,jj,jk) * mxlb(ji,jj,jk)**rnn487 psi(ji,jj,jk) = rc02 * eb(ji,jj,jk) * hmxl_b(ji,jj,jk)**rnn 467 488 END DO 468 489 END DO … … 475 496 ! Resolution of a tridiagonal linear system by a "methode de chasse" 476 497 ! computation from level 2 to jpkm1 (e(1) already computed and e(jpk)=0 ). 477 ! z _elem_b : diagonal z_elem_c : upper diagonal z_elem_a: lower diagonal498 ! zdiag : diagonal zd_up : upper diagonal zd_lw : lower diagonal 478 499 ! Warning : after this step, en : right hand side of the matrix 479 500 … … 485 506 zratio = psi(ji,jj,jk) / eb(ji,jj,jk) 486 507 ! 487 ! psi3+ : stable : B=-KhN²<0 => N²>0 if rn2>0 dir = 1 (stable) otherwisedir = 0 (unstable)488 dir = 0.5_wp + SIGN( 0.5_wp, rn2(ji,jj,jk) )489 ! 490 rpsi3 = dir * rpsi3m + ( 1._wp -dir ) * rpsi3p508 ! psi3+ : stable : B=-KhN²<0 => N²>0 if rn2>0 zdir = 1 (stable) otherwise zdir = 0 (unstable) 509 zdir = 0.5_wp + SIGN( 0.5_wp, rn2(ji,jj,jk) ) 510 ! 511 rpsi3 = zdir * rpsi3m + ( 1._wp - zdir ) * rpsi3p 491 512 ! 492 513 ! shear prod. - stratif. destruction 493 prod = rpsi1 * zratio * shear(ji,jj,jk)514 prod = rpsi1 * zratio * p_sh2(ji,jj,jk) 494 515 ! 495 516 ! stratif. destruction 496 buoy = rpsi3 * zratio * (- avt(ji,jj,jk) * rn2(ji,jj,jk) )517 buoy = rpsi3 * zratio * (- p_avt(ji,jj,jk) * rn2(ji,jj,jk) ) 497 518 ! 498 519 ! shear prod. - stratif. destruction 499 520 diss = rpsi2 * zratio * zwall(ji,jj,jk) * eps(ji,jj,jk) 500 521 ! 501 dir = 0.5_wp + SIGN( 0.5_wp, prod + buoy ) !dir =1(=0) if shear(ji,jj,jk)+buoy >0(<0)502 ! 503 zesh2 = dir * ( prod + buoy ) + (1._wp -dir ) * prod ! production term504 zdiss = dir * ( diss / psi(ji,jj,jk) ) + (1._wp -dir ) * (diss-buoy) / psi(ji,jj,jk) ! dissipation term522 zdir = 0.5_wp + SIGN( 0.5_wp, prod + buoy ) ! zdir =1(=0) if shear(ji,jj,jk)+buoy >0(<0) 523 ! 524 zesh2 = zdir * ( prod + buoy ) + (1._wp - zdir ) * prod ! production term 525 zdiss = zdir * ( diss / psi(ji,jj,jk) ) + (1._wp - zdir ) * (diss-buoy) / psi(ji,jj,jk) ! dissipation term 505 526 ! 506 527 ! building the matrix 507 528 zcof = rfact_psi * zwall_psi(ji,jj,jk) * tmask(ji,jj,jk) 508 ! lower diagonal 509 z_elem_a(ji,jj,jk) = zcof * ( avm (ji,jj,jk ) + avm (ji,jj,jk-1) ) & 510 & / ( e3t_n(ji,jj,jk-1) * e3w_n(ji,jj,jk ) ) 511 ! upper diagonal 512 z_elem_c(ji,jj,jk) = zcof * ( avm (ji,jj,jk+1) + avm (ji,jj,jk ) ) & 513 & / ( e3t_n(ji,jj,jk ) * e3w_n(ji,jj,jk) ) 514 ! diagonal 515 z_elem_b(ji,jj,jk) = 1._wp - z_elem_a(ji,jj,jk) - z_elem_c(ji,jj,jk) & 516 & + rdt * zdiss * tmask(ji,jj,jk) 517 ! 518 ! right hand side in psi 519 psi(ji,jj,jk) = psi(ji,jj,jk) + rdt * zesh2 * tmask(ji,jj,jk) 529 ! ! lower diagonal 530 zd_lw(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) ) 531 ! ! upper diagonal 532 zd_up(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) ) 533 ! ! diagonal 534 zdiag(ji,jj,jk) = 1._wp - zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) + rdt * zdiss * wmask(ji,jj,jk) 535 ! ! right hand side in psi 536 psi(ji,jj,jk) = psi(ji,jj,jk) + rdt * zesh2 * wmask(ji,jj,jk) 520 537 END DO 521 538 END DO 522 539 END DO 523 540 ! 524 z _elem_b(:,:,jpk) = 1._wp541 zdiag(:,:,jpk) = 1._wp 525 542 526 543 ! Surface boundary condition on psi … … 530 547 ! 531 548 CASE ( 0 ) ! Dirichlet boundary conditions 532 ! 533 ! Surface value 534 zdep(:,:) = zhsro(:,:) * rl_sf ! Cosmetic 535 psi (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 536 z_elem_a(:,:,1) = psi(:,:,1) 537 z_elem_c(:,:,1) = 0._wp 538 z_elem_b(:,:,1) = 1._wp 539 ! 540 ! One level below 541 zkar(:,:) = (rl_sf + (vkarmn-rl_sf)*(1._wp-exp(-rtrans*gdepw_n(:,:,2)/zhsro(:,:) ))) 542 zdep(:,:) = (zhsro(:,:) + gdepw_n(:,:,2)) * zkar(:,:) 543 psi (:,:,2) = rc0**rpp * en(:,:,2)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 544 z_elem_a(:,:,2) = 0._wp 545 z_elem_c(:,:,2) = 0._wp 546 z_elem_b(:,:,2) = 1._wp 547 ! 548 ! 549 ! 550 ! Surface value 551 zdep (:,:) = zhsro(:,:) * rl_sf ! Cosmetic 552 psi (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 553 zd_lw(:,:,1) = psi(:,:,1) 554 zd_up(:,:,1) = 0._wp 555 zdiag(:,:,1) = 1._wp 556 ! 557 ! One level below 558 zkar (:,:) = (rl_sf + (vkarmn-rl_sf)*(1._wp-EXP(-rtrans*gdepw_n(:,:,2)/zhsro(:,:) ))) 559 zdep (:,:) = (zhsro(:,:) + gdepw_n(:,:,2)) * zkar(:,:) 560 psi (:,:,2) = rc0**rpp * en(:,:,2)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 561 zd_lw(:,:,2) = 0._wp 562 zd_up(:,:,2) = 0._wp 563 zdiag(:,:,2) = 1._wp 564 ! 549 565 CASE ( 1 ) ! Neumann boundary condition on d(psi)/dz 550 ! 551 ! Surface value: Dirichlet 552 zdep(:,:) = zhsro(:,:) * rl_sf 553 psi (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 554 z_elem_a(:,:,1) = psi(:,:,1) 555 z_elem_c(:,:,1) = 0._wp 556 z_elem_b(:,:,1) = 1._wp 557 ! 558 ! Neumann condition at k=2 559 z_elem_b(:,:,2) = z_elem_b(:,:,2) + z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b 560 z_elem_a(:,:,2) = 0._wp 561 ! 562 ! Set psi vertical flux at the surface: 563 zkar(:,:) = rl_sf + (vkarmn-rl_sf)*(1._wp-exp(-rtrans*gdept_n(:,:,1)/zhsro(:,:) )) ! Lengh scale slope 564 zdep(:,:) = ((zhsro(:,:) + gdept_n(:,:,1)) / zhsro(:,:))**(rmm*ra_sf) 565 zflxs(:,:) = (rnn + rsbc_tke1 * (rnn + rmm*ra_sf) * zdep(:,:))*(1._wp + rsbc_tke1*zdep(:,:))**(2._wp*rmm/3._wp-1_wp) 566 zdep(:,:) = rsbc_psi1 * (zwall_psi(:,:,1)*avm(:,:,1)+zwall_psi(:,:,2)*avm(:,:,2)) * & 567 & ustars2(:,:)**rmm * zkar(:,:)**rnn * (zhsro(:,:) + gdept_n(:,:,1))**(rnn-1.) 568 zflxs(:,:) = zdep(:,:) * zflxs(:,:) 569 psi(:,:,2) = psi(:,:,2) + zflxs(:,:) / e3w_n(:,:,2) 570 571 ! 572 ! 566 ! 567 ! Surface value: Dirichlet 568 zdep (:,:) = zhsro(:,:) * rl_sf 569 psi (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 570 zd_lw(:,:,1) = psi(:,:,1) 571 zd_up(:,:,1) = 0._wp 572 zdiag(:,:,1) = 1._wp 573 ! 574 ! Neumann condition at k=2 575 zdiag(:,:,2) = zdiag(:,:,2) + zd_lw(:,:,2) ! Remove zd_lw from zdiag 576 zd_lw(:,:,2) = 0._wp 577 ! 578 ! Set psi vertical flux at the surface: 579 zkar (:,:) = rl_sf + (vkarmn-rl_sf)*(1._wp-EXP(-rtrans*gdept_n(:,:,1)/zhsro(:,:) )) ! Lengh scale slope 580 zdep (:,:) = ((zhsro(:,:) + gdept_n(:,:,1)) / zhsro(:,:))**(rmm*ra_sf) 581 zflxs(:,:) = (rnn + rsbc_tke1 * (rnn + rmm*ra_sf) * zdep(:,:))*(1._wp + rsbc_tke1*zdep(:,:))**(2._wp*rmm/3._wp-1_wp) 582 zdep (:,:) = rsbc_psi1 * (zwall_psi(:,:,1)*avm(:,:,1)+zwall_psi(:,:,2)*avm(:,:,2)) * & 583 & ustar2_surf(:,:)**rmm * zkar(:,:)**rnn * (zhsro(:,:) + gdept_n(:,:,1))**(rnn-1.) 584 zflxs(:,:) = zdep(:,:) * zflxs(:,:) 585 psi (:,:,2) = psi(:,:,2) + zflxs(:,:) / e3w_n(:,:,2) 586 ! 573 587 END SELECT 574 588 … … 576 590 ! -------------------------------- 577 591 ! 578 SELECT CASE ( nn_bc_bot ) 579 ! 592 !!gm should be done for ISF (top boundary cond.) 593 !!gm so, totally new staff needed ===>>> think about that ! 594 ! 595 SELECT CASE ( nn_bc_bot ) ! bottom boundary 580 596 ! 581 597 CASE ( 0 ) ! Dirichlet 582 ! ! en(ibot) = u*^2 / Co2 and mxln(ibot) = vkarmn * rn_bfrz0598 ! ! en(ibot) = u*^2 / Co2 and hmxl_n(ibot) = vkarmn * r_z0_bot 583 599 ! ! Balance between the production and the dissipation terms 584 600 DO jj = 2, jpjm1 … … 586 602 ibot = mbkt(ji,jj) + 1 ! k bottom level of w-point 587 603 ibotm1 = mbkt(ji,jj) ! k-1 bottom level of w-point but >=1 588 zdep(ji,jj) = vkarmn * r n_bfrz0604 zdep(ji,jj) = vkarmn * r_z0_bot 589 605 psi (ji,jj,ibot) = rc0**rpp * en(ji,jj,ibot)**rmm * zdep(ji,jj)**rnn 590 z _elem_a(ji,jj,ibot) = 0._wp591 z _elem_c(ji,jj,ibot) = 0._wp592 z _elem_b(ji,jj,ibot) = 1._wp606 zd_lw(ji,jj,ibot) = 0._wp 607 zd_up(ji,jj,ibot) = 0._wp 608 zdiag(ji,jj,ibot) = 1._wp 593 609 ! 594 610 ! Just above last level, Dirichlet condition again (GOTM like) 595 zdep(ji,jj) = vkarmn * ( r n_bfrz0+ e3t_n(ji,jj,ibotm1) )611 zdep(ji,jj) = vkarmn * ( r_z0_bot + e3t_n(ji,jj,ibotm1) ) 596 612 psi (ji,jj,ibotm1) = rc0**rpp * en(ji,jj,ibot )**rmm * zdep(ji,jj)**rnn 597 z _elem_a(ji,jj,ibotm1) = 0._wp598 z _elem_c(ji,jj,ibotm1) = 0._wp599 z _elem_b(ji,jj,ibotm1) = 1._wp613 zd_lw(ji,jj,ibotm1) = 0._wp 614 zd_up(ji,jj,ibotm1) = 0._wp 615 zdiag(ji,jj,ibotm1) = 1._wp 600 616 END DO 601 617 END DO … … 609 625 ! 610 626 ! Bottom level Dirichlet condition: 611 zdep(ji,jj) = vkarmn * r n_bfrz0627 zdep(ji,jj) = vkarmn * r_z0_bot 612 628 psi (ji,jj,ibot) = rc0**rpp * en(ji,jj,ibot)**rmm * zdep(ji,jj)**rnn 613 629 ! 614 z _elem_a(ji,jj,ibot) = 0._wp615 z _elem_c(ji,jj,ibot) = 0._wp616 z _elem_b(ji,jj,ibot) = 1._wp630 zd_lw(ji,jj,ibot) = 0._wp 631 zd_up(ji,jj,ibot) = 0._wp 632 zdiag(ji,jj,ibot) = 1._wp 617 633 ! 618 634 ! Just above last level: Neumann condition with flux injection 619 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_b620 z _elem_c(ji,jj,ibotm1) = 0.635 zdiag(ji,jj,ibotm1) = zdiag(ji,jj,ibotm1) + zd_up(ji,jj,ibotm1) ! Remove zd_up from zdiag 636 zd_up(ji,jj,ibotm1) = 0. 621 637 ! 622 638 ! Set psi vertical flux at the bottom: 623 zdep(ji,jj) = r n_bfrz0+ 0.5_wp*e3t_n(ji,jj,ibotm1)624 zflxb = rsbc_psi2 * ( avm(ji,jj,ibot) +avm(ji,jj,ibotm1) ) &639 zdep(ji,jj) = r_z0_bot + 0.5_wp*e3t_n(ji,jj,ibotm1) 640 zflxb = rsbc_psi2 * ( p_avm(ji,jj,ibot) + p_avm(ji,jj,ibotm1) ) & 625 641 & * (0.5_wp*(en(ji,jj,ibot)+en(ji,jj,ibotm1)))**rmm * zdep(ji,jj)**(rnn-1._wp) 626 642 psi(ji,jj,ibotm1) = psi(ji,jj,ibotm1) + zflxb / e3w_n(ji,jj,ibotm1) … … 636 652 DO jj = 2, jpjm1 637 653 DO ji = fs_2, fs_jpim1 ! vector opt. 638 z _elem_b(ji,jj,jk) = z_elem_b(ji,jj,jk) - z_elem_a(ji,jj,jk) * z_elem_c(ji,jj,jk-1) / z_elem_b(ji,jj,jk-1)654 zdiag(ji,jj,jk) = zdiag(ji,jj,jk) - zd_lw(ji,jj,jk) * zd_up(ji,jj,jk-1) / zdiag(ji,jj,jk-1) 639 655 END DO 640 656 END DO … … 643 659 DO jj = 2, jpjm1 644 660 DO ji = fs_2, fs_jpim1 ! vector opt. 645 z _elem_a(ji,jj,jk) = psi(ji,jj,jk) - z_elem_a(ji,jj,jk) / z_elem_b(ji,jj,jk-1) * z_elem_a(ji,jj,jk-1)661 zd_lw(ji,jj,jk) = psi(ji,jj,jk) - zd_lw(ji,jj,jk) / zdiag(ji,jj,jk-1) * zd_lw(ji,jj,jk-1) 646 662 END DO 647 663 END DO … … 650 666 DO jj = 2, jpjm1 651 667 DO ji = fs_2, fs_jpim1 ! vector opt. 652 psi(ji,jj,jk) = ( z _elem_a(ji,jj,jk) - z_elem_c(ji,jj,jk) * psi(ji,jj,jk+1) ) / z_elem_b(ji,jj,jk)668 psi(ji,jj,jk) = ( zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) * psi(ji,jj,jk+1) ) / zdiag(ji,jj,jk) 653 669 END DO 654 670 END DO … … 703 719 ! Limit dissipation rate under stable stratification 704 720 ! -------------------------------------------------- 705 DO jk = 1, jpkm1 ! Note that this set boundary conditions on mxln at the same time721 DO jk = 1, jpkm1 ! Note that this set boundary conditions on hmxl_n at the same time 706 722 DO jj = 2, jpjm1 707 723 DO ji = fs_2, fs_jpim1 ! vector opt. 708 724 ! limitation 709 eps (ji,jj,jk) = MAX( eps(ji,jj,jk), rn_epsmin )710 mxln(ji,jj,jk) = rc03 * en(ji,jj,jk) * SQRT( en(ji,jj,jk) ) / eps(ji,jj,jk)725 eps (ji,jj,jk) = MAX( eps(ji,jj,jk), rn_epsmin ) 726 hmxl_n(ji,jj,jk) = rc03 * en(ji,jj,jk) * SQRT( en(ji,jj,jk) ) / eps(ji,jj,jk) 711 727 ! Galperin criterium (NOTE : Not required if the proper value of C3 in stable cases is calculated) 712 728 zrn2 = MAX( rn2(ji,jj,jk), rsmall ) 713 IF (ln_length_lim) mxln(ji,jj,jk) = MIN( rn_clim_galp * SQRT( 2._wp * en(ji,jj,jk) / zrn2 ), mxln(ji,jj,jk) )729 IF( ln_length_lim ) hmxl_n(ji,jj,jk) = MIN( rn_clim_galp * SQRT( 2._wp * en(ji,jj,jk) / zrn2 ), hmxl_n(ji,jj,jk) ) 714 730 END DO 715 731 END DO … … 727 743 DO ji = fs_2, fs_jpim1 ! vector opt. 728 744 ! zcof = l²/q² 729 zcof = mxlb(ji,jj,jk) * mxlb(ji,jj,jk) / ( 2._wp*eb(ji,jj,jk) )745 zcof = hmxl_b(ji,jj,jk) * hmxl_b(ji,jj,jk) / ( 2._wp*eb(ji,jj,jk) ) 730 746 ! Gh = -N²l²/q² 731 747 gh = - rn2(ji,jj,jk) * zcof … … 736 752 sm = ( rb1**(-1._wp/3._wp) + ( 18._wp*ra1*ra1 + 9._wp*ra1*ra2*(1._wp-rc2) )*sh*gh ) / (1._wp-9._wp*ra1*ra2*gh) 737 753 ! 738 ! Store stability function in avmu and avmv739 avmu(ji,jj,jk) = rc_diff * sh * tmask(ji,jj,jk)740 avmv(ji,jj,jk) = rc_diff * sm * tmask(ji,jj,jk)754 ! Store stability function in zstt and zstm 755 zstt(ji,jj,jk) = rc_diff * sh * tmask(ji,jj,jk) 756 zstm(ji,jj,jk) = rc_diff * sm * tmask(ji,jj,jk) 741 757 END DO 742 758 END DO … … 748 764 DO ji = fs_2, fs_jpim1 ! vector opt. 749 765 ! zcof = l²/q² 750 zcof = mxlb(ji,jj,jk)*mxlb(ji,jj,jk) / ( 2._wp * eb(ji,jj,jk) )766 zcof = hmxl_b(ji,jj,jk)*hmxl_b(ji,jj,jk) / ( 2._wp * eb(ji,jj,jk) ) 751 767 ! Gh = -N²l²/q² 752 768 gh = - rn2(ji,jj,jk) * zcof … … 755 771 gh = gh * rf6 756 772 ! Gm = M²l²/q² Shear number 757 shr = shear(ji,jj,jk) / MAX(avm(ji,jj,jk), rsmall )773 shr = p_sh2(ji,jj,jk) / MAX( p_avm(ji,jj,jk), rsmall ) 758 774 gm = MAX( shr * zcof , 1.e-10 ) 759 775 gm = gm * rf6 … … 764 780 sh = (rs4 - rs5*gh + rs6*gm) / rcff 765 781 ! 766 ! Store stability function in avmu and avmv767 avmu(ji,jj,jk) = rc_diff * sh * tmask(ji,jj,jk)768 avmv(ji,jj,jk) = rc_diff * sm * tmask(ji,jj,jk)782 ! Store stability function in zstt and zstm 783 zstt(ji,jj,jk) = rc_diff * sh * tmask(ji,jj,jk) 784 zstm(ji,jj,jk) = rc_diff * sm * tmask(ji,jj,jk) 769 785 END DO 770 786 END DO … … 776 792 ! Lines below are useless if GOTM style Dirichlet conditions are used 777 793 778 avmv(:,:,1) = avmv(:,:,2)794 zstm(:,:,1) = zstm(:,:,2) 779 795 780 796 DO jj = 2, jpjm1 781 797 DO ji = fs_2, fs_jpim1 ! vector opt. 782 avmv(ji,jj,mbkt(ji,jj)+1) = avmv(ji,jj,mbkt(ji,jj))798 zstm(ji,jj,mbkt(ji,jj)+1) = zstm(ji,jj,mbkt(ji,jj)) 783 799 END DO 784 800 END DO 801 !!gm should be done for ISF (top boundary cond.) 802 !!gm so, totally new staff needed!!gm 785 803 786 804 ! Compute diffusivities/viscosities … … 789 807 DO jj = 2, jpjm1 790 808 DO ji = fs_2, fs_jpim1 ! vector opt. 791 zsqen = SQRT( 2._wp * en(ji,jj,jk) ) * mxln(ji,jj,jk)792 zav = zsqen * avmu(ji,jj,jk)793 avt(ji,jj,jk) = MAX( zav, avtb(jk) )*tmask(ji,jj,jk) ! apply mask for zdfmxl routine794 zav = zsqen * avmv(ji,jj,jk)795 avm(ji,jj,jk) = MAX( zav, avmb(jk) )! Note that avm is not masked at the surface and the bottom809 zsqen = SQRT( 2._wp * en(ji,jj,jk) ) * hmxl_n(ji,jj,jk) 810 zavt = zsqen * zstt(ji,jj,jk) 811 zavm = zsqen * zstm(ji,jj,jk) 812 p_avt(ji,jj,jk) = MAX( zavt, avtb(jk) ) * wmask(ji,jj,jk) ! apply mask for zdfmxl routine 813 p_avm(ji,jj,jk) = MAX( zavm, avmb(jk) ) ! Note that avm is not masked at the surface and the bottom 796 814 END DO 797 815 END DO 798 816 END DO 799 !800 ! Lateral boundary conditions (sign unchanged)801 817 avt(:,:,1) = 0._wp 802 CALL lbc_lnk( avm, 'W', 1. ) ; CALL lbc_lnk( avt, 'W', 1. ) 803 804 DO jk = 2, jpkm1 !* vertical eddy viscosity at u- and v-points 805 DO jj = 2, jpjm1 806 DO ji = fs_2, fs_jpim1 ! vector opt. 807 avmu(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji+1,jj ,jk) ) * umask(ji,jj,jk) 808 avmv(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji ,jj+1,jk) ) * vmask(ji,jj,jk) 809 END DO 810 END DO 811 END DO 812 avmu(:,:,1) = 0._wp ; avmv(:,:,1) = 0._wp ! set surface to zero 813 CALL lbc_lnk( avmu, 'U', 1. ) ; CALL lbc_lnk( avmv, 'V', 1. ) ! Lateral boundary conditions 814 818 ! 815 819 IF(ln_ctl) THEN 816 CALL prt_ctl( tab3d_1=en , clinfo1=' gls - e: ', tab3d_2=avt, clinfo2=' t: ', ovlap=1, kdim=jpk) 817 CALL prt_ctl( tab3d_1=avmu, clinfo1=' gls - u: ', mask1=umask, & 818 & tab3d_2=avmv, clinfo2= ' v: ', mask2=vmask, ovlap=1, kdim=jpk ) 820 CALL prt_ctl( tab3d_1=en , clinfo1=' gls - e: ', tab3d_2=avt, clinfo2=' t: ', ovlap=1, kdim=jpk) 821 CALL prt_ctl( tab3d_1=avm, clinfo1=' gls - m: ', ovlap=1, kdim=jpk ) 819 822 ENDIF 820 823 ! 821 avt_k (:,:,:) = avt (:,:,:) 822 avm_k (:,:,:) = avm (:,:,:) 823 avmu_k(:,:,:) = avmu(:,:,:) 824 avmv_k(:,:,:) = avmv(:,:,:) 825 ! 826 CALL wrk_dealloc( jpi,jpj, zdep, zkar, zflxs, zhsro ) 827 CALL wrk_dealloc( jpi,jpj,jpk, eb, mxlb, shear, eps, zwall_psi, z_elem_a, z_elem_b, z_elem_c, psi ) 828 ! 829 IF( nn_timing == 1 ) CALL timing_stop('zdf_gls') 830 ! 824 IF( ln_timing ) CALL timing_stop('zdf_gls') 831 825 ! 832 826 END SUBROUTINE zdf_gls … … 838 832 !! 839 833 !! ** Purpose : Initialization of the vertical eddy diffivity and 840 !! viscosity when using a glsturbulent closure scheme834 !! viscosity computed using a GLS turbulent closure scheme 841 835 !! 842 836 !! ** Method : Read the namzdf_gls namelist and check the parameters 843 !! called at the first timestep (nit000)844 837 !! 845 838 !! ** input : Namlist namzdf_gls … … 848 841 !! 849 842 !!---------------------------------------------------------------------- 850 USE dynzdf_exp851 USE trazdf_exp852 !853 843 INTEGER :: jk ! dummy loop indices 854 844 INTEGER :: ios ! Local integer output status for namelist read … … 862 852 !!---------------------------------------------------------- 863 853 ! 864 IF( nn_timing == 1 )CALL timing_start('zdf_gls_init')854 IF( ln_timing ) CALL timing_start('zdf_gls_init') 865 855 ! 866 856 REWIND( numnam_ref ) ! Namelist namzdf_gls in reference namelist : Vertical eddy diffivity and viscosity using gls turbulent closure scheme … … 875 865 IF(lwp) THEN !* Control print 876 866 WRITE(numout,*) 877 WRITE(numout,*) 'zdf_gls_init : glsturbulent closure scheme'867 WRITE(numout,*) 'zdf_gls_init : GLS turbulent closure scheme' 878 868 WRITE(numout,*) '~~~~~~~~~~~~' 879 869 WRITE(numout,*) ' Namelist namzdf_gls : set gls mixing parameters' … … 892 882 WRITE(numout,*) ' Type of closure nn_clos = ', nn_clos 893 883 WRITE(numout,*) ' Surface roughness (m) rn_hsro = ', rn_hsro 894 WRITE(numout,*) ' Bottom roughness (m) (nambfr namelist) rn_bfrz0 = ', rn_bfrz0 884 WRITE(numout,*) 885 WRITE(numout,*) ' Namelist namdrg_top/_bot: used values:' 886 WRITE(numout,*) ' top ocean cavity roughness (m) rn_z0(_top) = ', r_z0_top 887 WRITE(numout,*) ' Bottom seafloor roughness (m) rn_z0(_bot) = ', r_z0_bot 888 WRITE(numout,*) 895 889 ENDIF 896 890 897 ! !* allocate glsarrays891 ! !* allocate GLS arrays 898 892 IF( zdf_gls_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'zdf_gls_init : unable to allocate arrays' ) 899 893 900 894 ! !* Check of some namelist values 901 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' )902 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' )903 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' )904 IF( nn_z0_met == 3 .AND. .NOT.ln_wave )CALL ctl_stop( 'zdf_gls_init: nn_z0_met=3 requires ln_wave=T' )905 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' )906 IF( nn_clos < 0 .OR. nn_clos > 3 )CALL ctl_stop( 'zdf_gls_init: bad flag: nn_clos is 0, 1, 2 or 3' )895 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' ) 896 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' ) 897 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' ) 898 IF( nn_z0_met == 3 .AND. .NOT.ln_wave ) CALL ctl_stop( 'zdf_gls_init: nn_z0_met=3 requires ln_wave=T' ) 899 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' ) 900 IF( nn_clos < 0 .OR. nn_clos > 3 ) CALL ctl_stop( 'zdf_gls_init: bad flag: nn_clos is 0, 1, 2 or 3' ) 907 901 908 902 SELECT CASE ( nn_clos ) !* set the parameters for the chosen closure … … 910 904 CASE( 0 ) ! k-kl (Mellor-Yamada) 911 905 ! 912 IF(lwp) WRITE(numout,*) 'The choosen closure is k-kl closed to the classical Mellor-Yamada' 906 IF(lwp) WRITE(numout,*) ' ==>> k-kl closure chosen (i.e. closed to the classical Mellor-Yamada)' 907 IF(lwp) WRITE(numout,*) 913 908 rpp = 0._wp 914 909 rmm = 1._wp … … 928 923 CASE( 1 ) ! k-eps 929 924 ! 930 IF(lwp) WRITE(numout,*) 'The choosen closure is k-eps' 925 IF(lwp) WRITE(numout,*) ' ==>> k-eps closure chosen' 926 IF(lwp) WRITE(numout,*) 931 927 rpp = 3._wp 932 928 rmm = 1.5_wp … … 946 942 CASE( 2 ) ! k-omega 947 943 ! 948 IF(lwp) WRITE(numout,*) 'The choosen closure is k-omega' 944 IF(lwp) WRITE(numout,*) ' ==>> k-omega closure chosen' 945 IF(lwp) WRITE(numout,*) 949 946 rpp = -1._wp 950 947 rmm = 0.5_wp … … 964 961 CASE( 3 ) ! generic 965 962 ! 966 IF(lwp) WRITE(numout,*) 'The choosen closure is generic' 963 IF(lwp) WRITE(numout,*) ' ==>> generic closure chosen' 964 IF(lwp) WRITE(numout,*) 967 965 rpp = 2._wp 968 966 rmm = 1._wp … … 987 985 CASE ( 0 ) ! Galperin stability functions 988 986 ! 989 IF(lwp) WRITE(numout,*) ' Stability functions from Galperin'987 IF(lwp) WRITE(numout,*) ' ==>> Stability functions from Galperin' 990 988 rc2 = 0._wp 991 989 rc3 = 0._wp … … 999 997 CASE ( 1 ) ! Kantha-Clayson stability functions 1000 998 ! 1001 IF(lwp) WRITE(numout,*) ' Stability functions from Kantha-Clayson'999 IF(lwp) WRITE(numout,*) ' ==>> Stability functions from Kantha-Clayson' 1002 1000 rc2 = 0.7_wp 1003 1001 rc3 = 0.2_wp … … 1011 1009 CASE ( 2 ) ! Canuto A stability functions 1012 1010 ! 1013 IF(lwp) WRITE(numout,*) ' Stability functions from Canuto A'1011 IF(lwp) WRITE(numout,*) ' ==>> Stability functions from Canuto A' 1014 1012 rs0 = 1.5_wp * rl1 * rl5*rl5 1015 1013 rs1 = -rl4*(rl6+rl7) + 2._wp*rl4*rl5*(rl1-(1._wp/3._wp)*rl2-rl3) + 1.5_wp*rl1*rl5*rl8 … … 1035 1033 CASE ( 3 ) ! Canuto B stability functions 1036 1034 ! 1037 IF(lwp) WRITE(numout,*) ' Stability functions from Canuto B'1035 IF(lwp) WRITE(numout,*) ' ==>> Stability functions from Canuto B' 1038 1036 rs0 = 1.5_wp * rm1 * rm5*rm5 1039 1037 rs1 = -rm4 * (rm6+rm7) + 2._wp * rm4*rm5*(rm1-(1._wp/3._wp)*rm2-rm3) + 1.5_wp * rm1*rm5*rm8 … … 1079 1077 rl_sf = vkarmn 1080 1078 ELSE 1081 rl_sf = rc0 * SQRT(rc0/rcm_sf) * SQRT( ( (1._wp + 4._wp*rmm + 8._wp*rmm**2_wp)*rsc_tke & 1082 & + 12._wp * rsc_psi0*rpsi2 - (1._wp + 4._wp*rmm) & 1083 & *SQRT(rsc_tke*(rsc_tke & 1084 & + 24._wp*rsc_psi0*rpsi2)) ) & 1085 & /(12._wp*rnn**2.) & 1086 & ) 1079 rl_sf = rc0 * SQRT(rc0/rcm_sf) * SQRT( ( (1._wp + 4._wp*rmm + 8._wp*rmm**2_wp) * rsc_tke & 1080 & + 12._wp*rsc_psi0*rpsi2 - (1._wp + 4._wp*rmm) & 1081 & *SQRT(rsc_tke*(rsc_tke & 1082 & + 24._wp*rsc_psi0*rpsi2)) ) & 1083 & /(12._wp*rnn**2.) ) 1087 1084 ENDIF 1088 1085 … … 1090 1087 IF(lwp) THEN !* Control print 1091 1088 WRITE(numout,*) 1092 WRITE(numout,*) 'Limit values' 1093 WRITE(numout,*) '~~~~~~~~~~~~' 1094 WRITE(numout,*) 'Parameter m = ',rmm 1095 WRITE(numout,*) 'Parameter n = ',rnn 1096 WRITE(numout,*) 'Parameter p = ',rpp 1097 WRITE(numout,*) 'rpsi1 = ',rpsi1 1098 WRITE(numout,*) 'rpsi2 = ',rpsi2 1099 WRITE(numout,*) 'rpsi3m = ',rpsi3m 1100 WRITE(numout,*) 'rpsi3p = ',rpsi3p 1101 WRITE(numout,*) 'rsc_tke = ',rsc_tke 1102 WRITE(numout,*) 'rsc_psi = ',rsc_psi 1103 WRITE(numout,*) 'rsc_psi0 = ',rsc_psi0 1104 WRITE(numout,*) 'rc0 = ',rc0 1089 WRITE(numout,*) ' Limit values :' 1090 WRITE(numout,*) ' Parameter m = ', rmm 1091 WRITE(numout,*) ' Parameter n = ', rnn 1092 WRITE(numout,*) ' Parameter p = ', rpp 1093 WRITE(numout,*) ' rpsi1 = ', rpsi1 1094 WRITE(numout,*) ' rpsi2 = ', rpsi2 1095 WRITE(numout,*) ' rpsi3m = ', rpsi3m 1096 WRITE(numout,*) ' rpsi3p = ', rpsi3p 1097 WRITE(numout,*) ' rsc_tke = ', rsc_tke 1098 WRITE(numout,*) ' rsc_psi = ', rsc_psi 1099 WRITE(numout,*) ' rsc_psi0 = ', rsc_psi0 1100 WRITE(numout,*) ' rc0 = ', rc0 1105 1101 WRITE(numout,*) 1106 WRITE(numout,*) 'Shear free turbulence parameters:' 1107 WRITE(numout,*) 'rcm_sf = ',rcm_sf 1108 WRITE(numout,*) 'ra_sf = ',ra_sf 1109 WRITE(numout,*) 'rl_sf = ',rl_sf 1110 WRITE(numout,*) 1102 WRITE(numout,*) ' Shear free turbulence parameters:' 1103 WRITE(numout,*) ' rcm_sf = ', rcm_sf 1104 WRITE(numout,*) ' ra_sf = ', ra_sf 1105 WRITE(numout,*) ' rl_sf = ', rl_sf 1111 1106 ENDIF 1112 1107 … … 1123 1118 rsbc_psi1 = -0.5_wp * rdt * rc0**(rpp-2._wp*rmm) / rsc_psi 1124 1119 rsbc_psi2 = -0.5_wp * rdt * rc0**rpp * rnn * vkarmn**rnn / rsc_psi ! Neumann + NO Wave breaking 1125 1120 ! 1126 1121 rfact_tke = -0.5_wp / rsc_tke * rdt ! Cst used for the Diffusion term of tke 1127 1122 rfact_psi = -0.5_wp / rsc_psi * rdt ! Cst used for the Diffusion term of tke 1128 1123 ! 1129 1124 ! !* Wall proximity function 1130 zwall (:,:,:) = 1._wp * tmask(:,:,:) 1131 1132 ! !* set vertical eddy coef. to the background value 1133 DO jk = 1, jpk 1134 avt (:,:,jk) = avtb(jk) * tmask(:,:,jk) 1135 avm (:,:,jk) = avmb(jk) * tmask(:,:,jk) 1136 avmu(:,:,jk) = avmb(jk) * umask(:,:,jk) 1137 avmv(:,:,jk) = avmb(jk) * vmask(:,:,jk) 1138 END DO 1139 ! 1140 CALL gls_rst( nit000, 'READ' ) !* read or initialize all required files 1141 ! 1142 IF( nn_timing == 1 ) CALL timing_stop('zdf_gls_init') 1125 !!gm tmask or wmask ???? 1126 zwall(:,:,:) = 1._wp * tmask(:,:,:) 1127 1128 ! !* read or initialize all required files 1129 CALL gls_rst( nit000, 'READ' ) ! (en, avt_k, avm_k, hmxl_n) 1130 ! 1131 IF( ln_timing ) CALL timing_stop('zdf_gls_init') 1143 1132 ! 1144 1133 END SUBROUTINE zdf_gls_init … … 1155 1144 !! set to rn_emin or recomputed (nn_igls/=0) 1156 1145 !!---------------------------------------------------------------------- 1157 INTEGER , INTENT(in) :: kt 1158 CHARACTER(len=*), INTENT(in) :: cdrw 1146 INTEGER , INTENT(in) :: kt ! ocean time-step 1147 CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag 1159 1148 ! 1160 1149 INTEGER :: jit, jk ! dummy loop indices 1161 INTEGER :: id1, id2, id3, id4 , id5, id61150 INTEGER :: id1, id2, id3, id4 1162 1151 INTEGER :: ji, jj, ikbu, ikbv 1163 1152 REAL(wp):: cbx, cby … … 1167 1156 ! ! --------------- 1168 1157 IF( ln_rstart ) THEN !* Read the restart file 1169 id1 = iom_varid( numror, 'en' , ldstop = .FALSE. ) 1170 id2 = iom_varid( numror, 'avt' , ldstop = .FALSE. ) 1171 id3 = iom_varid( numror, 'avm' , ldstop = .FALSE. ) 1172 id4 = iom_varid( numror, 'avmu' , ldstop = .FALSE. ) 1173 id5 = iom_varid( numror, 'avmv' , ldstop = .FALSE. ) 1174 id6 = iom_varid( numror, 'mxln' , ldstop = .FALSE. ) 1158 id1 = iom_varid( numror, 'en' , ldstop = .FALSE. ) 1159 id2 = iom_varid( numror, 'avt_k' , ldstop = .FALSE. ) 1160 id3 = iom_varid( numror, 'avm_k' , ldstop = .FALSE. ) 1161 id4 = iom_varid( numror, 'hmxl_n', ldstop = .FALSE. ) 1175 1162 ! 1176 IF( MIN( id1, id2, id3, id4 , id5, id6) > 0 ) THEN ! all required arrays exist1163 IF( MIN( id1, id2, id3, id4 ) > 0 ) THEN ! all required arrays exist 1177 1164 CALL iom_get( numror, jpdom_autoglo, 'en' , en ) 1178 CALL iom_get( numror, jpdom_autoglo, 'avt' , avt ) 1179 CALL iom_get( numror, jpdom_autoglo, 'avm' , avm ) 1180 CALL iom_get( numror, jpdom_autoglo, 'avmu' , avmu ) 1181 CALL iom_get( numror, jpdom_autoglo, 'avmv' , avmv ) 1182 CALL iom_get( numror, jpdom_autoglo, 'mxln' , mxln ) 1165 CALL iom_get( numror, jpdom_autoglo, 'avt_k' , avt_k ) 1166 CALL iom_get( numror, jpdom_autoglo, 'avm_k' , avm_k ) 1167 CALL iom_get( numror, jpdom_autoglo, 'hmxl_n', hmxl_n ) 1183 1168 ELSE 1184 IF(lwp) WRITE(numout,*) ' ===>>>> : previous run without gls scheme, en and mxln computed by iterative loop' 1185 en (:,:,:) = rn_emin 1186 mxln(:,:,:) = 0.05 1187 avt_k (:,:,:) = avt (:,:,:) 1188 avm_k (:,:,:) = avm (:,:,:) 1189 avmu_k(:,:,:) = avmu(:,:,:) 1190 avmv_k(:,:,:) = avmv(:,:,:) 1191 DO jit = nit000 + 1, nit000 + 10 ; CALL zdf_gls( jit ) ; END DO 1169 IF(lwp) WRITE(numout,*) 1170 IF(lwp) WRITE(numout,*) ' ==>> previous run without GLS scheme, set en and hmxl_n to background values' 1171 en (:,:,:) = rn_emin 1172 hmxl_n(:,:,:) = 0.05_wp 1173 ! avt_k, avm_k already set to the background value in zdf_phy_init 1192 1174 ENDIF 1193 1175 ELSE !* Start from rest 1194 IF(lwp) WRITE(numout,*) ' ===>>>> : Initialisation of en and mxln by background values' 1195 en (:,:,:) = rn_emin 1196 mxln(:,:,:) = 0.05 1176 IF(lwp) WRITE(numout,*) 1177 IF(lwp) WRITE(numout,*) ' ==>> start from rest, set en and hmxl_n by background values' 1178 en (:,:,:) = rn_emin 1179 hmxl_n(:,:,:) = 0.05_wp 1180 ! avt_k, avm_k already set to the background value in zdf_phy_init 1197 1181 ENDIF 1198 1182 ! … … 1200 1184 ! ! ------------------- 1201 1185 IF(lwp) WRITE(numout,*) '---- gls-rst ----' 1202 CALL iom_rstput( kt, nitrst, numrow, 'en' , en ) 1203 CALL iom_rstput( kt, nitrst, numrow, 'avt' , avt_k ) 1204 CALL iom_rstput( kt, nitrst, numrow, 'avm' , avm_k ) 1205 CALL iom_rstput( kt, nitrst, numrow, 'avmu' , avmu_k ) 1206 CALL iom_rstput( kt, nitrst, numrow, 'avmv' , avmv_k ) 1207 CALL iom_rstput( kt, nitrst, numrow, 'mxln' , mxln ) 1186 CALL iom_rstput( kt, nitrst, numrow, 'en' , en ) 1187 CALL iom_rstput( kt, nitrst, numrow, 'avt_k' , avt_k ) 1188 CALL iom_rstput( kt, nitrst, numrow, 'avm_k' , avm_k ) 1189 CALL iom_rstput( kt, nitrst, numrow, 'hmxl_n', hmxl_n ) 1208 1190 ! 1209 1191 ENDIF 1210 1192 ! 1211 1193 END SUBROUTINE gls_rst 1212 1213 #else1214 !!----------------------------------------------------------------------1215 !! Dummy module : NO TKE scheme1216 !!----------------------------------------------------------------------1217 LOGICAL, PUBLIC, PARAMETER :: lk_zdfgls = .FALSE. !: TKE flag1218 CONTAINS1219 SUBROUTINE zdf_gls_init ! Empty routine1220 WRITE(*,*) 'zdf_gls_init: You should not have seen this print! error?'1221 END SUBROUTINE zdf_gls_init1222 SUBROUTINE zdf_gls( kt ) ! Empty routine1223 WRITE(*,*) 'zdf_gls: You should not have seen this print! error?', kt1224 END SUBROUTINE zdf_gls1225 SUBROUTINE gls_rst( kt, cdrw ) ! Empty routine1226 INTEGER , INTENT(in) :: kt ! ocean time-step1227 CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag1228 WRITE(*,*) 'gls_rst: You should not have seen this print! error?', kt, cdrw1229 END SUBROUTINE gls_rst1230 #endif1231 1194 1232 1195 !!======================================================================
Note: See TracChangeset
for help on using the changeset viewer.
