Changeset 256
- Timestamp:
- Aug 3, 2011, 4:49:20 PM (14 years ago)
- Location:
- trunk/LMDZ.MARS
- Files:
-
- 6 edited
-
README (modified) (1 diff)
-
libf/phymars/meso_inc/meso_inc_les.F (modified) (1 diff)
-
libf/phymars/physiq.F (modified) (6 diffs)
-
libf/phymars/thermcell_dqupdown.F90 (modified) (1 diff)
-
libf/phymars/vdif_cd.F (modified) (5 diffs)
-
libf/phymars/vdifc.F (modified) (3 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/LMDZ.MARS/README
r249 r256 836 836 == 01/08/2011 == EM 837 837 - Fixed bug in readtesassim (timeflag was not always set before being used) 838 839 == 03/08/2011 == AC 840 - Modified physiq.F to interface new SL parametrization 841 - New SL parametrization based on a bulk Richardson Monin-Obukhov theory formulation 842 Stability functions are taken from D.E. England et al. (95) 843 Similarity functions coefficients based on Dyer and Hicks (70). 844 Includes thermal roughness length computation, heat and momentum drag coefficient computation 845 Can be used to output hydrodynamic-related SL quantities (bulk Richardson, turbulent Prandtl number estimation, Reynolds number) 846 - Minor modification in thermcell_dqupdown.F to suit picky compilers 847 - vdifc.F Now takes into account sub-grid gustiness, evaluated from thermals activity (it's proxy being the maximum vertical velocity) 848 - Minor modification in meso_inc_les.F: u* is now taken from the new vdifc and not recomputed from a simple law -
trunk/LMDZ.MARS/libf/phymars/meso_inc/meso_inc_les.F
r226 r256 3 3 hfx(ig) = zflubid(ig)-capcal(ig)*zdtsdif(ig) 4 4 !! u star in similarity theory in m/s 5 ust(ig) = 0.4 6 . * sqrt( pu(ig,1)*pu(ig,1) + pv(ig,1)*pv(ig,1) ) 7 . / log( 1.E+0 + zzlay(ig,1)/z0_default ) 5 ! ust(ig) = 0.4 6 ! . * sqrt( pu(ig,1)*pu(ig,1) + pv(ig,1)*pv(ig,1) ) 7 ! . / log( 1.E+0 + zzlay(ig,1)/z0_default ) 8 9 ! New SL parametrization, ust is more accurately computed in vdif_cd : 10 ust(ig) = zcdv(ig) 11 8 12 ENDDO 9 10 13 ! write (*,*) 'PHYS HFX cp zdts', hfx(100), zflubid(100), 11 14 ! . capcal(100), -
trunk/LMDZ.MARS/libf/phymars/physiq.F
r234 r256 321 321 322 322 REAL lmax_th_out(ngridmx),zmax_th(ngridmx) 323 REAL wmax_th(ngridmx)323 REAL, SAVE :: wmax_th(ngridmx) 324 324 REAL, SAVE :: hfmax_th(ngridmx) 325 325 REAL pdu_th(ngridmx,nlayermx),pdv_th(ngridmx,nlayermx) … … 327 327 INTEGER lmax_th(ngridmx) 328 328 REAL dtke_th(ngridmx,nlayermx+1) 329 REAL dummycol(ngridmx) 330 329 REAL zcdv(ngridmx) 331 330 c======================================================================= 332 331 … … 697 696 $ zdum1,zdum2,zdh,pdq,zflubid, 698 697 $ zdudif,zdvdif,zdhdif,zdtsdif,q2, 699 & zdqdif,zdqsdif )698 & zdqdif,zdqsdif,wmax_th,zcdv) 700 699 701 700 #ifdef MESOSCALE … … 1641 1640 1642 1641 co2col(:)=0. 1643 dummycol(:)=0.1644 1642 if (tracer) then 1645 1643 do l=1,nlayermx … … 1647 1645 co2col(ig)=co2col(ig)+ 1648 1646 & zq(ig,l,1)*(pplev(ig,l)-pplev(ig,l+1))/g 1649 if (nqmx .gt. 1) then1650 iq=2 ! to avoid out-of-bounds spotting by picky compilers1651 ! when gcm is compiled with only one tracer1652 dummycol(ig)=dummycol(ig)+1653 & zq(ig,l,iq)*(pplev(ig,l)-pplev(ig,l+1))/g1654 endif1655 1647 enddo 1656 1648 enddo … … 1659 1651 call WRITEDIAGFI(ngrid,'co2col','integrated co2 mass' & 1660 1652 & ,'kg/m-2',0,co2col) 1661 call WRITEDIAGFI(ngrid,'dummycol','integrated dummy mass' &1662 & ,'kg/m-2',0,dummycol)1663 1653 endif 1664 1654 call WRITEDIAGFI(ngrid,'w','vertical velocity' & -
trunk/LMDZ.MARS/libf/phymars/thermcell_dqupdown.F90
r185 r256 61 61 ! of fmd in the equations, so it has to be positive 62 62 ! 63 !fmd(:,:)=-fm_down(:,:)63 fmd(:,:)=-fm_down(:,:) 64 64 ! 65 65 !! ========== Entrainment, Detrainement and Mass ================= -
trunk/LMDZ.MARS/libf/phymars/vdif_cd.F
r224 r256 1 SUBROUTINE vdif_cd( ngrid,nlay,pz0,pg,pz,pu,pv,pts,ph,pcdv,pcdh) 1 SUBROUTINE vdif_cd(ngrid,nlay,pz0, 2 & pg,pz,pu,pv,pts,ph,pcdv,pcdh) 2 3 IMPLICIT NONE 3 4 c======================================================================= … … 7 8 c 8 9 c Author: Frederic Hourdin 15 /10 /93 10 c Modified by : Arnaud Colaitis 03/08/11 9 11 c ------- 10 12 c … … 16 18 c ngrid size of the horizontal grid 17 19 c pg gravity (m s -2) 18 c pz(ngrid ) height of the first atmospheric layer19 c pu(ngrid ) u component of the wind in that layer20 c pv(ngrid ) v component of the wind in that layer21 c pts(ngrid) surfac te temperature20 c pz(ngrid,nlay) height of layers 21 c pu(ngrid,nlay) u component of the wind 22 c pv(ngrid,nlay) v component of the wind 23 c pts(ngrid) surface temperature 22 24 c ph(ngrid) potential temperature T*(p/ps)^kappa 23 25 c … … 33 35 c ------------- 34 36 37 #include "comcstfi.h" 38 35 39 c Arguments: 36 40 c ---------- 37 41 38 INTEGER ngrid,nlay39 REAL pz0(ngrid)40 REAL pg,pz(ngrid,nlay)41 REAL pu(ngrid,nlay),pv(ngrid,nlay)42 REAL pts(ngrid,nlay),ph(ngrid,nlay)43 REAL pcdv(ngrid),pcdh(ngrid)42 INTEGER, INTENT(IN) :: ngrid,nlay 43 REAL, INTENT(IN) :: pz0(ngrid) 44 REAL, INTENT(IN) :: pg,pz(ngrid,nlay) 45 REAL, INTENT(IN) :: pu(ngrid,nlay),pv(ngrid,nlay) 46 REAL, INTENT(IN) :: pts(ngrid),ph(ngrid,nlay) 47 REAL, INTENT(OUT) :: pcdv(ngrid),pcdh(ngrid) ! momentum and heat drag coefficient 44 48 45 49 c Local: … … 48 52 INTEGER ig 49 53 50 REAL zu2,z1,zri,zcd0,zz 51 52 REAL karman,b,c,d,c2b,c3bc,c3b,umin2 54 REAL karman 53 55 LOGICAL firstcal 54 DATA karman ,b,c,d,umin2/.4,5.,5.,5.,1.e-12/56 DATA karman/.41/ 55 57 DATA firstcal/.true./ 56 SAVE b,c,d,karman,c2b,c3bc,c3b,firstcal,umin2 57 58 SAVE karman 59 60 c Local(2): 61 c --------- 62 63 REAL rib(ngrid) ! Bulk Richardson number 64 REAL fm(ngrid) ! stability function for momentum 65 REAL fh(ngrid) ! stability function for heat 66 REAL z1z0,z1z0t ! ratios z1/z0 and z1/z0T 67 68 c phim = 1+betam*zeta or (1-bm*zeta)**am 69 c phih = alphah + betah*zeta or alphah(1.-bh*zeta)**ah 70 REAL betam, betah, alphah, bm, bh, lambda 71 c ah and am are assumed to be -0.25 and -0.5 respectively 72 73 REAL cdn(ngrid),chn(ngrid) ! neutral momentum and heat drag coefficient 74 REAL pz0t ! initial thermal roughness length. (local) 75 REAL ric ! critical richardson number 76 REAL prandtl(ngrid) ! prandtl number for UBL 77 REAL reynolds(ngrid) ! reynolds number for UBL 78 REAL pz0tcomp(ngrid) ! computed z0t 79 REAL ite 80 REAL residual 58 81 c----------------------------------------------------------------------- 59 82 c couche de surface: 60 83 c ------------------ 61 84 62 c DO ig=1,ngrid 63 c zu2=pu(ig)*pu(ig)+pv(ig)*pv(ig)+umin2 64 c pcdv(ig)=pz0(ig)*(1.+sqrt(zu2)) 65 c pcdh(ig)=pcdv(ig) 66 c ENDDO 67 c RETURN 68 c 69 c IF (firstcal) THEN 70 c c2b=2.*b 71 c c3bc=3.*b*c 72 c c3b=3.*b 73 c firstcal=.false. 74 c ENDIF 75 c 76 c!!!! WARNING, verifier la formule originale de Louis! 77 c DO ig=1,ngrid 78 c zu2=pu(ig)*pu(ig)+pv(ig)*pv(ig)+umin2 79 c zri=pg*pz(ig)*(ph(ig)-pts(ig))/(ph(ig)*zu2) 80 c._. 81 c zri=0.E+0 82 c._. 83 c z1=1.+pz(ig)/pz0(ig) 84 c zcd0=karman/log(z1) 85 c._. zcd0=zcd0*zcd0*sqrt(zu2) 86 c zcd0=zcd0*zcd0 87 c IF(zri.LT.0.) THEN 88 c._. z1=b*zri/(1.+c3bc*zcd0*sqrt(-z1*zri)) 89 c z1=b*zri/(1.+c3bc*zcd0*sqrt(-z1*zri*zu2)) 90 c pcdv(ig)=zcd0*(1.-2.*z1) 91 c pcdh(ig)=zcd0*(1.-3.*z1) 92 c ELSE 93 c zz=sqrt(1.+d*zri) 94 c pcdv(ig)=zcd0/(1.+c2b*zri/zz) 95 c pcdh(ig)=zcd0/(1.+c3b*zri*zz) 96 c ENDIF 97 c ENDDO 98 99 100 c On calcule un VRAI cdrag tout bete 85 reynolds(:)=0. 86 87 c Original formulation : 88 89 ! DO ig=1,ngrid 90 ! z1=1.E+0 + pz(ig,1)/pz0(ig) 91 ! zcd0=karman/log(z1) 92 ! zcd0=zcd0*zcd0 93 ! pcdv(ig)=zcd0 94 ! pcdh(ig)=zcd0 95 ! ENDDO 96 97 ! print*,'old : cd,ch; ',pcdv,pcdh 98 99 c New formulation (AC) : 100 101 c phim = 1+betam*zeta or (1-bm*zeta)**am 102 c phih = alphah + betah*zeta or alphah(1.-bh*zeta)**ah 103 c am=-0.5, ah=-0.25 104 105 pz0t = 0. ! for the sake of simplicity 106 pz0tcomp(:) = 0.1*pz0(:) 107 rib(:)=0. 108 pcdv(:)=0. 109 pcdh(:)=0. 110 111 c this formulation assumes alphah=1., implying betah=betam 112 c We use Dyer et al. parameters, as they cover a broad range of Richardson numbers : 113 bm=16. !UBL 114 bh=16. !UBL 115 alphah=1. 116 betam=5. !SBL 117 betah=5. !SBL 118 lambda=(sqrt(bh/bm))/alphah 119 ric=betah/(betam**2) 101 120 102 121 DO ig=1,ngrid 103 z1=1.E+0 + pz(ig,1)/pz0(ig) 104 zcd0=karman/log(z1) 105 zcd0=zcd0*zcd0 106 pcdv(ig)=zcd0 107 pcdh(ig)=zcd0 122 123 ite=0. 124 residual=abs(pz0tcomp(ig)-pz0t) 125 126 do while((residual .gt. 0.01*pz0(ig)) .and. (ite .lt. 10.)) 127 128 pz0t=pz0tcomp(ig) 129 130 if ((pu(ig,1) .ne. 0.) .or. (pv(ig,1) .ne. 0.)) then 131 132 c Classical Richardson number formulation 133 134 c rib(ig) = (pg/ph(ig,1))*pz(ig,1)*(ph(ig,1)-pts(ig)) 135 c & /(pu(ig,1)*pu(ig,1) + pv(ig,1)*pv(ig,1)) 136 137 c Richardson number formulation proposed by D.E. England et al. (1995) 138 139 rib(ig) = (pg/ph(ig,1)) 140 & *pz(ig,1)*pz0(ig)/sqrt(pz(ig,1)*pz0t) 141 & *(((log(pz(ig,1)/pz0(ig)))**2)/(log(pz(ig,1)/pz0t))) 142 & *(ph(ig,1)-pts(ig)) 143 & /(pu(ig,1)*pu(ig,1) + pv(ig,1)*pv(ig,1)) 144 145 else 146 print*,'warning, infinite Richardson at surface' 147 print*,pu(ig,1),pv(ig,1) 148 rib(ig) = ric ! traiter ce cas ! 149 endif 150 151 z1z0=pz(ig,1)/pz0(ig) 152 z1z0t=pz(ig,1)/pz0t 153 154 cdn(ig)=karman/log(z1z0) 155 cdn(ig)=cdn(ig)*cdn(ig) 156 chn(ig)=cdn(ig)*log(z1z0)/log(z1z0t) 157 158 c Stable case : 159 if (rib(ig) .gt. 0.) then 160 c From D.E. England et al. (95) 161 prandtl(ig)=1. 162 if(rib(ig) .lt. ric) then 163 c Assuming alphah=1. and bh=bm for stable conditions : 164 fm(ig)=((ric-rib(ig))/ric)**2 165 fh(ig)=((ric-rib(ig))/ric)**2 166 else 167 fm(ig)=0. 168 fh(ig)=0. 169 endif 170 c Unstable case : 171 else 172 c From D.E. England et al. (95) 173 fm(ig)=sqrt(1.-lambda*bm*rib(ig)) 174 fh(ig)=(1./alphah)*((1.-lambda*bh*rib(ig))**0.5)* 175 & (1.-lambda*bm*rib(ig))**0.25 176 prandtl(ig)=alphah*((1.-lambda*bm*rib(ig))**0.25)/ 177 & ((1.-lambda*bh*rib(ig))**0.5) 178 endif 179 180 reynolds(ig)=sqrt(fm(ig))*sqrt(pu(ig,1)**2 + pv(ig,1)**2)*pz0(ig) 181 & /(log(z1z0)) 182 pz0tcomp(ig)=pz0(ig)*exp(-karman*7.3* 183 & (reynolds(ig)**0.25)*(prandtl(ig)**0.5)) 184 185 186 residual = abs(pz0t-pz0tcomp(ig)) 187 ite = ite+1 188 ! print*, "iteration nnumber, residual",ite,residual 189 190 enddo ! of while 191 192 pz0t=0. 193 194 c Drag computation : 195 196 pcdv(ig)=cdn(ig)*fm(ig) 197 pcdh(ig)=chn(ig)*fh(ig) 198 108 199 ENDDO 109 200 110 111 112 113 c----------------------------------------------------------------------- 201 ! print*,'new : cd,ch; ',pcdv,pcdh 202 203 ! Some useful diagnostics : 204 205 ! call WRITEDIAGFI(ngrid,'Ri', 206 ! & 'Richardson nb','no units', 207 ! & 2,rib) 208 ! call WRITEDIAGFI(ngrid,'Pr', 209 ! & 'Prandtl nb','no units', 210 ! & 0,prandtl) 211 ! call WRITEDIAGFI(ngrid,'Re', 212 ! & 'Reynolds nb','no units', 213 ! & 0,reynolds) 214 ! call WRITEDIAGFI(ngrid,'z0tcomp', 215 ! & 'computed z0t','m', 216 ! & 2,pz0tcomp) 114 217 115 218 RETURN -
trunk/LMDZ.MARS/libf/phymars/vdifc.F
r224 r256 5 5 $ pdufi,pdvfi,pdhfi,pdqfi,pfluxsrf, 6 6 $ pdudif,pdvdif,pdhdif,pdtsrf,pq2, 7 $ pdqdif,pdqsdif )7 $ pdqdif,pdqsdif,wmax,zcdv) 8 8 IMPLICIT NONE 9 9 … … 59 59 REAL pz0(ngridmx) ! surface roughness length (m) 60 60 61 c Argument added to account for subgrid gustiness : 62 63 REAL wmax(ngridmx) 64 61 65 c Traceurs : 62 66 integer nq … … 228 232 c --------------------- 229 233 230 CALL vdif_cd( 234 CALL vdif_cd(ngrid,nlay,pz0,g,pzlay,pu,pv,ptsrf,ph 231 235 & ,zcdv_true,zcdh_true) 232 DO ig=1,ngrid 236 237 DO ig=1,ngrid 238 233 239 zu2=pu(ig,1)*pu(ig,1)+pv(ig,1)*pv(ig,1) 240 & +(1.2*wmax(ig))**2 241 234 242 zcdv(ig)=zcdv_true(ig)*sqrt(zu2) 235 243 zcdh(ig)=zcdh_true(ig)*sqrt(zu2) 236 ENDDO 244 245 ENDDO 246 247 ! Some usefull diagnostics for the new surface layer parametrization : 248 249 ! call WRITEDIAGFI(ngridmx,'Cd', 250 ! & 'momentum drag','no units', 251 ! & 2,zcdv_true) 252 ! call WRITEDIAGFI(ngridmx,'Ch', 253 ! & 'heat drag','no units', 254 ! & 2,zcdh_true) 255 ! call WRITEDIAGFI(ngridmx,'u*', 256 ! & 'friction velocity','m/s', 257 ! & 2,zcdv) 258 ! call WRITEDIAGFI(ngridmx,'T*', 259 ! & 'friction temperature','m/s', 260 ! & 2,zcdh) 261 237 262 238 263 c ** schema de diffusion turbulente dans la couche limite
Note: See TracChangeset
for help on using the changeset viewer.
