Changeset 1060 for trunk/NEMO/OPA_SRC/ZDF/zdftke.F90
- Timestamp:
- 2008-06-05T11:47:45+02:00 (16 years ago)
- File:
-
- 1 edited
-
trunk/NEMO/OPA_SRC/ZDF/zdftke.F90 (modified) (39 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMO/OPA_SRC/ZDF/zdftke.F90
r1050 r1060 2 2 !!====================================================================== 3 3 !! *** MODULE zdftke *** 4 !! Ocean physics: vertical mixing coefficient compute from the tke4 !! Ocean physics: vertical mixing coefficient computed from the tke 5 5 !! turbulent closure parameterization 6 6 !!===================================================================== 7 !! History : 6.0 ! 91-03 (b. blanke) Original code 8 !! 7.0 ! 91-11 (G. Madec) bug fix 9 !! 7.1 ! 92-10 (G. Madec) new mixing length and eav 10 !! 7.2 ! 93-03 (M. Guyon) symetrical conditions 11 !! 7.3 ! 94-08 (G. Madec, M. Imbard) npdl flag 12 !! 7.5 ! 96-01 (G. Madec) s-coordinates 13 !! 8.0 ! 97-07 (G. Madec) lbc 14 !! 8.1 ! 99-01 (E. Stretta) new option for the mixing length 15 !! 8.5 ! 02-06 (G. Madec) add zdf_tke_init routine 16 !! 8.5 ! 02-08 (G. Madec) ri_c and Free form, F90 17 !! 9.0 ! 04-10 (C. Ethe ) 1D configuration 18 !! 9.0 ! 02-08 (G. Madec) autotasking optimization 19 !! 9.0 ! 06-07 (S. Masson) distributed restart using iom 20 !! - ! 2005-07 (C. Ethe, G.Madec) : update TKE physics: 7 !! History : OPA ! 1991-03 (b. blanke) Original code 8 !! 7.0 ! 1991-11 (G. Madec) bug fix 9 !! 7.1 ! 1992-10 (G. Madec) new mixing length and eav 10 !! 7.2 ! 1993-03 (M. Guyon) symetrical conditions 11 !! 7.3 ! 1994-08 (G. Madec, M. Imbard) nn_pdl flag 12 !! 7.5 ! 1996-01 (G. Madec) s-coordinates 13 !! 8.0 ! 1997-07 (G. Madec) lbc 14 !! 8.1 ! 1999-01 (E. Stretta) new option for the mixing length 15 !! NEMO 1.0 ! 2002-06 (G. Madec) add zdf_tke_init routine 16 !! - ! 2002-08 (G. Madec) rn_cri and Free form, F90 17 !! - ! 2004-10 (C. Ethe ) 1D configuration 18 !! 2.0 ! 2006-07 (S. Masson) distributed restart using iom 19 !! 3.0 ! 2008-05 (C. Ethe, G.Madec) : update TKE physics: 21 20 !! - tke penetration (wind steering) 22 21 !! - suface condition for tke & mixing length 23 22 !! - Langmuir cells 24 !! 3.0 ! 2007-11 (G. Madec) avtb_2d, nn_avtb_2d 23 !! - ! 2008-05 (J.-M. Molines, G. Madec) 2D form of avtb 24 !! - ! 2008-06 (G. Madec) style + DOCTOR name for namelist parameters 25 25 !!---------------------------------------------------------------------- 26 26 #if defined key_zdftke || defined key_esopa … … 39 39 USE phycst ! physical constants 40 40 USE zdfmxl ! mixed layer 41 USE restart ! only for lrst_oce 41 42 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 42 43 USE prtctl ! Print control 43 44 USE in_out_manager ! I/O manager 44 USE iom 45 USE restart ! only for lrst_oce 45 USE iom ! I/O manager library 46 46 47 47 IMPLICIT NONE 48 48 PRIVATE 49 49 50 PUBLIC zdf_tke 50 PUBLIC zdf_tke ! routine called in step module 51 51 52 52 LOGICAL , PUBLIC, PARAMETER :: lk_zdftke = .TRUE. !: TKE vertical mixing flag … … 54 54 REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) :: en !: now turbulent kinetic energy 55 55 # if defined key_vectopt_memory 56 ! !!! key_vectopt_memory 56 57 REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) :: etmean !: coefficient used for horizontal smoothing 57 58 REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) :: eumean, evmean !: at t-, u- and v-points 58 59 # endif 59 60 !! * Namelist (namtke) 61 LOGICAL , PUBLIC :: ln_rstke = .FALSE. !: =T restart with tke from a run without tke with 62 ! ! a none zero initial value for en 63 INTEGER , PUBLIC :: nitke = 50 , & !: number of restart iterative loops 64 & nmxl = 2 , & !: = 0/1/2/3 flag for the type of mixing length used 65 & npdl = 1 , & !: = 0/1/2 flag on prandtl number on vert. eddy coeff. 66 & nave = 1 , & !: = 0/1 flag for horizontal average on avt, avmu, avmv 67 & navb = 0 !: = 0/1 flag for constant or profile background avt 68 REAL(wp), PUBLIC :: ediff = 0.1_wp , & !: coeff. for vertical eddy coef.; avt=ediff*mxl*sqrt(e) 69 & ediss = 0.7_wp , & !: coef. of the Kolmogoroff dissipation 70 & ebb = 3.75_wp , & !: coef. of the surface input of tke 71 & efave = 1._wp , & !: coef. for the tke vert. diff. coeff.; avtke=efave*avm 72 & emin = 0.7071e-6_wp , & !: minimum value of tke (m2/s2) 73 & emin0 = 1.e-4_wp , & !: surface minimum value of tke (m2/s2) 74 & ri_c = 2._wp / 9._wp !: critic Richardson number 75 ! !!! ** Namelist (namtke) ** 76 INTEGER :: nn_avtb_2d = 1 !: = 0/1 horizontal shape for avtb 77 INTEGER :: nn_etau = 0 !: = 0/1/2 tke depth penetration 78 INTEGER :: nn_htau = 0 !: = 0/1/2 type of tke profile of penetration 79 REAL(wp) :: rn_lmin = 0.4_wp !: surface min value of mixing turbulent length scale 80 REAL(wp) :: rn_efr = 1.0_wp !: fraction of TKE surface value which penetrates in the ocean 81 LOGICAL :: ln_lsfc = .FALSE. !: mixing length scale surface value as function of wind stress or not 82 LOGICAL :: ln_lc = .FALSE. !: Langmuir cells (LC) as a source term of TKE or not 83 REAL(wp) :: rn_lc = 0.15_wp !: coef to compute vertical velocity of LC 84 85 REAL(wp), DIMENSION (jpi,jpj) :: avtb_2d ! set in tke_init, for other modif than ice 86 87 # if defined key_c1d 88 ! ! 1D cfg only 89 REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) :: e_dis, e_mix, & ! dissipation and mixing turbulent lengh scales 90 & e_pdl, e_ric ! prandl and local Richardson numbers 60 #if defined key_c1d 61 ! !!! 1D cfg only 62 REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) :: e_dis, e_mix !: dissipation and mixing turbulent lengh scales 63 REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) :: e_pdl, e_ric !: prandl and local Richardson numbers 91 64 #endif 65 66 ! !!! ** Namelist namtke ** 67 LOGICAL :: ln_rstke = .FALSE. ! =T restart with tke from a run without tke 68 LOGICAL :: ln_mxl0 = .FALSE. ! mixing length scale surface value as function of wind stress or not 69 LOGICAL :: ln_lc = .FALSE. ! Langmuir cells (LC) as a source term of TKE or not 70 INTEGER :: nn_itke = 50 ! number of restart iterative loops 71 INTEGER :: nn_mxl = 2 ! type of mixing length (=0/1/2/3) 72 INTEGER :: nn_pdl = 1 ! Prandtl number or not (ratio avt/avm) (=0/1) 73 INTEGER :: nn_ave = 1 ! horizontal average or not on avt, avmu, avmv (=0/1) 74 INTEGER :: nn_avb = 0 ! constant or profile background on avt (=0/1) 75 REAL(wp) :: rn_ediff = 0.1_wp ! coefficient for avt: avt=rn_ediff*mxl*sqrt(e) 76 REAL(wp) :: rn_ediss = 0.7_wp ! coefficient of the Kolmogoroff dissipation 77 REAL(wp) :: rn_ebb = 3.75_wp ! coefficient of the surface input of tke 78 REAL(wp) :: rn_efave = 1._wp ! coefficient for ave : ave=rn_efave*avm 79 REAL(wp) :: rn_emin = 0.7071e-6_wp ! minimum value of tke (m2/s2) 80 REAL(wp) :: rn_emin0 = 1.e-4_wp ! surface minimum value of tke (m2/s2) 81 REAL(wp) :: rn_cri = 2._wp / 9._wp ! critic Richardson number 82 INTEGER :: nn_havtb = 1 ! horizontal shape or not for avtb (=0/1) 83 INTEGER :: nn_etau = 0 ! type of depth penetration of surface tke (=0/1/2) 84 INTEGER :: nn_htau = 0 ! type of tke profile of penetration (=0/1/2) 85 REAL(wp) :: rn_lmin0 = 0.4_wp ! surface min value of mixing length 86 REAL(wp) :: rn_lmin = 0.4_wp ! interior min value of mixing length 87 REAL(wp) :: rn_efr = 1.0_wp ! fraction of TKE surface value which penetrates in the ocean 88 REAL(wp) :: rn_lc = 0.15_wp ! coef to compute vertical velocity of Langmuir cells 89 90 REAL(wp), DIMENSION (jpi,jpj) :: avtb_2d ! set in tke_init, for other modif than ice 92 91 93 92 !! * Substitutions … … 95 94 # include "vectopt_loop_substitute.h90" 96 95 !!---------------------------------------------------------------------- 97 !! OPA 9.0 , LOCEAN-IPSL (2006)96 !! NEMO/OPA 3.0 , LOCEAN-IPSL (2008) 98 97 !! $Id$ 99 98 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) … … 111 110 !! ** Method : The time evolution of the turbulent kinetic energy 112 111 !! (tke) is computed from a prognostic equation : 113 !! d(en)/dt = eboost eav (d(u)/dz)**2 ! shear production114 !! + d( efave eav d(en)/dz )/dz! diffusion of tke115 !! + grav/rau0 pdl eav d(rau)/dz ! stratif. destruc.116 !! - ediss / emxl en**(2/3)! dissipation112 !! d(en)/dt = eboost eav (d(u)/dz)**2 ! shear production 113 !! + d( rn_efave eav d(en)/dz )/dz ! diffusion of tke 114 !! + grav/rau0 pdl eav d(rau)/dz ! stratif. destruc. 115 !! - rn_ediss / emxl en**(2/3) ! dissipation 117 116 !! with the boundary conditions: 118 !! surface: en = max( emin0,ebb sqrt(utau^2 + vtau^2) )119 !! bottom : en = emin117 !! surface: en = max( rn_emin0, rn_ebb sqrt(utau^2 + vtau^2) ) 118 !! bottom : en = rn_emin 120 119 !! -1- The dissipation and mixing turbulent lengh scales are computed 121 !! from the usual diagnostic buoyancy length scale: 122 !! mxl= 1/(sqrt(en)/N) WHERE N is the brunt-vaisala frequency 120 !! from the usual diagnostic buoyancy length scale: 121 !! mxl= sqrt(2*en)/N where N is the brunt-vaisala frequency 122 !! with mxl = rn_lmin at the bottom minimum value of 0.4 123 123 !! Four cases : 124 !! n mxl=0 : mxl bounded by the distance to surface and bottom.124 !! nn_mxl=0 : mxl bounded by the distance to surface and bottom. 125 125 !! zmxld = zmxlm = mxl 126 !! n mxl=1 : mxl bounded by the vertical scale factor.126 !! nn_mxl=1 : mxl bounded by the vertical scale factor. 127 127 !! zmxld = zmxlm = mxl 128 !! n mxl=2 : mxl bounded such that the vertical derivative of mxl128 !! nn_mxl=2 : mxl bounded such that the vertical derivative of mxl 129 129 !! is less than 1 (|d/dz(xml)|<1). 130 130 !! zmxld = zmxlm = mxl 131 !! n mxl=3 : lup = mxl bounded using |d/dz(xml)|<1 from the surface131 !! nn_mxl=3 : lup = mxl bounded using |d/dz(xml)|<1 from the surface 132 132 !! to the bottom 133 133 !! ldown = mxl bounded using |d/dz(xml)|<1 from the bottom … … 140 140 !! solved. 141 141 !! Note that - the shear production is multiplied by eboost in order 142 !! to set the critic richardson number to r i_c(namelist parameter)142 !! to set the critic richardson number to rn_cri (namelist parameter) 143 143 !! - the destruction by stratification term is multiplied 144 144 !! by the Prandtl number (defined by an empirical funtion of the local 145 !! Richardson number) if n pdl=1 (namelist parameter)145 !! Richardson number) if nn_pdl=1 (namelist parameter) 146 146 !! coefficient (zesh2): 147 147 !! -3- Compute the now vertical eddy vicosity and diffusivity 148 148 !! coefficients from en (before the time stepping) and zmxlm: 149 !! avm = max( avtb, ediff*zmxlm*en^1/2 )150 !! avt = max( avmb, pdl*avm ) (pdl=1 if n pdl=0)149 !! avm = max( avtb, rn_ediff*zmxlm*en^1/2 ) 150 !! avt = max( avmb, pdl*avm ) (pdl=1 if nn_pdl=0) 151 151 !! eav = max( avmb, avm ) 152 152 !! avt and avm are horizontally averaged to avoid numerical insta- … … 159 159 !! update avt, avmu, avmv (before vertical eddy coef.) 160 160 !! 161 !! References : Gaspar et al., jgr, 95, 1990, 162 !! Blanke and Delecluse, jpo, 1991 161 !! References : Gaspar et al., JGR, 1990, 162 !! Blanke and Delecluse, JPO, 1991 163 !! Mellor and Blumberg, JPO 2004 164 !! Axell, JGR, 2002 163 165 !!---------------------------------------------------------------------- 164 USE oce , zwd => ua, & ! use ua as workspace 165 & zmxlm => ta, & ! use ta as workspace 166 & zmxld => sa ! use sa as workspace 167 USE oce , ztkelc => va ! use va as workspace 168 ! 169 INTEGER, INTENT(in) :: kt ! ocean time step 170 ! 171 INTEGER :: ji, jj, jk ! dummy loop arguments 172 REAL(wp) :: zmlmin, zbbrau, & ! temporary scalars 173 & zfact1, zfact2, zfact3, & ! 174 & zrn2, zesurf, & ! 175 & ztx2, zty2, zav, & ! 176 & zcoef, zcof, zsh2, & ! 177 & zdku, zdkv, zpdl, zri, & ! 178 & zsqen, zesh2, & ! 179 & zemxl, zemlm, zemlp 166 USE oce, zwd => ua ! use ua as workspace 167 USE oce, zmxlm => va ! use va as workspace 168 USE oce, zmxld => ta ! use ta as workspace 169 USE oce, ztkelc => sa ! use sa as workspace 170 ! 171 INTEGER, INTENT(in) :: kt ! ocean time step 172 ! 173 INTEGER :: ji, jj, jk ! dummy loop arguments 174 REAL(wp) :: zbbrau, zrn2, zesurf ! temporary scalars 175 REAL(wp) :: zfact1, ztx2, zdku ! - - 176 REAL(wp) :: zfact2, zty2, zdkv ! - - 177 REAL(wp) :: zfact3, zcoef, zcof, zav ! - - 178 REAL(wp) :: zsh2, zpdl, zri, zsqen, zesh2 ! - - 179 REAL(wp) :: zemxl, zemlm, zemlp ! - - 180 REAL(wp) :: zraug, zus, zwlc, zind ! - - 180 181 INTEGER , DIMENSION(jpi,jpj) :: imlc ! 2D workspace 181 REAL(wp) :: zraug, zus, zwlc, zind ! temporary scalar 182 REAL(wp), DIMENSION(jpi,jpj) :: zhtau ! 2D workspace 183 REAL(wp), DIMENSION(jpi,jpj) :: zhlc ! 2D workspace 182 REAL(wp), DIMENSION(jpi,jpj) :: zhtau ! - - 183 REAL(wp), DIMENSION(jpi,jpj) :: zhlc ! - - 184 184 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zpelc ! 3D workspace 185 185 !!-------------------------------------------------------------------- … … 188 188 189 189 ! ! Local constant initialization 190 zmlmin = 0.4 191 zbbrau = .5 * ebb / rau0 192 zfact1 = -.5 * rdt * efave 193 zfact2 = 1.5 * rdt * ediss 194 zfact3 = 0.5 * rdt * ediss 190 zbbrau = .5 * rn_ebb / rau0 191 zfact1 = -.5 * rdt * rn_efave 192 zfact2 = 1.5 * rdt * rn_ediss 193 zfact3 = 0.5 * rdt * rn_ediss 195 194 196 195 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> … … 200 199 ! Buoyancy length scale: l=sqrt(2*e/n**2) 201 200 ! --------------------- 202 IF( ln_lsfc ) THEN ! lsurf is function of wind stress : l=2.e5*sqrt(utau^2 + vtau^2)/(rau0*g) 201 IF( ln_mxl0 ) THEN ! surface mixing length = F(stress) : l=2.e5*sqrt(utau^2 + vtau^2)/(rau0*g) 202 !!gm this should be useless 203 203 zmxlm(:,:,1) = 0.e0 204 !!gm end 204 205 zraug = 0.5 * 2.e5 / ( rau0 * grav ) 205 206 DO jj = 2, jpjm1 … … 208 209 ztx2 = utau(ji-1,jj ) + utau(ji,jj) 209 210 zty2 = vtau(ji ,jj-1) + vtau(ji,jj) 210 zmxlm(ji,jj,1 ) = zraug * SQRT( ztx2 * ztx2 + zty2 * zty2 ) 211 ! set the surface minimum value to rn_lmin 212 zmxlm(ji,jj,1 ) = MAX( zmxlm(ji,jj,1) , rn_lmin ) 213 END DO 214 END DO 215 ELSE ! surface set to the minimum value 216 zmxlm(:,:,1) = rn_lmin 211 zmxlm(ji,jj,1) = MAX( rn_lmin0, zraug * SQRT( ztx2 * ztx2 + zty2 * zty2 ) ) 212 END DO 213 END DO 214 ELSE ! surface set to the minimum value 215 zmxlm(:,:,1) = rn_lmin0 217 216 ENDIF 218 zmxlm(:,:,jpk) = rn_lmin ! bottom set to the minimum value 219 !CDIR NOVERRCHK 220 DO jk = 2, jpkm1 217 zmxlm(:,:,jpk) = rn_lmin ! bottom set to the interior minium value 218 ! 219 !CDIR NOVERRCHK 220 DO jk = 2, jpkm1 ! interior value : l=sqrt(2*e/n**2) 221 221 !CDIR NOVERRCHK 222 222 DO jj = 2, jpjm1 … … 224 224 DO ji = fs_2, fs_jpim1 ! vector opt. 225 225 zrn2 = MAX( rn2(ji,jj,jk), rsmall ) 226 zmxlm(ji,jj,jk) = MAX( SQRT( 2. * en(ji,jj,jk) / zrn2 ), zmlmin)226 zmxlm(ji,jj,jk) = MAX( rn_lmin, SQRT( 2. * en(ji,jj,jk) / zrn2 ) ) 227 227 END DO 228 228 END DO … … 232 232 ! ------------------------------------- 233 233 zmxld(:,:, 1 ) = zmxlm(:,:,1) ! surface set to the minimum value 234 zmxld(:,:,jpk) = rn_lmin 235 236 SELECT CASE ( n mxl )237 234 zmxld(:,:,jpk) = rn_lmin ! bottom set to the minimum value 235 236 SELECT CASE ( nn_mxl ) 237 ! 238 238 CASE ( 0 ) ! bounded by the distance to surface and bottom 239 239 DO jk = 2, jpkm1 … … 247 247 END DO 248 248 END DO 249 249 ! 250 250 CASE ( 1 ) ! bounded by the vertical scale factor 251 251 DO jk = 2, jpkm1 … … 258 258 END DO 259 259 END DO 260 260 ! 261 261 CASE ( 2 ) ! |dk[xml]| bounded by e3t : 262 262 DO jk = 2, jpkm1 ! from the surface to the bottom : … … 276 276 END DO 277 277 END DO 278 278 ! 279 279 CASE ( 3 ) ! lup and ldown, |dk[xml]| bounded by e3t : 280 280 DO jk = 2, jpkm1 ! from the surface to the bottom : lup … … 305 305 END DO 306 306 END DO 307 307 ! 308 308 END SELECT 309 309 310 310 # if defined key_c1d 311 ! save mixing and dissipation turbulent length scales311 ! c1d configuration : save mixing and dissipation turbulent length scales 312 312 e_dis(:,:,:) = zmxld(:,:,:) 313 313 e_mix(:,:,:) = zmxlm(:,:,:) … … 346 346 ! 347 347 # if defined key_c1d 348 hlc(:,:) = zhlc(:,:) * tmask(:,:,1) ! save finite Langmuir Circulationdepth348 hlc(:,:) = zhlc(:,:) * tmask(:,:,1) ! c1d configuration: save finite Langmuir Circulation depth 349 349 # endif 350 350 ! … … 382 382 DO ji = fs_2, fs_jpim1 ! vector opt. 383 383 zsqen = SQRT( en(ji,jj,jk) ) 384 zav = ediff * zmxlm(ji,jj,jk) * zsqen384 zav = rn_ediff * zmxlm(ji,jj,jk) * zsqen 385 385 avt (ji,jj,jk) = MAX( zav, avtb_2d(ji,jj) * avtb(jk) ) * tmask(ji,jj,jk) 386 386 zmxlm(ji,jj,jk) = MAX( zav, avmb(jk) ) * tmask(ji,jj,jk) … … 392 392 ! 2. Surface boundary condition on tke and its eddy viscosity (zmxlm) 393 393 ! ------------------------------------------------- 394 ! en(1) = ebb sqrt(utau^2+vtau^2) / rau0 (min valueemin0)394 ! en(1) = rn_ebb sqrt(utau^2+vtau^2) / rau0 (min value rn_emin0) 395 395 ! zmxlm(1) = avmb(1) and zmxlm(jpk) = 0. 396 396 !CDIR NOVERRCHK … … 401 401 zty2 = vtau(ji ,jj-1) + vtau(ji,jj) 402 402 zesurf = zbbrau * SQRT( ztx2 * ztx2 + zty2 * zty2 ) 403 en (ji,jj,1) = MAX( zesurf, emin0 ) * tmask(ji,jj,1)404 zav = ediff * zmxlm(ji,jj,1) * SQRT( en(ji,jj,1) )* tmask(ji,jj,1)405 zmxlm(ji,jj,1 ) = MAX( zav, avmb(1)) * tmask(ji,jj,1)406 avt (ji,jj,1 ) = MAX( zav, avtb(1) * avtb_2d(ji,jj) ) * tmask(ji,jj,1)403 en (ji,jj,1) = MAX( zesurf, rn_emin0 ) * tmask(ji,jj,1) 404 zav = rn_ediff * zmxlm(ji,jj,1) * SQRT( en(ji,jj,1) ) 405 zmxlm(ji,jj,1 ) = MAX( zav, avmb(1) ) * tmask(ji,jj,1) 406 avt (ji,jj,1 ) = MAX( zav, avtb(1) * avtb_2d(ji,jj) ) * tmask(ji,jj,1) 407 407 zmxlm(ji,jj,jpk) = 0.e0 408 408 END DO … … 412 412 ! ------------------------------- 413 413 ! Resolution of a tridiagonal linear system by a "methode de chasse" 414 ! computation from level 2 to jpkm1 (e(1) already computed and 415 ! e(jpk)=0 ). 416 417 SELECT CASE ( npdl ) 418 414 ! computation from level 2 to jpkm1 (e(1) already computed and e(jpk)=0 ). 415 416 SELECT CASE ( nn_pdl ) 417 ! 419 418 CASE ( 0 ) ! No Prandtl number 420 419 DO jk = 2, jpkm1 421 420 DO jj = 2, jpjm1 422 421 DO ji = fs_2, fs_jpim1 ! vector opt. 423 ! zesh2 = eboost * (du/dz)^2 - N^2422 ! ! shear prod. - stratif. destruction 424 423 zcoef = 0.5 / fse3w(ji,jj,jk) 425 ! shear 426 zdku = zcoef * ( ub(ji-1, jj ,jk-1) + ub(ji,jj,jk-1) & 427 & - ub(ji-1, jj ,jk ) - ub(ji,jj,jk ) ) 428 zdkv = zcoef * ( vb( ji ,jj-1,jk-1) + vb(ji,jj,jk-1) & 429 & - vb( ji ,jj-1,jk ) - vb(ji,jj,jk ) ) 430 ! coefficient (zesh2) 431 zesh2 = eboost * ( zdku*zdku + zdkv*zdkv ) - rn2(ji,jj,jk) 432 433 ! Matrix 424 zdku = zcoef * ( ub(ji-1, jj ,jk-1) + ub(ji,jj,jk-1) & ! shear 425 & - ub(ji-1, jj ,jk ) - ub(ji,jj,jk ) ) 426 zdkv = zcoef * ( vb( ji ,jj-1,jk-1) + vb(ji,jj,jk-1) & 427 & - vb( ji ,jj-1,jk ) - vb(ji,jj,jk ) ) 428 zesh2 = eboost * ( zdku*zdku + zdkv*zdkv ) - rn2(ji,jj,jk) ! coefficient (zesh2) 429 ! 430 ! ! Matrix 434 431 zcof = zfact1 * tmask(ji,jj,jk) 435 ! lower diagonal432 ! ! lower diagonal 436 433 avmv(ji,jj,jk) = zcof * ( zmxlm(ji,jj,jk ) + zmxlm(ji,jj,jk-1) ) & 437 &/ ( fse3t(ji,jj,jk-1) * fse3w(ji,jj,jk ) )438 ! upper diagonal434 & / ( fse3t(ji,jj,jk-1) * fse3w(ji,jj,jk ) ) 435 ! ! upper diagonal 439 436 avmu(ji,jj,jk) = zcof * ( zmxlm(ji,jj,jk+1) + zmxlm(ji,jj,jk ) ) & 440 &/ ( fse3t(ji,jj,jk ) * fse3w(ji,jj,jk) )441 ! diagonal437 & / ( fse3t(ji,jj,jk ) * fse3w(ji,jj,jk) ) 438 ! ! diagonal 442 439 zwd(ji,jj,jk) = 1. - avmv(ji,jj,jk) - avmu(ji,jj,jk) + zfact2 * zmxld(ji,jj,jk) 443 ! right hand side in en 444 IF( .NOT. ln_lc ) THEN 440 ! 441 ! ! right hand side in en 442 IF( .NOT. ln_lc ) THEN ! No Langmuir cells 445 443 en(ji,jj,jk) = en(ji,jj,jk) + zfact3 * zmxld (ji,jj,jk) * en(ji,jj,jk) & 446 444 & + rdt * zmxlm (ji,jj,jk) * zesh2 447 ELSE 445 ELSE ! Langmuir cell source term 448 446 en(ji,jj,jk) = en(ji,jj,jk) + zfact3 * zmxld (ji,jj,jk) * en(ji,jj,jk) & 449 447 & + rdt * zmxlm (ji,jj,jk) * zesh2 & … … 453 451 END DO 454 452 END DO 455 453 ! 456 454 CASE ( 1 ) ! Prandtl number 457 455 DO jk = 2, jpkm1 458 456 DO jj = 2, jpjm1 459 457 DO ji = fs_2, fs_jpim1 ! vector opt. 460 ! zesh2 = eboost * (du/dz)^2 - pdl * N^2 461 zcoef = 0.5 / fse3w(ji,jj,jk) 462 ! shear 463 zdku = zcoef * ( ub(ji-1,jj ,jk-1) + ub(ji,jj,jk-1) & 464 & - ub(ji-1,jj ,jk ) - ub(ji,jj,jk ) ) 465 zdkv = zcoef * ( vb(ji ,jj-1,jk-1) + vb(ji,jj,jk-1) & 466 & - vb(ji ,jj-1,jk ) - vb(ji,jj,jk ) ) 467 ! square of vertical shear 468 zsh2 = zdku * zdku + zdkv * zdkv 469 ! local Richardson number 470 zri = MAX( rn2(ji,jj,jk), 0. ) / ( zsh2 + 1.e-20 ) 458 ! ! shear prod. - stratif. destruction 459 zcoef = 0.5 / fse3w(ji,jj,jk) 460 zdku = zcoef * ( ub(ji-1,jj ,jk-1) + ub(ji,jj,jk-1) & ! shear 461 & - ub(ji-1,jj ,jk ) - ub(ji,jj,jk ) ) 462 zdkv = zcoef * ( vb(ji ,jj-1,jk-1) + vb(ji,jj,jk-1) & 463 & - vb(ji ,jj-1,jk ) - vb(ji,jj,jk ) ) 464 zsh2 = zdku * zdku + zdkv * zdkv ! square of shear 465 zri = MAX( rn2(ji,jj,jk), 0. ) / ( zsh2 + 1.e-20 ) ! local Richardson number 471 466 # if defined key_c1d 472 ! save masked local Richardson number in zmxlm array 473 e_ric(ji,jj,jk) = zri * tmask(ji,jj,jk) 467 e_ric(ji,jj,jk) = zri * tmask(ji,jj,jk) ! c1d config. : save Ri 474 468 # endif 475 ! Prandtl number 476 zpdl = 1.0 477 IF( zri >= 0.2 ) zpdl = 0.2 / zri 469 zpdl = 1.0 ! Prandtl number 470 IF( zri >= 0.2 ) zpdl = 0.2 / zri 478 471 zpdl = MAX( 0.1, zpdl ) 479 ! coefficient (esh2) 480 zesh2 = eboost * zsh2 - zpdl * rn2(ji,jj,jk) 481 482 ! Matrix 472 zesh2 = eboost * zsh2 - zpdl * rn2(ji,jj,jk) ! coefficient (esh2) 473 ! 474 ! ! Matrix 483 475 zcof = zfact1 * tmask(ji,jj,jk) 484 ! lower diagonal476 ! ! lower diagonal 485 477 avmv(ji,jj,jk) = zcof * ( zmxlm(ji,jj,jk ) + zmxlm(ji,jj,jk-1) ) & 486 478 & / ( fse3t(ji,jj,jk-1) * fse3w(ji,jj,jk ) ) 487 ! upper diagonal479 ! ! upper diagonal 488 480 avmu(ji,jj,jk) = zcof * ( zmxlm(ji,jj,jk+1) + zmxlm(ji,jj,jk ) ) & 489 481 & / ( fse3t(ji,jj,jk ) * fse3w(ji,jj,jk) ) 490 ! diagonal482 ! ! diagonal 491 483 zwd(ji,jj,jk) = 1. - avmv(ji,jj,jk) - avmu(ji,jj,jk) + zfact2 * zmxld(ji,jj,jk) 492 ! right hand side in en 493 IF( .NOT. ln_lc ) THEN 484 ! 485 ! ! right hand side in en 486 IF( .NOT. ln_lc ) THEN ! No Langmuir cells 494 487 en(ji,jj,jk) = en(ji,jj,jk) + zfact3 * zmxld (ji,jj,jk) * en (ji,jj,jk) & 495 488 & + rdt * zmxlm (ji,jj,jk) * zesh2 496 ELSE 489 ELSE ! Langmuir cell source term 497 490 en(ji,jj,jk) = en(ji,jj,jk) + zfact3 * zmxld (ji,jj,jk) * en (ji,jj,jk) & 498 491 & + rdt * zmxlm (ji,jj,jk) * zesh2 & 499 492 & + rdt * ztkelc(ji,jj,jk) 500 493 ENDIF 501 ! save masked Prandlt number in zmxld array 502 zmxld(ji,jj,jk) = zpdl * tmask(ji,jj,jk) 503 END DO 504 END DO 505 END DO 506 494 zmxld(ji,jj,jk) = zpdl * tmask(ji,jj,jk) ! store masked Prandlt number in zmxld array 495 END DO 496 END DO 497 END DO 498 ! 507 499 END SELECT 508 500 509 501 # if defined key_c1d 510 ! save masked Prandlt number 511 e_pdl(:,:,2:jpkm1) = zmxld(:,:,2:jpkm1) 502 e_pdl(:,:,2:jpkm1) = zmxld(:,:,2:jpkm1) ! c1d configuration : save masked Prandlt number 512 503 e_pdl(:,:, 1) = e_pdl(:,:, 2) 513 504 e_pdl(:,:, jpk) = e_pdl(:,:, jpkm1) … … 516 507 ! 4. Matrix inversion from level 2 (tke prescribed at level 1) 517 508 !!-------------------------------- 518 ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 519 DO jk = 3, jpkm1 509 DO jk = 3, jpkm1 ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 520 510 DO jj = 2, jpjm1 521 511 DO ji = fs_2, fs_jpim1 ! vector opt. … … 524 514 END DO 525 515 END DO 526 527 ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1 528 DO jj = 2, jpjm1 516 DO jj = 2, jpjm1 ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1 529 517 DO ji = fs_2, fs_jpim1 ! vector opt. 530 518 avmv(ji,jj,2) = en(ji,jj,2) - avmv(ji,jj,2) * en(ji,jj,1) ! Surface boudary conditions on tke … … 538 526 END DO 539 527 END DO 540 541 ! thrid recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk 542 DO jj = 2, jpjm1 528 DO jj = 2, jpjm1 ! thrid recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk 543 529 DO ji = fs_2, fs_jpim1 ! vector opt. 544 530 en(ji,jj,jpkm1) = avmv(ji,jj,jpkm1) / zwd(ji,jj,jpkm1) … … 552 538 END DO 553 539 END DO 554 555 ! Save the result in en and set minimum value of tke : emin 556 DO jk = 2, jpkm1 540 DO jk = 2, jpkm1 ! set the minimum value of tke 557 541 DO jj = 2, jpjm1 558 542 DO ji = fs_2, fs_jpim1 ! vector opt. 559 en(ji,jj,jk) = MAX( en(ji,jj,jk), emin ) * tmask(ji,jj,jk)543 en(ji,jj,jk) = MAX( en(ji,jj,jk), rn_emin ) * tmask(ji,jj,jk) 560 544 END DO 561 545 END DO 562 546 END DO 563 547 564 ! Modify the value of TKE by an exponential assumption 565 ! en(ji,jj,jk)=en(ji,jj,jk)+en(ji,jj,1)*EXP(-fsdepw(ji,jj,jk)/ zhtau(ji,jj) ) 566 567 SELECT CASE( nn_htau ) ! Choice of H value 568 ! 569 CASE( 0 ) 570 DO jj = 2, jpjm1 571 DO ji = fs_2, fs_jpim1 ! vector opt. 572 zhtau(ji,jj) = 10. 573 END DO 574 END DO 575 ! 576 CASE( 1 ) 577 DO jj = 2, jpjm1 578 DO ji = fs_2, fs_jpim1 ! vector opt. 579 zhtau(ji,jj) = MAX( 5., MIN( 40., 60.*ABS( SIN( rpi/180. * gphit(ji,jj) ) ) ) ) 580 END DO 581 END DO 582 ! 583 CASE( 2 ) 584 DO jj = 2, jpjm1 585 DO ji = fs_2, fs_jpim1 ! vector opt. 586 zhtau(ji,jj) = MAX( 5.,6./4.* MIN( 40., 60.*ABS( SIN( rpi/180. * gphit(ji,jj) ) ) ) ) 587 END DO 588 END DO 589 ! 590 END SELECT 591 592 IF( nn_etau == 1 ) THEN 593 DO jk = 2, jpkm1 594 DO jj = 2, jpjm1 595 DO ji = fs_2, fs_jpim1 ! vector opt. 596 en(ji,jj,jk) = en(ji,jj,jk) & 548 ! 5. Add extra TKE due to surface and internal wave breaking (nn_etau /= 0) 549 !!---------------------------------------------------------- 550 IF( nn_etau /= 0 ) THEN ! extra tke : en = en + rn_efr * en(1) * exp( -z/zhtau ) 551 ! 552 SELECT CASE( nn_htau ) ! Choice of the depth of penetration 553 CASE( 0 ) ! constant depth penetration (here 10 meters) 554 DO jj = 2, jpjm1 555 DO ji = fs_2, fs_jpim1 ! vector opt. 556 zhtau(ji,jj) = 10. 557 END DO 558 END DO 559 CASE( 1 ) ! meridional profile 1 560 DO jj = 2, jpjm1 ! ( 5m in the tropics to a maximum of 40 m at high lat.) 561 DO ji = fs_2, fs_jpim1 ! vector opt. 562 zhtau(ji,jj) = MAX( 5., MIN( 40., 60.*ABS( SIN( rpi/180. * gphit(ji,jj) ) ) ) ) 563 END DO 564 END DO 565 CASE( 2 ) ! meridional profile 2 566 DO jj = 2, jpjm1 ! ( 5m in the tropics to a maximum of 60 m at high lat.) 567 DO ji = fs_2, fs_jpim1 ! vector opt. 568 zhtau(ji,jj) = MAX( 5.,6./4.* MIN( 40., 60.*ABS( SIN( rpi/180. * gphit(ji,jj) ) ) ) ) 569 END DO 570 END DO 571 END SELECT 572 ! 573 IF( nn_etau == 1 ) THEN ! extra term throughout the water column 574 DO jk = 2, jpkm1 575 DO jj = 2, jpjm1 576 DO ji = fs_2, fs_jpim1 ! vector opt. 577 en(ji,jj,jk) = en(ji,jj,jk) & 578 & + rn_efr * en(ji,jj,1)*EXP( -fsdepw(ji,jj,jk) / zhtau(ji,jj) ) & 579 & * ( 1.e0 - fr_i(ji,jj) ) * tmask(ji,jj,jk) 580 END DO 581 END DO 582 END DO 583 ELSEIF( nn_etau == 2 ) THEN ! extra term only at the base of the mixed layer 584 DO jj = 2, jpjm1 585 DO ji = fs_2, fs_jpim1 ! vector opt. 586 jk = nmln(ji,jj) 587 en(ji,jj,jk) = en(ji,jj,jk) & 597 588 & + rn_efr * en(ji,jj,1)*EXP( -fsdepw(ji,jj,jk) / zhtau(ji,jj) ) & 598 & * ( 1.e0 - fr_i(ji,jj) ) * tmask(ji,jj,jk) 599 END DO 600 END DO 601 END DO 602 ELSEIF( nn_etau == 2 ) THEN ! only at the base of the mixed layer 603 DO jj = 2, jpjm1 604 DO ji = fs_2, fs_jpim1 ! vector opt. 605 jk = nmln(ji,jj) 606 en(ji,jj,jk) = en(ji,jj,jk) & 607 & + rn_efr * en(ji,jj,1)*EXP( -fsdepw(ji,jj,jk) / zhtau(ji,jj) ) & 608 & * ( 1.e0 - fr_i(ji,jj) ) * tmask(ji,jj,jk) 609 END DO 610 END DO 589 & * ( 1.e0 - fr_i(ji,jj) ) * tmask(ji,jj,jk) 590 END DO 591 END DO 592 ENDIF 593 ! 611 594 ENDIF 612 595 613 596 614 ! Lateral boundary conditions on ( avt, en ) (sign unchanged) 615 CALL lbc_lnk( en , 'W', 1. ) ; CALL lbc_lnk( avt, 'W', 1. ) 597 ! Lateral boundary conditions (sign unchanged) 598 CALL lbc_lnk( en, 'W', 1. ) ; CALL lbc_lnk( avt, 'W', 1. ) 599 616 600 617 601 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 618 602 ! IV. Before vertical eddy vicosity and diffusivity coefficients 619 603 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 620 621 SELECT CASE ( nave ) 622 623 CASE ( 0 ) ! no horizontal average 624 625 ! Vertical eddy viscosity 626 627 DO jk = 2, jpkm1 ! Horizontal slab 604 ! 605 SELECT CASE ( nn_ave ) 606 CASE ( 0 ) ! no horizontal average 607 DO jk = 2, jpkm1 ! only vertical eddy viscosity 628 608 DO jj = 2, jpjm1 629 609 DO ji = fs_2, fs_jpim1 ! vector opt. … … 635 615 END DO 636 616 END DO 637 638 ! Lateral boundary conditions (avmu,avmv) (U- and V- points, sign unchanged) 639 CALL lbc_lnk( avmu, 'U', 1. ) ; CALL lbc_lnk( avmv, 'V', 1. ) 640 641 CASE ( 1 ) ! horizontal average 642 643 ! ( 1/2 1/2 ) 644 ! Eddy viscosity: horizontal average: avmu = 1/4 ( 1 1 ) 645 ! ( 1/2 1 1/2 ) ( 1/2 1/2 ) 646 ! avmv = 1/4 ( 1/2 1 1/2 ) 647 648 !! caution vectopt_memory change the solution (last digit of the solver stat) 617 CALL lbc_lnk( avmu, 'U', 1. ) ; CALL lbc_lnk( avmv, 'V', 1. ) ! Lateral boundary conditions 618 ! 619 CASE ( 1 ) ! horizontal average ( 1/2 1/2 ) 620 ! ! Vertical eddy viscosity avmu = 1/4 ( 1 1 ) 621 ! ! ( 1/2 1/2 ) 622 ! ! 623 ! ! ( 1/2 1 1/2 ) 624 ! ! avmv = 1/4 ( 1/2 1 1/2 ) 625 DO jk = 2, jpkm1 626 DO jj = 2, jpjm1 627 DO ji = fs_2, fs_jpim1 ! vector opt. 628 # if defined key_vectopt_memory 629 ! ! caution : vectopt_memory change the solution 630 ! ! (last digit of the solver stat) 631 avmu(ji,jj,jk) = ( avt(ji,jj ,jk) + avt(ji+1,jj ,jk) & 632 & +.5*( avt(ji,jj-1,jk) + avt(ji+1,jj-1,jk) & 633 & +avt(ji,jj+1,jk) + avt(ji+1,jj+1,jk) ) ) * eumean(ji,jj,jk) 634 635 avmv(ji,jj,jk) = ( avt(ji ,jj,jk) + avt(ji ,jj+1,jk) & 636 & +.5*( avt(ji-1,jj,jk) + avt(ji-1,jj+1,jk) & 637 & +avt(ji+1,jj,jk) + avt(ji+1,jj+1,jk) ) ) * evmean(ji,jj,jk) 638 # else 639 avmu(ji,jj,jk) = ( avt (ji,jj ,jk) + avt (ji+1,jj ,jk) & 640 & +.5*( avt (ji,jj-1,jk) + avt (ji+1,jj-1,jk) & 641 & +avt (ji,jj+1,jk) + avt (ji+1,jj+1,jk) ) ) * umask(ji,jj,jk) & 642 & / MAX( 1., tmask(ji,jj ,jk) + tmask(ji+1,jj ,jk) & 643 & +.5*( tmask(ji,jj-1,jk) + tmask(ji+1,jj-1,jk) & 644 & +tmask(ji,jj+1,jk) + tmask(ji+1,jj+1,jk) ) ) 645 646 avmv(ji,jj,jk) = ( avt (ji ,jj,jk) + avt (ji ,jj+1,jk) & 647 & +.5*( avt (ji-1,jj,jk) + avt (ji-1,jj+1,jk) & 648 & +avt (ji+1,jj,jk) + avt (ji+1,jj+1,jk) ) ) * vmask(ji,jj,jk) & 649 & / MAX( 1., tmask(ji ,jj,jk) + tmask(ji ,jj+1,jk) & 650 & +.5*( tmask(ji-1,jj,jk) + tmask(ji-1,jj+1,jk) & 651 & +tmask(ji+1,jj,jk) + tmask(ji+1,jj+1,jk) ) ) 652 # endif 653 END DO 654 END DO 655 END DO 656 CALL lbc_lnk( avmu, 'U', 1. ) ; CALL lbc_lnk( avmv, 'V', 1. ) ! Lateral boundary conditions 657 ! 658 ! ! Vertical eddy diffusivity (1 2 1) 659 ! ! avt = 1/16 (2 4 2) 660 ! ! (1 2 1) 661 DO jk = 2, jpkm1 662 DO jj = 2, jpjm1 663 DO ji = fs_2, fs_jpim1 ! vector opt. 649 664 # if defined key_vectopt_memory 650 DO jk = 2, jpkm1 ! Horizontal slab 651 DO jj = 2, jpjm1 652 DO ji = fs_2, fs_jpim1 ! vector opt. 653 avmu(ji,jj,jk) = ( avt(ji,jj ,jk) + avt(ji+1,jj ,jk) & 654 & +.5*( avt(ji,jj-1,jk) + avt(ji+1,jj-1,jk) & 655 & +avt(ji,jj+1,jk) + avt(ji+1,jj+1,jk) ) ) * eumean(ji,jj,jk) 656 657 avmv(ji,jj,jk) = ( avt(ji ,jj,jk) + avt(ji ,jj+1,jk) & 658 & +.5*( avt(ji-1,jj,jk) + avt(ji-1,jj+1,jk) & 659 & +avt(ji+1,jj,jk) + avt(ji+1,jj+1,jk) ) ) * evmean(ji,jj,jk) 660 END DO 661 END DO 662 END DO 665 avt(ji,jj,jk) = ( avmu(ji,jj,jk) + avmu(ji-1,jj ,jk) & 666 & + avmv(ji,jj,jk) + avmv(ji ,jj-1,jk) ) * etmean(ji,jj,jk) 663 667 # else 664 DO jk = 2, jpkm1 ! Horizontal slab 665 DO jj = 2, jpjm1 666 DO ji = fs_2, fs_jpim1 ! vector opt. 667 avmu(ji,jj,jk) = ( avt (ji,jj ,jk) + avt (ji+1,jj ,jk) & 668 & +.5*( avt (ji,jj-1,jk) + avt (ji+1,jj-1,jk) & 669 & +avt (ji,jj+1,jk) + avt (ji+1,jj+1,jk) ) ) * umask(ji,jj,jk) & 670 & / MAX( 1., tmask(ji,jj ,jk) + tmask(ji+1,jj ,jk) & 671 & +.5*( tmask(ji,jj-1,jk) + tmask(ji+1,jj-1,jk) & 672 & +tmask(ji,jj+1,jk) + tmask(ji+1,jj+1,jk) ) ) 673 674 avmv(ji,jj,jk) = ( avt (ji ,jj,jk) + avt (ji ,jj+1,jk) & 675 & +.5*( avt (ji-1,jj,jk) + avt (ji-1,jj+1,jk) & 676 & +avt (ji+1,jj,jk) + avt (ji+1,jj+1,jk) ) ) * vmask(ji,jj,jk) & 677 & / MAX( 1., tmask(ji ,jj,jk) + tmask(ji ,jj+1,jk) & 678 & +.5*( tmask(ji-1,jj,jk) + tmask(ji-1,jj+1,jk) & 679 & +tmask(ji+1,jj,jk) + tmask(ji+1,jj+1,jk) ) ) 680 END DO 681 END DO 682 END DO 668 avt(ji,jj,jk) = ( avmu (ji,jj,jk) + avmu (ji-1,jj ,jk) & 669 & + avmv (ji,jj,jk) + avmv (ji ,jj-1,jk) ) * tmask(ji,jj,jk) & 670 & / MAX( 1., umask(ji,jj,jk) + umask(ji-1,jj ,jk) & 671 & + vmask(ji,jj,jk) + vmask(ji ,jj-1,jk) ) 683 672 # endif 684 685 ! Lateral boundary conditions (avmu,avmv) (sign unchanged) 686 CALL lbc_lnk( avmu, 'U', 1. ) ; CALL lbc_lnk( avmv, 'V', 1. ) 687 688 ! Vertical eddy diffusivity 689 ! ------------------------------ 690 ! (1 2 1) 691 ! horizontal average avt = 1/16 (2 4 2) 692 ! (1 2 1) 693 DO jk = 2, jpkm1 ! Horizontal slab 694 # if defined key_vectopt_memory 695 DO jj = 2, jpjm1 696 DO ji = fs_2, fs_jpim1 ! vector opt. 697 avt(ji,jj,jk) = ( avmu(ji,jj,jk) + avmu(ji-1,jj ,jk) & 698 & + avmv(ji,jj,jk) + avmv(ji ,jj-1,jk) ) * etmean(ji,jj,jk) 699 END DO 700 END DO 701 # else 702 DO jj = 2, jpjm1 703 DO ji = fs_2, fs_jpim1 ! vector opt. 704 avt(ji,jj,jk) = ( avmu (ji,jj,jk) + avmu (ji-1,jj ,jk) & 705 & + avmv (ji,jj,jk) + avmv (ji ,jj-1,jk) ) * tmask(ji,jj,jk) & 706 & / MAX( 1., umask(ji,jj,jk) + umask(ji-1,jj ,jk) & 707 & + vmask(ji,jj,jk) + vmask(ji ,jj-1,jk) ) 708 END DO 709 END DO 710 # endif 711 END DO 712 673 END DO 674 END DO 675 END DO 676 ! 713 677 END SELECT 714 715 ! multiplied by the Prandtl number (npdl>1) 716 ! ---------------------------------------- 717 IF( npdl == 1 ) THEN 718 DO jk = 2, jpkm1 ! Horizontal slab 678 ! 679 IF( nn_pdl == 1 ) THEN ! Ponderation by the Prandtl number (nn_pdl=1) 680 DO jk = 2, jpkm1 719 681 DO jj = 2, jpjm1 720 682 DO ji = fs_2, fs_jpim1 ! vector opt. … … 725 687 END DO 726 688 ENDIF 727 728 ! Minimum value on the eddy viscosity 729 ! ---------------------------------------- 730 DO jk = 2, jpkm1 ! Horizontal slab 731 DO jj = 1, jpj 732 DO ji = 1, jpi 733 avmu(ji,jj,jk) = MAX( avmu(ji,jj,jk), avmb(jk) ) * umask(ji,jj,jk) 734 avmv(ji,jj,jk) = MAX( avmv(ji,jj,jk), avmb(jk) ) * vmask(ji,jj,jk) 735 END DO 736 END DO 737 END DO 738 739 ! Lateral boundary conditions on avt (sign unchanged) 740 ! ------------------------------===== 741 CALL lbc_lnk( avt, 'W', 1. ) 742 743 ! write en in restart file 744 ! ------------------------ 745 IF( lrst_oce ) CALL tke_rst( kt, 'WRITE' ) 689 ! 690 DO jk = 2, jpkm1 ! Minimum value on the eddy viscosity 691 avmu(:,:,jk) = MAX( avmu(:,:,jk), avmb(jk) ) * umask(:,:,jk) 692 avmv(:,:,jk) = MAX( avmv(:,:,jk), avmb(jk) ) * vmask(:,:,jk) 693 END DO 694 695 CALL lbc_lnk( avt, 'W', 1. ) ! Lateral boundary conditions on avt (sign unchanged) 696 697 IF( lrst_oce ) CALL tke_rst( kt, 'WRITE' ) ! write en in restart file 746 698 747 699 IF(ln_ctl) THEN 748 CALL prt_ctl( tab3d_1=en , clinfo1=' tke - e: ', tab3d_2=avt, clinfo2=' t: ', ovlap=1, kdim=jpk)749 CALL prt_ctl( tab3d_1=avmu, clinfo1=' tke - u: ', mask1=umask,&750 & tab3d_2=avmv, clinfo2= ' v: ', mask2=vmask, ovlap=1, kdim=jpk)700 CALL prt_ctl( tab3d_1=en , clinfo1=' tke - e: ', tab3d_2=avt, clinfo2=' t: ', ovlap=1, kdim=jpk) 701 CALL prt_ctl( tab3d_1=avmu, clinfo1=' tke - u: ', mask1=umask, & 702 & tab3d_2=avmv, clinfo2= ' v: ', mask2=vmask, ovlap=1, kdim=jpk ) 751 703 ENDIF 752 704 ! 753 705 END SUBROUTINE zdf_tke 754 706 … … 773 725 ! 774 726 # if defined key_vectopt_memory 775 ! caution vectopt_memory change the solution (last digit of the solver stat)776 727 INTEGER :: ji, jj, jk ! dummy loop indices 777 728 # else … … 779 730 # endif 780 731 !! 781 NAMELIST/namtke/ ln_rstke, ediff, ediss, ebb, efave, emin, emin0, &782 ri_c, nitke, nmxl, npdl, nave, navb,&783 ln_lsfc, rn_lmin, nn_avtb_2d, nn_etau, nn_htau, rn_efr,&784 ln_lc, rn_lc732 NAMELIST/namtke/ ln_rstke, rn_ediff, rn_ediss, rn_ebb , rn_efave, rn_emin, & 733 & rn_emin0, rn_cri , nn_itke , nn_mxl , nn_pdl , nn_ave , & 734 & nn_avb , ln_mxl0 , rn_lmin , rn_lmin0, nn_havtb, nn_etau, & 735 & nn_htau , rn_efr , ln_lc , rn_lc 785 736 !!---------------------------------------------------------------------- 786 737 … … 791 742 792 743 ! Compute boost associated with the Richardson critic 793 ! (control values: ri_c = 0.3 ==> eboost=1.25 for npdl=1 or 2) 794 ! ( ri_c = 0.222 ==> eboost=1. ) 795 eboost = ri_c * ( 2. + ediss / ediff ) / 2. 744 ! (control values: rn_cri = 0.3 ==> eboost=1.25 for nn_pdl=1) 745 ! ( rn_cri = 0.222 ==> eboost=1. ) 746 eboost = rn_cri * ( 2. + rn_ediss / rn_ediff ) / 2. 747 796 748 797 749 … … 804 756 WRITE(numout,*) '~~~~~~~~~~~~' 805 757 WRITE(numout,*) ' Namelist namtke : set tke mixing parameters' 806 WRITE(numout,*) ' restart with tke from no tke ln_rstke = ', ln_rstke 807 WRITE(numout,*) ' coef. to compute avt ediff = ', ediff 808 WRITE(numout,*) ' Kolmogoroff dissipation coef. ediss = ', ediss 809 WRITE(numout,*) ' tke surface input coef. ebb = ', ebb 810 WRITE(numout,*) ' tke diffusion coef. efave = ', efave 811 WRITE(numout,*) ' minimum value of tke emin = ', emin 812 WRITE(numout,*) ' surface minimum value of tke emin0 = ', emin0 813 WRITE(numout,*) ' number of restart iter loops nitke = ', nitke 814 WRITE(numout,*) ' mixing length type nmxl = ', nmxl 815 WRITE(numout,*) ' prandl number flag npdl = ', npdl 816 WRITE(numout,*) ' horizontal average flag nave = ', nave 817 WRITE(numout,*) ' critic Richardson nb ri_c = ', ri_c 818 WRITE(numout,*) ' and its associated coeff. eboost = ', eboost 819 WRITE(numout,*) ' constant background or profile navb = ', navb 820 WRITE(numout,*) ' flag for compu.of bls as fun. of wind ln_lsfc = ', ln_lsfc 821 WRITE(numout,*) ' buoyancy lenght scale surface min value rn_lmin = ', rn_lmin 822 WRITE(numout,*) ' horizontal variation for avtb nn_avtb_2d = ', nn_avtb_2d 823 WRITE(numout,*) ' test param. to add tke induced by wind nn_etau = ', nn_etau 824 WRITE(numout,*) ' flag for computation of exp. tke profile nn_htau = ', nn_htau 825 WRITE(numout,*) ' % of emin which pene. the thermocline rn_efr = ', rn_efr 826 WRITE(numout,*) ' flag to take into acc. Langmuir circ. ln_lc = ', ln_lc 827 WRITE(numout,*) ' coef to compute verticla velocity of LC rn_lc = ', rn_lc 758 WRITE(numout,*) ' restart with tke from no tke ln_rstke = ', ln_rstke 759 WRITE(numout,*) ' coef. to compute avt rn_ediff = ', rn_ediff 760 WRITE(numout,*) ' Kolmogoroff dissipation coef. rn_ediss = ', rn_ediss 761 WRITE(numout,*) ' tke surface input coef. rn_ebb = ', rn_ebb 762 WRITE(numout,*) ' tke diffusion coef. rn_efave = ', rn_efave 763 WRITE(numout,*) ' minimum value of tke rn_emin = ', rn_emin 764 WRITE(numout,*) ' surface minimum value of tke rn_emin0 = ', rn_emin0 765 WRITE(numout,*) ' number of restart iter loops nn_itke = ', nn_itke 766 WRITE(numout,*) ' mixing length type nn_mxl = ', nn_mxl 767 WRITE(numout,*) ' prandl number flag nn_pdl = ', nn_pdl 768 WRITE(numout,*) ' horizontal average flag nn_ave = ', nn_ave 769 WRITE(numout,*) ' critic Richardson nb rn_cri = ', rn_cri 770 WRITE(numout,*) ' and its associated coeff. eboost = ', eboost 771 WRITE(numout,*) ' constant background or profile nn_avb = ', nn_avb 772 WRITE(numout,*) ' surface mixing length = F(stress) or not ln_mxl0 = ', ln_mxl0 773 WRITE(numout,*) ' surface mixing length minimum value rn_lmin0 = ', rn_lmin0 774 WRITE(numout,*) ' interior mixing length minimum value rn_lmin0 = ', rn_lmin0 775 WRITE(numout,*) ' horizontal variation for avtb nn_havtb = ', nn_havtb 776 WRITE(numout,*) ' test param. to add tke induced by wind nn_etau = ', nn_etau 777 WRITE(numout,*) ' flag for computation of exp. tke profile nn_htau = ', nn_htau 778 WRITE(numout,*) ' % of rn_emin0 which pene. the thermocline rn_efr = ', rn_efr 779 WRITE(numout,*) ' flag to take into acc. Langmuir circ. ln_lc = ', ln_lc 780 WRITE(numout,*) ' coef to compute verticla velocity of LC rn_lc = ', rn_lc 828 781 WRITE(numout,*) 829 782 ENDIF 830 783 831 ! Check nmxl and npdl values 832 IF( nmxl < 0 .OR. nmxl > 3 ) CALL ctl_stop( ' bad flag: nmxl is < 0 or > 3 ' ) 833 IF( npdl < 0 .OR. npdl > 1 ) CALL ctl_stop( ' bad flag: npdl is < 0 or > 1 ' ) 834 IF( nn_htau < 0 .OR. nn_htau > 2 ) CALL ctl_stop( 'bad flag: nn_htau is < 0 or > 3 ' ) 835 IF( rn_lc < 0.15 .OR. rn_lc > 0.2 ) CALL ctl_stop( 'bad value: rn_lc must be > 0.15 and < 0.2 ' ) 784 ! Check of some namelist values 785 IF( nn_mxl < 0 .OR. nn_mxl > 3 ) CALL ctl_stop( 'bad flag: nn_mxl is 0, 1 or 2 ' ) 786 IF( nn_pdl < 0 .OR. nn_pdl > 1 ) CALL ctl_stop( 'bad flag: nn_pdl is 0 or 1 ' ) 787 IF( nn_ave < 0 .OR. nn_ave > 1 ) CALL ctl_stop( 'bad flag: nn_ave is 0 or 1 ' ) 788 IF( nn_htau < 0 .OR. nn_htau > 2 ) CALL ctl_stop( 'bad flag: nn_htau is 0, 1 or 2 ' ) 789 IF( rn_lc < 0.15 .OR. rn_lc > 0.2 ) CALL ctl_stop( 'bad value: rn_lc must be between 0.15 and 0.2 ' ) 836 790 837 791 IF( nn_etau == 2 ) CALL zdf_mxl( nit000 ) ! Initialization of nmln … … 839 793 ! Horizontal average : initialization of weighting arrays 840 794 ! ------------------- 841 842 SELECT CASE ( nave ) 843 795 ! 796 SELECT CASE ( nn_ave ) 797 ! Notice: mean arrays have not to by defined at local domain boundaries due to the vertical nature of TKE 798 ! 844 799 CASE ( 0 ) ! no horizontal average 845 800 IF(lwp) WRITE(numout,*) ' no horizontal average on avt, avmu, avmv' … … 854 809 eumean(:,:,:) = 0.e0 855 810 evmean(:,:,:) = 0.e0 856 811 ! 857 812 DO jk = 1, jpkm1 858 813 DO jj = 2, jpjm1 859 814 DO ji = fs_2, fs_jpim1 ! vector opt. 860 etmean(ji,jj,jk) = tmask(ji,jj,jk) & 861 & / MAX( 1., umask(ji-1,jj ,jk) + umask(ji,jj,jk) & 862 & + vmask(ji ,jj-1,jk) + vmask(ji,jj,jk) ) 863 864 eumean(ji,jj,jk) = umask(ji,jj,jk) & 865 & / MAX( 1., tmask(ji,jj,jk) + tmask(ji+1,jj ,jk) ) 866 867 evmean(ji,jj,jk) = vmask(ji,jj,jk) & 868 & / MAX( 1., tmask(ji,jj,jk) + tmask(ji ,jj+1,jk) ) 815 etmean(ji,jj,jk) = tmask(ji,jj,jk) / MAX( 1., umask(ji,jj,jk) + umask(ji-1,jj ,jk) & 816 & + vmask(ji,jj,jk) + vmask(ji ,jj-1,jk) ) 817 eumean(ji,jj,jk) = umask(ji,jj,jk) / MAX( 1., tmask(ji,jj,jk) + tmask(ji+1,jj ,jk) ) 818 evmean(ji,jj,jk) = vmask(ji,jj,jk) / MAX( 1., tmask(ji,jj,jk) + tmask(ji ,jj+1,jk) ) 869 819 END DO 870 820 END DO 871 821 END DO 872 822 # endif 873 823 ! 874 824 CASE ( 1 ) ! horizontal average 875 825 IF(lwp) WRITE(numout,*) ' horizontal average on avt, avmu, avmv' … … 883 833 eumean(:,:,:) = 0.e0 884 834 evmean(:,:,:) = 0.e0 885 835 ! 886 836 DO jk = 1, jpkm1 887 837 DO jj = 2, jpjm1 888 838 DO ji = fs_2, fs_jpim1 ! vector opt. 889 etmean(ji,jj,jk) = tmask(ji,jj,jk) & 890 & / MAX( 1., umask(ji-1,jj ,jk) + umask(ji,jj,jk) & 891 & + vmask(ji ,jj-1,jk) + vmask(ji,jj,jk) ) 892 893 eumean(ji,jj,jk) = umask(ji,jj,jk) & 894 & / MAX( 1., tmask(ji,jj ,jk) + tmask(ji+1,jj ,jk) & 895 & +.5 * ( tmask(ji,jj-1,jk) + tmask(ji+1,jj-1,jk) & 896 & +tmask(ji,jj+1,jk) + tmask(ji+1,jj+1,jk) ) ) 897 898 evmean(ji,jj,jk) = vmask(ji,jj,jk) & 899 & / MAX( 1., tmask(ji ,jj,jk) + tmask(ji ,jj+1,jk) & 900 & +.5 * ( tmask(ji-1,jj,jk) + tmask(ji-1,jj+1,jk) & 901 & +tmask(ji+1,jj,jk) + tmask(ji+1,jj+1,jk) ) ) 839 etmean(ji,jj,jk) = tmask(ji,jj,jk) / MAX( 1., umask(ji-1,jj ,jk) + umask(ji ,jj ,jk) & 840 & + vmask(ji ,jj-1,jk) + vmask(ji ,jj ,jk) ) 841 eumean(ji,jj,jk) = umask(ji,jj,jk) / MAX( 1., tmask(ji ,jj ,jk) + tmask(ji+1,jj ,jk) & 842 & +.5 * ( tmask(ji ,jj-1,jk) + tmask(ji+1,jj-1,jk) & 843 & + tmask(ji ,jj+1,jk) + tmask(ji+1,jj+1,jk) ) ) 844 evmean(ji,jj,jk) = vmask(ji,jj,jk) / MAX( 1., tmask(ji ,jj ,jk) + tmask(ji ,jj+1,jk) & 845 & +.5 * ( tmask(ji-1,jj ,jk) + tmask(ji-1,jj+1,jk) & 846 & + tmask(ji+1,jj ,jk) + tmask(ji+1,jj+1,jk) ) ) 902 847 END DO 903 848 END DO 904 849 END DO 905 850 # endif 906 907 CASE DEFAULT 908 WRITE(ctmp1,*) ' bad flag value for nave = ', nave 909 CALL ctl_stop( ctmp1 ) 910 851 ! 911 852 END SELECT 912 853 … … 914 855 ! Background eddy viscosity and diffusivity profil 915 856 ! ------------------------------------------------ 916 IF( n avb == 0 ) THEN ! Define avmb, avtb from namelist parameter857 IF( nn_avb == 0 ) THEN ! Define avmb, avtb from namelist parameter 917 858 avmb(:) = avm0 918 859 avtb(:) = avt0 919 860 ELSE ! Background profile of avt (fit a theoretical/observational profile (Krauss 1990) 920 861 avmb(:) = avm0 921 !!bug this is not valide neither in scoord922 IF(ln_sco .AND. lwp) WRITE(numout,cform_war)923 IF(ln_sco .AND. lwp) WRITE(numout,*) ' avtb profile not valid in sco'924 925 862 avtb(:) = avt0 + ( 3.0e-4 - 2 * avt0 ) * 1.0e-4 * gdepw_0(:) ! m2/s 863 IF(ln_sco .AND. lwp) CALL ctl_warn( ' avtb profile not valid in sco' ) 926 864 ENDIF 927 865 ! 928 866 ! ! 2D shape of the avtb 929 867 avtb_2d(:,:) = 1.e0 ! uniform 930 868 ! 931 IF( nn_ avtb_2d== 1 ) THEN ! decrease avtb in the equatorial band869 IF( nn_havtb == 1 ) THEN ! decrease avtb in the equatorial band 932 870 ! -15S -5S : linear decrease from avt0 to avt0/10. 933 871 ! -5S +5N : cst value avt0/10. … … 938 876 ENDIF 939 877 878 !!gm the increase is only required when using cen2 advection scheme 879 !!gm for the equatorial upwelling. ===> use of TVD in ORCA so this have to be suppressed 940 880 ! Increase the background in the surface layers 941 !!gm the increase is only required when using cen2 advection scheme942 !!gm for the equatorial upwelling.943 881 avmb(1) = 10. * avmb(1) ; avtb(1) = 10. * avtb(1) 944 882 avmb(2) = 10. * avmb(2) ; avtb(2) = 10. * avtb(2) 945 883 avmb(3) = 5. * avmb(3) ; avtb(3) = 5. * avtb(3) 946 884 avmb(4) = 2.5 * avmb(4) ; avtb(4) = 2.5 * avtb(4) 885 !!gm end 947 886 948 887 … … 967 906 !! *** ROUTINE ts_rst *** 968 907 !! 969 !! ** Purpose : Read or write filtered free surface arraysin restart file908 !! ** Purpose : Read or write TKE file (en) in restart file 970 909 !! 971 !! ** Method : 972 !! 910 !! ** Method : use of IOM library 911 !! if the restart does not contain TKE, en is either 912 !! set to rn_emin or recomputed (nn_itke/=0) 973 913 !!---------------------------------------------------------------------- 974 914 INTEGER , INTENT(in) :: kt ! ocean time-step … … 986 926 & WRITE(numout,*) ' ===>>>> : previous run without tke scheme' 987 927 IF( lwp .AND. ln_rstke ) WRITE(numout,*) ' ===>>>> : We do not use en from the restart file' 988 IF( lwp ) WRITE(numout,*) ' ===>>>> : en set by iterative loop' 989 IF( lwp ) WRITE(numout,*) ' ======= =========' 990 en (:,:,:) = emin * tmask(:,:,:) 991 DO jit = 2, nitke+1 928 IF( lwp ) WRITE(numout,*) ' ===>>>> : en is set by iterative loop' 929 en (:,:,:) = rn_emin * tmask(:,:,:) 930 DO jit = 2, nn_itke + 1 992 931 CALL zdf_tke( jit ) 993 932 END DO 994 933 ENDIF 995 934 ELSE 996 en(:,:,:) = emin * tmask(:,:,:) ! no restart: en set to emin935 en(:,:,:) = rn_emin * tmask(:,:,:) ! no restart: en set to its minimum value 997 936 ENDIF 998 937 ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN
Note: See TracChangeset
for help on using the changeset viewer.
