- Timestamp:
- 2017-05-30T10:13:14+02:00 (7 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfgls.F90
r7991 r8093 19 19 USE domvvl ! ocean space and time domain : variable volume layer 20 20 USE zdf_oce ! ocean vertical physics 21 USE zdfsh2 ! vertical physics: shear production term of TKE 22 USE zdfbfr ! bottom friction (only for rn_bfrz0) 21 !!gm old 22 USE zdfbfr , ONLY : rn_tfrz0, rn_bfrz0 ! top/bottom roughness 23 !!gm new 24 USE zdfdrg , ONLY : r_z0_top , r_z0_bot ! top/bottom roughness 25 USE zdfdrg , ONLY : rCdU_top , rCdU_bot ! top/bottom friction 26 !!gm 23 27 USE sbc_oce ! surface boundary condition: ocean 24 28 USE phycst ! physical constants 25 29 USE zdfmxl ! mixed layer 26 USE sbcwave , ONLY: hsw ! significant wave height30 USE sbcwave , ONLY : hsw ! significant wave height 27 31 ! 28 32 USE lbclnk ! ocean lateral boundary conditions (or mpp link) … … 45 49 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: hmxl_n !: now mixing length 46 50 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: zwall !: wall function 47 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ustars2 !: Squared surface velocity scale at T-points 48 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ustarb2 !: Squared bottom velocity scale at T-points 51 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ustar2_surf !: Squared surface velocity scale at T-points 52 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ustar2_top !: Squared top velocity scale at T-points 53 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ustar2_bot !: Squared bottom velocity scale at T-points 49 54 50 55 ! !! ** Namelist namzdf_gls ** … … 101 106 REAL(wp) :: rsc_tke, rsc_psi, rpsi1, rpsi2, rpsi3, rsc_psi0 ! - - - - 102 107 REAL(wp) :: rpsi3m, rpsi3p, rpp, rmm, rnn ! - - - - 108 ! 109 REAL(wp) :: r2_3 = 2._wp/3._wp ! constant=2/3 103 110 104 111 !! * Substitutions … … 115 122 !! *** FUNCTION zdf_gls_alloc *** 116 123 !!---------------------------------------------------------------------- 117 ALLOCATE( hmxl_n(jpi,jpj,jpk) , ustar s2(jpi,jpj) ,&118 & zwall (jpi,jpj,jpk) , ustar b2(jpi,jpj) ,STAT= zdf_gls_alloc )124 ALLOCATE( hmxl_n(jpi,jpj,jpk) , ustar2_surf(jpi,jpj) , & 125 & zwall (jpi,jpj,jpk) , ustar2_top (jpi,jpj) , ustar2_bot(jpi,jpj) , STAT= zdf_gls_alloc ) 119 126 ! 120 127 IF( lk_mpp ) CALL mpp_sum ( zdf_gls_alloc ) … … 123 130 124 131 125 SUBROUTINE zdf_gls( kt )132 SUBROUTINE zdf_gls( kt, p_sh2, p_avm, p_avt ) 126 133 !!---------------------------------------------------------------------- 127 134 !! *** ROUTINE zdf_gls *** … … 130 137 !! coefficients using the GLS turbulent closure scheme. 131 138 !!---------------------------------------------------------------------- 132 INTEGER, INTENT(in) :: kt ! ocean time step 133 INTEGER :: ji, jj, jk, ibot, ibotm1, dir ! dummy loop arguments 134 REAL(wp) :: zesh2, zsigpsi, zcoef, zex1, zex2 ! local scalars 135 REAL(wp) :: ztx2, zty2, zup, zdown, zcof ! - - 136 REAL(wp) :: zratio, zrn2, zflxb, sh ! - - 139 INTEGER , INTENT(in ) :: kt ! ocean time step 140 REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: p_sh2 ! shear production term 141 REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: p_avm, p_avt ! momentum and tracer Kz (w-points) 142 ! 143 INTEGER :: ji, jj, jk ! dummy loop arguments 144 INTEGER :: ibot, ibotm1 ! local integers 145 INTEGER :: itop, itopp1 ! - - 146 REAL(wp) :: zesh2, zsigpsi, zcoef, zex1 , zex2 ! local scalars 147 REAL(wp) :: ztx2, zty2, zup, zdown, zcof, zdir ! - - 148 REAL(wp) :: zratio, zrn2, zflxb, sh , z_en ! - - 137 149 REAL(wp) :: prod, buoy, diss, zdiss, sm ! - - 138 REAL(wp) :: gh, gm, shr, dif, zsqen, zav ! - - 139 REAL(wp), POINTER, DIMENSION(:,: ) :: zdep 140 REAL(wp), POINTER, DIMENSION(:,: ) :: zkar 141 REAL(wp), POINTER, DIMENSION(:,: ) :: zflxs ! Turbulence fluxed induced by internal waves 142 REAL(wp), POINTER, DIMENSION(:,: ) :: zhsro ! Surface roughness (surface waves) 143 REAL(wp), POINTER, DIMENSION(:,:,:) :: eb ! tke at time before 144 REAL(wp), POINTER, DIMENSION(:,:,:) :: hmxl_b ! mixing length at time before 145 REAL(wp), POINTER, DIMENSION(:,:,:) :: shear ! vertical shear 146 REAL(wp), POINTER, DIMENSION(:,:,:) :: eps ! dissipation rate 147 REAL(wp), POINTER, DIMENSION(:,:,:) :: zwall_psi ! Wall function use in the wb case (ln_sigpsi) 148 REAL(wp), POINTER, DIMENSION(:,:,:) :: psi ! psi at time now 149 REAL(wp), POINTER, DIMENSION(:,:,:) :: z_elem_a ! element of the first matrix diagonal 150 REAL(wp), POINTER, DIMENSION(:,:,:) :: z_elem_b ! element of the second matrix diagonal 151 REAL(wp), POINTER, DIMENSION(:,:,:) :: z_elem_c ! element of the third matrix diagonal 152 REAL(wp), POINTER, DIMENSION(:,:,:) :: zstt, zstm ! stability function on tracer and momentum 150 REAL(wp) :: gh, gm, shr, dif, zsqen, zavt, zavm ! - - 151 REAL(wp), DIMENSION(jpi,jpj) :: zdep 152 REAL(wp), DIMENSION(jpi,jpj) :: zkar 153 REAL(wp), DIMENSION(jpi,jpj) :: zflxs ! Turbulence fluxed induced by internal waves 154 REAL(wp), DIMENSION(jpi,jpj) :: zhsro ! Surface roughness (surface waves) 155 REAL(wp), DIMENSION(jpi,jpj,jpk) :: eb ! tke at time before 156 REAL(wp), DIMENSION(jpi,jpj,jpk) :: hmxl_b ! mixing length at time before 157 REAL(wp), DIMENSION(jpi,jpj,jpk) :: eps ! dissipation rate 158 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zwall_psi ! Wall function use in the wb case (ln_sigpsi) 159 REAL(wp), DIMENSION(jpi,jpj,jpk) :: psi ! psi at time now 160 REAL(wp), DIMENSION(jpi,jpj,jpk) :: z_elem_a ! element of the first matrix diagonal 161 REAL(wp), DIMENSION(jpi,jpj,jpk) :: z_elem_b ! element of the second matrix diagonal 162 REAL(wp), DIMENSION(jpi,jpj,jpk) :: z_elem_c ! element of the third matrix diagonal 163 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zstt, zstm ! stability function on tracer and momentum 153 164 !!-------------------------------------------------------------------- 154 165 ! 155 166 IF( nn_timing == 1 ) CALL timing_start('zdf_gls') 156 167 ! 157 CALL wrk_alloc( jpi,jpj, zdep, zkar, zflxs, zhsro )158 CALL wrk_alloc( jpi,jpj,jpk, eb, hmxl_b, shear, eps, zwall_psi, z_elem_a, z_elem_b, z_elem_c, psi )159 CALL wrk_alloc( jpi,jpj,jpk, zstt, zstm )160 161 168 ! Preliminary computing 162 169 163 ustar s2 = 0._wp ; ustarb2 = 0._wp ; psi = 0._wp ; zwall_psi = 0._wp164 165 ! restore before values computed GLS alone166 avt(:,:,:) = avt_k(:,:,:)167 avm(:,:,:) = avm_k(:,:,:) 168 169 ! Compute surface and bottom friction at T-points170 ustar2_surf(:,:) = 0._wp ; psi(:,:,:) = 0._wp 171 ustar2_top (:,:) = 0._wp ; zwall_psi(:,:,:) = 0._wp 172 ustar2_bot (:,:) = 0._wp 173 174 175 176 ! Compute surface, top and bottom friction at T-points 170 177 DO jj = 2, jpjm1 171 178 DO ji = fs_2, fs_jpim1 ! vector opt. 172 179 ! 173 180 ! surface friction 174 ustar s2(ji,jj) = r1_rau0 * taum(ji,jj) * tmask(ji,jj,1)181 ustar2_surf(ji,jj) = r1_rau0 * taum(ji,jj) * tmask(ji,jj,1) 175 182 ! 183 !!gm old 176 184 ! bottom friction (explicit before friction) 177 185 ! Note that we chose here not to bound the friction as in dynbfr) … … 180 188 zty2 = ( bfrva(ji,jj) * vb(ji,jj,mbkv(ji,jj)) + bfrva(ji,jj-1) * vb(ji,jj-1,mbkv(ji,jj-1)) ) & 181 189 & * ( 1._wp - 0.5_wp * vmask(ji,jj,1) * vmask(ji,jj-1,1) ) 182 ustar b2(ji,jj) = SQRT( ztx2 * ztx2 + zty2 * zty2 ) * tmask(ji,jj,1)190 ustar2_bot(ji,jj) = SQRT( ztx2 * ztx2 + zty2 * zty2 ) * tmask(ji,jj,1) 183 191 END DO 184 192 END DO 185 193 !!gm new 194 !!gm Rq we may add here r_ke0(_top/_bot) ? ==>> think about that... 195 ! ! bottom friction (explicit before friction) 196 ! zmsku = ( 2._wp - umask(ji-1,jj,mbkt(ji,jj)) * umask(ji,jj,mbkt(ji,jj)) ) 197 ! zmskv = ( 2._wp - vmask(ji,jj-1,mbkt(ji,jj)) * vmask(ji,jj,mbkt(ji,jj)) ) ! (CAUTION: CdU<0) 198 ! ustar2_bot(ji,jj) = - rCdU_bot(ji,jj) * SQRT( ( zmsku*( ub(ji,jj,mbkt(ji,jj))+ub(ji-1,jj,mbkt(ji,jj)) ) )**2 199 ! & + ( zmskv*( vb(ji,jj,mbkt(ji,jj))+vb(ji,jj-1,mbkt(ji,jj)) ) )**2 ) 200 ! END DO 201 ! END DO 202 ! IF( ln_isfcav ) THEN !top friction 203 ! DO jj = 2, jpjm1 204 ! DO ji = fs_2, fs_jpim1 ! vector opt. 205 ! zmsku = ( 2. - umask(ji-1,jj,mikt(ji,jj)) * umask(ji,jj,mikt(ji,jj)) ) 206 ! zmskv = ( 2. - vmask(ji,jj-1,mikt(ji,jj)) * vmask(ji,jj,mikt(ji,jj)) ) ! (CAUTION: CdU<0) 207 ! ustar2_top(ji,jj) = - rCdU_top(ji,jj) * SQRT( ( zmsku*( ub(ji,jj,mikt(ji,jj))+ub(ji-1,jj,mikt(ji,jj)) ) )**2 208 ! & + ( zmskv*( vb(ji,jj,mikt(ji,jj))+vb(ji,jj-1,mikt(ji,jj)) ) )**2 ) 209 ! END DO 210 ! END DO 211 ! ENDIF 212 !!gm 186 213 SELECT CASE ( nn_z0_met ) !== Set surface roughness length ==! 187 214 CASE ( 0 ) ! Constant roughness 188 215 zhsro(:,:) = rn_hsro 189 216 CASE ( 1 ) ! Standard Charnock formula 190 zhsro(:,:) = MAX( rsbc_zs1 * ustars2(:,:), rn_hsro)217 zhsro(:,:) = MAX( rsbc_zs1 * ustar2_surf(:,:) , rn_hsro ) 191 218 CASE ( 2 ) ! Roughness formulae according to Rascle et al., Ocean Modelling (2008) 192 219 !!gm zcof = 2._wp * 0.6_wp / 28._wp 193 !!gm zdep(:,:) = 30._wp * TANH( zcof/ SQRT( MAX(ustar s2(:,:),rsmall) ) ) ! Wave age (eq. 10)194 zdep (:,:) = 30.*TANH(2.*0.3/(28.*SQRT(MAX(ustars2(:,:),rsmall))))! Wave age (eq. 10)195 zhsro(:,:) = MAX(rsbc_zs2 * ustar s2(:,:) * zdep(:,:)**1.5, rn_hsro) ! zhsro = rn_frac_hs * Hsw (eq. 11)220 !!gm zdep(:,:) = 30._wp * TANH( zcof/ SQRT( MAX(ustar2_surf(:,:),rsmall) ) ) ! Wave age (eq. 10) 221 zdep (:,:) = 30.*TANH( 2.*0.3/(28.*SQRT(MAX(ustar2_surf(:,:),rsmall))) ) ! Wave age (eq. 10) 222 zhsro(:,:) = MAX(rsbc_zs2 * ustar2_surf(:,:) * zdep(:,:)**1.5, rn_hsro) ! zhsro = rn_frac_hs * Hsw (eq. 11) 196 223 CASE ( 3 ) ! Roughness given by the wave model (coupled or read in file) 197 224 !!gm BUG missing a multiplicative coefficient.... 198 225 zhsro(:,:) = hsw(:,:) 199 226 END SELECT 200 201 ! !== Compute shear and dissipation rate ==! 202 CALL zdf_sh2( shear ) 203 204 227 ! 205 228 DO jk = 2, jpkm1 !== Compute dissipation rate ==! 206 229 DO jj = 1, jpjm1 … … 245 268 DO ji = 2, jpim1 246 269 ! 247 buoy = - avt(ji,jj,jk) * rn2(ji,jj,jk)! stratif. destruction270 buoy = - p_avt(ji,jj,jk) * rn2(ji,jj,jk) ! stratif. destruction 248 271 ! 249 272 diss = eps(ji,jj,jk) ! dissipation 250 273 ! 251 dir = 0.5_wp + SIGN( 0.5_wp, shear(ji,jj,jk) + buoy ) ! dir =1(=0) if shear(ji,jj,jk)+buoy >0(<0) 252 ! 253 zesh2 = dir*(shear(ji,jj,jk)+buoy)+(1._wp-dir)*shear(ji,jj,jk) ! production term 254 zdiss = dir*(diss/en(ji,jj,jk)) +(1._wp-dir)*(diss-buoy)/en(ji,jj,jk) ! dissipation term 274 zdir = 0.5_wp + SIGN( 0.5_wp, p_sh2(ji,jj,jk) + buoy ) ! zdir =1(=0) if shear(ji,jj,jk)+buoy >0(<0) 275 ! 276 zesh2 = zdir*(p_sh2(ji,jj,jk)+buoy)+(1._wp-zdir)*p_sh2(ji,jj,jk) ! production term 277 zdiss = zdir*(diss/en(ji,jj,jk)) +(1._wp-zdir)*(diss-buoy)/en(ji,jj,jk) ! dissipation term 278 !!gm better coding, identical results 279 ! zesh2 = p_sh2(ji,jj,jk) + zdir*buoy ! production term 280 ! zdiss = ( diss - (1._wp-zdir)*buoy ) / en(ji,jj,jk) ! dissipation term 281 !!gm 255 282 ! 256 283 ! Compute a wall function from 1. to rsc_psi*zwall/rsc_psi0 … … 269 296 zcof = rfact_tke * tmask(ji,jj,jk) 270 297 ! ! lower diagonal 271 z_elem_a(ji,jj,jk) = zcof * ( avm(ji,jj,jk ) +avm(ji,jj,jk-1) ) / ( e3t_n(ji,jj,jk-1) * e3w_n(ji,jj,jk) )298 z_elem_a(ji,jj,jk) = zcof * ( p_avm(ji,jj,jk ) + p_avm(ji,jj,jk-1) ) / ( e3t_n(ji,jj,jk-1) * e3w_n(ji,jj,jk) ) 272 299 ! ! upper diagonal 273 z_elem_c(ji,jj,jk) = zcof * ( avm(ji,jj,jk+1) +avm(ji,jj,jk ) ) / ( e3t_n(ji,jj,jk ) * e3w_n(ji,jj,jk) )300 z_elem_c(ji,jj,jk) = zcof * ( p_avm(ji,jj,jk+1) + p_avm(ji,jj,jk ) ) / ( e3t_n(ji,jj,jk ) * e3w_n(ji,jj,jk) ) 274 301 ! ! diagonal 275 302 z_elem_b(ji,jj,jk) = 1._wp - z_elem_a(ji,jj,jk) - z_elem_c(ji,jj,jk) + rdt * zdiss * wmask(ji,jj,jk) … … 293 320 CASE ( 0 ) ! Dirichlet case 294 321 ! First level 295 en(:,:,1) = rc02r * ustar s2(:,:) * (1._wp + rsbc_tke1)**(2._wp/3._wp)322 en(:,:,1) = rc02r * ustar2_surf(:,:) * (1._wp + rsbc_tke1)**r2_3 296 323 en(:,:,1) = MAX(en(:,:,1), rn_emin) 297 324 z_elem_a(:,:,1) = en(:,:,1) … … 300 327 ! 301 328 ! One level below 302 en(:,:,2) = rc02r * ustar s2(:,:) * (1._wp + rsbc_tke1 * ((zhsro(:,:)+gdepw_n(:,:,2)) &329 en(:,:,2) = rc02r * ustar2_surf(:,:) * (1._wp + rsbc_tke1 * ((zhsro(:,:)+gdepw_n(:,:,2)) & 303 330 & / zhsro(:,:) )**(1.5_wp*ra_sf))**(2._wp/3._wp) 304 331 en(:,:,2) = MAX(en(:,:,2), rn_emin ) … … 311 338 ! 312 339 ! Dirichlet conditions at k=1 313 en(:,:,1) = rc02r * ustar s2(:,:) * (1._wp + rsbc_tke1)**(2._wp/3._wp)340 en(:,:,1) = rc02r * ustar2_surf(:,:) * (1._wp + rsbc_tke1)**r2_3 314 341 en(:,:,1) = MAX(en(:,:,1), rn_emin) 315 342 z_elem_a(:,:,1) = en(:,:,1) … … 322 349 z_elem_a(:,:,2) = 0._wp 323 350 zkar(:,:) = (rl_sf + (vkarmn-rl_sf)*(1.-exp(-rtrans*gdept_n(:,:,1)/zhsro(:,:)) )) 324 zflxs(:,:) = rsbc_tke2 * ustar s2(:,:)**1.5_wp * zkar(:,:) &351 zflxs(:,:) = rsbc_tke2 * ustar2_surf(:,:)**1.5_wp * zkar(:,:) & 325 352 & * ((zhsro(:,:)+gdept_n(:,:,1)) / zhsro(:,:) )**(1.5_wp*ra_sf) 326 353 … … 340 367 DO jj = 2, jpjm1 341 368 DO ji = fs_2, fs_jpim1 ! vector opt. 369 !!gm This means that bottom and ocean w-level above have a specified "en" value. Sure ???? 370 !! With thick deep ocean level thickness, this may be quite large, no ??? 371 !! in particular in ocean cavities where top stratification can be large... 342 372 ibot = mbkt(ji,jj) + 1 ! k bottom level of w-point 343 373 ibotm1 = mbkt(ji,jj) ! k-1 bottom level of w-point but >=1 344 374 ! 345 ! Bottom level Dirichlet condition: 346 z_elem_a(ji,jj,ibot ) = 0._wp 347 z_elem_c(ji,jj,ibot ) = 0._wp 348 z_elem_b(ji,jj,ibot ) = 1._wp 349 en(ji,jj,ibot ) = MAX( rc02r * ustarb2(ji,jj), rn_emin ) 350 ! 351 ! Just above last level, Dirichlet condition again 352 z_elem_a(ji,jj,ibotm1) = 0._wp 353 z_elem_c(ji,jj,ibotm1) = 0._wp 354 z_elem_b(ji,jj,ibotm1) = 1._wp 355 en(ji,jj,ibotm1) = MAX( rc02r * ustarb2(ji,jj), rn_emin ) 356 END DO 357 END DO 375 z_en = MAX( rc02r * ustar2_bot(ji,jj), rn_emin ) 376 ! 377 ! Dirichlet condition applied at: 378 ! Bottom level (ibot) & Just above it (ibotm1) 379 z_elem_a(ji,jj,ibot) = 0._wp ; z_elem_a(ji,jj,ibotm1) = 0._wp 380 z_elem_c(ji,jj,ibot) = 0._wp ; z_elem_c(ji,jj,ibotm1) = 0._wp 381 z_elem_b(ji,jj,ibot) = 1._wp ; z_elem_b(ji,jj,ibotm1) = 1._wp 382 en (ji,jj,ibot) = z_en ; en (ji,jj,ibotm1) = z_en 383 END DO 384 END DO 385 !!gm new 386 ! IF( ln_isfcav) THEN ! top boundary (ocean cavity) 387 ! DO jj = 2, jpjm1 388 ! DO ji = fs_2, fs_jpim1 ! vector opt. 389 ! itop = mikt(ji,jj) ! k top w-point 390 ! itopp1 = mikt(ji,jj) + 1 ! k+1 1st w-point below the top one 391 ! ! ! mask at the ocean surface points 392 ! z_en = MAX( rc02r * ustar2_top(ji,jj), rn_emin ) * ( 1._wp - tmask(ji,jj,1) ) 393 ! ! 394 ! ! Dirichlet condition applied at: 395 ! ! top level (itop) & Just below it (itopp1) 396 ! z_elem_a(ji,jj,itop) = 0._wp ; z_elem_a(ji,jj,ipotm1) = 0._wp 397 ! z_elem_c(ji,jj,itop) = 0._wp ; z_elem_c(ji,jj,itopp1) = 0._wp 398 ! z_elem_b(ji,jj,itop) = 1._wp ; z_elem_b(ji,jj,itopp1) = 1._wp 399 ! en (ji,jj,itop) = zen ; en (ji,jj,itopp1) = z_en 400 ! END DO 401 ! END DO 402 ! ENDIF 403 !!gm 358 404 ! 359 405 CASE ( 1 ) ! Neumman boundary condition … … 364 410 ibotm1 = mbkt(ji,jj) ! k-1 bottom level of w-point but >=1 365 411 ! 412 z_en = MAX( rc02r * ustar2_bot(ji,jj), rn_emin ) 413 ! 366 414 ! Bottom level Dirichlet condition: 367 z_elem_a(ji,jj,ibot) = 0._wp 368 z_elem_c(ji,jj,ibot) = 0._wp 369 z_elem_b(ji,jj,ibot) = 1._wp 370 en(ji,jj,ibot) = MAX( rc02r * ustarb2(ji,jj), rn_emin ) 371 ! 372 ! Just above last level: Neumann condition 373 z_elem_b(ji,jj,ibotm1) = z_elem_b(ji,jj,ibotm1) + z_elem_c(ji,jj,ibotm1) ! Remove z_elem_c from z_elem_b 374 z_elem_c(ji,jj,ibotm1) = 0._wp 375 END DO 376 END DO 415 ! Bottom level (ibot) & Just above it (ibotm1) 416 ! Dirichlet ! Neumann 417 z_elem_a(ji,jj,ibot) = 0._wp ! ! Remove z_elem_c from z_elem_b 418 z_elem_b(ji,jj,ibot) = 1._wp ; z_elem_b(ji,jj,ibotm1) = z_elem_b(ji,jj,ibotm1) + z_elem_c(ji,jj,ibotm1) 419 z_elem_c(ji,jj,ibot) = 0._wp ; z_elem_c(ji,jj,ibotm1) = 0._wp 420 END DO 421 END DO 422 !!gm new 423 ! IF( ln_isfcav) THEN ! top boundary (ocean cavity) 424 ! DO jj = 2, jpjm1 425 ! DO ji = fs_2, fs_jpim1 ! vector opt. 426 ! itop = mikt(ji,jj) ! k top w-point 427 ! itopp1 = mikt(ji,jj) + 1 ! k+1 1st w-point below the top one 428 ! ! ! mask at the ocean surface points 429 ! z_en = MAX( rc02r * ustar2_top(ji,jj), rn_emin ) * ( 1._wp - tmask(ji,jj,1) ) 430 ! ! 431 ! ! Bottom level Dirichlet condition: 432 ! ! Bottom level (ibot) & Just above it (ibotm1) 433 ! ! Dirichlet ! Neumann 434 ! z_elem_a(ji,jj,itop) = 0._wp ! ! Remove z_elem_c from z_elem_b 435 ! z_elem_b(ji,jj,itop) = 1._wp ; z_elem_b(ji,jj,itopp1) = z_elem_b(ji,jj,itopp1) + z_elem_c(ji,jj,itopp1) 436 ! z_elem_c(ji,jj,itop) = 0._wp ; z_elem_c(ji,jj,itopp1) = 0._wp 437 ! END DO 438 ! END DO 439 ! ENDIF 440 !!gm 377 441 ! 378 442 END SELECT … … 465 529 zratio = psi(ji,jj,jk) / eb(ji,jj,jk) 466 530 ! 467 ! psi3+ : stable : B=-KhN²<0 => N²>0 if rn2>0 dir = 1 (stable) otherwisedir = 0 (unstable)468 dir = 0.5_wp + SIGN( 0.5_wp, rn2(ji,jj,jk) )469 ! 470 rpsi3 = dir * rpsi3m + ( 1._wp -dir ) * rpsi3p531 ! psi3+ : stable : B=-KhN²<0 => N²>0 if rn2>0 zdir = 1 (stable) otherwise zdir = 0 (unstable) 532 zdir = 0.5_wp + SIGN( 0.5_wp, rn2(ji,jj,jk) ) 533 ! 534 rpsi3 = zdir * rpsi3m + ( 1._wp - zdir ) * rpsi3p 471 535 ! 472 536 ! shear prod. - stratif. destruction 473 prod = rpsi1 * zratio * shear(ji,jj,jk)537 prod = rpsi1 * zratio * p_sh2(ji,jj,jk) 474 538 ! 475 539 ! stratif. destruction 476 buoy = rpsi3 * zratio * (- avt(ji,jj,jk) * rn2(ji,jj,jk) )540 buoy = rpsi3 * zratio * (- p_avt(ji,jj,jk) * rn2(ji,jj,jk) ) 477 541 ! 478 542 ! shear prod. - stratif. destruction 479 543 diss = rpsi2 * zratio * zwall(ji,jj,jk) * eps(ji,jj,jk) 480 544 ! 481 dir = 0.5_wp + SIGN( 0.5_wp, prod + buoy ) !dir =1(=0) if shear(ji,jj,jk)+buoy >0(<0)482 ! 483 zesh2 = dir * ( prod + buoy ) + (1._wp -dir ) * prod ! production term484 zdiss = dir * ( diss / psi(ji,jj,jk) ) + (1._wp -dir ) * (diss-buoy) / psi(ji,jj,jk) ! dissipation term545 zdir = 0.5_wp + SIGN( 0.5_wp, prod + buoy ) ! zdir =1(=0) if shear(ji,jj,jk)+buoy >0(<0) 546 ! 547 zesh2 = zdir * ( prod + buoy ) + (1._wp - zdir ) * prod ! production term 548 zdiss = zdir * ( diss / psi(ji,jj,jk) ) + (1._wp - zdir ) * (diss-buoy) / psi(ji,jj,jk) ! dissipation term 485 549 ! 486 550 ! building the matrix 487 551 zcof = rfact_psi * zwall_psi(ji,jj,jk) * tmask(ji,jj,jk) 488 552 ! ! lower diagonal 489 z_elem_a(ji,jj,jk) = zcof * ( avm(ji,jj,jk ) +avm(ji,jj,jk-1) ) / ( e3t_n(ji,jj,jk-1) * e3w_n(ji,jj,jk) )553 z_elem_a(ji,jj,jk) = zcof * ( p_avm(ji,jj,jk ) + p_avm(ji,jj,jk-1) ) / ( e3t_n(ji,jj,jk-1) * e3w_n(ji,jj,jk) ) 490 554 ! ! upper diagonal 491 z_elem_c(ji,jj,jk) = zcof * ( avm(ji,jj,jk+1) +avm(ji,jj,jk ) ) / ( e3t_n(ji,jj,jk ) * e3w_n(ji,jj,jk) )555 z_elem_c(ji,jj,jk) = zcof * ( p_avm(ji,jj,jk+1) + p_avm(ji,jj,jk ) ) / ( e3t_n(ji,jj,jk ) * e3w_n(ji,jj,jk) ) 492 556 ! ! diagonal 493 557 z_elem_b(ji,jj,jk) = 1._wp - z_elem_a(ji,jj,jk) - z_elem_c(ji,jj,jk) + rdt * zdiss * wmask(ji,jj,jk) … … 506 570 ! 507 571 CASE ( 0 ) ! Dirichlet boundary conditions 508 ! 509 ! Surface value 510 zdep(:,:) = zhsro(:,:) * rl_sf ! Cosmetic 511 psi (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 512 z_elem_a(:,:,1) = psi(:,:,1) 513 z_elem_c(:,:,1) = 0._wp 514 z_elem_b(:,:,1) = 1._wp 515 ! 516 ! One level below 517 zkar(:,:) = (rl_sf + (vkarmn-rl_sf)*(1._wp-exp(-rtrans*gdepw_n(:,:,2)/zhsro(:,:) ))) 518 zdep(:,:) = (zhsro(:,:) + gdepw_n(:,:,2)) * zkar(:,:) 519 psi (:,:,2) = rc0**rpp * en(:,:,2)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 520 z_elem_a(:,:,2) = 0._wp 521 z_elem_c(:,:,2) = 0._wp 522 z_elem_b(:,:,2) = 1._wp 523 ! 524 ! 572 ! 573 ! Surface value 574 zdep (:,:) = zhsro(:,:) * rl_sf ! Cosmetic 575 psi (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 576 z_elem_a(:,:,1) = psi(:,:,1) 577 z_elem_c(:,:,1) = 0._wp 578 z_elem_b(:,:,1) = 1._wp 579 ! 580 ! One level below 581 zkar (:,:) = (rl_sf + (vkarmn-rl_sf)*(1._wp-exp(-rtrans*gdepw_n(:,:,2)/zhsro(:,:) ))) 582 zdep (:,:) = (zhsro(:,:) + gdepw_n(:,:,2)) * zkar(:,:) 583 psi (:,:,2) = rc0**rpp * en(:,:,2)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 584 z_elem_a(:,:,2) = 0._wp 585 z_elem_c(:,:,2) = 0._wp 586 z_elem_b(:,:,2) = 1._wp 587 ! 525 588 CASE ( 1 ) ! Neumann boundary condition on d(psi)/dz 526 ! 527 ! Surface value: Dirichlet 528 zdep(:,:) = zhsro(:,:) * rl_sf 529 psi (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 530 z_elem_a(:,:,1) = psi(:,:,1) 531 z_elem_c(:,:,1) = 0._wp 532 z_elem_b(:,:,1) = 1._wp 533 ! 534 ! Neumann condition at k=2 535 z_elem_b(:,:,2) = z_elem_b(:,:,2) + z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b 536 z_elem_a(:,:,2) = 0._wp 537 ! 538 ! Set psi vertical flux at the surface: 539 zkar(:,:) = rl_sf + (vkarmn-rl_sf)*(1._wp-exp(-rtrans*gdept_n(:,:,1)/zhsro(:,:) )) ! Lengh scale slope 540 zdep(:,:) = ((zhsro(:,:) + gdept_n(:,:,1)) / zhsro(:,:))**(rmm*ra_sf) 541 zflxs(:,:) = (rnn + rsbc_tke1 * (rnn + rmm*ra_sf) * zdep(:,:))*(1._wp + rsbc_tke1*zdep(:,:))**(2._wp*rmm/3._wp-1_wp) 542 zdep(:,:) = rsbc_psi1 * (zwall_psi(:,:,1)*avm(:,:,1)+zwall_psi(:,:,2)*avm(:,:,2)) * & 543 & ustars2(:,:)**rmm * zkar(:,:)**rnn * (zhsro(:,:) + gdept_n(:,:,1))**(rnn-1.) 544 zflxs(:,:) = zdep(:,:) * zflxs(:,:) 545 psi(:,:,2) = psi(:,:,2) + zflxs(:,:) / e3w_n(:,:,2) 546 ! 547 ! 589 ! 590 ! Surface value: Dirichlet 591 zdep (:,:) = zhsro(:,:) * rl_sf 592 psi (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 593 z_elem_a(:,:,1) = psi(:,:,1) 594 z_elem_c(:,:,1) = 0._wp 595 z_elem_b(:,:,1) = 1._wp 596 ! 597 ! Neumann condition at k=2 598 z_elem_b(:,:,2) = z_elem_b(:,:,2) + z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b 599 z_elem_a(:,:,2) = 0._wp 600 ! 601 ! Set psi vertical flux at the surface: 602 zkar (:,:) = rl_sf + (vkarmn-rl_sf)*(1._wp-exp(-rtrans*gdept_n(:,:,1)/zhsro(:,:) )) ! Lengh scale slope 603 zdep (:,:) = ((zhsro(:,:) + gdept_n(:,:,1)) / zhsro(:,:))**(rmm*ra_sf) 604 zflxs(:,:) = (rnn + rsbc_tke1 * (rnn + rmm*ra_sf) * zdep(:,:))*(1._wp + rsbc_tke1*zdep(:,:))**(2._wp*rmm/3._wp-1_wp) 605 zdep (:,:) = rsbc_psi1 * (zwall_psi(:,:,1)*avm(:,:,1)+zwall_psi(:,:,2)*avm(:,:,2)) * & 606 & ustar2_surf(:,:)**rmm * zkar(:,:)**rnn * (zhsro(:,:) + gdept_n(:,:,1))**(rnn-1.) 607 zflxs(:,:) = zdep(:,:) * zflxs(:,:) 608 psi (:,:,2) = psi(:,:,2) + zflxs(:,:) / e3w_n(:,:,2) 609 ! 548 610 END SELECT 549 611 … … 552 614 ! 553 615 SELECT CASE ( nn_bc_bot ) 554 !555 616 ! 556 617 CASE ( 0 ) ! Dirichlet … … 597 658 ! Set psi vertical flux at the bottom: 598 659 zdep(ji,jj) = rn_bfrz0 + 0.5_wp*e3t_n(ji,jj,ibotm1) 599 zflxb = rsbc_psi2 * ( avm(ji,jj,ibot) +avm(ji,jj,ibotm1) ) &660 zflxb = rsbc_psi2 * ( p_avm(ji,jj,ibot) + p_avm(ji,jj,ibotm1) ) & 600 661 & * (0.5_wp*(en(ji,jj,ibot)+en(ji,jj,ibotm1)))**rmm * zdep(ji,jj)**(rnn-1._wp) 601 662 psi(ji,jj,ibotm1) = psi(ji,jj,ibotm1) + zflxb / e3w_n(ji,jj,ibotm1) … … 730 791 gh = gh * rf6 731 792 ! Gm = M²l²/q² Shear number 732 shr = shear(ji,jj,jk) / MAX(avm(ji,jj,jk), rsmall )793 shr = p_sh2(ji,jj,jk) / MAX( p_avm(ji,jj,jk), rsmall ) 733 794 gm = MAX( shr * zcof , 1.e-10 ) 734 795 gm = gm * rf6 … … 765 826 DO jj = 2, jpjm1 766 827 DO ji = fs_2, fs_jpim1 ! vector opt. 767 zsqen 768 zav 769 avt(ji,jj,jk) = MAX( zav, avtb(jk) ) * wmask(ji,jj,jk) ! apply mask for zdfmxl routine770 zav = zsqen * zstm(ji,jj,jk)771 avm(ji,jj,jk) = MAX( zav, avmb(jk) ) ! Note that avm is not masked at the surface and the bottom828 zsqen = SQRT( 2._wp * en(ji,jj,jk) ) * hmxl_n(ji,jj,jk) 829 zavt = zsqen * zstt(ji,jj,jk) 830 zavm = zsqen * zstm(ji,jj,jk) 831 p_avt(ji,jj,jk) = MAX( zavt, avtb(jk) ) * wmask(ji,jj,jk) ! apply mask for zdfmxl routine 832 p_avm(ji,jj,jk) = MAX( zavm, avmb(jk) ) ! Note that avm is not masked at the surface and the bottom 772 833 END DO 773 834 END DO 774 835 END DO 775 836 avt(:,:,1) = 0._wp 776 !!gm I'm sure this is needed to compute the shear term at next time-step777 CALL lbc_lnk( avm, 'W', 1. ) ; CALL lbc_lnk( avt, 'W', 1. )778 837 ! 779 838 IF(ln_ctl) THEN … … 781 840 CALL prt_ctl( tab3d_1=avm, clinfo1=' gls - m: ', ovlap=1, kdim=jpk ) 782 841 ENDIF 783 !784 avt_k(:,:,:) = avt(:,:,:) !== store avt, avm values computed by GLS only ==!785 avm_k(:,:,:) = avm(:,:,:)786 !787 IF( lrst_oce ) CALL gls_rst( kt, 'WRITE' ) ! write the GLS info in the restart file788 !789 CALL wrk_dealloc( jpi,jpj, zdep, zkar, zflxs, zhsro )790 CALL wrk_dealloc( jpi,jpj,jpk, eb, hmxl_b, shear, eps, zwall_psi, z_elem_a, z_elem_b, z_elem_c, psi )791 CALL wrk_dealloc( jpi,jpj,jpk, zstt, zstm )792 842 ! 793 843 IF( nn_timing == 1 ) CALL timing_stop('zdf_gls') … … 837 887 IF(lwp) THEN !* Control print 838 888 WRITE(numout,*) 839 WRITE(numout,*) 'zdf_gls_init : glsturbulent closure scheme'889 WRITE(numout,*) 'zdf_gls_init : GLS turbulent closure scheme' 840 890 WRITE(numout,*) '~~~~~~~~~~~~' 841 891 WRITE(numout,*) ' Namelist namzdf_gls : set gls mixing parameters' … … 854 904 WRITE(numout,*) ' Type of closure nn_clos = ', nn_clos 855 905 WRITE(numout,*) ' Surface roughness (m) rn_hsro = ', rn_hsro 906 !!gm old 907 WRITE(numout,*) ' top roughness (m) (nambfr namelist) rn_tfrz0 = ', rn_tfrz0 856 908 WRITE(numout,*) ' Bottom roughness (m) (nambfr namelist) rn_bfrz0 = ', rn_bfrz0 909 !!gm new 910 ! WRITE(numout,*) ' Namelist namdrg_top/_bot used values:' 911 ! WRITE(numout,*) ' top ocean cavity roughness (m) rn_z0(_top) = ', r_z0_top 912 ! WRITE(numout,*) ' Bottom seafloor roughness (m) rn_z0(_bot) = ', r_z0_bot 913 !!gm 914 WRITE(numout,*) 857 915 ENDIF 858 916 … … 861 919 862 920 ! !* Check of some namelist values 863 IF( nn_bc_surf < 0 .OR. nn_bc_surf > 1 )CALL ctl_stop( 'zdf_gls_init: bad flag: nn_bc_surf is 0 or 1' )864 IF( nn_bc_surf < 0 .OR. nn_bc_surf > 1 )CALL ctl_stop( 'zdf_gls_init: bad flag: nn_bc_surf is 0 or 1' )865 IF( nn_z0_met < 0 .OR. nn_z0_met > 3 )CALL ctl_stop( 'zdf_gls_init: bad flag: nn_z0_met is 0, 1, 2 or 3' )866 IF( nn_z0_met == 3 .AND. .NOT.ln_wave )CALL ctl_stop( 'zdf_gls_init: nn_z0_met=3 requires ln_wave=T' )867 IF( nn_stab_func < 0 .OR. nn_stab_func > 3 )CALL ctl_stop( 'zdf_gls_init: bad flag: nn_stab_func is 0, 1, 2 and 3' )868 IF( nn_clos < 0 .OR. nn_clos > 3 )CALL ctl_stop( 'zdf_gls_init: bad flag: nn_clos is 0, 1, 2 or 3' )921 IF( nn_bc_surf < 0 .OR. nn_bc_surf > 1 ) CALL ctl_stop( 'zdf_gls_init: bad flag: nn_bc_surf is 0 or 1' ) 922 IF( nn_bc_surf < 0 .OR. nn_bc_surf > 1 ) CALL ctl_stop( 'zdf_gls_init: bad flag: nn_bc_surf is 0 or 1' ) 923 IF( nn_z0_met < 0 .OR. nn_z0_met > 3 ) CALL ctl_stop( 'zdf_gls_init: bad flag: nn_z0_met is 0, 1, 2 or 3' ) 924 IF( nn_z0_met == 3 .AND. .NOT.ln_wave ) CALL ctl_stop( 'zdf_gls_init: nn_z0_met=3 requires ln_wave=T' ) 925 IF( nn_stab_func < 0 .OR. nn_stab_func > 3 ) CALL ctl_stop( 'zdf_gls_init: bad flag: nn_stab_func is 0, 1, 2 and 3' ) 926 IF( nn_clos < 0 .OR. nn_clos > 3 ) CALL ctl_stop( 'zdf_gls_init: bad flag: nn_clos is 0, 1, 2 or 3' ) 869 927 870 928 SELECT CASE ( nn_clos ) !* set the parameters for the chosen closure … … 872 930 CASE( 0 ) ! k-kl (Mellor-Yamada) 873 931 ! 874 IF(lwp) WRITE(numout,*) 'The choosen closure is k-kl closed to the classical Mellor-Yamada' 932 IF(lwp) WRITE(numout,*) ' ==>> k-kl closure chosen (i.e. closed to the classical Mellor-Yamada)' 933 IF(lwp) WRITE(numout,*) 875 934 rpp = 0._wp 876 935 rmm = 1._wp … … 890 949 CASE( 1 ) ! k-eps 891 950 ! 892 IF(lwp) WRITE(numout,*) 'The choosen closure is k-eps' 951 IF(lwp) WRITE(numout,*) ' ==>> k-eps closure chosen' 952 IF(lwp) WRITE(numout,*) 893 953 rpp = 3._wp 894 954 rmm = 1.5_wp … … 908 968 CASE( 2 ) ! k-omega 909 969 ! 910 IF(lwp) WRITE(numout,*) 'The choosen closure is k-omega' 970 IF(lwp) WRITE(numout,*) ' ==>> k-omega closure chosen' 971 IF(lwp) WRITE(numout,*) 911 972 rpp = -1._wp 912 973 rmm = 0.5_wp … … 926 987 CASE( 3 ) ! generic 927 988 ! 928 IF(lwp) WRITE(numout,*) 'The choosen closure is generic' 989 IF(lwp) WRITE(numout,*) ' ==>> generic closure chosen' 990 IF(lwp) WRITE(numout,*) 929 991 rpp = 2._wp 930 992 rmm = 1._wp … … 949 1011 CASE ( 0 ) ! Galperin stability functions 950 1012 ! 951 IF(lwp) WRITE(numout,*) ' Stability functions from Galperin'1013 IF(lwp) WRITE(numout,*) ' ==>> Stability functions from Galperin' 952 1014 rc2 = 0._wp 953 1015 rc3 = 0._wp … … 961 1023 CASE ( 1 ) ! Kantha-Clayson stability functions 962 1024 ! 963 IF(lwp) WRITE(numout,*) ' Stability functions from Kantha-Clayson'1025 IF(lwp) WRITE(numout,*) ' ==>> Stability functions from Kantha-Clayson' 964 1026 rc2 = 0.7_wp 965 1027 rc3 = 0.2_wp … … 973 1035 CASE ( 2 ) ! Canuto A stability functions 974 1036 ! 975 IF(lwp) WRITE(numout,*) ' Stability functions from Canuto A'1037 IF(lwp) WRITE(numout,*) ' ==>> Stability functions from Canuto A' 976 1038 rs0 = 1.5_wp * rl1 * rl5*rl5 977 1039 rs1 = -rl4*(rl6+rl7) + 2._wp*rl4*rl5*(rl1-(1._wp/3._wp)*rl2-rl3) + 1.5_wp*rl1*rl5*rl8 … … 997 1059 CASE ( 3 ) ! Canuto B stability functions 998 1060 ! 999 IF(lwp) WRITE(numout,*) ' Stability functions from Canuto B'1061 IF(lwp) WRITE(numout,*) ' ==>> Stability functions from Canuto B' 1000 1062 rs0 = 1.5_wp * rm1 * rm5*rm5 1001 1063 rs1 = -rm4 * (rm6+rm7) + 2._wp * rm4*rm5*(rm1-(1._wp/3._wp)*rm2-rm3) + 1.5_wp * rm1*rm5*rm8 … … 1052 1114 IF(lwp) THEN !* Control print 1053 1115 WRITE(numout,*) 1054 WRITE(numout,*) 'Limit values' 1055 WRITE(numout,*) '~~~~~~~~~~~~' 1056 WRITE(numout,*) 'Parameter m = ',rmm 1057 WRITE(numout,*) 'Parameter n = ',rnn 1058 WRITE(numout,*) 'Parameter p = ',rpp 1059 WRITE(numout,*) 'rpsi1 = ',rpsi1 1060 WRITE(numout,*) 'rpsi2 = ',rpsi2 1061 WRITE(numout,*) 'rpsi3m = ',rpsi3m 1062 WRITE(numout,*) 'rpsi3p = ',rpsi3p 1063 WRITE(numout,*) 'rsc_tke = ',rsc_tke 1064 WRITE(numout,*) 'rsc_psi = ',rsc_psi 1065 WRITE(numout,*) 'rsc_psi0 = ',rsc_psi0 1066 WRITE(numout,*) 'rc0 = ',rc0 1116 WRITE(numout,*) ' Limit values :' 1117 WRITE(numout,*) ' Parameter m = ', rmm 1118 WRITE(numout,*) ' Parameter n = ', rnn 1119 WRITE(numout,*) ' Parameter p = ', rpp 1120 WRITE(numout,*) ' rpsi1 = ', rpsi1 1121 WRITE(numout,*) ' rpsi2 = ', rpsi2 1122 WRITE(numout,*) ' rpsi3m = ', rpsi3m 1123 WRITE(numout,*) ' rpsi3p = ', rpsi3p 1124 WRITE(numout,*) ' rsc_tke = ', rsc_tke 1125 WRITE(numout,*) ' rsc_psi = ', rsc_psi 1126 WRITE(numout,*) ' rsc_psi0 = ', rsc_psi0 1127 WRITE(numout,*) ' rc0 = ', rc0 1067 1128 WRITE(numout,*) 1068 WRITE(numout,*) 'Shear free turbulence parameters:' 1069 WRITE(numout,*) 'rcm_sf = ',rcm_sf 1070 WRITE(numout,*) 'ra_sf = ',ra_sf 1071 WRITE(numout,*) 'rl_sf = ',rl_sf 1072 WRITE(numout,*) 1129 WRITE(numout,*) ' Shear free turbulence parameters:' 1130 WRITE(numout,*) ' rcm_sf = ', rcm_sf 1131 WRITE(numout,*) ' ra_sf = ', ra_sf 1132 WRITE(numout,*) ' rl_sf = ', rl_sf 1073 1133 ENDIF 1074 1134 … … 1090 1150 ! 1091 1151 ! !* Wall proximity function 1092 zwall (:,:,:) = 1._wp * tmask(:,:,:) 1152 !!gm tmask or wmask ???? 1153 zwall(:,:,:) = 1._wp * tmask(:,:,:) 1093 1154 1094 1155 ! !* read or initialize all required files … … 1133 1194 CALL iom_get( numror, jpdom_autoglo, 'hmxl_n', hmxl_n ) 1134 1195 ELSE 1135 IF(lwp) WRITE(numout,*) ' ===>>>> : previous run without gls scheme, en and hmxl_n computed by iterative loop' 1196 IF(lwp) WRITE(numout,*) 1197 IF(lwp) WRITE(numout,*) ' ==>> previous run without GLS scheme, set en and hmxl_n to background values' 1136 1198 en (:,:,:) = rn_emin 1137 1199 hmxl_n(:,:,:) = 0.05_wp 1138 avt_k (:,:,:) = avt(:,:,:) 1139 avm_k (:,:,:) = avm(:,:,:) 1140 DO jit = nit000 + 1, nit000 + 10 ; CALL zdf_gls( jit ) ; END DO 1200 ! avt_k, avm_k already set to the background value in zdf_phy_init 1141 1201 ENDIF 1142 1202 ELSE !* Start from rest 1143 IF(lwp) WRITE(numout,*) ' ===>>>> : Initialisation of en and hmxl_n by background values' 1203 IF(lwp) WRITE(numout,*) 1204 IF(lwp) WRITE(numout,*) ' ==>> start from rest, set en and hmxl_n by background values' 1144 1205 en (:,:,:) = rn_emin 1145 1206 hmxl_n(:,:,:) = 0.05_wp 1146 DO jk = 1, jpk 1147 avt_k(:,:,jk) = avtb(jk) * wmask(:,:,jk) 1148 avm_k(:,:,jk) = avmb(jk) * wmask(:,:,jk) 1149 END DO 1207 ! avt_k, avm_k already set to the background value in zdf_phy_init 1150 1208 ENDIF 1151 1209 !
Note: See TracChangeset
for help on using the changeset viewer.
