Changeset 2102 for trunk/LMDZ.GENERIC
- Timestamp:
- Feb 18, 2019, 11:40:47 AM (6 years ago)
- Location:
- trunk/LMDZ.GENERIC
- Files:
-
- 3 edited
-
README (modified) (1 diff)
-
libf/phystd/thermcell_dq.F90 (modified) (16 diffs)
-
libf/phystd/thermcell_main.F90 (modified) (21 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/LMDZ.GENERIC/README
r2101 r2102 1455 1455 - Remove useless variable zmax0 in thermcell_main, thermcell_height and physiq_mod. 1456 1456 - Some minor changes in thermcell_plume and thermcell_main. 1457 1458 ==18/02/2018 == AB 1459 - use detr as thermcell_dq argument (called in thermcell_main) and clean up thermcell_dq -
trunk/LMDZ.GENERIC/libf/phystd/thermcell_dq.F90
r2060 r2102 2 2 ! 3 3 ! 4 SUBROUTINE thermcell_dq(ngrid,nlay,impl,ptimestep,fm,entr, & 5 & masse,q,dq,qa,lmin,lev_out) 6 7 8 USE print_control_mod, ONLY: prt_level 9 10 IMPLICIT NONE 4 SUBROUTINE thermcell_dq(ngrid,nlay,impl,ptimestep,fm,entr,detr,masse, & 5 q,dq,qa,lmin,lev_out) 6 11 7 12 8 !============================================================================== … … 25 21 !============================================================================== 26 22 23 USE print_control_mod, ONLY: prt_level 24 25 IMPLICIT NONE 26 27 27 28 !============================================================================== 28 29 ! Declaration … … 32 33 ! ------- 33 34 34 INTEGER ngrid,nlay,impl 35 INTEGER ngrid, nlay 36 INTEGER impl 35 37 INTEGER lmin(ngrid) 36 38 INTEGER lev_out ! niveau pour les print 37 39 38 40 REAL ptimestep 39 REAL masse(ngrid,nlay),fm(ngrid,nlay+1) 41 REAL masse(ngrid,nlay) 42 REAL fm(ngrid,nlay+1) 40 43 REAL entr(ngrid,nlay) 44 REAL detr(ngrid,nlay) 41 45 REAL q(ngrid,nlay) 42 46 REAL qa(ngrid,nlay) 43 REAL detr(ngrid,nlay)44 47 45 48 ! outputs: … … 51 54 ! ------ 52 55 53 INTEGER ig, k54 INTEGER niter, iter56 INTEGER ig, l 57 INTEGER niter, iter 55 58 56 59 REAL cfl 57 REAL qold(ngrid,nlay),fqa(ngrid,nlay+1) 60 REAL qold(ngrid,nlay) 61 REAL fqa(ngrid,nlay+1) 58 62 REAL wqd(ngrid,nlay+1) 59 63 REAL zzm … … 78 82 79 83 !------------------------------------------------------------------------------ 80 ! C alcul du critere CFL pour l'advection dans la subsidence84 ! CFL criterion computation for advection in downdraft 81 85 !------------------------------------------------------------------------------ 82 86 83 87 cfl = 0. 84 88 85 DO k=1,nlay89 DO l=1,nlay 86 90 DO ig=1,ngrid 87 zzm = masse(ig, k) / ptimestep88 cfl = max(cfl, fm(ig, k) / zzm)91 zzm = masse(ig,l) / ptimestep 92 cfl = max(cfl, fm(ig,l) / zzm) 89 93 90 IF (entr(ig, k).gt.zzm) THEN94 IF (entr(ig,l).gt.zzm) THEN 91 95 print *, 'ERROR: entrainment is greater than the layer mass!' 92 print *, 'entr*dt,mass', entr(ig,k)*ptimestep, masse(ig,k) 93 94 print *, 'ig,l', ig, k 95 print *, 'fm-', fm(ig,k) 96 print *, 'entr,detr', entr(ig,k), detr(ig,k) 97 print *, 'fm+', fm(ig,k+1) 98 96 print *, 'ig,l,entr', ig, l, entr(ig,l) 97 print *, '-------------------------------' 98 print *, 'entr*dt,mass', entr(ig,l)*ptimestep, masse(ig,l) 99 print *, '-------------------------------' 100 print *, 'fm+', fm(ig,l+1) 101 print *, 'entr,detr', entr(ig,l), detr(ig,l) 102 print *, 'fm ', fm(ig,l) 103 print *, 'entr,detr', entr(ig,l-1), detr(ig,l-1) 104 print *, 'fm-', fm(ig,l-1) 99 105 abort_message = 'entr dt > m, 1st' 100 106 CALL abort_physic(modname,abort_message,1) … … 104 110 105 111 !------------------------------------------------------------------------------ 106 ! Detrainement computation 107 !------------------------------------------------------------------------------ 108 109 DO k=1,nlay 110 DO ig=1,ngrid 111 detr(ig,k) = fm(ig,k) - fm(ig,k+1) + entr(ig,k) 112 ! Computation of tracer concentrations in the ascending plume 113 !------------------------------------------------------------------------------ 114 115 DO ig=1,ngrid 116 DO l=1,lmin(ig) 117 qa(ig,l) = q(ig,l) 118 ENDDO 119 ENDDO 120 121 DO ig=1,ngrid 122 DO l=lmin(ig)+1,nlay 123 if ((fm(ig,l+1)+detr(ig,l))*ptimestep.gt.1.e-5*masse(ig,l)) then 124 qa(ig,l) = (fm(ig,l) * qa(ig,l-1) + entr(ig,l) * q(ig,l)) & 125 & / (fm(ig,l+1) + detr(ig,l)) 126 else 127 qa(ig,l) = q(ig,l) 128 endif 112 129 113 IF (detr(ig,k).lt.0.) THEN 114 entr(ig,k) = entr(ig,k) - detr(ig,k) 115 detr(ig,k) = 0. 116 ! print *, 'WARNING: detr2 is negative!' 117 ! print *, 'detr', detr(ig,k) 118 ENDIF 119 120 ! IF (qa(ig,k).lt.0.) THEN 121 ! print *, 'WARNING: fm2 is negative!' 122 ! print *, 'fm', fm(ig,k) 130 ! IF (qa(ig,l).lt.0.) THEN 131 ! print *, 'WARNING: qa is negative!' 132 ! print *, 'qa', qa(ig,l) 123 133 ! ENDIF 124 134 125 ! IF (q(ig, k).lt.0.) THEN126 ! print *, 'WARNING: entr2is negative!'127 ! print *, ' entr', entr(ig,k)135 ! IF (q(ig,l).lt.0.) THEN 136 ! print *, 'WARNING: q is negative!' 137 ! print *, 'q', q(ig,l) 128 138 ! ENDIF 129 139 ENDDO … … 131 141 132 142 !------------------------------------------------------------------------------ 133 ! Computation of tracer concentrations in the ascending plume134 !------------------------------------------------------------------------------135 136 DO ig=1,ngrid137 DO k=1,lmin(ig)138 qa(ig,k) = q(ig,k)139 ENDDO140 ENDDO141 142 DO ig=1,ngrid143 DO k=lmin(ig)+1,nlay144 if ((fm(ig,k+1)+detr(ig,k))*ptimestep.gt.1.e-5*masse(ig,k)) then145 qa(ig,k) = (fm(ig,k) * qa(ig,k-1) + entr(ig,k) * q(ig,k)) &146 & / (fm(ig,k+1) + detr(ig,k))147 else148 qa(ig,k) = q(ig,k)149 endif150 151 ! IF (qa(ig,k).lt.0.) THEN152 ! print *, 'WARNING: qa is negative!'153 ! print *, 'qa', qa(ig,k)154 ! ENDIF155 156 ! IF (q(ig,k).lt.0.) THEN157 ! print *, 'WARNING: q is negative!'158 ! print *, 'q', q(ig,k)159 ! ENDIF160 ENDDO161 ENDDO162 163 !------------------------------------------------------------------------------164 143 ! Plume vertical flux of q 165 144 !------------------------------------------------------------------------------ 166 145 167 DO k=2,nlay-1168 fqa(:, k) = fm(:,k) * qa(:,k-1)146 DO l=2,nlay-1 147 fqa(:,l) = fm(:,l) * qa(:,l-1) 169 148 ENDDO 170 149 … … 177 156 178 157 IF (impl==0) THEN 179 do k=1,nlay-1180 q(:, k) = q(:,k) + (fqa(:,k) - fqa(:,k+1) - fm(:,k) * q(:,k) &181 & + fm(:, k+1) * q(:,k+1)) * ptimestep / masse(:,k)158 do l=1,nlay-1 159 q(:,l) = q(:,l) + (fqa(:,l) - fqa(:,l+1) - fm(:,l) * q(:,l) & 160 & + fm(:,l+1) * q(:,l+1)) * ptimestep / masse(:,l) 182 161 enddo 183 162 ELSE 184 do k=nlay-1,1,-1185 q(:, k) = ( q(:,k) + ptimestep / masse(:,k) &186 & * ( fqa(:, k) - fqa(:,k+1) + fm(:,k+1) * q(:,k+1) ) ) &187 & / ( 1. + fm(:, k) * ptimestep / masse(:,k) )163 do l=nlay-1,1,-1 164 q(:,l) = ( q(:,l) + ptimestep / masse(:,l) & 165 & * ( fqa(:,l) - fqa(:,l+1) + fm(:,l+1) * q(:,l+1) ) ) & 166 & / ( 1. + fm(:,l) * ptimestep / masse(:,l) ) 188 167 ENDDO 189 168 ENDIF … … 193 172 !============================================================================== 194 173 195 DO k=1,nlay174 DO l=1,nlay 196 175 DO ig=1,ngrid 197 dq(ig, k) = (q(ig,k) - qold(ig,k)) / ptimestep198 q(ig, k) = qold(ig,k)176 dq(ig,l) = (q(ig,l) - qold(ig,l)) / ptimestep 177 q(ig,l) = qold(ig,l) 199 178 ENDDO 200 179 ENDDO … … 205 184 206 185 207 208 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 209 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 210 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 211 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 212 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 213 214 215 216 ! Obsolete version kept for convergence with Cmip5 NPv3.1 simulations 217 218 219 subroutine thermcell_dq_o(ngrid,nlay,impl,ptimestep,fm,entr, & 220 & masse,q,dq,qa,lev_out) 221 222 USE print_control_mod, ONLY: prt_level 223 224 implicit none 225 226 !======================================================================= 186 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 187 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 188 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 189 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 190 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 191 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 192 193 194 Subroutine thermcell_dq_o(ngrid,nlay,impl,ptimestep,fm,entr,masse, & 195 q,dq,qa,lev_out) 196 197 198 !============================================================================== 227 199 ! 228 200 ! Calcul du transport verticale dans la couche limite en presence … … 230 202 ! calcul du dq/dt une fois qu'on connait les ascendances 231 203 ! 232 !======================================================================= 233 204 !============================================================================== 205 206 USE print_control_mod, ONLY: prt_level 207 208 implicit none 209 210 211 !============================================================================== 212 ! Declaration 213 !============================================================================== 214 234 215 integer ngrid,nlay,impl 235 216 … … 245 226 real zzm 246 227 247 integer ig, k228 integer ig,l 248 229 real cfl 249 230 … … 258 239 ! Calcul du critere CFL pour l'advection dans la subsidence 259 240 cfl = 0. 260 do k=1,nlay241 do l=1,nlay 261 242 do ig=1,ngrid 262 zzm=masse(ig, k)/ptimestep263 cfl=max(cfl,fm(ig, k)/zzm)264 if (entr(ig, k).gt.zzm) then265 print*,'entr*dt>m,2', k,entr(ig,k)*ptimestep,masse(ig,k)243 zzm=masse(ig,l)/ptimestep 244 cfl=max(cfl,fm(ig,l)/zzm) 245 if (entr(ig,l).gt.zzm) then 246 print*,'entr*dt>m,2',l,entr(ig,l)*ptimestep,masse(ig,l) 266 247 abort_message = 'entr dt > m, 2nd' 267 248 CALL abort_physic (modname,abort_message,1) … … 269 250 enddo 270 251 enddo 271 252 272 253 !IM 090508 print*,'CFL CFL CFL CFL ',cfl 273 254 274 255 #undef CFL 275 256 #ifdef CFL … … 279 260 niter=1 280 261 #endif 281 262 282 263 ztimestep=ptimestep/niter 283 264 qold=q 284 285 286 do iter=1,niter 287 if (prt_level.ge.1) print*,'Q2 THERMCEL_DQ 0' 288 265 266 do iter=1,niter 267 if (prt_level.ge.1) print*,'Q2 THERMCEL_DQ 0' 268 289 269 ! calcul du detrainement 290 do k=1,nlay291 do ig=1,ngrid292 detr(ig,k)=fm(ig,k)-fm(ig,k+1)+entr(ig,k)293 ! print*,'Q2 DQ ',detr(ig,k),fm(ig,k),entr(ig,k)270 do l=1,nlay 271 do ig=1,ngrid 272 detr(ig,l)=fm(ig,l)-fm(ig,l+1)+entr(ig,l) 273 ! print*,'Q2 DQ ',detr(ig,l),fm(ig,l),entr(ig,l) 294 274 !test 295 if (detr(ig,k).lt.0.) then 296 entr(ig,k)=entr(ig,k)-detr(ig,k) 297 detr(ig,k)=0. 298 ! print*,'detr2<0!!!','ig=',ig,'k=',k,'f=',fm(ig,k), 299 ! s 'f+1=',fm(ig,k+1),'e=',entr(ig,k),'d=',detr(ig,k) 300 endif 301 if (fm(ig,k+1).lt.0.) then 302 ! print*,'fm2<0!!!' 303 endif 304 if (entr(ig,k).lt.0.) then 305 ! print*,'entr2<0!!!' 306 endif 307 enddo 308 enddo 309 310 ! calcul de la valeur dans les ascendances 311 do ig=1,ngrid 312 qa(ig,1)=q(ig,1) 313 enddo 314 315 do k=2,nlay 316 do ig=1,ngrid 317 if ((fm(ig,k+1)+detr(ig,k))*ztimestep.gt. & 318 & 1.e-5*masse(ig,k)) then 319 qa(ig,k)=(fm(ig,k)*qa(ig,k-1)+entr(ig,k)*q(ig,k)) & 320 & /(fm(ig,k+1)+detr(ig,k)) 321 else 322 qa(ig,k)=q(ig,k) 323 endif 324 if (qa(ig,k).lt.0.) then 325 ! print*,'qa<0!!!' 326 endif 327 if (q(ig,k).lt.0.) then 328 ! print*,'q<0!!!' 329 endif 330 enddo 331 enddo 332 275 if (detr(ig,l).lt.0.) then 276 entr(ig,l)=entr(ig,l)-detr(ig,l) 277 detr(ig,l)=0. 278 ! print*,'detr2<0!!!','ig=',ig,'l=',l,'f=',fm(ig,l), 279 ! s 'f+1=',fm(ig,l+1),'e=',entr(ig,l),'d=',detr(ig,l) 280 endif 281 282 ! if (fm(ig,l+1).lt.0.) then 283 ! print*,'fm2<0!!!' 284 ! endif 285 286 ! if (entr(ig,l).lt.0.) then 287 ! print*,'entr2<0!!!' 288 ! endif 289 enddo 290 enddo 291 292 ! calcul de la valeur dans les ascendances 293 do ig=1,ngrid 294 qa(ig,1)=q(ig,1) 295 enddo 296 297 do l=2,nlay 298 do ig=1,ngrid 299 if ((fm(ig,l+1)+detr(ig,l))*ztimestep.gt.1.e-5*masse(ig,l)) then 300 qa(ig,l) = (fm(ig,l)*qa(ig,l-1)+entr(ig,l)*q(ig,l)) & 301 & / (fm(ig,l+1)+detr(ig,l)) 302 else 303 qa(ig,l)=q(ig,l) 304 endif 305 306 if (qa(ig,l).lt.0.) then 307 ! print*,'qa<0!!!' 308 endif 309 310 if (q(ig,l).lt.0.) then 311 ! print*,'q<0!!!' 312 endif 313 enddo 314 enddo 315 333 316 ! Calcul du flux subsident 334 335 do k=2,nlay336 do ig=1,ngrid317 318 do l=2,nlay 319 do ig=1,ngrid 337 320 #undef centre 338 321 #ifdef centre 339 wqd(ig,k)=fm(ig,k)*0.5*(q(ig,k-1)+q(ig,k))322 wqd(ig,l)=fm(ig,l)*0.5*(q(ig,l-1)+q(ig,l)) 340 323 #else 341 324 … … 343 326 #ifdef plusqueun 344 327 ! Schema avec advection sur plus qu'une maille. 345 zzm=masse(ig,k)/ztimestep346 if (fm(ig,k)>zzm) then347 wqd(ig,k)=zzm*q(ig,k)+(fm(ig,k)-zzm)*q(ig,k+1)348 else349 wqd(ig,k)=fm(ig,k)*q(ig,k)350 endif328 zzm=masse(ig,l)/ztimestep 329 if (fm(ig,l)>zzm) then 330 wqd(ig,l)=zzm*q(ig,l)+(fm(ig,l)-zzm)*q(ig,l+1) 331 else 332 wqd(ig,l)=fm(ig,l)*q(ig,l) 333 endif 351 334 #else 352 wqd(ig,k)=fm(ig,k)*q(ig,k)335 wqd(ig,l)=fm(ig,l)*q(ig,l) 353 336 #endif 354 337 #endif 355 338 356 if (wqd(ig,k).lt.0.) then 357 ! print*,'wqd<0!!!' 358 endif 339 if (wqd(ig,l).lt.0.) then 340 ! print*,'wqd<0!!!' 341 endif 342 enddo 343 enddo 344 345 do ig=1,ngrid 346 wqd(ig,1)=0. 347 wqd(ig,nlay+1)=0. 348 enddo 349 350 ! Calcul des tendances 351 do l=1,nlay 352 do ig=1,ngrid 353 q(ig,l)=q(ig,l)+(detr(ig,l)*qa(ig,l)-entr(ig,l)*q(ig,l) & 354 & -wqd(ig,l)+wqd(ig,l+1)) & 355 & *ztimestep/masse(ig,l) 356 ! if (dq(ig,l).lt.0.) then 357 ! print*,'dq<0!!!' 358 ! endif 359 enddo 359 360 enddo 360 361 enddo 361 do ig=1,ngrid 362 wqd(ig,1)=0. 363 wqd(ig,nlay+1)=0. 362 363 ! Calcul des tendances 364 do l=1,nlay 365 do ig=1,ngrid 366 dq(ig,l) = (q(ig,l) - qold(ig,l)) / ptimestep 367 q(ig,l) = qold(ig,l) 368 enddo 364 369 enddo 365 366 367 ! Calcul des tendances 368 do k=1,nlay 369 do ig=1,ngrid 370 q(ig,k)=q(ig,k)+(detr(ig,k)*qa(ig,k)-entr(ig,k)*q(ig,k) & 371 & -wqd(ig,k)+wqd(ig,k+1)) & 372 & *ztimestep/masse(ig,k) 373 ! if (dq(ig,k).lt.0.) then 374 ! print*,'dq<0!!!' 375 ! endif 376 enddo 377 enddo 378 379 380 enddo 381 382 383 ! Calcul des tendances 384 do k=1,nlay 385 do ig=1,ngrid 386 dq(ig,k)=(q(ig,k)-qold(ig,k))/ptimestep 387 q(ig,k)=qold(ig,k) 388 enddo 389 enddo 390 391 return 392 end 370 371 372 return 373 end -
trunk/LMDZ.GENERIC/libf/phystd/thermcell_main.F90
r2101 r2102 2 2 ! 3 3 ! 4 SUBROUTINE thermcell_main(itap,ngrid,nlay,ptimestep 5 ,pplay,pplev,pphi,firstcall&6 ,pu,pv,pt,po&7 ,pduadj,pdvadj,pdtadj,pdoadj&8 ,f0,fm0,entr0,detr0&9 ,zqta,zqla,ztv,ztva,ztla,zthl,zqsatth&10 ,zw2,fraca&11 ,lmin,lmix,lalim,lmax&12 ,zpspsk,ratqscth,ratqsdiff&13 ,Ale_bl,Alp_bl,lalim_conv,wght_th&4 SUBROUTINE thermcell_main(itap,ngrid,nlay,ptimestep, & 5 pplay,pplev,pphi,firstcall, & 6 pu,pv,pt,po, & 7 pduadj,pdvadj,pdtadj,pdoadj, & 8 f0,fm0,entr0,detr0, & 9 zqta,zqla,ztv,ztva,ztla,zthl,zqsatth, & 10 zw2,fraca, & 11 lmin,lmix,lalim,lmax, & 12 zpspsk,ratqscth,ratqsdiff, & 13 Ale_bl,Alp_bl,lalim_conv,wght_th, & 14 14 !!! nrlmd le 10/04/2012 15 ,pbl_tke,pctsrf,omega,airephy&16 ,zlcl,fraca0,w0,w_conv,therm_tke_max0,env_tke_max0&17 ,n2,s2,ale_bl_stat&18 ,therm_tke_max,env_tke_max&19 ,alp_bl_det,alp_bl_fluct_m,alp_bl_fluct_tke&20 ,alp_bl_conv,alp_bl_stat)15 pbl_tke,pctsrf,omega,airephy, & 16 zlcl,fraca0,w0,w_conv,therm_tke_max0,env_tke_max0, & 17 n2,s2,ale_bl_stat, & 18 therm_tke_max,env_tke_max, & 19 alp_bl_det,alp_bl_fluct_m,alp_bl_fluct_tke, & 20 alp_bl_conv,alp_bl_stat) 21 21 !!! fin nrlmd le 10/04/2012 22 22 … … 213 213 !------Sorties 214 214 real zlcl(ngrid),fraca0(ngrid),w0(ngrid),w_conv(ngrid) 215 real therm_tke_max0(ngrid),env_tke_max0(ngrid) 215 real therm_tke_max0(ngrid) 216 real env_tke_max0(ngrid) 216 217 real n2(ngrid),s2(ngrid) 217 218 real ale_bl_stat(ngrid) 218 real therm_tke_max(ngrid,nlay),env_tke_max(ngrid,nlay) 219 real alp_bl_det(ngrid),alp_bl_fluct_m(ngrid),alp_bl_fluct_tke(ngrid),alp_bl_conv(ngrid),alp_bl_stat(ngrid) 219 real therm_tke_max(ngrid,nlay) 220 real env_tke_max(ngrid,nlay) 221 real alp_bl_det(ngrid) 222 real alp_bl_fluct_m(ngrid) 223 real alp_bl_fluct_tke(ngrid) 224 real alp_bl_conv(ngrid) 225 real alp_bl_stat(ngrid) 220 226 !------Local 221 227 integer nsrf … … 225 231 real interp(ngrid) ! Coef d'interpolation pour le LCL 226 232 !------Triggering 227 real Su ! Surface unite: celle d'un updraft elementaire 228 parameter(Su=4e4) 229 real hcoef ! Coefficient directeur pour le calcul de s2 230 parameter(hcoef=1) 231 real hmincoef ! Coefficient directeur pour l'ordonnee a l'origine pour le calcul de s2 232 parameter(hmincoef=0.3) 233 real eps1 ! Fraction de surface occupee par la population 1 : eps1=n1*s1/(fraca0*Sd) 234 parameter(eps1=0.3) 233 real,parameter :: Su=4e4 ! Surface unite: celle d'un updraft elementaire 234 real,parameter :: hcoef ! Coefficient directeur pour le calcul de s2 235 real,parameter :: hmincoef=0.3 ! Coefficient directeur pour l'ordonnee a l'origine pour le calcul de s2 (hmincoef=0.3) 236 real,parameter :: eps1=0.3 ! Fraction de surface occupee par la population 1 : eps1=n1*s1/(fraca0*Sd) 235 237 real hmin(ngrid) ! Ordonnee a l'origine pour le calcul de s2 236 238 real zmax_moy(ngrid) ! Hauteur moyenne des thermiques : zmax_moy = zlcl + 0.33 (zmax-zlcl) 237 real zmax_moy_coef 238 parameter(zmax_moy_coef=0.33) 239 real,parameter :: zmax_moy_coef=0.33 239 240 real depth(ngrid) ! Epaisseur moyenne du cumulus 240 241 real w_max(ngrid) ! Vitesse max statistique … … 244 245 real pbl_tke_max0(ngrid) ! TKE moyenne au LCL 245 246 real w_ls(ngrid,nlay) ! Vitesse verticale grande echelle (m/s) 246 real coef_m ! On considere un rendement pour alp_bl_fluct_m 247 parameter(coef_m=1.) 248 real coef_tke ! On considere un rendement pour alp_bl_fluct_tke 249 parameter(coef_tke=1.) 247 real,parameter :: coef_m=1. ! On considere un rendement pour alp_bl_fluct_m 248 real,parameter :: coef_tke=1. ! On considere un rendement pour alp_bl_fluct_tke 250 249 251 250 !!! fin nrlmd le 10/04/2012 252 251 253 252 ! Nouvelles variables pour la convection 253 integer lalim_conv(ngrid) 254 integer n_int(ngrid) 254 255 real Ale_bl(ngrid) 255 256 real Alp_bl(ngrid) 256 real alp_int(ngrid),dp_int(ngrid),zdp 257 real alp_int(ngrid) 258 real dp_int(ngrid),zdp 257 259 real ale_int(ngrid) 258 integer n_int(ngrid)259 260 real fm_tot(ngrid) 260 261 real wght_th(ngrid,nlay) 261 integer lalim_conv(ngrid)262 262 263 263 CHARACTER*2 str2 … … 515 515 IF (tau_thermals>1.) THEN 516 516 lambda = exp(-ptimestep/tau_thermals) 517 f0 = (1.-lambda) * f + lambda * f0517 f0(:) = (1.-lambda) * f(:) + lambda * f0(:) 518 518 ELSE 519 519 f0(:) = f(:) … … 527 527 print *, 'f0 =', f0(1) 528 528 CALL abort_physic(modname,abort_message,1) 529 ELSE530 print *, 'f0 =', f0(1)531 529 ENDIF 532 530 … … 587 585 !------------------------------------------------------------------------------ 588 586 589 CALL thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0, masse,&587 CALL thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,detr0,masse, & 590 588 & zthl,zdthladj,zta,lmin,lev_out) 591 589 592 CALL thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0, masse,&590 CALL thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,detr0,masse, & 593 591 & po,pdoadj,zoa,lmin,lev_out) 594 592 … … 614 612 ! Calcul purement conservatif pour le transport de V=(u,v) 615 613 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 616 CALL thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0, masse,&614 CALL thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,detr0,masse, & 617 615 & zu,pduadj,zua,lmin,lev_out) 618 616 619 CALL thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0, masse,&617 CALL thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,detr0,masse, & 620 618 & zv,pdvadj,zva,lmin,lev_out) 621 619 ENDIF … … 631 629 !------------------------------------------------------------------------------ 632 630 633 if (prt_level.ge.1) print*,'14a OK convect8' 634 635 do ig=1,ngrid 631 DO ig=1,ngrid 636 632 nivcon(ig) = 0 637 633 zcon(ig) = 0. 638 enddo634 ENDDO 639 635 !nouveau calcul 640 636 do ig=1,ngrid 641 637 CHI = zh(ig,1)/(1669.0-122.0*zo(ig,1)/zqsat(ig,1)-zh(ig,1)) 642 638 pcon(ig) = pplay(ig,1)*(zo(ig,1)/zqsat(ig,1))**CHI 643 enddo639 ENDDO 644 640 645 641 !IM do k=1,nlay 646 642 do k=1,nlay-1 647 643 do ig=1,ngrid 648 if ((pcon(ig).le.pplay(ig,k)).and.(pcon(ig).gt.pplay(ig,k+1))) then 649 zcon2(ig) = zlay(ig,k)-(pcon(ig)-pplay(ig,k))/(RG*rho(ig,k))/100. 650 endif 651 enddo 652 enddo 644 IF ((pcon(ig).le.pplay(ig,k)).and.(pcon(ig).gt.pplay(ig,k+1))) then 645 zcon2(ig) = zlay(ig,k) - (pcon(ig) - pplay(ig,k)) & 646 / (RG * rho(ig,k)) / 100. 647 ENDIF 648 ENDDO 649 ENDDO 653 650 654 651 ierr = 0 655 652 656 653 do ig=1,ngrid 657 if (pcon(ig).le.pplay(ig,nlay)) then 658 zcon2(ig) = zlay(ig,nlay)-(pcon(ig)-pplay(ig,nlay))/(RG*rho(ig,nlay))/100. 654 IF (pcon(ig).le.pplay(ig,nlay)) then 655 zcon2(ig) = zlay(ig,nlay) - (pcon(ig) - pplay(ig,nlay)) & 656 / (RG * rho(ig,nlay)) / 100. 659 657 ierr = 1 660 endif661 enddo662 663 if(ierr==1) then658 ENDIF 659 ENDDO 660 661 IF (ierr==1) then 664 662 write(*,*) 'ERROR: thermal plumes rise the model top' 665 663 CALL abort 666 endif667 668 if(prt_level.ge.1) print*,'14b OK convect8'664 ENDIF 665 666 IF (prt_level.ge.1) print*,'14b OK convect8' 669 667 670 668 do k=nlay,1,-1 671 669 do ig=1,ngrid 672 if(zqla(ig,k).gt.1e-10) then670 IF (zqla(ig,k).gt.1e-10) then 673 671 nivcon(ig) = k 674 672 zcon(ig) = zlev(ig,k) 675 endif676 enddo677 enddo678 679 if(prt_level.ge.1) print*,'14c OK convect8'673 ENDIF 674 ENDDO 675 ENDDO 676 677 IF (prt_level.ge.1) print*,'14c OK convect8' 680 678 681 679 !------------------------------------------------------------------------------ … … 690 688 ratqscth(ig,l) = 0. 691 689 ratqsdiff(ig,l) = 0. 692 enddo693 enddo694 695 if(prt_level.ge.1) print*,'14d OK convect8'690 ENDDO 691 ENDDO 692 693 IF (prt_level.ge.1) print*,'14d OK convect8' 696 694 697 695 do l=1,nlay … … 701 699 thetath2(ig,l) = zf2*(ztla(ig,l)-zthl(ig,l))**2 702 700 703 if(zw2(ig,l).gt.1.e-10) then701 IF (zw2(ig,l).gt.1.e-10) then 704 702 wth2(ig,l) = zf2*(zw2(ig,l))**2 705 703 else 706 704 wth2(ig,l) = 0. 707 endif705 ENDIF 708 706 709 wth3(ig,l) = zf2*(1-2.*fraca(ig,l))/(1-fraca(ig,l)) 707 wth3(ig,l) = zf2*(1-2.*fraca(ig,l))/(1-fraca(ig,l)) & 710 708 & *zw2(ig,l)*zw2(ig,l)*zw2(ig,l) 711 709 q2(ig,l) = zf2*(zqta(ig,l)*1000.-po(ig,l)*1000.)**2 712 710 !test: on calcul q2/po=ratqsc 713 711 ratqscth(ig,l) = sqrt(max(q2(ig,l),1.e-6)/(po(ig,l)*1000.)) 714 enddo715 enddo712 ENDDO 713 ENDDO 716 714 717 715 !------------------------------------------------------------------------------ … … 724 722 wthl(ig,l)=fraca(ig,l)*zw2(ig,l)*(ztla(ig,l)-zthl(ig,l)) 725 723 wthv(ig,l)=fraca(ig,l)*zw2(ig,l)*(ztva(ig,l)-ztv(ig,l)) 726 enddo727 enddo728 729 CALL thermcell_alp(ngrid,nlay,ptimestep, 730 & pplay,pplev, 731 & fm0,entr0,lmax, 732 & Ale_bl,Alp_bl,lalim_conv,wght_th, 733 & zw2,fraca, 734 & pcon,rhobarz,wth3,wmax_sec,lalim,fm,alim_star,zmax,&735 & pbl_tke,pctsrf,omega,airephy, 736 & zlcl,fraca0,w0,w_conv,therm_tke_max0,env_tke_max0, 737 & n2,s2,ale_bl_stat, 738 & therm_tke_max,env_tke_max, 739 & alp_bl_det,alp_bl_fluct_m,alp_bl_fluct_tke, 724 ENDDO 725 ENDDO 726 727 CALL thermcell_alp(ngrid,nlay,ptimestep, & 728 & pplay,pplev, & 729 & fm0,entr0,lmax, & 730 & Ale_bl,Alp_bl,lalim_conv,wght_th, & 731 & zw2,fraca,pcon, & 732 & rhobarz,wth3,wmax_sec,lalim,fm,alim_star,zmax, & 733 & pbl_tke,pctsrf,omega,airephy, & 734 & zlcl,fraca0,w0,w_conv,therm_tke_max0,env_tke_max0,& 735 & n2,s2,ale_bl_stat, & 736 & therm_tke_max,env_tke_max, & 737 & alp_bl_det,alp_bl_fluct_m,alp_bl_fluct_tke, & 740 738 & alp_bl_conv,alp_bl_stat) 741 739 … … 756 754 ENDDO 757 755 758 if(prt_level.ge.1) print*,'14f OK convect8'756 IF (prt_level.ge.1) print*,'14f OK convect8' 759 757 760 758 DO l=1,nlay … … 763 761 zf = fraca(ig,l) 764 762 zf2 = zf / (1.-zf) 765 vardiff = vardiff + alim_star(ig,l) * (zqta(ig,l) * 1000. - var)**2 763 vardiff = vardiff + alim_star(ig,l) & 764 * (zqta(ig,l) * 1000. - var)**2 766 765 ENDIF 767 766 ENDDO … … 840 839 ! test sur la hauteur des thermiques ... 841 840 do i=1,ngrid 842 !IMtemp if(pplay(i,long(i)).lt.seuil*pplev(i,1)) then843 if(prt_level.ge.10) then841 !IMtemp IF (pplay(i,long(i)).lt.seuil*pplev(i,1)) then 842 IF (prt_level.ge.10) then 844 843 print *, 'WARNING ',comment,' au point ',i,' K= ',long(i) 845 844 print *, ' K P(MB) THV(K) Qenv(g/kg)THVA QLA(g/kg) F* W2' 846 845 do k=1,nlay 847 846 write(6,'(i3,7f10.3)') k,pplay(i,k),ztv(i,k),1000*po(i,k),ztva(i,k),1000*zqla(i,k),f_star(i,k),zw2(i,k) 848 enddo849 endif850 enddo847 ENDDO 848 ENDIF 849 ENDDO 851 850 852 851 … … 911 910 detr0(:,k) = fm0(:,k) - fm0(:,k+1) + entr0(:,k) 912 911 masse0(:,k) = (pplev(:,k) - pplev(:,k+1)) / RG 913 enddo912 ENDDO 914 913 915 914 !------------------------------------------------------------------------------ … … 927 926 detr(:,k+1)=0.5*(detr0(:,k)+detr0(:,k+1)) 928 927 fm(:,k+1)=fm(:,k)+entr(:,k)-detr(:,k) 929 enddo928 ENDDO 930 929 931 930 fm(:,nlay+1)=0. … … 935 934 ! do ig=1,ngrid 936 935 ! qa(ig,1) = q(ig,1) 937 ! enddo936 ! ENDDO 938 937 939 938 q(:,:)=therm_tke_max(:,:) … … 941 940 do ig=1,ngrid 942 941 qa(ig,1)=q(ig,1) 943 enddo944 945 if(1==1) then942 ENDDO 943 944 IF (1==1) then 946 945 do k=2,nlay 947 946 do ig=1,ngrid 948 if((fm(ig,k+1)+detr(ig,k))*ptimestep.gt.1.e-5*masse(ig,k)) then947 IF ((fm(ig,k+1)+detr(ig,k))*ptimestep.gt.1.e-5*masse(ig,k)) then 949 948 qa(ig,k) = (fm(ig,k)*qa(ig,k-1)+entr(ig,k)*q(ig,k)) & 950 949 & / (fm(ig,k+1)+detr(ig,k)) 951 950 else 952 951 qa(ig,k)=q(ig,k) 953 endif952 ENDIF 954 953 955 ! if(qa(ig,k).lt.0.) print *, 'WARNING: qa is negative!'956 ! if(q(ig,k).lt.0.) print *, 'WARNING: q is negative!'957 enddo958 enddo954 ! IF (qa(ig,k).lt.0.) print *, 'WARNING: qa is negative!' 955 ! IF (q(ig,k).lt.0.) print *, 'WARNING: q is negative!' 956 ENDDO 957 ENDDO 959 958 960 959 !------------------------------------------------------------------------------ … … 965 964 do ig=1,ngrid 966 965 wqd(ig,k) = fm(ig,k) * q(ig,k) 967 if(wqd(ig,k).lt.0.) print*,'WARNING: wqd is negative!'968 enddo969 enddo966 IF (wqd(ig,k).lt.0.) print*,'WARNING: wqd is negative!' 967 ENDDO 968 ENDDO 970 969 971 970 do ig=1,ngrid 972 971 wqd(ig,1) = 0. 973 972 wqd(ig,nlay+1) = 0. 974 enddo973 ENDDO 975 974 976 975 !------------------------------------------------------------------------------ … … 983 982 & * (detr(ig,k) * qa(ig,k) - entr(ig,k) * q(ig,k) & 984 983 & - wqd(ig,k) + wqd(ig,k+1)) 985 enddo986 enddo987 endif984 ENDDO 985 ENDDO 986 ENDIF 988 987 989 988 therm_tke_max(:,:) = q(:,:)
Note: See TracChangeset
for help on using the changeset viewer.
