Changeset 1492 for trunk/NEMO/OPA_SRC/ZDF/zdftke2.F90
- Timestamp:
- 2009-07-17T11:48:07+02:00 (15 years ago)
- File:
-
- 1 edited
-
trunk/NEMO/OPA_SRC/ZDF/zdftke2.F90 (modified) (26 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMO/OPA_SRC/ZDF/zdftke2.F90
r1481 r1492 5 5 !! turbulent closure parameterization 6 6 !!===================================================================== 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_tke2_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: 20 !! - tke penetration (wind steering) 21 !! - suface condition for tke & mixing length 22 !! - Langmuir cells 23 !! - ! 2008-05 (J.-M. Molines, G. Madec) 2D form of avtb 24 !! - ! 2008-06 (G. Madec) style + DOCTOR name for namelist parameters 25 !! - ! 2008-12 (G. Reffray) stable discretization of the production term 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 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: 20 !! ! - tke penetration (wind steering) 21 !! ! - suface condition for tke & mixing length 22 !! ! - Langmuir cells 23 !! - ! 2008-05 (J.-M. Molines, G. Madec) 2D form of avtb 24 !! - ! 2008-06 (G. Madec) style + DOCTOR name for namelist parameters 25 !! - ! 2008-12 (G. Reffray) stable discretization of the production term 26 !! 3.2 ! 2009-06 (G. Madec, S. Masson) TKE restart compatible with key_cpl 27 !! ! + cleaning of the parameters + bugs correction 26 28 !!---------------------------------------------------------------------- 27 29 #if defined key_zdftke2 || defined key_esopa … … 29 31 !! 'key_zdftke2' TKE vertical physics 30 32 !!---------------------------------------------------------------------- 31 !!----------------------------------------------------------------------32 33 !! zdf_tke2 : update momentum and tracer Kz from a tke scheme 33 !! zdf_tke2_init : initialization, namelist read, and parameters control 34 !! tke_tke : tke time stepping: update tke at now time step (en) 35 !! tke_avn : compute mixing length scale and deduce avm and avt 36 !! tke_init : initialization, namelist read, and parameters control 34 37 !! tke2_rst : read/write tke restart in ocean restart file 35 38 !!---------------------------------------------------------------------- 36 USE oce ! ocean dynamics and active tracers 37 USE dom_oce ! ocean space and time domain 38 USE zdf_oce ! ocean vertical physics 39 USE sbc_oce ! surface boundary condition: ocean 40 USE phycst ! physical constants 41 USE zdfmxl ! mixed layer 42 USE restart ! only for lrst_oce 43 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 44 USE prtctl ! Print control 45 USE in_out_manager ! I/O manager 46 USE iom ! I/O manager library 39 USE oce ! ocean dynamics and active tracers 40 USE dom_oce ! ocean space and time domain 41 USE domvvl ! ocean space and time domain : variable volume layer 42 USE zdf_oce ! ocean vertical physics 43 USE sbc_oce ! surface boundary condition: ocean 44 USE phycst ! physical constants 45 USE zdfmxl ! mixed layer 46 USE restart ! only for lrst_oce 47 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 48 USE prtctl ! Print control 49 USE in_out_manager ! I/O manager 50 USE iom ! I/O manager library 47 51 48 52 IMPLICIT NONE 49 53 PRIVATE 50 54 51 PUBLIC zdf_tke2 52 PUBLIC tke2_rst 55 PUBLIC zdf_tke2 ! routine called in step module 56 PUBLIC tke2_rst ! routine called in step module 53 57 54 58 LOGICAL , PUBLIC, PARAMETER :: lk_zdftke2 = .TRUE. !: TKE vertical mixing flag 55 REAL(wp), PUBLIC :: eboost !: multiplicative coeff of the shear product.56 REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) :: en !: now turbulent kinetic energy57 REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) :: avm !: now mixing coefficient of diffusion for TKE58 REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) :: dissl !: now mixing lenght of dissipation59 59 60 60 #if defined key_c1d 61 ! !! ! 1D cfg only61 ! !!* 1D cfg only * 62 62 REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) :: e_dis, e_mix !: dissipation and mixing turbulent lengh scales 63 63 REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) :: e_pdl, e_ric !: prandl and local Richardson numbers 64 64 #endif 65 65 66 !!gm - variables to be suppressed from the namelist: 67 LOGICAL :: ln_rstke = .FALSE. ! =T restart with tke from a run without tke 68 INTEGER :: nn_itke = 50 ! number of restart iterative loops 69 INTEGER :: nn_ave = 1 ! horizontal average or not on avt, avmu, avmv (=0/1) 70 REAL(wp) :: rn_cri = 2._wp / 9._wp ! critic Richardson number 71 REAL(wp) :: rn_efave = 1._wp ! coefficient for ave : ave=rn_efave*avm 72 !!gm end 73 74 !!gm - variables to be added in the namelist (with a larger default value) 75 REAL(wp) :: rn_bshear = 1.e-20 ! backgrounf shear (>0) 76 !!gm 77 66 78 ! !!! ** Namelist namtke ** 67 LOGICAL :: ln_rstke = .FALSE. ! =T restart with tke from a run without tke68 79 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 not70 INTEGER :: nn_itke = 50 ! number of restart iterative loops71 80 INTEGER :: nn_mxl = 2 ! type of mixing length (=0/1/2/3) 81 !!gm to be change: the ref value of lmin0 and lmin 82 REAL(wp) :: rn_lmin0 = 0.4_wp ! surface min value of mixing length [m] 83 REAL(wp) :: rn_lmin = 0.1_wp ! interior min value of mixing length [m] 72 84 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 85 INTEGER :: nn_avb = 0 ! constant or profile background on avt (=0/1) 75 86 REAL(wp) :: rn_ediff = 0.1_wp ! coefficient for avt: avt=rn_ediff*mxl*sqrt(e) 76 87 REAL(wp) :: rn_ediss = 0.7_wp ! coefficient of the Kolmogoroff dissipation 77 88 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 89 REAL(wp) :: rn_emin = 0.7071e-6_wp ! minimum value of tke [m2/s2] 90 REAL(wp) :: rn_emin0 = 1.e-4_wp ! surface minimum value of tke [m2/s2] 82 91 INTEGER :: nn_havtb = 1 ! horizontal shape or not for avtb (=0/1) 83 92 INTEGER :: nn_etau = 0 ! type of depth penetration of surface tke (=0/1/2) 84 93 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 length86 REAL(wp) :: rn_lmin = 0.1_wp ! interior min value of mixing length87 94 REAL(wp) :: rn_efr = 1.0_wp ! fraction of TKE surface value which penetrates in the ocean 95 LOGICAL :: ln_lc = .FALSE. ! Langmuir cells (LC) as a source term of TKE or not 88 96 REAL(wp) :: rn_lc = 0.15_wp ! coef to compute vertical velocity of Langmuir cells 89 97 90 REAL(wp), DIMENSION (jpi,jpj) :: avtb_2d ! set in tke2_init, for other modif than ice 98 REAL(wp) :: ri_cri ! critic Richardson number (deduced from rn_ediff and rn_ediss values) 99 100 REAL(wp), DIMENSION(jpi,jpj) :: htau ! depth of tke penetration (nn_htau) 101 REAL(wp), DIMENSION(jpi,jpj) :: avtb_2d ! set in tke2_init, for other modif than ice 102 REAL(wp), DIMENSION(jpi,jpj,jpk) :: en ! now turbulent kinetic energy [m2/s2] 103 REAL(wp), DIMENSION(jpi,jpj,jpk) :: dissl ! now mixing lenght of dissipation 91 104 92 105 !! * Substitutions … … 94 107 # include "vectopt_loop_substitute.h90" 95 108 !!---------------------------------------------------------------------- 96 !! NEMO/OPA 3. 0 , LOCEAN-IPSL (2008)109 !! NEMO/OPA 3.2 , LOCEAN-IPSL (2009) 97 110 !! $Id: zdftke2.F90 1201 2008-09-24 13:24:21Z rblod $ 98 111 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) … … 106 119 !! 107 120 !! ** Purpose : Compute the vertical eddy viscosity and diffusivity 108 !! coefficients using a 1.5 turbulent closure scheme.109 !! 110 !! ** Method : The time evolution of the turbulent kinetic energy 111 !! (tke)is computed from a prognostic equation :112 !! d(en)/dt = eboost eav (d(u)/dz)**2! shear production113 !! + d( rn_efave eav d(en)/dz )/dz! diffusion of tke114 !! + grav/rau0 pdl eav d(rau)/dz! stratif. destruc.115 !! - rn_ediss / emxl en**(2/3) !dissipation121 !! coefficients using a turbulent closure scheme (TKE). 122 !! 123 !! ** Method : The time evolution of the turbulent kinetic energy (tke) 124 !! is computed from a prognostic equation : 125 !! d(en)/dt = avm (d(u)/dz)**2 ! shear production 126 !! + d( avm d(en)/dz )/dz ! diffusion of tke 127 !! + avt N^2 ! stratif. destruc. 128 !! - rn_ediss / emxl en**(2/3) ! Kolmogoroff dissipation 116 129 !! with the boundary conditions: 117 130 !! surface: en = max( rn_emin0, rn_ebb sqrt(utau^2 + vtau^2) ) 118 131 !! bottom : en = rn_emin 119 !! -1- The dissipation and mixing turbulent lengh scales are computed 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 !! Four cases : 124 !! nn_mxl=0 : mxl bounded by the distance to surface and bottom. 125 !! zmxld = zmxlm = mxl 126 !! nn_mxl=1 : mxl bounded by the vertical scale factor. 127 !! zmxld = zmxlm = mxl 128 !! nn_mxl=2 : mxl bounded such that the vertical derivative of mxl 129 !! is less than 1 (|d/dz(xml)|<1). 130 !! zmxld = zmxlm = mxl 131 !! nn_mxl=3 : lup = mxl bounded using |d/dz(xml)|<1 from the surface 132 !! to the bottom 133 !! ldown = mxl bounded using |d/dz(xml)|<1 from the bottom 134 !! to the surface 135 !! zmxld = sqrt (lup*ldown) ; zmxlm = min(lup,ldown) 136 !! -2- Compute the now Turbulent kinetic energy. The time differencing 137 !! is implicit for vertical diffusion term, linearized for kolmo- 138 !! goroff dissipation term, and explicit forward for both buoyancy 139 !! and dynamic production terms. Thus a tridiagonal linear system is 140 !! solved. 141 !! Note that - the shear production is multiplied by eboost in order 142 !! to set the critic richardson number to rn_cri (namelist parameter) 143 !! - the destruction by stratification term is multiplied 144 !! by the Prandtl number (defined by an empirical funtion of the local 145 !! Richardson number) if nn_pdl=1 (namelist parameter) 146 !! coefficient (zesh2): 147 !! -3- Compute the now vertical eddy vicosity and diffusivity 148 !! coefficients from en (before the time stepping) and zmxlm: 149 !! avm = max( avtb, rn_ediff*zmxlm*en^1/2 ) 150 !! avt = max( avmb, pdl*avm ) (pdl=1 if nn_pdl=0) 132 !! The associated critical Richardson number is: ri_cri = 2/(2+rn_ediss/rn_ediff) 133 !! 134 !! The now Turbulent kinetic energy is computed using the following 135 !! time stepping: implicit for vertical diffusion term, linearized semi 136 !! implicit for kolmogoroff dissipation term, and explicit forward for 137 !! both buoyancy and shear production terms. Therefore a tridiagonal 138 !! linear system is solved. Note that buoyancy and shear terms are 139 !! discretized in a energy conserving form (Bruchard 2002). 140 !! 141 !! The dissipative and mixing length scale are computed from en and 142 !! the stratification (see tke_avn) 143 !! 144 !! The now vertical eddy vicosity and diffusivity coefficients are 145 !! given by: 146 !! avm = max( avtb, rn_ediff * zmxlm * en^1/2 ) 147 !! avt = max( avmb, pdl * avm ) 151 148 !! eav = max( avmb, avm ) 152 !! avt and avm are horizontally averaged to avoid numerical insta- 153 !! bilities. 154 !! N.B. The computation is done from jk=2 to jpkm1 except for 155 !! en. Surface value of avt avmu avmv are set once a time to 156 !! their background value in routine zdf_tke2_init. 149 !! where pdl, the inverse of the Prandtl number is 1 if nn_pdl=0 and 150 !! given by an empirical funtion of the localRichardson number if nn_pdl=1 157 151 !! 158 152 !! ** Action : compute en (now turbulent kinetic energy) … … 163 157 !! Mellor and Blumberg, JPO 2004 164 158 !! Axell, JGR, 2002 159 !! Bruchard OM 2002 165 160 !!---------------------------------------------------------------------- 166 167 INTEGER, INTENT(in) :: kt ! ocean time step 168 169 IF( kt == nit000 ) CALL zdf_tke2_init ! Initialization (first time-step only) 170 ! 171 CALL tke_tke ! now tke (en) 172 CALL tke_mlmc ! now avmu, avmv, avt 173 161 INTEGER, INTENT(in) :: kt ! ocean time step 162 !!---------------------------------------------------------------------- 163 ! 164 IF( kt == nit000 ) CALL tke_init ! initialisation 165 ! 166 CALL tke_tke ! now tke (en) 167 ! 168 CALL tke_avn ! now avt, avm, avmu, avmv 169 ! 174 170 END SUBROUTINE zdf_tke2 171 175 172 176 173 SUBROUTINE tke_tke 177 174 !!---------------------------------------------------------------------- 178 !! Now Turbulent kinetic energy (output in en) 175 !! *** ROUTINE tke_tke *** 176 !! 177 !! ** Purpose : Compute the now Turbulente Kinetic Energy (TKE) 178 !! 179 !! ** Method : - TKE surface boundary condition 180 !! - source term due to Langmuir cells (ln_lc=T) 181 !! - source term due to shear (saved in avmu, avmv arrays) 182 !! - Now TKE : resolution of the TKE equation by inverting 183 !! a tridiagonal linear system by a "methode de chasse" 184 !! - increase TKE due to surface and internal wave breaking 185 !! 186 !! ** Action : - en : now turbulent kinetic energy) 187 !! - avmu, avmv : production of TKE by shear at u and v-points 188 !! (= Kz dz[Ub] * dz[Un] ) 179 189 !! --------------------------------------------------------------------- 180 !! Surface boundary condition for TKE 181 !! Resolution of a tridiagonal linear system by a "methode de chasse" 182 !! computation from level 2 to jpkm1 (e(1) computed and e(jpk)=0 ). 183 !! avm represents the turbulent coefficient of the TKE 184 !!---------------------------------------------------------------------- 185 USE oce, zwd => ua ! use ua as workspace 186 USE oce, zdiag1 => va ! use va as workspace 187 USE oce, zdiag2 => ta ! use ta as workspace 188 USE oce, ztkelc => sa ! use sa as workspace 189 ! 190 INTEGER :: ji, jj, jk ! dummy loop arguments 191 REAL(wp) :: zbbrau, zesurf, zesh2 ! temporary scalars 192 REAL(wp) :: zfact1, ztx2, zdku ! - - 193 REAL(wp) :: zfact2, zty2, zdkv ! - - 194 REAL(wp) :: zfact3, zcoef, zcof ! - - 195 REAL(wp) :: zus, zwlc, zind ! - - 190 USE oce, zdiag => ua ! use ua as workspace 191 USE oce, zd_up => va ! use va as workspace 192 USE oce, zd_lw => ta ! use ta as workspace 193 !! 194 INTEGER :: ji, jj, jk ! dummy loop arguments 195 REAL(wp) :: zbbrau, zesurf, zesh2 ! temporary scalars 196 REAL(wp) :: zfact1, zfact2, zfact3 ! - - 197 REAL(wp) :: ztx2 , zty2 , zcof ! - - 198 REAL(wp) :: zus , zwlc , zind ! - - 199 REAL(wp) :: zzd_up, zzd_lw ! - - 196 200 INTEGER , DIMENSION(jpi,jpj) :: imlc ! 2D workspace 197 REAL(wp), DIMENSION(jpi,jpj) :: zhtau ! - -198 201 REAL(wp), DIMENSION(jpi,jpj) :: zhlc ! - - 199 202 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zpelc ! 3D workspace 200 203 !!-------------------------------------------------------------------- 201 ! ! Local constant initialization202 zbbrau = .5 * rn_ebb / rau0 203 zfact1 = -.5 * rdt * rn_efave204 ! 205 zbbrau = .5 * rn_ebb / rau0 ! Local constant initialisation 206 zfact1 = -.5 * rdt 204 207 zfact2 = 1.5 * rdt * rn_ediss 205 208 zfact3 = 0.5 * rn_ediss 206 207 ! Surface boundary condition on tke208 ! -------------------------------------------------209 ! en(1) = rn_ebb sqrt(utau^2+vtau^2) / rau0 (min value rn_emin0)210 !CDIR NOVERRCHK 211 DO jj = 2, jpjm1 209 ! 210 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 211 ! ! Surface boundary condition on tke 212 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 213 !CDIR NOVERRCHK 214 DO jj = 2, jpjm1 ! en(1) = rn_ebb sqrt(utau^2+vtau^2) / rau0 (min value rn_emin0) 212 215 !CDIR NOVERRCHK 213 216 DO ji = fs_2, fs_jpim1 ! vector opt. 214 ztx2 = utau(ji-1,jj ) + utau(ji,jj)215 zty2 = vtau(ji ,jj-1) + vtau(ji,jj)217 ztx2 = utau(ji-1,jj ) + utau(ji,jj) 218 zty2 = vtau(ji ,jj-1) + vtau(ji,jj) 216 219 zesurf = zbbrau * SQRT( ztx2 * ztx2 + zty2 * zty2 ) 217 220 en(ji,jj,1) = MAX( zesurf, rn_emin0 ) * tmask(ji,jj,1) 218 221 END DO 219 222 END DO 220 221 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 222 ! Langmuir circulation source term 223 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 224 IF( ln_lc ) THEN 223 ! 224 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 225 ! ! Bottom boundary condition on tke 226 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 227 !!gm to be added soon 228 ! 229 ! 230 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 231 IF( ln_lc ) THEN ! Langmuir circulation source term added to tke 232 ! !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 225 233 ! 226 ! Computation oftotal energy produce by LC : cumulative sum over jk234 ! !* total energy produce by LC : cumulative sum over jk 227 235 zpelc(:,:,1) = MAX( rn2b(:,:,1), 0. ) * fsdepw(:,:,1) * fse3w(:,:,1) 228 236 DO jk = 2, jpk 229 237 zpelc(:,:,jk) = zpelc(:,:,jk-1) + MAX( rn2b(:,:,jk), 0. ) * fsdepw(:,:,jk) * fse3w(:,:,jk) 230 238 END DO 231 ! 232 ! Computation of finite Langmuir Circulation depth 233 ! Initialization to the number of w ocean point mbathy 234 imlc(:,:) = mbathy(:,:) 239 ! !* finite Langmuir Circulation depth 240 imlc(:,:) = mbathy(:,:) ! Initialization to the number of w ocean point mbathy 235 241 DO jk = jpkm1, 2, -1 236 DO jj = 1, jpj 237 DO ji = 1, jpi 238 ! Last w-level at which zpelc>=0.5*us*us with us=0.016*wind(starting from jpk-1) 242 DO jj = 1, jpj ! Last w-level at which zpelc>=0.5*us*us 243 DO ji = 1, jpi ! with us=0.016*wind(starting from jpk-1) 239 244 zus = 0.000128 * wndm(ji,jj) * wndm(ji,jj) 240 245 IF( zpelc(ji,jj,jk) > zus ) imlc(ji,jj) = jk … … 242 247 END DO 243 248 END DO 244 ! 245 ! finite LC depth 246 DO jj = 1, jpj 249 ! ! finite LC depth 250 # if defined key_vectopt_loop 251 DO jj = 1, 1 252 DO ji = 1, jpij ! vector opt. (forced unrolling) 253 # else 254 DO jj = 1, jpj 247 255 DO ji = 1, jpi 256 # endif 248 257 zhlc(ji,jj) = fsdepw(ji,jj,imlc(ji,jj)) 249 258 END DO 250 259 END DO 251 !252 260 # if defined key_c1d 253 261 hlc(:,:) = zhlc(:,:) * tmask(:,:,1) ! c1d configuration: save finite Langmuir Circulation depth 254 262 # endif 255 ! 256 ! TKE Langmuir circulation source term added to en 257 !CDIR NOVERRCHK 258 DO jk = 2, jpkm1 259 !CDIR NOVERRCHK 260 DO jj = 2, jpjm1 261 !CDIR NOVERRCHK 262 DO ji = fs_2, fs_jpim1 ! vector opt. 263 ! Stokes drift 264 zus = 0.016 * wndm(ji,jj) 265 ! computation of vertical velocity due to LC 263 !CDIR NOVERRCHK 264 DO jk = 2, jpkm1 !* TKE Langmuir circulation source term added to en 265 !CDIR NOVERRCHK 266 DO jj = 2, jpjm1 267 !CDIR NOVERRCHK 268 DO ji = fs_2, fs_jpim1 ! vector opt. 269 !!gm replace here wndn by a formulation with the stress module 270 zus = 0.016 * wndm(ji,jj) ! Stokes drift 271 ! ! vertical velocity due to LC 266 272 zind = 0.5 - SIGN( 0.5, fsdepw(ji,jj,jk) - zhlc(ji,jj) ) 267 273 zwlc = zind * rn_lc * zus * SIN( rpi * fsdepw(ji,jj,jk) / zhlc(ji,jj) ) 268 ! TKE Langmuir circulation source term269 en(ji,jj,jk) = en(ji,jj,jk) + rdt * ( ( zwlc * zwlc * zwlc ) / zhlc(ji,jj)) * tmask(ji,jj,jk)274 ! ! TKE Langmuir circulation source term 275 en(ji,jj,jk) = en(ji,jj,jk) + rdt * ( zwlc * zwlc * zwlc ) / zhlc(ji,jj) * tmask(ji,jj,jk) 270 276 END DO 271 277 END DO … … 273 279 ! 274 280 ENDIF 275 276 277 ! Shear production at uw- and vw-points (energy conserving form) 278 DO jk = 2, jpkm1 281 ! 282 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 283 ! ! Now Turbulent kinetic energy (output in en) 284 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 285 ! ! Resolution of a tridiagonal linear system by a "methode de chasse" 286 ! ! computation from level 2 to jpkm1 (e(1) already computed and e(jpk)=0 ). 287 ! ! zdiag : diagonal zd_up : upper diagonal zd_lw : lower diagonal 288 ! 289 DO jk = 2, jpkm1 !* Shear production at uw- and vw-points (energy conserving form) 290 DO jj = 1, jpj ! here avmu, avmv used as workspace 291 DO ji = 1, jpi 292 avmu(ji,jj,jk) = avmu(ji,jj,jk) * ( un(ji,jj,jk-1) - un(ji,jj,jk) ) & 293 & * ( ub(ji,jj,jk-1) - ub(ji,jj,jk) ) & 294 & / ( fse3uw_n(ji,jj,jk) & 295 & * fse3uw_b(ji,jj,jk) ) 296 avmv(ji,jj,jk) = avmv(ji,jj,jk) * ( vn(ji,jj,jk-1) - vn(ji,jj,jk) ) & 297 & * ( vb(ji,jj,jk-1) - vb(ji,jj,jk) ) & 298 & / ( fse3vw_n(ji,jj,jk) & 299 & * fse3vw_b(ji,jj,jk) ) 300 END DO 301 END DO 302 END DO 303 ! 304 DO jk = 2, jpkm1 !* Matrix and right hand side in en 279 305 DO jj = 2, jpjm1 280 306 DO ji = fs_2, fs_jpim1 ! vector opt. 281 avmu(ji,jj,jk) = avmu(ji,jj,jk) * ( un(ji,jj,jk-1) - un(ji,jj,jk) ) * & 282 & ( ub(ji,jj,jk-1) - ub(ji,jj,jk) ) / & 283 & ( fse3uw_n(ji,jj,jk) * fse3uw_b(ji,jj,jk) ) 284 avmv(ji,jj,jk) = avmv(ji,jj,jk) * ( vn(ji,jj,jk-1) - vn(ji,jj,jk) ) * & 285 & ( vb(ji,jj,jk-1) - vb(ji,jj,jk) ) / & 286 & ( fse3vw_n(ji,jj,jk) * fse3vw_b(ji,jj,jk) ) 287 ENDDO 288 ENDDO 289 ENDDO 290 291 ! Lateral boundary conditions (avmu,avmv) (sign unchanged) 292 CALL lbc_lnk( avmu, 'U', 1. ) ; CALL lbc_lnk( avmv, 'V', 1. ) 293 294 ! Now Turbulent kinetic energy (output in en) 295 ! ------------------------------- 296 ! Resolution of a tridiagonal linear system by a "methode de chasse" 297 ! computation from level 2 to jpkm1 (e(1) already computed and e(jpk)=0 ). 298 ! zwd : diagonal zdiag1 : upper diagonal zdiag2 : lower diagonal 299 ! Warning : after this step, en : right hand side of the matrix 300 DO jk = 2, jpkm1 301 DO jj = 2, jpjm1 302 DO ji = fs_2, fs_jpim1 ! vector opt. 303 ! ! shear prod. at w-point weightened by mask 307 zcof = zfact1 * tmask(ji,jj,jk) 308 zzd_up = zcof * ( avm (ji,jj,jk+1) + avm (ji,jj,jk ) ) & ! upper diagonal 309 & / ( fse3t(ji,jj,jk ) * fse3w(ji,jj,jk ) ) 310 zzd_lw = zcof * ( avm (ji,jj,jk ) + avm (ji,jj,jk-1) ) & ! lower diagonal 311 & / ( fse3t(ji,jj,jk-1) * fse3w(ji,jj,jk ) ) 312 ! ! shear prod. at w-point weightened by mask 304 313 zesh2 = ( avmu(ji-1,jj,jk) + avmu(ji,jj,jk) ) / MAX( 1.e0 , umask(ji-1,jj,jk) + umask(ji,jj,jk) ) & 305 314 & + ( avmv(ji,jj-1,jk) + avmv(ji,jj,jk) ) / MAX( 1.e0 , vmask(ji,jj-1,jk) + vmask(ji,jj,jk) ) 306 ! ! Matrix 307 zcof = zfact1 * tmask(ji,jj,jk) 308 ! ! lower diagonal 309 zdiag2(ji,jj,jk) = zcof * ( avm (ji,jj,jk ) + avm (ji,jj,jk-1) ) & 310 & / ( fse3t(ji,jj,jk-1) * fse3w(ji,jj,jk ) ) 311 ! ! upper diagonal 312 zdiag1(ji,jj,jk) = zcof * ( avm (ji,jj,jk+1) + avm (ji,jj,jk ) ) & 313 & / ( fse3t(ji,jj,jk ) * fse3w(ji,jj,jk) ) 314 ! ! diagonal 315 zwd(ji,jj,jk) = 1. - zdiag2(ji,jj,jk) - zdiag1(ji,jj,jk) + zfact2 * dissl(ji,jj,jk) * tmask(ji,jj,jk) 315 ! 316 zd_up(ji,jj,jk) = zzd_up ! Matrix (zdiag, zd_up, zd_lw) 317 zd_lw(ji,jj,jk) = zzd_lw 318 zdiag(ji,jj,jk) = 1.e0 - zzd_lw - zzd_up + zfact2 * dissl(ji,jj,jk) * tmask(ji,jj,jk) 316 319 ! 317 ! 320 ! ! right hand side in en 318 321 en(ji,jj,jk) = en(ji,jj,jk) + rdt * ( zesh2 - avt(ji,jj,jk) * rn2(ji,jj,jk) & 319 322 & + zfact3 * dissl(ji,jj,jk) * en (ji,jj,jk) ) * tmask(ji,jj,jk) … … 321 324 END DO 322 325 END DO 323 ! Matrix inversion from level 2 (tke prescribed at level 1) 324 !!-------------------------------- 326 ! !* Matrix inversion from level 2 (tke prescribed at level 1) 325 327 DO jk = 3, jpkm1 ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 326 328 DO jj = 2, jpjm1 327 329 DO ji = fs_2, fs_jpim1 ! vector opt. 328 z wd(ji,jj,jk) = zwd(ji,jj,jk) - zdiag2(ji,jj,jk) * zdiag1(ji,jj,jk-1) / zwd(ji,jj,jk-1)330 zdiag(ji,jj,jk) = zdiag(ji,jj,jk) - zd_lw(ji,jj,jk) * zd_up(ji,jj,jk-1) / zdiag(ji,jj,jk-1) 329 331 END DO 330 332 END DO … … 332 334 DO jj = 2, jpjm1 ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1 333 335 DO ji = fs_2, fs_jpim1 ! vector opt. 334 zd iag2(ji,jj,2) = en(ji,jj,2) - zdiag2(ji,jj,2) * en(ji,jj,1) ! Surface boudary conditions on tke336 zd_lw(ji,jj,2) = en(ji,jj,2) - zd_lw(ji,jj,2) * en(ji,jj,1) ! Surface boudary conditions on tke 335 337 END DO 336 338 END DO … … 338 340 DO jj = 2, jpjm1 339 341 DO ji = fs_2, fs_jpim1 ! vector opt. 340 zd iag2(ji,jj,jk) = en(ji,jj,jk) - zdiag2(ji,jj,jk) / zwd(ji,jj,jk-1) *zdiag2(ji,jj,jk-1)342 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) 341 343 END DO 342 344 END DO … … 344 346 DO jj = 2, jpjm1 ! thrid recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk 345 347 DO ji = fs_2, fs_jpim1 ! vector opt. 346 en(ji,jj,jpkm1) = zd iag2(ji,jj,jpkm1) / zwd(ji,jj,jpkm1)348 en(ji,jj,jpkm1) = zd_lw(ji,jj,jpkm1) / zdiag(ji,jj,jpkm1) 347 349 END DO 348 350 END DO … … 350 352 DO jj = 2, jpjm1 351 353 DO ji = fs_2, fs_jpim1 ! vector opt. 352 en(ji,jj,jk) = ( zd iag2(ji,jj,jk) - zdiag1(ji,jj,jk) * en(ji,jj,jk+1) ) / zwd(ji,jj,jk)354 en(ji,jj,jk) = ( zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) * en(ji,jj,jk+1) ) / zdiag(ji,jj,jk) 353 355 END DO 354 356 END DO … … 362 364 END DO 363 365 364 ! 5. Add extra TKE due to surface and internal wave breaking (nn_etau /= 0) 365 !!---------------------------------------------------------- 366 IF( nn_etau /= 0 ) THEN ! extra tke : en = en + rn_efr * en(1) * exp( -z/zhtau ) 367 ! 368 SELECT CASE( nn_htau ) ! Choice of the depth of penetration 369 CASE( 0 ) ! constant depth penetration (here 10 meters) 370 DO jj = 2, jpjm1 371 DO ji = fs_2, fs_jpim1 ! vector opt. 372 zhtau(ji,jj) = 10. 373 END DO 374 END DO 375 CASE( 1 ) ! meridional profile 1 376 DO jj = 2, jpjm1 ! ( 5m in the tropics to a maximum of 40 m at high lat.) 377 DO ji = fs_2, fs_jpim1 ! vector opt. 378 zhtau(ji,jj) = MAX( 5., MIN( 40., 60.*ABS( SIN( rpi/180. * gphit(ji,jj) ) ) ) ) 379 END DO 380 END DO 381 CASE( 2 ) ! meridional profile 2 382 DO jj = 2, jpjm1 ! ( 5m in the tropics to a maximum of 60 m at high lat.) 383 DO ji = fs_2, fs_jpim1 ! vector opt. 384 zhtau(ji,jj) = MAX( 5.,6./4.* MIN( 40., 60.*ABS( SIN( rpi/180. * gphit(ji,jj) ) ) ) ) 385 END DO 386 END DO 387 END SELECT 388 ! 389 IF( nn_etau == 1 ) THEN ! extra term throughout the water column 390 DO jk = 2, jpkm1 391 DO jj = 2, jpjm1 392 DO ji = fs_2, fs_jpim1 ! vector opt. 393 en(ji,jj,jk) = en(ji,jj,jk) & 394 & + rn_efr * en(ji,jj,1)*EXP( -fsdepw(ji,jj,jk) / zhtau(ji,jj) ) & 395 & * ( 1.e0 - fr_i(ji,jj) ) * tmask(ji,jj,jk) 396 END DO 397 END DO 398 END DO 399 ELSEIF( nn_etau == 2 ) THEN ! extra term only at the base of the mixed layer 400 DO jj = 2, jpjm1 401 DO ji = fs_2, fs_jpim1 ! vector opt. 402 jk = nmln(ji,jj) 403 en(ji,jj,jk) = en(ji,jj,jk) & 404 & + rn_efr * en(ji,jj,1)*EXP( -fsdepw(ji,jj,jk) / zhtau(ji,jj) ) & 405 & * ( 1.e0 - fr_i(ji,jj) ) * tmask(ji,jj,jk) 406 END DO 407 END DO 408 ENDIF 409 ! 366 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 367 ! ! TKE due to surface and internal wave breaking 368 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 369 IF( nn_etau == 1 ) THEN !* penetration throughout the water column 370 DO jk = 2, jpkm1 371 DO jj = 2, jpjm1 372 DO ji = fs_2, fs_jpim1 ! vector opt. 373 en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -fsdepw(ji,jj,jk) / htau(ji,jj) ) & 374 & * ( 1.e0 - fr_i(ji,jj) ) * tmask(ji,jj,jk) 375 END DO 376 END DO 377 END DO 378 ELSEIF( nn_etau == 2 ) THEN !* act only at the base of the mixed layer (jk=nmln) 379 DO jj = 2, jpjm1 380 DO ji = fs_2, fs_jpim1 ! vector opt. 381 jk = nmln(ji,jj) 382 en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -fsdepw(ji,jj,jk) / htau(ji,jj) ) & 383 & * ( 1.e0 - fr_i(ji,jj) ) * tmask(ji,jj,jk) 384 END DO 385 END DO 410 386 ENDIF 411 ! Lateral boundary conditions (sign unchanged)412 CALL lbc_lnk( en, 'W', 1. ) 413 387 ! 388 CALL lbc_lnk( en, 'W', 1. ) ! Lateral boundary conditions (sign unchanged) 389 ! 414 390 END SUBROUTINE tke_tke 415 391 416 SUBROUTINE tke_mlmc 392 393 SUBROUTINE tke_avn 417 394 !!---------------------------------------------------------------------- 418 !! 395 !! *** ROUTINE tke_avn *** 396 !! 397 !! ** Purpose : Compute the vertical eddy viscosity and diffusivity 398 !! 399 !! ** Method : At this stage, en, the now TKE, is known (computed in 400 !! the tke_tke routine). First, the now mixing lenth is 401 !! computed from en and the strafification (N^2), then the mixings 402 !! coefficients are computed. 403 !! - Mixing length : a first evaluation of the mixing lengh 404 !! scales is: 405 !! mxl = sqrt(2*en) / N 406 !! where N is the brunt-vaisala frequency, with a minimum value set 407 !! to rn_lmin (rn_lmin0) in the interior (surface) ocean. 408 !! The mixing and dissipative length scale are bound as follow : 409 !! nn_mxl=0 : mxl bounded by the distance to surface and bottom. 410 !! zmxld = zmxlm = mxl 411 !! nn_mxl=1 : mxl bounded by the e3w and zmxld = zmxlm = mxl 412 !! nn_mxl=2 : mxl bounded such that the vertical derivative of mxl is 413 !! less than 1 (|d/dz(mxl)|<1) and zmxld = zmxlm = mxl 414 !! nn_mxl=3 : mxl is bounded from the surface to the bottom usings 415 !! |d/dz(xml)|<1 to obtain lup, and from the bottom to 416 !! the surface to obtain ldown. the resulting length 417 !! scales are: 418 !! zmxld = sqrt( lup * ldown ) 419 !! zmxlm = min ( lup , ldown ) 420 !! - Vertical eddy viscosity and diffusivity: 421 !! avm = max( avtb, rn_ediff * zmxlm * en^1/2 ) 422 !! avt = max( avmb, pdlr * avm ) 423 !! with pdlr=1 if nn_pdl=0, pdlr=1/pdl=F(Ri) otherwise. 424 !! 425 !! ** Action : - avt : now vertical eddy diffusivity (w-point) 426 !! - avmu, avmv : now vertical eddy viscosity at uw- and vw-points 419 427 !!---------------------------------------------------------------------- 420 428 USE oce, zmpdl => ua ! use ua as workspace 421 429 USE oce, zmxlm => va ! use va as workspace 422 430 USE oce, zmxld => ta ! use ta as workspace 423 ! 424 INTEGER :: ji, jj, jk 425 REAL(wp) :: zrn2, zraug 426 REAL(wp) :: ztx2, zdku 427 REAL(wp) :: zty2, zdkv 428 REAL(wp) :: zcoef, zav 429 REAL(wp) :: zpdl , zri, zsqen! - -430 REAL(wp) :: zemxl, zemlm, zemlp 431 !! 432 INTEGER :: ji, jj, jk ! dummy loop arguments 433 REAL(wp) :: zrn2, zraug ! temporary scalars 434 REAL(wp) :: ztx2, zdku ! - - 435 REAL(wp) :: zty2, zdkv ! - - 436 REAL(wp) :: zcoef, zav ! - - 437 REAL(wp) :: zpdlr, zri, zsqen ! - - 438 REAL(wp) :: zemxl, zemlm, zemlp ! - - 431 439 !!-------------------------------------------------------------------- 432 440 433 ! >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>434 ! Mixing length435 ! <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<436 437 ! Buoyancy length scale: l=sqrt(2*e/n**2)438 ! ---------------------441 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 442 ! ! Mixing length 443 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 444 ! 445 ! !* Buoyancy length scale: l=sqrt(2*e/n**2) 446 ! 439 447 IF( ln_mxl0 ) THEN ! surface mixing length = F(stress) : l=vkarmn*2.e5*sqrt(utau^2 + vtau^2)/(rau0*g) 440 448 !!gm this should be useless … … 466 474 END DO 467 475 END DO 468 469 ! Physical limits for the mixing length 470 ! ------------------------------------- 476 ! 477 !!gm CAUTION: to be added here the bottom boundary condition on zmxlm 478 ! 479 ! !* Physical limits for the mixing length 480 ! 471 481 zmxld(:,:, 1 ) = zmxlm(:,:,1) ! surface set to the minimum value 472 482 zmxld(:,:,jpk) = rn_lmin ! bottom set to the minimum value 473 483 ! 474 484 SELECT CASE ( nn_mxl ) 475 485 ! … … 545 555 ! 546 556 END SELECT 547 557 ! 548 558 # if defined key_c1d 549 ! c1d configuration : save mixing and dissipation turbulent length scales 550 e_dis(:,:,:) = zmxld(:,:,:) 559 e_dis(:,:,:) = zmxld(:,:,:) ! c1d configuration : save mixing and dissipation turbulent length scales 551 560 e_mix(:,:,:) = zmxlm(:,:,:) 552 561 # endif 553 562 554 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>! 555 ! II Vertical eddy viscosity on tke (put in zmxlm) and first estimate of avt ! 556 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<! 557 ! Surface boundary condition on avt and avm : jk = 1 558 ! ------------------------------------------------- 559 ! avt(1) = avmb(1) and avt(jpk) = 0. 560 ! avm(1) = avmb(1) and avm(jpk) = 0. 561 ! 562 !CDIR NOVERRCHK 563 DO jk = 1, jpkm1 563 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 564 ! ! Vertical eddy viscosity and diffusivity (avmu, avmv, avt) 565 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 566 !CDIR NOVERRCHK 567 DO jk = 1, jpkm1 !* vertical eddy viscosity & diffivity at w-points 564 568 !CDIR NOVERRCHK 565 569 DO jj = 2, jpjm1 … … 568 572 zsqen = SQRT( en(ji,jj,jk) ) 569 573 zav = rn_ediff * zmxlm(ji,jj,jk) * zsqen 570 avm (ji,jj,jk) = MAX( zav, avmb(jk) ) * tmask(ji,jj,jk)571 avt (ji,jj,jk) = MAX( zav, avtb_2d(ji,jj) * avtb(jk) ) * tmask(ji,jj,jk)574 avm (ji,jj,jk) = MAX( zav, avmb(jk) ) * tmask(ji,jj,jk) 575 avt (ji,jj,jk) = MAX( zav, avtb_2d(ji,jj) * avtb(jk) ) * tmask(ji,jj,jk) 572 576 dissl(ji,jj,jk) = zsqen / zmxld(ji,jj,jk) 573 577 END DO 574 578 END DO 575 579 END DO 576 577 ! Lateral boundary conditions (sign unchanged) 578 CALL lbc_lnk( avm, 'W', 1. ) 579 580 DO jk = 2, jpkm1 ! only vertical eddy viscosity 580 CALL lbc_lnk( avm, 'W', 1. ) ! Lateral boundary conditions (sign unchanged) 581 ! 582 DO jk = 2, jpkm1 !* vertical eddy viscosity at u- and v-points 581 583 DO jj = 2, jpjm1 582 584 DO ji = fs_2, fs_jpim1 ! vector opt. … … 587 589 END DO 588 590 CALL lbc_lnk( avmu, 'U', 1. ) ; CALL lbc_lnk( avmv, 'V', 1. ) ! Lateral boundary conditions 589 590 591 ! Prandtl number 592 ! ---------------------------------------- 593 IF( nn_pdl == 1 ) THEN 591 ! 592 IF( nn_pdl == 1 ) THEN !* Prandtl number case: update avt 594 593 DO jk = 2, jpkm1 595 594 DO jj = 2, jpjm1 596 595 DO ji = fs_2, fs_jpim1 ! vector opt. 597 596 zcoef = 0.5 / ( fse3w(ji,jj,jk) * fse3w(ji,jj,jk) ) 598 ! shear 599 zdku = avmu(ji-1,jj,jk) * (un(ji-1,jj,jk-1)-un(ji-1,jj,jk)) * (ub(ji-1,jj,jk-1)-ub(ji-1,jj,jk)) + & 600 & avmu(ji ,jj,jk) * (un(ji ,jj,jk-1)-un(ji ,jj,jk)) * (ub(ji ,jj,jk-1)-ub(ji ,jj,jk)) 601 zdkv = avmv(ji,jj-1,jk) * (vn(ji,jj-1,jk-1)-vn(ji,jj-1,jk)) * (vb(ji,jj-1,jk-1)-vb(ji,jj-1,jk)) + & 602 & avmv(ji,jj ,jk) * (vn(ji,jj ,jk-1)-vn(ji,jj ,jk)) * (vb(ji,jj ,jk-1)-vb(ji,jj ,jk)) 603 zri = MAX( rn2b(ji,jj,jk), 0. ) * avm(ji,jj,jk)/ (zcoef * (zdku + zdkv + 1.e-20) ) ! local Richardson number 604 # if defined key_cfg_1d 597 ! ! shear 598 zdku = avmu(ji-1,jj,jk) * ( un(ji-1,jj,jk-1) - un(ji-1,jj,jk) ) * ( ub(ji-1,jj,jk-1) - ub(ji-1,jj,jk) ) & 599 & + avmu(ji ,jj,jk) * ( un(ji ,jj,jk-1) - un(ji ,jj,jk) ) * ( ub(ji ,jj,jk-1) - ub(ji ,jj,jk) ) 600 zdkv = avmv(ji,jj-1,jk) * ( vn(ji,jj-1,jk-1) - vn(ji,jj-1,jk) ) * ( vb(ji,jj-1,jk-1) - vb(ji,jj-1,jk) ) & 601 & + avmv(ji,jj ,jk) * ( vn(ji,jj ,jk-1) - vn(ji,jj ,jk) ) * ( vb(ji,jj ,jk-1) - vb(ji,jj ,jk) ) 602 ! ! local Richardson number 603 zri = MAX( rn2b(ji,jj,jk), 0. ) * avm(ji,jj,jk) / (zcoef * (zdku + zdkv + rn_bshear ) ) 604 zpdlr = MAX( 0.1, 0.2 / MAX( 0.2 , zri ) ) 605 !!gm and even better with the use of the "true" ri_crit=0.22222... (this change the results!) 606 !!gm zpdlr = MAX( 0.1, ri_crit / MAX( ri_crit , zri ) ) 607 avt(ji,jj,jk) = MAX( zpdlr * avt(ji,jj,jk), avtb_2d(ji,jj) * avtb(jk) ) * tmask(ji,jj,jk) 608 # if defined key_c1d 609 e_pdl(ji,jj,jk) = zpdlr * tmask(ji,jj,jk) ! c1d configuration : save masked Prandlt number 605 610 e_ric(ji,jj,jk) = zri * tmask(ji,jj,jk) ! c1d config. : save Ri 606 611 # endif 607 zpdl = 1.0 ! Prandtl number608 IF( zri >= 0.2 ) zpdl = 0.2 / zri609 zpdl = MAX( 0.1, zpdl )610 avt(ji,jj,jk) = MAX( zpdl * avt(ji,jj,jk), avtb_2d(ji,jj) * avtb(jk) ) * tmask(ji,jj,jk)611 zmpdl(ji,jj,jk) = zpdl * tmask(ji,jj,jk)612 612 END DO 613 613 END DO 614 614 END DO 615 615 ENDIF 616 617 # if defined key_c1d618 e_pdl(:,:,2:jpkm1) = zmxld(:,:,2:jpkm1) ! c1d configuration : save masked Prandlt number619 e_pdl(:,:, 1) = e_pdl(:,:, 2)620 e_pdl(:,:, jpk) = e_pdl(:,:, jpkm1)621 # endif622 623 616 CALL lbc_lnk( avt, 'W', 1. ) ! Lateral boundary conditions on avt (sign unchanged) 624 617 … … 629 622 ENDIF 630 623 ! 631 END SUBROUTINE tke_mlmc 632 633 SUBROUTINE zdf_tke2_init 624 END SUBROUTINE tke_avn 625 626 627 SUBROUTINE tke_init 634 628 !!---------------------------------------------------------------------- 635 629 !! *** ROUTINE zdf_tke2_init *** 636 630 !! 637 631 !! ** Purpose : Initialization of the vertical eddy diffivity and 638 !! viscosity when using a tke turbulent closure scheme632 !! viscosity when using a tke turbulent closure scheme 639 633 !! 640 634 !! ** Method : Read the namtke namelist and check the parameters 641 !! called at the first timestep (nit000)635 !! called at the first timestep (nit000) 642 636 !! 643 637 !! ** input : Namlist namtke 644 638 !! 645 639 !! ** Action : Increase by 1 the nstop flag is setting problem encounter 646 !!647 640 !!---------------------------------------------------------------------- 648 USE dynzdf_exp649 USE trazdf_exp650 !651 # if defined key_vectopt_memory652 641 INTEGER :: ji, jj, jk ! dummy loop indices 653 # else654 INTEGER :: jk ! dummy loop indices655 # endif656 642 !! 657 643 NAMELIST/namtke/ ln_rstke, rn_ediff, rn_ediss, rn_ebb , rn_efave, rn_emin, & … … 661 647 !!---------------------------------------------------------------------- 662 648 663 ! Read Namelist namtke : Turbulente Kinetic Energy 664 ! -------------------- 665 REWIND ( numnam ) 649 REWIND ( numnam ) !* Read Namelist namtke : Turbulente Kinetic Energy 666 650 READ ( numnam, namtke ) 667 668 ! Compute boost associated with the Richardson critic 669 ! (control values: rn_cri = 0.3 ==> eboost=1.25 for nn_pdl=1) 670 ! ( rn_cri = 0.222 ==> eboost=1. ) 671 eboost = rn_cri * ( 2. + rn_ediss / rn_ediff ) / 2. 672 673 674 675 ! Parameter control and print 676 ! --------------------------- 677 ! Control print 678 IF(lwp) THEN 651 652 ri_cri = 2. / ( 2. + rn_ediss / rn_ediff ) ! resulting critical Richardson number 653 654 IF(lwp) THEN !* Control print 679 655 WRITE(numout,*) 680 WRITE(numout,*) 'zdf_tke2 _init : tke turbulent closure scheme'681 WRITE(numout,*) '~~~~~~~~ ~~~~'656 WRITE(numout,*) 'zdf_tke2 : tke turbulent closure scheme - initialisation' 657 WRITE(numout,*) '~~~~~~~~' 682 658 WRITE(numout,*) ' Namelist namtke : set tke mixing parameters' 683 659 WRITE(numout,*) ' restart with tke from no tke ln_rstke = ', ln_rstke … … 685 661 WRITE(numout,*) ' Kolmogoroff dissipation coef. rn_ediss = ', rn_ediss 686 662 WRITE(numout,*) ' tke surface input coef. rn_ebb = ', rn_ebb 687 WRITE(numout,*) ' tke diffusion coef. rn_efave = ', rn_efave688 663 WRITE(numout,*) ' minimum value of tke rn_emin = ', rn_emin 689 664 WRITE(numout,*) ' surface minimum value of tke rn_emin0 = ', rn_emin0 690 WRITE(numout,*) ' number of restart iter loops nn_itke = ', nn_itke691 665 WRITE(numout,*) ' mixing length type nn_mxl = ', nn_mxl 692 666 WRITE(numout,*) ' prandl number flag nn_pdl = ', nn_pdl 693 WRITE(numout,*) ' horizontal average flag nn_ave = ', nn_ave694 WRITE(numout,*) ' critic Richardson nb rn_cri = ', rn_cri695 WRITE(numout,*) ' and its associated coeff. eboost = ', eboost696 667 WRITE(numout,*) ' constant background or profile nn_avb = ', nn_avb 697 668 WRITE(numout,*) ' surface mixing length = F(stress) or not ln_mxl0 = ', ln_mxl0 … … 705 676 WRITE(numout,*) ' coef to compute verticla velocity of LC rn_lc = ', rn_lc 706 677 WRITE(numout,*) 678 WRITE(numout,*) ' critical Richardson nb with your choice of coefs. = ', ri_cri 707 679 ENDIF 708 680 709 ! Check of some namelist values681 ! !* Check of some namelist values 710 682 IF( nn_mxl < 0 .OR. nn_mxl > 3 ) CALL ctl_stop( 'bad flag: nn_mxl is 0, 1 or 2 ' ) 711 683 IF( nn_pdl < 0 .OR. nn_pdl > 1 ) CALL ctl_stop( 'bad flag: nn_pdl is 0 or 1 ' ) 712 IF( nn_ave < 0 .OR. nn_ave > 1 ) CALL ctl_stop( 'bad flag: nn_ave is 0 or 1 ' )713 684 IF( nn_htau < 0 .OR. nn_htau > 2 ) CALL ctl_stop( 'bad flag: nn_htau is 0, 1 or 2 ' ) 714 685 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 ' ) 715 686 716 IF( nn_etau == 2 ) CALL zdf_mxl( nit000 ) ! Initialization of nmln 717 718 ! Horizontal average : initialization of weighting arrays 719 ! ------------------- 720 ! 721 IF(lwp) WRITE(numout,*) ' no horizontal average on avt, avmu, avmv' 722 723 ! Background eddy viscosity and diffusivity profil 724 ! ------------------------------------------------ 725 IF( nn_avb == 0 ) THEN ! Define avmb, avtb from namelist parameter 687 IF( nn_etau == 2 ) CALL zdf_mxl( nit000 ) ! Initialization of nmln 688 689 ! !* Background eddy viscosity and diffusivity profil 690 IF( nn_avb == 0 ) THEN ! Define avmb, avtb from namelist parameter 726 691 avmb(:) = avm0 727 692 avtb(:) = avt0 728 ELSE ! Background profile of avt (fit a theoretical/observational profile (Krauss 1990)693 ELSE ! Background profile of avt (fit a theoretical/observational profile (Krauss 1990) 729 694 avmb(:) = avm0 730 695 avtb(:) = avt0 + ( 3.0e-4 - 2 * avt0 ) * 1.0e-4 * gdepw_0(:) ! m2/s … … 732 697 ENDIF 733 698 ! 734 ! ! 2D shape of the avtb735 avtb_2d(:,:) = 1.e0 ! uniform736 ! 737 IF( nn_havtb == 1 ) THEN ! decrease avtb in the equatorial band699 ! ! 2D shape of the avtb 700 avtb_2d(:,:) = 1.e0 ! uniform 701 ! 702 IF( nn_havtb == 1 ) THEN ! decrease avtb in the equatorial band 738 703 ! -15S -5S : linear decrease from avt0 to avt0/10. 739 704 ! -5S +5N : cst value avt0/10. … … 744 709 ENDIF 745 710 746 ! Initialization of vertical eddy coef. to the background value 747 ! ------------------------------------------------------------- 711 ! !* depth of penetration of surface tke 712 IF( nn_etau /= 0 ) THEN 713 SELECT CASE( nn_htau ) ! Choice of the depth of penetration 714 CASE( 0 ) ! constant depth penetration (here 10 meters) 715 htau(:,:) = 10.e0 716 CASE( 1 ) ! F(latitude) : 5m to 40m at high lat. 717 DO jj = 1, jpj 718 DO ji = 1, jpi 719 htau(ji,jj) = MAX( 5., MIN( 40., 60.*ABS( SIN( rpi/180. * gphit(ji,jj) ) ) ) ) 720 END DO 721 END DO 722 CASE( 2 ) ! F(latitude) : 5m to 60m at high lat. 723 DO jj = 1, jpj 724 DO ji = 1, jpi 725 htau(ji,jj) = MAX( 5.,6./4.* MIN( 40., 60.*ABS( SIN( rpi/180. * gphit(ji,jj) ) ) ) ) 726 END DO 727 END DO 728 END SELECT 729 ENDIF 730 731 ! !* set vertical eddy coef. to the background value 748 732 DO jk = 1, jpk 749 733 avt (:,:,jk) = avtb(jk) * tmask(:,:,jk) … … 753 737 END DO 754 738 dissl(:,:,:) = 1.e-12 755 756 ! read or initialize turbulent kinetic energy ( en ) 757 ! ------------------------------------------------- 739 ! !* read or initialize all required files 758 740 CALL tke2_rst( nit000, 'READ' ) 759 741 ! 760 END SUBROUTINE zdf_tke2_init742 END SUBROUTINE tke_init 761 743 762 744 … … 796 778 CALL iom_get( numror, jpdom_autoglo, 'avmv' , avmv ) 797 779 CALL iom_get( numror, jpdom_autoglo, 'dissl', dissl ) 798 ELSE ! one at least array is missing799 CALL tke_ mlmc ! recompute avt, avm, avmu, avmv and dissl (approximation)780 ELSE ! one at least array is missing 781 CALL tke_avn ! compute avt, avm, avmu, avmv and dissl (approximation) 800 782 ENDIF 801 783 ELSE ! No TKE array found: initialisation 802 784 IF(lwp) WRITE(numout,*) ' ===>>>> : previous run without tke scheme, en computed by iterative loop' 803 785 en (:,:,:) = rn_emin * tmask(:,:,:) 804 CALL tke_ mlmc! recompute avt, avm, avmu, avmv and dissl (approximation)786 CALL tke_avn ! recompute avt, avm, avmu, avmv and dissl (approximation) 805 787 DO jit = nit000 + 1, nit000 + 10 ; CALL zdf_tke2( jit ) ; END DO 806 788 ENDIF … … 838 820 WRITE(*,*) 'zdf_tke2: You should not have seen this print! error?', kt 839 821 END SUBROUTINE zdf_tke2 822 SUBROUTINE tke2_rst( kt, cdrw ) 823 CHARACTER(len=*) :: cdrw 824 WRITE(*,*) 'tke2_rst: You should not have seen this print! error?', kt, cdwr 825 END SUBROUTINE tke2_rst 840 826 #endif 841 827
Note: See TracChangeset
for help on using the changeset viewer.
