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