Changeset 530
- Timestamp:
- Feb 15, 2012, 2:16:43 AM (13 years ago)
- Location:
- trunk/LMDZ.MARS
- Files:
-
- 5 edited
-
README (modified) (1 diff)
-
libf/phymars/growthrate.F (modified) (3 diffs)
-
libf/phymars/improvedclouds.F (modified) (8 diffs)
-
libf/phymars/newsedim.F (modified) (2 diffs)
-
libf/phymars/nuclea.F (modified) (7 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/LMDZ.MARS/README
r520 r530 1399 1399 >> Update in deftank: callphys.def.outliers,run.def.1d; added traceur.def.scavenging 1400 1400 1401 1401 == 15/02/12 == AS 1402 >> Optimization of watercycle new capabilities, e.g.: 1403 ... in nuclea, removed 'if' from function fshape, modified **, decrease number of overall variables 1404 ... in growthrate and newsedim, modified ** 1405 ... in improvedclouds, modified **, improved recursive calls to erf (now loop w/ 2 erf + 1 log instead of 4 erf + 4 log) 1406 >> Checked consistency with previous formulations 1407 >> Code 35% faster with a Valles Marineris mesoscale run in springtime (set imicro=25) 1408 ... probably even better at cloudier seasons -
trunk/LMDZ.MARS/libf/phymars/growthrate.F
r358 r530 58 58 k = (0.17913 * temp - 13.9789) * 4.184e-4 59 59 c - Latent heat of h2o (J.kg-1) 60 Lv = (2834.3 - 0.28 * (temp-To) - 0.004 * (temp-To)**2 ) * 1.e+3 60 Lv = (2834.3 61 & - 0.28 * (temp-To) 62 & - 0.004 * (temp-To) * (temp-To) ) * 1.e+3 61 63 62 64 c - Constant to compute gas mean free path 63 65 c l= (T/P)*a, with a = ( 0.707*8.31/(4*pi*molrad**2 * avogadro)) 64 afactor = 0.707*rgp/(4 * pi * molco2**2* nav)66 afactor = 0.707*rgp/(4 * pi * molco2 * molco2 * nav) 65 67 66 68 c - Compute Dv, water vapor diffusion coefficient … … 69 71 70 72 Dv = 1./3. * sqrt( 8*kbz*temp/(pi*mh2o/nav) )* kbz * temp / 71 & ( pi * pmid * (molco2+molh2o)**2 * sqrt(1.+mh2o/mco2) ) 73 & ( pi * pmid * (molco2+molh2o)*(molco2+molh2o) 74 & * sqrt(1.+mh2o/mco2) ) 72 75 73 76 knudsen = temp / pmid * afactor / rcrystal … … 76 79 77 80 c - Compute Rk 78 Rk = Lv**2 * rho_ice * mh2o / (k*rgp*temp* *2.)81 Rk = Lv**2 * rho_ice * mh2o / (k*rgp*temp*temp) 79 82 c - Compute Rd 80 83 Rd = rgp * temp *rho_ice / (Dv*psat*mh2o) -
trunk/LMDZ.MARS/libf/phymars/improvedclouds.F
r520 r530 25 25 c (October 2011) 26 26 c T. Navarro, debug & correction (October-December 2011) 27 c A. Spiga, optimization (February 2012) 27 28 c------------------------------------------------------------------ 28 29 #include "dimensions.h" … … 94 95 REAL*8 Mo,No 95 96 REAL*8 dN,dM,newvap 96 REAL*8 Rn, Rm, dev2 97 REAL*8 Rn, Rm, dev2, yeah, yeahn, yeahm 97 98 REAL*8 n_aer(nbin_cld) ! number conc. of particle/each size bin 98 99 REAL*8 m_aer(nbin_cld) ! mass mixing ratio of particle/each size bin … … 166 167 rb_cld(1) = rbmin_cld 167 168 rad_cld(1) = rmin_cld 168 vol_cld(1) = 4./3. * dble(pi) * rmin_cld* *3.169 vol_cld(1) = 4./3. * dble(pi) * rmin_cld*rmin_cld*rmin_cld 169 170 170 171 do i=1,nbin_cld-1 171 rad_cld(i+1) = rad_cld(i) * vrat_cld**(1./3.)172 rad_cld(i+1) = rad_cld(i) * CBRT(vrat_cld) !!**(1./3.) 172 173 vol_cld(i+1) = vol_cld(i) * vrat_cld 173 174 enddo 174 175 175 176 do i=1,nbin_cld 176 rb_cld(i+1)= ( (2.*vrat_cld) / (vrat_cld+1.) )**(1./3.) *177 rb_cld(i+1)= CBRT( (2.*vrat_cld) / (vrat_cld+1.) ) * 177 178 & rad_cld(i) 178 179 dr_cld(i) = rb_cld(i+1) - rb_cld(i) … … 285 286 Mo = zq(ig,l,igcm_dust_mass)* tauscaling(ig) 286 287 No = zq(ig,l,igcm_dust_number)* tauscaling(ig)+ 1.e-30 288 dev2 = 1. / ( sqrt(2.) * sigma_ice ) 287 289 Rn = rdust(ig,l) 288 Rm = Rn * exp( 3. * sigma_ice**2. ) 289 Rn = 1. / Rn 290 Rm = 1. / Rm 291 dev2 = 1. / ( sqrt(2.) * sigma_ice ) 290 Rn = -dlog(Rn) 291 Rm = Rn - 3. * sigma_ice*sigma_ice 292 yeah = dlog(rb_cld(1)) 293 yeahn = derf( (yeah+Rn) *dev2) 294 yeahm = derf( (yeah+Rm) *dev2) 292 295 do i = 1, nbin_cld 293 n_aer(i) = 0.5 * No * ( derf( dlog(rb_cld(i+1)*Rn) * dev2 ) 294 & -derf( dlog(rb_cld(i) * Rn) * dev2 ) ) 295 m_aer(i) = 0.5 * Mo * ( derf( dlog(rb_cld(i+1)*Rm) * dev2 ) 296 & -derf( dlog(rb_cld(i) * Rm) * dev2 ) ) 296 n_aer(i) = -0.5 * No * yeahn !! this ith previously computed 297 m_aer(i) = -0.5 * Mo * yeahm !! this ith previously computed 298 yeah = dlog(rb_cld(i+1)) !! this (i+1)th now computed 299 yeahn = derf( (yeah+Rn) *dev2) 300 yeahm = derf( (yeah+Rm) *dev2) 301 n_aer(i) = n_aer(i) + 0.5 * No * yeahn 302 m_aer(i) = m_aer(i) + 0.5 * Mo * yeahm 297 303 enddo 298 304 !!! MORE EFFICIENT COMPUTATIONALLY THAN 305 ! Rn = rdust(ig,l) 306 ! Rm = Rn * exp( 3. * sigma_ice**2. ) 307 ! Rn = 1. / Rn 308 ! Rm = 1. / Rm 309 ! do i = 1, nbin_cld 310 ! n_aer(i) = 0.5 * No * ( derf( dlog(rb_cld(i+1)*Rn) * dev2 ) 311 ! & -derf( dlog(rb_cld(i) * Rn) * dev2 ) ) 312 ! m_aer(i) = 0.5 * Mo * ( derf( dlog(rb_cld(i+1)*Rm) * dev2 ) 313 ! & -derf( dlog(rb_cld(i) * Rm) * dev2 ) ) 314 ! enddo 315 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 316 299 317 ! sumcheck = 0 300 318 ! do i = 1, nbin_cld … … 346 364 rhocloud(ig,l) = min(max(rhocloud(ig,l),rho_ice),rho_dust) 347 365 rice(ig,l) = 348 & ( Mo / No * 0.75 / pi / rhocloud(ig,l) )**(1./3.)366 & CBRT( Mo / No * 0.75 / pi / rhocloud(ig,l) ) !**(1./3.) 349 367 c nuice(ig,l)=nuice_ref ! used for rad. transfer calculations 350 368 if ((Mo.lt.1.e-20) .or. (No.le.1)) rice(ig,l) = 1.e-8 … … 466 484 rice_out(ig,l)=rice(ig,l) 467 485 rice(ig,l) = 468 & ( Mo / No * 0.75 / pi / rhocloud(ig,l) )**(1./3.)486 & CBRT( Mo / No * 0.75 / pi / rhocloud(ig,l) ) !**(1./3.) 469 487 if ((Mo.lt.1.e-20) .or. (No.le.1)) rice(ig,l) = 1.e-8 470 488 rice_out(ig,l)=rice(ig,l)-rice_out(ig,l) … … 478 496 rdust(ig,l)=min(max(rdust(ig,l),1.e-10),500.e-6) 479 497 480 rsedcloud(ig,l)=max( rice(ig,l)*(1.+nuice_sed)**3., 498 rsedcloud(ig,l)=max( rice(ig,l)*(1.+nuice_sed)* 499 & (1.+nuice_sed)*(1.+nuice_sed), 481 500 & rdust(ig,l) ) 482 501 rsedcloud(ig,l)=min(rsedcloud(ig,l),1.e-4) … … 540 559 541 560 drice_col(ig) = drice_col(ig)/icetot(ig) 542 543 561 544 562 IF (ngrid.ne.1) THEN ! 3D -
trunk/LMDZ.MARS/libf/phymars/newsedim.F
r358 r530 91 91 c - Constant to compute gas mean free path 92 92 c l= (T/P)*a, avec a = ( 0.707*8.31/(4*pi*molrad**2 * avogadro)) 93 a = 0.707*8.31/(4*3.1416* molrad* *2* 6.023e23)93 a = 0.707*8.31/(4*3.1416* molrad*molrad * 6.023e23) 94 94 95 95 c - Correction to account for non-spherical shape (Murphy et al. 1990) … … 124 124 endif 125 125 c vstokes 126 vstokes(ig,l) = b * rhofall * rfall* *2*126 vstokes(ig,l) = b * rhofall * rfall*rfall * 127 127 & (1 + 1.333* ( a*pt(ig,l)/pplev(ig,l) )/rfall) 128 128 -
trunk/LMDZ.MARS/libf/phymars/nuclea.F
r520 r530 11 11 * Adapted for the LMD/GCM by J.-B. Madeleine * 12 12 * (October 2011) * 13 * Optimisation by A. Spiga (February 2012) * 13 14 ******************************************************* 14 15 … … 34 35 DOUBLE PRECISION gstar ! # of molecules forming a critical embryo 35 36 DOUBLE PRECISION fistar ! Activation energy required to form a critical embryo (J) 36 DOUBLE PRECISION zeldov ! Zeldovitch factor (no dim)37 ! DOUBLE PRECISION zeldov ! Zeldovitch factor (no dim) 37 38 DOUBLE PRECISION fshape ! function defined at the end of the file 38 39 DOUBLE PRECISION deltaf 39 40 40 41 c Ratio rstar/radius of the nucleating dust particle 41 42 c double precision xratio 42 43 43 44 double precision mtetalocal ! local mteta in double precision 45 46 double precision fshapesimple,zefshape 47 44 48 45 49 integer i … … 51 55 c ************************************************* 52 56 53 mtetalocal = mteta 57 mtetalocal = mteta !! use mtetalocal for better performance 54 58 55 59 cccccccccccccccccccccccccccccccccccccccccccccccccc … … 86 90 nh2o = ph2o / kbz / temp 87 91 rstar = 2. * sig(temp) * vo1 / (rgp*temp*dlog(sat)) 88 gstar = 4. * nav * pi * (rstar**3) / (3.*vo1) 92 gstar = 4. * nav * pi * (rstar * rstar * rstar) / (3.*vo1) 93 94 fshapesimple = (2.+mtetalocal)*(1.-mtetalocal)*(1.-mtetalocal) 95 & / 4. 89 96 90 97 c Loop over size bins … … 97 104 endif 98 105 99 xratio = rad_cld(i) / rstar 100 fistar = (4./3.*pi) * sig(temp) * (rstar**2.) * 101 & fshape(mtetalocal,xratio) 106 if (rad_cld(i).gt.3000.*rstar) then 107 zefshape = fshapesimple 108 else 109 zefshape = fshape(mtetalocal,rad_cld(i)/rstar) 110 endif 111 112 fistar = (4./3.*pi) * sig(temp) * (rstar * rstar) * 113 & zefshape 102 114 deltaf = (2.*desorp-surfdif-fistar)/ 103 115 & (kbz*temp) … … 107 119 nucrate(i) = 0. 108 120 else 109 zeldov= sqrt ( fistar /110 & (3.*pi*kbz*temp*(gstar* *2.)) )111 nucrate(i)= zeldov* kbz * temp * rstar121 nucrate(i)= sqrt ( fistar / 122 & (3.*pi*kbz*temp*(gstar*gstar)) ) 123 & * kbz * temp * rstar 112 124 & * rstar * 4. * pi 113 & * ( nh2o*rad_cld(i) )**2. 114 & / ( fshape(mtetalocal,xratio) * nus * m0 ) 125 & * ( nh2o*rad_cld(i) ) 126 & * ( nh2o*rad_cld(i) ) 127 & / ( zefshape * nus * m0 ) 115 128 & * dexp (deltaf) 116 129 endif … … 137 150 138 151 double precision cost,rap 139 double precision phi 140 double precision a,b,c 152 double precision yeah 141 153 142 phi = sqrt( 1. - 2.*cost*rap + rap**2 ) 143 a = 1. + ( (1.-cost*rap)/phi )**3 144 b = (rap**3) * (2.-3.*(rap-cost)/phi+((rap-cost)/phi)**3) 145 c = 3. * cost * (rap**2) * ((rap-cost)/phi-1.) 154 !! PHI 155 yeah = sqrt( 1. - 2.*cost*rap + rap*rap ) 156 !! FSHAPE = TERM A 157 fshape = (1.-cost*rap) / yeah 158 fshape = fshape * fshape * fshape 159 fshape = 1. + fshape 160 !! ... + TERM B 161 yeah = (rap-cost)/yeah 162 fshape = fshape + rap*rap*rap*(2.-3.*yeah+yeah*yeah*yeah) 163 !! ... + TERM C 164 fshape = fshape + 3. * cost * rap * rap * (yeah-1.) 165 !! FACTOR 1/2 166 fshape = 0.5*fshape 146 167 147 fshape = 0.5*(a+b+c) 148 149 if (rap.gt.3000.) fshape = ((2.+cost)*(1.-cost)**2)/4. 150 151 return 168 return 152 169 end
Note: See TracChangeset
for help on using the changeset viewer.
