Changeset 1720 for trunk/LMDZ.MARS/libf/phymars
- Timestamp:
- Jul 6, 2017, 5:12:58 PM (8 years ago)
- Location:
- trunk/LMDZ.MARS/libf/phymars
- Files:
-
- 11 edited
-
callkeys.h (modified) (2 diffs)
-
callsedim.F (modified) (2 diffs)
-
co2cloud.F (modified) (32 diffs)
-
conf_phys.F (modified) (4 diffs)
-
improvedCO2clouds.F (modified) (30 diffs)
-
initracer.F (modified) (2 diffs)
-
massflowrateCO2.F (modified) (2 diffs)
-
microphys.h (modified) (1 diff)
-
nucleaCO2.F (modified) (1 diff)
-
physiq_mod.F (modified) (11 diffs)
-
tracer_mod.F90 (modified) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
trunk/LMDZ.MARS/libf/phymars/callkeys.h
r1711 r1720 14 14 & ,activice,water,tifeedback,microphys,supersat,caps,photochem & 15 15 & ,calltherm,callrichsl,callslope,tituscap,callyamada4,co2clouds,& 16 & microphysco2,CLFvarying16 & co2useh2o,meteo_flux,CLFvaryingCO2,spantCO2,CLFvarying 17 17 18 18 COMMON/callkeys_i/iradia,iaervar,iddist,ilwd,ilwb,ilwn,ncouche & … … 60 60 logical sedimentation 61 61 logical activice,tifeedback,supersat,caps 62 logical co2clouds,co2useh2o,meteo_flux,CLFvaryingCO2 63 integer spantCO2 62 64 logical CLFvarying 63 logical co2clouds64 65 logical water 65 66 logical microphys 66 67 logical photochem 67 logical microphysco268 68 integer nltemodel 69 69 integer nircorr -
trunk/LMDZ.MARS/libf/phymars/callsedim.F
r1617 r1720 205 205 ENDIF !of if (microphys) 206 206 207 IF ( microphysco2) THEN207 IF (co2clouds) THEN 208 208 iccnco2_mass=0 209 209 iccnco2_number=0 … … 227 227 write(*,*) 'callsedim: error! could not identify' 228 228 write(*,*) ' tracers for ccn co2 mass and number mixing' 229 write(*,*) ' ratio and microphysco2 isactivated!'229 write(*,*) ' ratio and co2clouds are activated!' 230 230 stop 231 231 endif 232 ENDIF !of if ( microphysco2)232 ENDIF !of if (co2clouds) 233 233 234 234 -
trunk/LMDZ.MARS/libf/phymars/co2cloud.F
r1685 r1720 4 4 & nq,tau,tauscaling,rdust,rice,riceco2,nuice, 5 5 & rsedcloudco2,rhocloudco2, 6 & rsedcloud,rhocloud, zlev,pdqs_sedco2)6 & rsedcloud,rhocloud,pzlev,pdqs_sedco2) 7 7 ! to use 'getin' 8 8 use dimradmars_mod, only: naerkind … … 16 16 & igcm_ccnco2_mass, igcm_ccnco2_number, 17 17 & rho_dust, nuiceco2_sed, nuiceco2_ref, 18 & rho_ice_co2,r3n_q,rho_ice,nuice_sed, 19 & meteo_flux_mass,meteo_flux_number, 20 & meteo_alt 18 & rho_ice_co2,r3n_q,rho_ice,nuice_sed 19 21 20 22 21 IMPLICIT NONE 22 23 24 #include "datafile.h" 25 #include "callkeys.h" 26 #include "microphys.h" 23 27 24 28 … … 38 42 c of Montmessin, Navarro & al. 39 43 c 40 c 44 c 07/2017 J.Audouard 45 c Several logicals can be set to .true. in callphys.def 46 c co2clouds=.true. call this routine 47 c co2useh2o=.true. allow the use of water ice particles as CCN for CO2 48 c meteo_flux=.true. supply meteoritic particles 49 c CLFvaryingCO2=.true. allows a subgrid temperature distribution 50 c of amplitude spantCO2(=integer in callphys.def) 51 c 41 52 c======================================================================= 42 53 … … 44 55 c declarations: 45 56 c ------------- 46 47 !#include "dimensions.h"48 !#include "dimphys.h"49 #include "callkeys.h"50 !#include "tracer.h"51 !#include "comgeomfi.h"52 !#include "dimradmars.h"53 ! naerkind is set in scatterers.h (built when compiling with makegcm -s #)54 !#include"scatterers.h"55 #include "microphys.h"56 57 57 58 58 c Inputs: … … 68 68 REAL pt(ngrid,nlay) ! temperature at the middle of the layers (K) 69 69 REAL pdt(ngrid,nlay) ! tendence temperature des autres param. 70 real,intent(in) :: zlev(ngrid,nlay+1) ! altitude at the boundaries of the layers70 real,intent(in) :: pzlev(ngrid,nlay+1) ! altitude at the boundaries of the layers 71 71 72 72 real pq(ngrid,nlay,nq) ! traceur (kg/kg) … … 101 101 real rhocloud(ngrid,nlay) ! Cloud density (kg.m-3) 102 102 ! for ice radius computation 103 REAL Mo,No103 104 104 REAl ccntyp 105 105 … … 123 123 REAL zqsatco2(ngrid,nlay) ! saturation co2 124 124 125 INTEGER iq,ig,l 125 INTEGER iq,ig,l,i 126 126 LOGICAL,SAVE :: firstcall=.true. 127 DOUBLE PRECISION Nccnco2, Niceco2, mdustJA,ndustJA127 DOUBLE PRECISION Nccnco2, Niceco2,Nco2,mdustJA,ndustJA 128 128 DOUBLE PRECISION Qccnco2 129 129 real :: beta … … 132 132 real masse (ngrid,nlay) ! Layer mass (kg.m-2) 133 133 double precision diff,diff0 134 integer meteo_lvl135 134 real tempo_traceur_t(ngrid,nlay) 136 135 real tempo_traceurs(ngrid,nlay,nq) 137 136 real sav_trac(ngrid,nlay,nq) 138 137 real pdqsed(ngrid,nlay,nq) 139 140 DOUBLE PRECISION,allocatable,save :: memdMMccn(:,:) 138 REAL lw !Latent heat of sublimation (J.kg-1) 139 REAL,save :: l0,l1,l2,l3,l4 !Coeffs for lw 140 141 DOUBLE PRECISION,allocatable,save :: memdMMccn(:,:) ! Nb particules intégré 141 142 DOUBLE PRECISION,allocatable,save :: memdMMh2o(:,:) 142 143 DOUBLE PRECISION,allocatable,save :: memdNNccn(:,:) 144 DOUBLE PRECISION :: myT 145 ! What we need for Qext reading and tau computation 146 DOUBLE PRECISION vrat_cld ! Volume ratio 147 DOUBLE PRECISION rb_cldco2(nbinco2_cld+1) ! boundary values of each rad_cldco2 bin (m) 148 SAVE rb_cldco2 149 DOUBLE PRECISION, PARAMETER :: rmin_cld = 1.e-11 ! Minimum radius (m) 150 DOUBLE PRECISION, PARAMETER :: rmax_cld = 3.e-6 ! Maximum radius (m) 151 DOUBLE PRECISION, PARAMETER :: rbmin_cld =1.e-12 152 ! Minimum boundary radius (m) 153 DOUBLE PRECISION, PARAMETER :: rbmax_cld = 5.e-6 ! Maximum boundary radius (m) 154 155 DOUBLE PRECISION dr_cld(nbinco2_cld) ! width of each rad_cldco2 bin (m) 156 DOUBLE PRECISION vol_cld(nbinco2_cld) ! particle volume for each bin (m3) 157 REAL sigma_iceco2 ! Variance of the ice and CCN distributions 158 159 logical :: file_ok 160 double precision :: radv(10000),Qextv1mic(10000) 161 double precision :: Qext1bins(100),Qtemp 162 double precision :: ltemp1(10000),ltemp2(10000) 163 integer :: nelem,lebon1,lebon2,uQext 164 save Qext1bins 165 character(len=100) scanline 166 DOUBLE PRECISION n_aer(nbinco2_cld),Rn,No,n_derf,dev2 167 DOUBLE PRECISION Qext1bins2(ngrid,nlay) 168 DOUBLE PRECISION tau1mic(ngrid) !co2 ice column opacity at 1µm 169 170 ! For sub grid T distribution 171 172 REAL zt(ngrid,nlay) ! local value of temperature 173 REAL :: zq(ngrid, nlay,nq) 174 175 REAL :: tcond(ngrid,nlay) 176 REAL :: zqvap(ngrid,nlay) 177 REAL :: zqice(ngrid,nlay) 178 REAL :: spant,zdelt ! delta T for the temperature distribution 179 REAL :: zteff(ngrid, nlay)! effective temperature in the cloud,neb 180 REAL :: pqeff(ngrid, nlay, nq)! effective tracers quantities in the cloud 181 REAL :: cloudfrac(ngrid,nlay) ! cloud fraction 182 REAL :: mincloud ! min cloud frac 183 c logical :: CLFvaryingCO2 143 184 c ** un petit test de coherence 144 185 c -------------------------- 145 186 146 187 IF (firstcall) THEN 147 148 188 if (nq.gt.nqmx) then 149 189 write(*,*) 'stop in co2cloud (nq.gt.nqmx)!' … … 152 192 endif 153 193 write(*,*) "co2cloud.F: rho_ice_co2 = ",rho_ice_co2 154 155 194 write(*,*) "co2cloud: igcm_co2=",igcm_co2 156 195 write(*,*) " igcm_co2_ice=",igcm_co2_ice … … 169 208 write(*,*)"CO2 Microphysics timestep is",microtimestep 170 209 171 172 if (meteo_flux_number .ne. 0) then 173 174 write(*,*) "Warning ! Meteoritic flux of dust is turned on" 175 write(*,*) "Dust mass = ",meteo_flux_mass 176 write(*,*) "Dust number = ",meteo_flux_number 177 write(*,*) "Are added at the z-level (km)",meteo_alt 178 write(*,*) "Every timestep in co2cloud.F" 179 endif 210 ! Values for latent heat computation 211 l0=595594d0 212 l1=903.111d0 213 l2=-11.5959d0 214 l3=0.0528288d0 215 l4=-0.000103183d0 216 c$$$ 217 c$$$ if (meteo_flux_number .ne. 0) then 218 c$$$ 219 c$$$ write(*,*) "Warning ! Meteoritic flux of dust is turned on" 220 c$$$ write(*,*) "Dust mass = ",meteo_flux_mass 221 c$$$ write(*,*) "Dust number = ",meteo_flux_number 222 c$$$ write(*,*) "Are added at the z-level (km)",meteo_alt 223 c$$$ write(*,*) "Every timestep in co2cloud.F" 224 c$$$ endif 225 180 226 if (.not. allocated(memdMMccn)) allocate(memdMMccn(ngrid,nlay)) 181 227 if (.not. allocated(memdNNccn)) allocate(memdNNccn(ngrid,nlay)) … … 185 231 memdMMh2o(:,:)=0. 186 232 memdNNccn(:,:)=0. 233 234 c Compute the size bins of the distribution of CO2 ice particles 235 c --> used for opacity calculations 236 237 c rad_cldco2 is the primary radius grid used for microphysics computation. 238 c The grid spacing is computed assuming a constant volume ratio 239 c between two consecutive bins; i.e. vrat_cld. 240 c vrat_cld is determined from the boundary values of the size grid: 241 c rmin_cld and rmax_cld. 242 c The rb_cldco2 array contains the boundary values of each rad_cldco2 bin. 243 c dr_cld is the width of each rad_cldco2 bin. 244 sigma_iceco2 = sqrt(log(1.+nuiceco2_sed)) 245 c Volume ratio between two adjacent bins 246 ! vrat_cld 247 vrat_cld = log(rmax_cld/rmin_cld) / float(nbinco2_cld-1) *3. 248 vrat_cld = exp(vrat_cld) 249 rb_cldco2(1) = rbmin_cld 250 rad_cldco2(1) = rmin_cld 251 vol_cld(1) = 4./3. * dble(pi) * rmin_cld*rmin_cld*rmin_cld 252 do i=1,nbinco2_cld-1 253 rad_cldco2(i+1) = rad_cldco2(i) * vrat_cld**(1./3.) 254 vol_cld(i+1) = vol_cld(i) * vrat_cld 255 enddo 256 do i=1,nbinco2_cld 257 rb_cldco2(i+1)= ( (2.*vrat_cld) / (vrat_cld+1.) )**(1./3.) * 258 & rad_cldco2(i) 259 dr_cld(i) = rb_cldco2(i+1) - rb_cldco2(i) 260 enddo 261 rb_cldco2(nbinco2_cld+1) = rbmax_cld 262 dr_cld(nbinco2_cld) = rb_cldco2(nbinco2_cld+1) - 263 & rb_cldco2(nbinco2_cld) 264 265 c read the Qext values 266 INQUIRE(FILE=datafile(1:LEN_TRIM(datafile))// 267 & '/optprop_co2ice_1mic.dat', EXIST=file_ok) 268 IF (.not. file_ok) THEN 269 write(*,*) 'file optprop_co2ice_1mic.dat should be in ' 270 & ,datafile 271 STOP 272 endif 273 open(newunit=uQext,file=trim(datafile)// 274 & '/optprop_co2ice_1mic.dat' 275 & ,FORM='formatted') 276 read(uQext,*) !skip 1 line 277 do i=1,10000 278 read(uQext,'(E11.5)') radv(i) 279 enddo 280 read(uQext,*) !skip 1 line 281 do i=1,10000 282 read(uQext,'(E11.5)') Qextv1mic(i) 283 enddo 284 close(uQext) 285 c innterpol the Qext values 286 !rice_out=rad_cldco2 287 do i=1,nbinco2_cld 288 ltemp1=abs(radv(:)-rb_cldco2(i)) 289 ltemp2=abs(radv(:)-rb_cldco2(i+1)) 290 lebon1=minloc(ltemp1,DIM=1) 291 lebon2=minloc(ltemp2,DIM=1) 292 nelem=lebon2-lebon1+1. 293 Qtemp=0d0 294 do l=0,nelem 295 Qtemp=Qtemp+Qextv1mic(lebon1+l) !mean value in the interval 296 enddo 297 Qtemp=Qtemp/nelem 298 Qext1bins(i)=Qtemp 299 enddo 300 Qext1bins(:)=Qext1bins(:)*rad_cldco2(:)*rad_cldco2(:)*pi 301 ! The actuall tau computation and output is performed in co2cloud.F 302 303 print*,'--------------------------------------------' 304 print*,'Microphysics co2: size bin-Qext information:' 305 print*,' i, rad_cldco2(i), Qext1bins(i)' 306 do i=1,nbinco2_cld 307 write(*,'(i3,3x,3(e12.6,4x))') i, rad_cldco2(i), 308 & Qext1bins(i) 309 enddo 310 print*,'--------------------------------------------' 311 312 313 do i=1,nbinco2_cld+1 314 rb_cldco2(i) = log(rb_cldco2(i)) !! we save that so that it is not computed at each timestep and gridpoint 315 enddo 316 if (CLFvaryingCO2) then 317 write(*,*) 318 write(*,*) "CLFvaryingCO2 is set to true is callphys.def" 319 write(*,*) "The temperature field is enlarged to +/-",spantCO2 320 write(*,*) "for the CO2 microphysics " 321 endif 187 322 firstcall=.false. 188 323 ENDIF ! of IF (firstcall) 189 324 190 325 c-----Initialization 191 326 dev2 = 1. / ( sqrt(2.) * sigma_iceco2 ) 192 327 beta=0.85 193 328 subpdq(1:ngrid,1:nlay,1:nq) = 0 … … 195 330 subpdqcloudco2(1:ngrid,1:nlay,1:nq) = 0 196 331 subpdtcloudco2(1:ngrid,1:nlay) = 0 197 198 c add meteoritic flux of dust 199 !convert meteo_alt (in km) to z-level 200 !pzlay altitudes of the layers 201 if (meteo_flux_number .ne. 0) then 202 do ig=1, ngrid 203 diff0=1000 204 meteo_lvl=1 205 do l=1,nlay 206 diff=pzlay(ig,l)/1000-meteo_alt 207 if (abs(diff) .lt. diff0) then 208 diff0=abs(diff) 209 meteo_lvl=l 210 endif 211 enddo 212 c write(*,*) "meteo_lvl=",meteo_lvl 213 pdq(ig,meteo_lvl,igcm_dust_mass)= 214 & pdq(ig,meteo_lvl,igcm_dust_mass) 215 & +dble(meteo_flux_mass)/tauscaling(ig) 216 pdq(ig,meteo_lvl,igcm_dust_number)= 217 & pdq(ig,meteo_lvl,igcm_dust_number) 218 & +dble(meteo_flux_number)/tauscaling(ig) 219 enddo 220 endif 221 c call WRITEDIAGFI(ngrid,"pzlay","altitude","km",3, 222 c & pzlay) 223 332 c$$$ 333 c$$$ call WRITEDIAGFI(ngrid,"pzlay","pzlay","km",3, 334 c$$$ & pzlay) 335 c$$$ call WRITEDIAGFI(ngrid,"pplay","pplay","Pa",3, 336 c$$$ & pplay) 224 337 225 338 wq(:,:)=0 … … 230 343 masse(1:ngrid,1:nlay)=0 231 344 232 tempo_traceur_t(1:ngrid,1:nlay)=0233 tempo_traceurs(1:ngrid,1:nlay,1:nq)=0234 345 sav_trac(1:ngrid,1:nlay,1:nq)=0 235 346 pdqsed(1:ngrid,1:nlay,1:nq)=0 … … 238 349 do ig=1, ngrid 239 350 masse(ig,l)=(pplev(ig,l) - pplev(ig,l+1)) /g 240 epaisseur(ig,l)= zlev(ig,l+1) -zlev(ig,l)351 epaisseur(ig,l)= pzlev(ig,l+1) - pzlev(ig,l) 241 352 enddo 242 353 enddo 243 354 244 355 !CLFvaryingCO2=.true. 245 356 246 357 c------------------------------------------------------------------- 358 c 0. Representation of sub-grid water ice clouds 359 c------------------ 360 IF (CLFvaryingCO2) THEN 361 spant=spantCO2 ! delta T for the temprature distribution 362 mincloud=0.1 ! min cloudfrac when there is ice 363 ! IF (flagcloudco2) THEN 364 ! write(*,*) "CLFCO2 Delta T", spant 365 ! write(*,*) "CLFCO2 mincloud", mincloud 366 ! flagcloudco2=.false. 367 ! END IF 368 369 c-----Tendencies 370 DO l=1,nlay 371 DO ig=1,ngrid 372 zt(ig,l)=pt(ig,l)+ pdt(ig,l)*ptimestep 373 ENDDO 374 ENDDO 375 DO l=1,nlay 376 DO ig=1,ngrid 377 DO iq=1,nq 378 zq(ig,l,iq)=pq(ig,l,iq)+pdq(ig,l,iq)*ptimestep 379 ENDDO 380 ENDDO 381 ENDDO 382 zqvap=zq(:,:,igcm_co2) 383 zqice=zq(:,:,igcm_co2_ice) 384 CALL tcondco2(ngrid,nlay,pplay,zqvap,tcond) 385 ! A tester: CALL tcondco2(ngrid,nlay,pplay,zqvap,tcond) 386 zdelt=spant 387 DO l=1,nlay 388 DO ig=1,ngrid 389 IF (tcond(ig,l) .ge. (zt(ig,l)+zdelt)) THEN !Toute la fraction est saturée 390 zteff(ig,l)=zt(ig,l) 391 cloudfrac(ig,l)=1. 392 ELSE IF (tcond(ig,l) .le. (zt(ig,l)-zdelt)) THEN !Rien n'est saturé 393 zteff(ig,l)=zt(ig,l)-zdelt 394 cloudfrac(ig,l)=mincloud 395 ELSE 396 cloudfrac(ig,l)=(tcond(ig,l)-zt(ig,l)+zdelt)/ 397 & (2.0*zdelt) 398 zteff(ig,l)=(tcond(ig,l)+zt(ig,l)-zdelt)/2. !Ja Temperature moyenne de la fraction nuageuse 399 400 END IF 401 zteff(ig,l)=zteff(ig,l)-pdt(ig,l)*ptimestep 402 IF (cloudfrac(ig,l).le. mincloud) THEN 403 cloudfrac(ig,l)=mincloud 404 ELSE IF (cloudfrac(ig,l).ge. 1) THEN 405 cloudfrac(ig,l)=1. 406 END IF 407 ENDDO 408 ENDDO 409 ! Totalcloud frac of the column missing here 410 c----------------------- 411 c-----No sub-grid cloud representation (CLFvarying=false) 412 ELSE 413 DO l=1,nlay 414 DO ig=1,ngrid 415 zteff(ig,l)=pt(ig,l) 416 END DO 417 END DO 418 END IF ! end if (CLFvaryingco2)-------------------------------------------- 247 419 c 1. Tendencies: 248 420 c------------------ 249 421 250 422 c------ Effective tracers quantities in the cloud fraction 423 IF (CLFvaryingCO2) THEN 424 pqeff(:,:,:)=pq(:,:,:) ! prevent from buggs (A. Pottier) 425 426 c pqeff(:,:,igcm_ccn_mass) =pq(:,:,igcm_ccn_mass)/ 427 c & cloudfrac(:,:) 428 c pqeff(:,:,igcm_ccn_number)= 429 c & pq(:,:,igcm_ccn_number)/cloudfrac(:,:) 430 c pqeff(:,:,igcm_h2o_ice)= pq(:,:,igcm_h2o_ice)/ 431 c & cloudfrac(:,:) 432 pqeff(:,:,igcm_ccnco2_mass) =pq(:,:,igcm_ccnco2_mass)/ 433 & cloudfrac(:,:) 434 pqeff(:,:,igcm_ccnco2_number)= 435 & pq(:,:,igcm_ccnco2_number)/cloudfrac(:,:) 436 pqeff(:,:,igcm_co2_ice)= pq(:,:,igcm_co2_ice)/ 437 & cloudfrac(:,:) 438 439 ELSE 440 pqeff(:,:,:)=pq(:,:,:) 441 c pqeff(:,:,igcm_ccn_mass)= pq(:,:,igcm_ccn_mass) 442 c pqeff(:,:,igcm_ccn_number)= pq(:,:,igcm_ccn_number) 443 c pqeff(:,:,igcm_h2o_ice)= pq(:,:,igcm_h2o_ice) 444 c pqeff(:,:,igcm_ccnco2_mass)= pq(:,:,igcm_ccnco2_mass) 445 c pqeff(:,:,igcm_ccnco2_number)= pq(:,:,igcm_ccnco2_number) 446 c pqeff(:,:,igcm_co2_ice)= pq(:,:,igcm_co2_ice) 447 END IF 448 tempo_traceur_t(1:ngrid,1:nlay)=zteff(1:ngrid,1:nlay) 449 tempo_traceurs(1:ngrid,1:nlay,1:nq)=pqeff(1:ngrid,1:nlay,1:nq) 251 450 252 451 c------------------------------------------------------------------ … … 259 458 DO l=1,nlay 260 459 DO ig=1,ngrid 261 c tempo_traceur_t(ig,l)=tempo_traceur_t(ig,l)262 c & + subpdtcloudco2(ig,l)263 !write(*,*) 'T micro= ', tempo_traceur_t(ig,l)264 c tempo_traceurs(ig,l,:)=tempo_traceurs(ig,l,:)265 c & +subpdqcloudco2(ig,l,:)266 460 267 461 subpdt(ig,l) = subpdt(ig,l) 268 462 & + pdt(ig,l) ! At each micro timestep we add pdt in order to have a stepped entry 269 270 463 subpdq(ig,l,igcm_dust_mass) = 271 464 & subpdq(ig,l,igcm_dust_mass) 272 465 & + pdq(ig,l,igcm_dust_mass) 273 274 466 subpdq(ig,l,igcm_dust_number) = 275 467 & subpdq(ig,l,igcm_dust_number) 276 468 & + pdq(ig,l,igcm_dust_number) 277 278 469 subpdq(ig,l,igcm_ccnco2_mass) = 279 470 & subpdq(ig,l,igcm_ccnco2_mass) 280 471 & + pdq(ig,l,igcm_ccnco2_mass) 281 282 472 subpdq(ig,l,igcm_ccnco2_number) = 283 473 & subpdq(ig,l,igcm_ccnco2_number) 284 474 & + pdq(ig,l,igcm_ccnco2_number) 285 286 475 subpdq(ig,l,igcm_co2_ice) = 287 476 & subpdq(ig,l,igcm_co2_ice) 288 477 & + pdq(ig,l,igcm_co2_ice) 289 290 478 subpdq(ig,l,igcm_co2) = 291 479 & subpdq(ig,l,igcm_co2) 292 480 & + pdq(ig,l,igcm_co2) 293 294 481 subpdq(ig,l,igcm_h2o_ice) = 295 482 & subpdq(ig,l,igcm_h2o_ice) 296 483 & + pdq(ig,l,igcm_h2o_ice) 297 298 484 subpdq(ig,l,igcm_ccn_mass) = 299 485 & subpdq(ig,l,igcm_ccn_mass) 300 486 & + pdq(ig,l,igcm_ccn_mass) 301 302 487 subpdq(ig,l,igcm_ccn_number) = 303 488 & subpdq(ig,l,igcm_ccn_number) 304 489 & + pdq(ig,l,igcm_ccn_number) 305 306 tempo_traceur_t(ig,l)= pt(ig,l)+subpdt(ig,l)*microtimestep307 tempo_traceurs(ig,l,:)= pq(ig,l,:)+subpdq(ig,l,:)308 & *microtimestep309 !Stepped entry for sedimentation310 490 ENDDO 311 491 ENDDO 312 313 !RSEDCLOUD AND RICECO2 HERE 492 493 c add meteoritic flux of dust (old version) 494 cNew version (John Plane values) are added in improvedCO2clouds 495 !convert meteo_alt (in km) to z-level 496 !pzlay altitudes of the layers 497 c$$$!zlev altitudes (nlayl+1) of the boundaries 498 c$$$ if (meteo_flux_number .ge. 0) then 499 c$$$ do ig=1, ngrid 500 c$$$ l=1 501 c$$$ do while ( (((meteo_alt .ge. pplev(ig,l)) .and. 502 c$$$ & (meteo_alt .le. pplev(ig,l+1))) .eq. .false.) 503 c$$$ & .and. (l .lt. nlay) ) 504 c$$$ l=l+1 505 c$$$ enddo 506 c$$$ meteo_lvl=l 507 c$$$ subpdq(ig,meteo_lvl,igcm_dust_mass)= 508 c$$$ & subpdq(ig,meteo_lvl,igcm_dust_mass) 509 c$$$ & +meteo_flux_mass 510 c$$$ subpdq(ig,meteo_lvl,igcm_dust_number)= 511 c$$$ & subpdq(ig,meteo_lvl,igcm_dust_number) 512 c$$$ & +meteo_flux_number 513 c$$$ enddo 514 c$$$ endif 515 c------------------------------------------------------------------- 516 c 2. Main call to the cloud schemes: 517 c------------------------------------------------ 518 CALL improvedCO2clouds(ngrid,nlay,microtimestep, 519 & pplay,pzlev,zteff,subpdt, 520 & pqeff,subpdq,subpdqcloudco2,subpdtcloudco2, 521 & nq,tauscaling,memdMMccn,memdMMh2o,memdNNccn) 522 c------------------------------------------------------------------- 523 c 3. Updating tendencies after cloud scheme: 524 c----------------------------------------------- 525 DO l=1,nlay 526 DO ig=1,ngrid 527 subpdt(ig,l) = 528 & subpdt(ig,l) + subpdtcloudco2(ig,l) 529 subpdq(ig,l,igcm_dust_mass) = 530 & subpdq(ig,l,igcm_dust_mass) 531 & + subpdqcloudco2(ig,l,igcm_dust_mass) 532 subpdq(ig,l,igcm_dust_number) = 533 & subpdq(ig,l,igcm_dust_number) 534 & + subpdqcloudco2(ig,l,igcm_dust_number) 535 subpdq(ig,l,igcm_ccnco2_mass) = 536 & subpdq(ig,l,igcm_ccnco2_mass) 537 & + subpdqcloudco2(ig,l,igcm_ccnco2_mass) 538 subpdq(ig,l,igcm_ccnco2_number) = 539 & subpdq(ig,l,igcm_ccnco2_number) 540 & + subpdqcloudco2(ig,l,igcm_ccnco2_number) 541 subpdq(ig,l,igcm_co2_ice) = 542 & subpdq(ig,l,igcm_co2_ice) 543 & + subpdqcloudco2(ig,l,igcm_co2_ice) 544 subpdq(ig,l,igcm_co2) = 545 & subpdq(ig,l,igcm_co2) 546 & + subpdqcloudco2(ig,l,igcm_co2) 547 subpdq(ig,l,igcm_h2o_ice) = 548 & subpdq(ig,l,igcm_h2o_ice) 549 & + subpdqcloudco2(ig,l,igcm_h2o_ice) 550 subpdq(ig,l,igcm_ccn_mass) = 551 & subpdq(ig,l,igcm_ccn_mass) 552 & + subpdqcloudco2(ig,l,igcm_ccn_mass) 553 subpdq(ig,l,igcm_ccn_number) = 554 & subpdq(ig,l,igcm_ccn_number) 555 & + subpdqcloudco2(ig,l,igcm_ccn_number) 556 ENDDO 557 ENDDO 558 559 !Sedimentation 560 !Update traceurs, compute rsed 314 561 315 562 DO l=1, nlay 316 DO ig=1,ngrid 317 563 DO ig=1,ngrid 564 tempo_traceur_t(ig,l)=zteff(ig,l)+subpdt(ig,l) 565 & *microtimestep 566 tempo_traceurs(ig,l,:)=pqeff(ig,l,:) 567 & +subpdq(ig,l,:)*microtimestep 568 318 569 rho_ice_co2T(ig,l)=1000.*(1.72391-2.53e-4* 319 570 & tempo_traceur_t(ig,l)-2.87e-6* 320 571 & tempo_traceur_t(ig,l)*tempo_traceur_t(ig,l)) 321 322 572 573 rho_ice_co2=rho_ice_co2T(ig,l) 323 574 Niceco2=max(tempo_traceurs(ig,l,igcm_co2_ice),1.e-30) 324 575 Nccnco2=max(tempo_traceurs(ig,l,igcm_ccnco2_number), … … 334 585 rdust(ig,l)=(3./4./pi/2500.*Qccnco2/Nccnco2)**(1./3.) 335 586 rdust(ig,l)=max(rdust(ig,l),1.e-10) 336 !rdust(ig,l)=min(rdust(ig,l),5.e-6)587 c rdust(ig,l)=min(rdust(ig,l),5.e-6) 337 588 end if 338 589 rhocloudco2t(ig,l) = (Niceco2 *rho_ice_co2 339 590 & + Qccnco2*tauscaling(ig)*rho_dust) 340 591 & / (Niceco2 + Qccnco2*tauscaling(ig)) 341 c riceco2(ig,l)= Niceco2*3.0/ 342 c & (4.0*rho_ice_co2*pi*Nccnco2) 343 c & +rdust(ig,l)*rdust(ig,l)*rdust(ig,l) 344 c riceco2(ig,l)=riceco2(ig,l)**(1.0/3.0) 345 c write(*,*) "in co2clouds, rice = ",riceco2(ig,l) 346 c write(*,*) "in co2clouds, rho = ",rhocloudco2t(ig,l) 347 592 593 rhocloudco2t(ig,l)=min(max(rhocloudco2t(ig,l) 594 & ,rho_ice_co2),rho_dust) 348 595 riceco2(ig,l)=(Niceco2*3.0/ 349 596 & (4.0*rho_ice_co2*pi*Nccnco2 350 597 & *tauscaling(ig)) +rdust(ig,l)*rdust(ig,l) 351 598 & *rdust(ig,l))**(1.0/3.0) 352 353 599 riceco2(ig,l)=max(1.e-10,riceco2(ig,l)) 354 !riceco2(ig,l)=min(5.e-5,riceco2(ig,l)) 355 356 c call updaterice_microCO2(Niceco2,Qccnco2,Nccnco2, 357 c & tauscaling(ig),riceco2(ig,l),rhocloudco2t(ig,l)) 358 c write(*,*) "in co2clouds, rice update = ",riceco2(ig,l) 359 c write(*,*) "in co2clouds, rho update = " 360 c & ,rhocloudco2t(ig,l) 361 362 rsedcloudco2(ig,l)=max(riceco2(ig,l)* 600 601 rsedcloudco2(ig,l)=max(riceco2(ig,l)* 363 602 & (1.+nuiceco2_sed)*(1.+nuiceco2_sed)*(1.+nuiceco2_sed), 364 603 & rdust(ig,l)) 365 rsedcloudco2(ig,l)=min(rsedcloudco2(ig,l),1.e-5) 366 c write(*,*) 'Rsedcloud = ',rsedcloudco2(ig,l) 367 !write(*,*) 'Rhocloudco2 = ',rhocloudco2t(ig,l) 368 604 369 605 ENDDO 370 606 ENDDO 371 607 372 608 ! Gravitational sedimentation 373 374 ! sedimentation computed from radius computed from q in module radii_mod 609 375 610 sav_trac(:,:,igcm_co2_ice)=tempo_traceurs(:,:,igcm_co2_ice) 376 611 sav_trac(:,:,igcm_ccnco2_mass)= … … 384 619 & tempo_traceurs(:,:,igcm_co2_ice),wq,beta) ! 3 traceurs 385 620 386 ! sedim at the surface of co2 ice621 ! sedim at the surface of co2 ice : keep track of it for physiq_mod 387 622 do ig=1,ngrid 388 pdqs_sedco2(ig)=pdqs_sedco2(ig)+ wq(ig,1)623 pdqs_sedco2(ig)=pdqs_sedco2(ig)+ wq(ig,1)/microtimestep 389 624 end do 390 625 … … 398 633 & rsedcloudco2,rhocloudco2t, 399 634 & tempo_traceurs(:,:,igcm_ccnco2_number),wq,beta) 400 401 402 DO l = 1, nlay 635 636 DO l = 1, nlay !Compute tendencies 403 637 DO ig=1,ngrid 404 638 pdqsed(ig,l,igcm_ccnco2_mass)= … … 415 649 !pdqsed est la tendance due a la sedimentation 416 650 417 DO l = 1, nlay418 DO ig=1,ngrid419 pdqsed(ig,l,igcm_ccnco2_mass)=420 & (tempo_traceurs(ig,l,igcm_ccnco2_mass)-421 & sav_trac(ig,l,igcm_ccnco2_mass))/microtimestep422 pdqsed(ig,l,igcm_ccnco2_number)=423 & (tempo_traceurs(ig,l,igcm_ccnco2_number)-424 & sav_trac(ig,l,igcm_ccnco2_number))/microtimestep425 pdqsed(ig,l,igcm_co2_ice)=426 & (tempo_traceurs(ig,l,igcm_co2_ice)-427 & sav_trac(ig,l,igcm_co2_ice))/microtimestep428 ENDDO429 ENDDO430 !pdqsed est la tendance due a la sedimentation431 651 DO l=1,nlay 432 652 DO ig=1,ngrid … … 434 654 & subpdq(ig,l,igcm_ccnco2_mass) 435 655 & +pdqsed(ig,l,igcm_ccnco2_mass) 436 437 656 subpdq(ig,l,igcm_ccnco2_number) = 438 657 & subpdq(ig,l,igcm_ccnco2_number) 439 658 & +pdqsed(ig,l,igcm_ccnco2_number) 440 441 659 subpdq(ig,l,igcm_co2_ice) = 442 660 & subpdq(ig,l,igcm_co2_ice) … … 444 662 ENDDO 445 663 ENDDO 446 c-------------------------------------------------------------------447 c 2. Main call to the different cloud schemes:448 c------------------------------------------------449 IF (microphysco2) THEN450 CALL improvedCO2clouds(ngrid,nlay,microtimestep,451 & pplay,pt,subpdt,452 & pq,subpdq,subpdqcloudco2,subpdtcloudco2,453 & nq,tauscaling,memdMMccn,memdMMh2o,memdNNccn)454 455 ELSE456 457 write(*,*) ' no simpleCO2clouds procedure: STOP' ! listo458 write(*,*) ' microphysco2 and co2clouds must be .true.' ! listo459 STOP460 461 ENDIF462 463 464 c-------------------------------------------------------------------465 c 3. Updating tendencies after cloud scheme:466 c-----------------------------------------------467 468 c IF (microphysco2) THEN469 DO l=1,nlay470 DO ig=1,ngrid471 subpdq(ig,l,igcm_dust_mass) =472 & subpdq(ig,l,igcm_dust_mass)473 & + subpdqcloudco2(ig,l,igcm_dust_mass)474 475 subpdq(ig,l,igcm_dust_number) =476 & subpdq(ig,l,igcm_dust_number)477 & + subpdqcloudco2(ig,l,igcm_dust_number)478 479 subpdq(ig,l,igcm_ccnco2_mass) =480 & subpdq(ig,l,igcm_ccnco2_mass)481 & + subpdqcloudco2(ig,l,igcm_ccnco2_mass)482 c & +pdqsed(ig,l,igcm_ccnco2_mass)483 484 subpdq(ig,l,igcm_ccnco2_number) =485 & subpdq(ig,l,igcm_ccnco2_number)486 & + subpdqcloudco2(ig,l,igcm_ccnco2_number)487 c & +pdqsed(ig,l,igcm_ccnco2_number)488 489 subpdq(ig,l,igcm_co2_ice) =490 & subpdq(ig,l,igcm_co2_ice)491 & + subpdqcloudco2(ig,l,igcm_co2_ice)492 c & +pdqsed(ig,l,igcm_co2_ice)493 494 subpdq(ig,l,igcm_co2) =495 & subpdq(ig,l,igcm_co2)496 & + subpdqcloudco2(ig,l,igcm_co2)497 498 subpdq(ig,l,igcm_h2o_ice) =499 & subpdq(ig,l,igcm_h2o_ice)500 & + subpdqcloudco2(ig,l,igcm_h2o_ice)501 502 subpdq(ig,l,igcm_ccn_mass) =503 & subpdq(ig,l,igcm_ccn_mass)504 & + subpdqcloudco2(ig,l,igcm_ccn_mass)505 506 subpdq(ig,l,igcm_ccn_number) =507 & subpdq(ig,l,igcm_ccn_number)508 & + subpdqcloudco2(ig,l,igcm_ccn_number)509 ENDDO510 ENDDO511 512 513 514 515 IF (activice) THEN516 DO l=1,nlay517 DO ig=1,ngrid518 subpdt(ig,l) =519 & subpdt(ig,l) + subpdtcloudco2(ig,l)520 ENDDO521 ENDDO522 ENDIF523 524 664 525 ENDDO ! of DO microstep=1,imicro665 ENDDO ! of DO microstep=1,imicro 526 666 527 667 c------------------------------------------------------------------- … … 530 670 c CO2 flux at surface (kg.m-2.s-1) 531 671 do ig=1,ngrid 532 pdqs_sedco2(ig)=pdqs_sedco2(ig)/ ptimestep672 pdqs_sedco2(ig)=pdqs_sedco2(ig)/real(imicro) 533 673 enddo 534 674 … … 546 686 547 687 pdqcloudco2(ig,l,igcm_co2_ice) = 548 & subpdq(ig,l,igcm_co2_ice)/real(imicro) 688 & subpdq(ig,l,igcm_co2_ice)/real(imicro) 549 689 & - pdq(ig,l,igcm_co2_ice) 550 551 690 pdqcloudco2(ig,l,igcm_co2) = 552 & subpdq(ig,l,igcm_co2)/real(imicro) 691 & subpdq(ig,l,igcm_co2)/real(imicro) 553 692 & - pdq(ig,l,igcm_co2) 554 555 693 pdqcloudco2(ig,l,igcm_h2o_ice) = 556 & subpdq(ig,l,igcm_h2o_ice)/real(imicro) 694 & subpdq(ig,l,igcm_h2o_ice)/real(imicro) 557 695 & - pdq(ig,l,igcm_h2o_ice) 558 696 ENDDO 559 697 ENDDO 560 698 561 562 ! call WRITEdiagfi(ngrid,"co2cloud00","co2 traceur","kg/kg",1,563 ! & pq(1,:,igcm_co2_ice) + ptimestep*564 ! & (pdq(1,:,igcm_co2_ice) + pdqcloudco2(1,:,igcm_co2_ice)))565 566 567 c IF(microphysco2) THEN568 699 DO l=1,nlay 569 700 DO ig=1,ngrid … … 582 713 ENDDO 583 714 ENDDO 584 c ENDIF585 715 586 716 587 IF(scavenging) THEN588 717 DO l=1,nlay 589 718 DO ig=1,ngrid 590 719 pdqcloudco2(ig,l,igcm_dust_mass) = 591 & subpdq(ig,l,igcm_dust_mass)/real(imicro) 720 & subpdq(ig,l,igcm_dust_mass)/real(imicro) 592 721 & - pdq(ig,l,igcm_dust_mass) 593 722 pdqcloudco2(ig,l,igcm_dust_number) = … … 596 725 ENDDO 597 726 ENDDO 598 ENDIF 599 600 c ENDIF 727 601 728 c------- Due to stepped entry, other processes tendencies can add up to negative values 602 729 c------- Therefore, enforce positive values and conserve mass 603 730 604 731 605 IF(microphysco2) THEN 732 606 733 DO l=1,nlay 607 DO ig=1,ngrid 608 IF ((pq(ig,l,igcm_ccnco2_number) + 609 & ptimestep* (pdq(ig,l,igcm_ccnco2_number) + 610 & pdqcloudco2(ig,l,igcm_ccnco2_number)) 611 & .lt. 1.e-15) 612 & .or. (pq(ig,l,igcm_ccnco2_mass) + 613 & ptimestep* (pdq(ig,l,igcm_ccnco2_mass) + 614 & pdqcloudco2(ig,l,igcm_ccnco2_mass)) 615 & .lt. 1.e-30)) THEN 616 617 pdqcloudco2(ig,l,igcm_ccnco2_number) = 618 & - pq(ig,l,igcm_ccnco2_number)/ptimestep 619 & - pdq(ig,l,igcm_ccnco2_number)+1.e-15 620 621 pdqcloudco2(ig,l,igcm_dust_number) = 622 & -pdqcloudco2(ig,l,igcm_ccnco2_number) 623 624 pdqcloudco2(ig,l,igcm_ccnco2_mass) = 625 & - pq(ig,l,igcm_ccnco2_mass)/ptimestep 626 & - pdq(ig,l,igcm_ccnco2_mass)+1.e-30 627 628 pdqcloudco2(ig,l,igcm_dust_mass) = 629 & -pdqcloudco2(ig,l,igcm_ccnco2_mass) 630 631 ENDIF 632 ENDDO 734 DO ig=1,ngrid 735 IF ((pqeff(ig,l,igcm_ccnco2_number) + 736 & ptimestep* (pdq(ig,l,igcm_ccnco2_number) + 737 & pdqcloudco2(ig,l,igcm_ccnco2_number)) 738 & .lt. 1.e-20) 739 & .or. (pqeff(ig,l,igcm_ccnco2_mass) + 740 & ptimestep* (pdq(ig,l,igcm_ccnco2_mass) + 741 & pdqcloudco2(ig,l,igcm_ccnco2_mass)) 742 & .lt. 1.e-30)) THEN 743 744 pdqcloudco2(ig,l,igcm_ccnco2_number) = 745 & - pqeff(ig,l,igcm_ccnco2_number)/ptimestep 746 & - pdq(ig,l,igcm_ccnco2_number)+1.e-20 747 pdqcloudco2(ig,l,igcm_dust_number) = 748 & -pdqcloudco2(ig,l,igcm_ccnco2_number) 749 pdqcloudco2(ig,l,igcm_ccnco2_mass) = 750 & - pqeff(ig,l,igcm_ccnco2_mass)/ptimestep 751 & - pdq(ig,l,igcm_ccnco2_mass)+1.e-30 752 pdqcloudco2(ig,l,igcm_dust_mass) = 753 & -pdqcloudco2(ig,l,igcm_ccnco2_mass) 754 755 ENDIF 756 ENDDO 633 757 ENDDO 634 ENDIF758 635 759 636 760 637 IF(scavenging) THEN 638 DO l=1,nlay 639 DO ig=1,ngrid 640 IF ( (pq(ig,l,igcm_dust_number) + 641 & ptimestep* (pdq(ig,l,igcm_dust_number) + 642 & pdqcloudco2(ig,l,igcm_dust_number)) .le. 1.e-15) 643 & .or. (pq(ig,l,igcm_dust_mass)+ 644 & ptimestep* (pdq(ig,l,igcm_dust_mass) + 645 & pdqcloudco2(ig,l,igcm_dust_mass)) 646 & .le. 1.e-30)) then 647 648 pdqcloudco2(ig,l,igcm_dust_number) = 649 & - pq(ig,l,igcm_dust_number)/ptimestep 650 & - pdq(ig,l,igcm_dust_number)+1.e-15 651 652 pdqcloudco2(ig,l,igcm_ccnco2_number) = 761 DO l=1,nlay 762 DO ig=1,ngrid 763 IF ( (pqeff(ig,l,igcm_dust_number) + 764 & ptimestep* (pdq(ig,l,igcm_dust_number) + 765 & pdqcloudco2(ig,l,igcm_dust_number)) .le. 1.e-30) 766 & .or. (pqeff(ig,l,igcm_dust_mass)+ 767 & ptimestep* (pdq(ig,l,igcm_dust_mass) + 768 & pdqcloudco2(ig,l,igcm_dust_mass)) 769 & .le. 1.e-30)) then 770 771 pdqcloudco2(ig,l,igcm_dust_number) = 772 & - pqeff(ig,l,igcm_dust_number)/ptimestep 773 & - pdq(ig,l,igcm_dust_number)+1.e-30 774 pdqcloudco2(ig,l,igcm_ccnco2_number) = 653 775 & -pdqcloudco2(ig,l,igcm_dust_number) 654 655 pdqcloudco2(ig,l,igcm_dust_mass) = 656 & - pq(ig,l,igcm_dust_mass)/ptimestep 776 pdqcloudco2(ig,l,igcm_dust_mass) = 777 & - pqeff(ig,l,igcm_dust_mass)/ptimestep 657 778 & - pdq(ig,l,igcm_dust_mass) +1.e-30 658 659 pdqcloudco2(ig,l,igcm_ccnco2_mass) = 660 & -pdqcloudco2(ig,l,igcm_dust_mass) 661 ENDIF 662 ENDDO 663 ENDDO 664 ENDIF !pq+ptime*(pdq+pdqc)=1 ! pdqc=1-pq/ptime-pdq 665 666 667 c DO l=1,nlay 668 c DO ig=1,ngrid 669 c IF (pq(ig,l,igcm_co2_ice) + ptimestep* 670 c & (pdq(ig,l,igcm_co2_ice) + pdqcloudco2(ig,l,igcm_co2_ice)) 671 c & .lt. 1.e-25) THEN 672 c pdqcloudco2(ig,l,igcm_co2_ice) = 673 c & - pq(ig,l,igcm_co2_ice)/ptimestep - pdq(ig,l,igcm_co2_ice) 674 c pdqcloudco2(ig,l,igcm_co2) = -pdqcloudco2(ig,l,igcm_co2_ice) 675 c write(*,*) "WARNING CO2 ICE in co2cloud.F" 676 677 c ENDIF 678 c IF (pq(ig,l,igcm_co2) + ptimestep* 679 c & (pdq(ig,l,igcm_co2) + pdqcloudco2(ig,l,igcm_co2)) 680 c & .lt. 1.e-10) THEN 681 c pdqcloudco2(ig,l,igcm_co2) = 682 c & - pq(ig,l,igcm_co2)/ptimestep - pdq(ig,l,igcm_co2) 683 c pdqcloudco2(ig,l,igcm_co2_ice) = -pdqcloudco2(ig,l,igcm_co2) 684 c write(*,*) "WARNING CO2 in co2cloud.F" 685 c ENDIF 686 c ENDDO 687 c ENDDO 688 689 690 691 692 c------Update the ice and dust particle size "riceco2" for output or photochemistry 693 c------Only rsedcloudco2 is used for the co2 (cloud) cycle 694 695 c IF(scavenging) THEN 696 c DO l=1, nlay 697 c DO ig=1,ngrid 698 699 c call updaterdust( 700 c & pq(ig,l,igcm_dust_mass) + ! dust mass 701 c & (pdq(ig,l,igcm_dust_mass) + ! dust mass 702 c & pdqcloudco2(ig,l,igcm_dust_mass))*ptimestep, ! dust mass 703 c & pq(ig,l,igcm_dust_number) + ! dust number 704 c & (pdq(ig,l,igcm_dust_number) + ! dust number 705 c & pdqcloudco2(ig,l,igcm_dust_number))*ptimestep, ! dust number 706 c & rdust(ig,l)) 707 c write(*,*) "in co2clouds, rdust(ig,l)= ",rdust(ig,l) 708 c mdustJA= pq(ig,l,igcm_dust_mass) + 709 c & (pdq(ig,l,igcm_dust_mass) + 710 c & pdqcloudco2(ig,l,igcm_dust_mass))*ptimestep 711 c ndustJA=pq(ig,l,igcm_dust_number) + 712 c & (pdq(ig,l,igcm_dust_number) + 713 c & pdqcloudco2(ig,l,igcm_dust_number))*ptimestep 714 c if ((ndustJA .lt. tauscaling(ig)) .or. (mdustJA .lt. 715 c & 1.e-30 *tauscaling(ig))) then 716 c rdust(ig,l)=1.e-10 717 c else 718 c rdust(ig,l)=(3./4./pi/2500.*mdustJA/ndustJA)**(1./3.) 719 c rdust(ig,l)=min(rdust(ig,l),5.e-6) 720 c rdust(ig,l)=max(rdust(ig,l),1.e-10) 721 c endif 722 c ENDDO 723 c ENDDO 724 c ENDIF 779 pdqcloudco2(ig,l,igcm_ccnco2_mass) = 780 & -pdqcloudco2(ig,l,igcm_dust_mass) 781 782 ENDIF 783 ENDDO 784 ENDDO 785 !pq+ptime*(pdq+pdqc)=1 ! pdqc=1-pq/ptime-pdq 786 c$$$ 787 c$$$ 788 c$$$ DO l=1,nlay 789 c$$$ DO ig=1,ngrid 790 c$$$ IF (pq(ig,l,igcm_co2_ice) + ptimestep* 791 c$$$ & (pdq(ig,l,igcm_co2_ice) + pdqcloudco2(ig,l,igcm_co2_ice)) 792 c$$$ & .lt. 1.e-30) THEN 793 c$$$ pdqcloudco2(ig,l,igcm_co2_ice) = 794 c$$$ & - pq(ig,l,igcm_co2_ice)/ptimestep - pdq(ig,l,igcm_co2_ice) 795 c$$$ pdqcloudco2(ig,l,igcm_co2) = -pdqcloudco2(ig,l,igcm_co2_ice) 796 c$$$ !write(*,*) "WARNING CO2 ICE in co2cloud.F" 797 c$$$ 798 c$$$ ENDIF 799 c$$$ IF (pq(ig,l,igcm_co2) + ptimestep* 800 c$$$ & (pdq(ig,l,igcm_co2) + pdqcloudco2(ig,l,igcm_co2)) 801 c$$$ & .lt. 0.5) THEN 802 c$$$ pdqcloudco2(ig,l,igcm_co2) = 803 c$$$ & - pdq(ig,l,igcm_co2_ice) !- pdq(ig,l,igcm_co2) 804 c$$$c pdqcloudco2(ig,l,igcm_co2_ice) = -pdqcloudco2(ig,l,igcm_co2) 805 c$$$ pdqcloudco2(ig,l,igcm_co2_ice) = 806 c$$$ & - pq(ig,l,igcm_co2_ice)/ptimestep - pdq(ig,l,igcm_co2_ice) 807 c$$$ ! write(*,*) "WARNING CO2 in co2cloud.F" 808 c$$$ ENDIF 809 c$$$ ENDDO 810 c$$$ ENDDO 811 c$$$ 725 812 726 727 IF(microphysco2) THEN 728 729 DO l=1, nlay 730 DO ig=1,ngrid 813 DO l=1, nlay 814 DO ig=1,ngrid 731 815 732 816 c call updaterice_microco2( … … 742 826 c & tauscaling(ig),riceco2(ig,l),rhocloudco2(ig,l)) 743 827 c write(*,*) "in co2clouds, riceco2(ig,l)= ",riceco2(ig,l) 828 Niceco2=pqeff(ig,l,igcm_co2_ice) + 829 & (pdq(ig,l,igcm_co2_ice) + 830 & pdqcloudco2(ig,l,igcm_co2_ice))*ptimestep 831 Nco2=pqeff(ig,l,igcm_co2) + 832 & (pdq(ig,l,igcm_co2) + 833 & pdqcloudco2(ig,l,igcm_co2))*ptimestep 834 Nccnco2=max((pqeff(ig,l,igcm_ccnco2_number) + 835 & (pdq(ig,l,igcm_ccnco2_number) + 836 & pdqcloudco2(ig,l,igcm_ccnco2_number))*ptimestep)* 837 & tauscaling(ig),1.e-30) 838 Qccnco2=max((pqeff(ig,l,igcm_ccnco2_mass) + 839 & (pdq(ig,l,igcm_ccnco2_mass) + 840 & pdqcloudco2(ig,l,igcm_ccnco2_mass))*ptimestep)* 841 & tauscaling(ig),1.e-30) 842 843 if (Nccnco2 .lt. 0.1) then 844 rdust(ig,l)=1.e-10 845 else 744 846 745 746 Niceco2=pq(ig,l,igcm_co2_ice) + 747 & (pdq(ig,l,igcm_co2_ice) + 748 & pdqcloudco2(ig,l,igcm_co2_ice))*ptimestep 749 Nccnco2=max((pq(ig,l,igcm_ccnco2_number) + 750 & (pdq(ig,l,igcm_ccnco2_number) + 751 & pdqcloudco2(ig,l,igcm_ccnco2_number))*ptimestep)* 752 & tauscaling(ig),1.e-30) 753 Qccnco2=max((pq(ig,l,igcm_ccnco2_mass) + 754 & (pdq(ig,l,igcm_ccnco2_mass) + 755 & pdqcloudco2(ig,l,igcm_ccnco2_mass))*ptimestep)* 756 & tauscaling(ig),1.e-30) 757 758 759 if ((Nccnco2 .lt. 1) .or. (Qccnco2 .lt. 760 & 1.e-30)) then 761 rdust(ig,l)=1.e-10 762 else 763 764 rdust(ig,l)= Qccnco2 765 & *0.75/pi/rho_dust 766 & / Nccnco2 767 rdust(ig,l)= rdust(ig,l)**(1./3.) 768 rdust(ig,l)=max(1.e-10,rdust(ig,l)) 769 !rdust(ig,l)=min(5.e-5,rdust(ig,l)) 770 endif 771 772 rho_ice_co2T(ig,l)=1000.*(1.72391-2.53e-4* 773 & (pt(ig,l)+(subpdt(ig,l)+pdtcloudco2(ig,l))*ptimestep) 774 & -2.87e-6* 775 & (pt(ig,l)+(subpdt(ig,l)+pdtcloudco2(ig,l))*ptimestep)* 776 & (pt(ig,l)+(subpdt(ig,l)+pdtcloudco2(ig,l))*ptimestep)) 777 rho_ice_co2=rho_ice_co2T(ig,l) 778 779 riceco2(ig,l)=(Niceco2*3.0/ 780 & (4.0*rho_ice_co2*pi*Nccnco2) 781 & +rdust(ig,l)*rdust(ig,l) 782 & *rdust(ig,l) )**(1.0/3.0) 783 784 if ( (Niceco2 .le. 1.e-23) .or. 785 & (riceco2(ig,l) .le. rdust(ig,l)) 786 & .or. (Nccnco2 .le. 0.01))then 787 788 c write(*,*) "Riceco2 co2cloud before reset=",riceco2(ig,l) 789 c write(*,*) "Niceco2 co2cloud before reset=",Niceco2 790 c write(*,*) "Nccnco2 co2cloud before reset=",Nccnco2 791 c write(*,*) "Rdust co2cloud before reset=",rdust(ig,l) 792 793 c$$$ !NO CLOUD : RESET TRACERS AND CONSERVE MASS 794 if (pq(ig,l,igcm_co2_ice)+(pdq(ig,l,igcm_co2_ice)+ 795 & pdqcloudco2(ig,l,igcm_co2_ice))*ptimestep .le. 0.) then 796 pdqcloudco2(ig,l,igcm_co2)=0. 797 pdqcloudco2(ig,l,igcm_co2_ice)=0. 798 else 799 pdqcloudco2(ig,l,igcm_co2)= pq(ig,l,igcm_co2_ice) 800 & /ptimestep+pdq(ig,l,igcm_co2_ice) 801 pdqcloudco2(ig,l,igcm_co2_ice)=-pq(ig,l,igcm_co2_ice) 847 rdust(ig,l)= Qccnco2 848 & *0.75/pi/rho_dust 849 & / Nccnco2 850 rdust(ig,l)= rdust(ig,l)**(1./3.) 851 rdust(ig,l)=max(1.e-10,rdust(ig,l)) 852 c rdust(ig,l)=min(5.e-6,rdust(ig,l)) 853 endif 854 myT=zteff(ig,l)+(pdt(ig,l)+pdtcloudco2(ig,l))*ptimestep 855 rho_ice_co2T(ig,l)=1000.*(1.72391-2.53e-4* 856 & myT-2.87e-6* myT* myT) 857 rho_ice_co2=rho_ice_co2T(ig,l) 858 859 lw = l0 + l1 * myT + l2 *myT**2 + 860 & l3 * myT**3 + l4 * myT**4 !J.kg-1 861 862 riceco2(ig,l)=(Niceco2*3.0/ 863 & (4.0*rho_ice_co2*pi*Nccnco2) 864 & +rdust(ig,l)*rdust(ig,l) 865 & *rdust(ig,l) )**(1.0/3.0) 866 c & .or. (riceco2(ig,l) .le. rdust(ig,l)) 867 if ( (Niceco2 .le. 1.e-25).or. (Nccnco2 .le. 0.1) )THEN !anciennement 200 868 c riceco2(ig,l)=0. 869 870 c & .or. (Nccnco2 .le. 1.) 871 c endif 872 !Flux chaleur latente <0 quand sublimation 873 874 pdtcloudco2(ig,l)= pdtcloudco2(ig,l)-Niceco2*lw/cpp/ptimestep 875 c$$$ !NO CLOUD : RESET TRACERS AND CONSERVE MASS 876 c if (pq(ig,l,igcm_co2_ice)+(pdq(ig,l,igcm_co2_ice)+ 877 c & pdqcloudco2(ig,l,igcm_co2_ice))*ptimestep .le. 0.) then 878 c pdqcloudco2(ig,l,igcm_co2)=0. 879 c pdqcloudco2(ig,l,igcm_co2_ice)=0. 880 c else 881 pdqcloudco2(ig,l,igcm_co2_ice)=-1.*pqeff(ig,l,igcm_co2_ice) 802 882 & /ptimestep-pdq(ig,l,igcm_co2_ice) 803 endif 804 !Reverse h2o ccn and ices into h2o tracers 883 pdqcloudco2(ig,l,igcm_co2)=-1.* 884 & pdqcloudco2(ig,l,igcm_co2_ice) 885 c endif 886 ! Reverse h2o ccn and ices into h2o tracers 805 887 if (memdMMccn(ig,l) .gt. 0) then 806 pdqcloudco2(ig,l,igcm_ccn_mass)=memdMMccn(ig,l) 888 pdqcloudco2(ig,l,igcm_ccn_mass)=memdMMccn(ig,l)/ptimestep 807 889 else 808 890 memdMMccn(ig,l)=0. … … 810 892 endif 811 893 if (memdNNccn(ig,l) .gt. 0) then 812 pdqcloudco2(ig,l,igcm_ccn_number)=memdNNccn(ig,l)894 pdqcloudco2(ig,l,igcm_ccn_number)=memdNNccn(ig,l)/ptimestep 813 895 else 814 896 memdNNccn(ig,l)=0. … … 816 898 endif 817 899 if (memdMMh2o(ig,l) .gt. 0) then 818 pdqcloudco2(ig,l,igcm_h2o_ice)=memdMMh2o(ig,l) 900 pdqcloudco2(ig,l,igcm_h2o_ice)=memdMMh2o(ig,l)/ptimestep 819 901 else 820 902 memdMMh2o(ig,l)=0. 821 903 pdqcloudco2(ig,l,igcm_h2o_ice)=0. 822 904 endif 823 824 if (pq(ig,l,igcm_ccnco2_mass)+(pdq(ig,l,igcm_ccnco2_mass)+ 825 & pdqcloudco2(ig,l,igcm_ccnco2_mass))*ptimestep 826 & .le. 1.e-30) then 827 pdqcloudco2(ig,l,igcm_ccnco2_mass)=0. 828 pdqcloudco2(ig,l,igcm_ccnco2_number)=0. 829 pdqcloudco2(ig,l,igcm_co2)=0. 830 pdqcloudco2(ig,l,igcm_co2_ice)=0. 831 else 905 c if (pq(ig,l,igcm_ccnco2_mass)+(pdq(ig,l,igcm_ccnco2_mass)+ 906 c & pdqcloudco2(ig,l,igcm_ccnco2_mass))*ptimestep 907 c & .le. 1.e-30) then 908 c pdqcloudco2(ig,l,igcm_ccnco2_mass)=0. 909 c pdqcloudco2(ig,l,igcm_ccnco2_number)=0. 910 c pdqcloudco2(ig,l,igcm_co2)=0. 911 c pdqcloudco2(ig,l,igcm_co2_ice)=0. 912 c else 832 913 pdqcloudco2(ig,l,igcm_ccnco2_mass)= 833 & - pq(ig,l,igcm_ccnco2_mass)914 & -1.*pqeff(ig,l,igcm_ccnco2_mass) 834 915 & /ptimestep-pdq(ig,l,igcm_ccnco2_mass) 835 endif 836 837 if (pq(ig,l,igcm_ccnco2_number)+ 838 & (pdq(ig,l,igcm_ccnco2_number)+ 839 & pdqcloudco2(ig,l,igcm_ccnco2_number)) 840 & *ptimestep.le. 1.e-30) 841 & then 842 pdqcloudco2(ig,l,igcm_ccnco2_mass)=0. 843 pdqcloudco2(ig,l,igcm_ccnco2_number)=0. 844 pdqcloudco2(ig,l,igcm_co2)=0. 845 pdqcloudco2(ig,l,igcm_co2_ice)=0. 846 else 916 c endif 917 c if (pq(ig,l,igcm_ccnco2_number)+ 918 c & (pdq(ig,l,igcm_ccnco2_number)+ 919 c & pdqcloudco2(ig,l,igcm_ccnco2_number)) 920 c & *ptimestep.le. 1.e-30) 921 c & then 922 c pdqcloudco2(ig,l,igcm_ccnco2_mass)=0. 923 c pdqcloudco2(ig,l,igcm_ccnco2_number)=0. 924 c pdqcloudco2(ig,l,igcm_co2)=0. 925 c pdqcloudco2(ig,l,igcm_co2_ice)=0. 926 c else 847 927 pdqcloudco2(ig,l,igcm_ccnco2_number)= 848 & - pq(ig,l,igcm_ccnco2_number)928 & -1.*pqeff(ig,l,igcm_ccnco2_number) 849 929 & /ptimestep-pdq(ig,l,igcm_ccnco2_number) 850 endif 851 852 if (pq(ig,l,igcm_dust_number)+ 853 & (pdq(ig,l,igcm_dust_number)+ 854 & pdqcloudco2(ig,l,igcm_dust_number)) 855 & *ptimestep.le. 1.e-30) 856 & then 857 pdqcloudco2(ig,l,igcm_dust_number)=0. 858 pdqcloudco2(ig,l,igcm_dust_mass)=0. 859 else 930 c endif 931 c if (pq(ig,l,igcm_dust_number)+ 932 c & (pdq(ig,l,igcm_dust_number)+ 933 c & pdqcloudco2(ig,l,igcm_dust_number)) 934 c & *ptimestep.le. 1.e-30) 935 c & then 936 c pdqcloudco2(ig,l,igcm_dust_number)=0. 937 c pdqcloudco2(ig,l,igcm_dust_mass)=0. 938 c else 860 939 pdqcloudco2(ig,l,igcm_dust_number)= 861 & pq (ig,l,igcm_ccnco2_number)940 & pqeff(ig,l,igcm_ccnco2_number) 862 941 & /ptimestep+pdq(ig,l,igcm_ccnco2_number) 863 & -memdNNccn(ig,l) 864 endif 865 866 if (pq(ig,l,igcm_dust_mass)+ 867 & (pdq(ig,l,igcm_dust_mass)+ 868 & pdqcloudco2(ig,l,igcm_dust_mass)) 869 & *ptimestep .le. 1.e-30) 870 & then 871 pdqcloudco2(ig,l,igcm_dust_number)=0. 872 pdqcloudco2(ig,l,igcm_dust_mass)=0. 873 else 942 & -memdNNccn(ig,l)/ptimestep 943 c endif 944 c if (pq(ig,l,igcm_dust_mass)+ 945 c & (pdq(ig,l,igcm_dust_mass)+ 946 c & pdqcloudco2(ig,l,igcm_dust_mass)) 947 c & *ptimestep .le. 1.e-30) 948 c & then 949 c pdqcloudco2(ig,l,igcm_dust_number)=0. 950 c pdqcloudco2(ig,l,igcm_dust_mass)=0. 951 c else 874 952 pdqcloudco2(ig,l,igcm_dust_mass)= 875 & pq (ig,l,igcm_ccnco2_mass)953 & pqeff(ig,l,igcm_ccnco2_mass) 876 954 & /ptimestep+pdq(ig,l,igcm_ccnco2_mass) 877 & -(memdMMh2o(ig,l)+memdMMccn(ig,l)) 878 endif955 & -(memdMMh2o(ig,l)+memdMMccn(ig,l))/ptimestep 956 c endif 879 957 memdMMccn(ig,l)=0. 880 958 memdMMh2o(ig,l)=0. 881 959 memdNNccn(ig,l)=0. 882 960 riceco2(ig,l)=0. 883 pdtcloudco2(ig,l)=0.884 961 endif 885 886 887 962 c Compute opacities 963 No=Nccnco2 964 Rn=-log(riceco2(ig,l)) 965 n_derf = erf( (rb_cldco2(1)+Rn) *dev2) 966 Qext1bins2(ig,l)=0. 967 do i = 1, nbinco2_cld !Qext below 50 is negligible 968 n_aer(i) = -0.5 * No * n_derf !! this ith previously computed 969 n_derf = derf((rb_cldco2(i+1)+Rn) *dev2) 970 n_aer(i) = n_aer(i) + 0.5 * No * n_derf 971 Qext1bins2(ig,l)=Qext1bins2(ig,l)+Qext1bins(i)*n_aer(i) 972 enddo 973 !l'opacité de la case ig est la somme sur l de Qext1bins2 974 975 888 976 !update rice water 889 977 call updaterice_micro( 890 & pq (ig,l,igcm_h2o_ice) + ! ice mass978 & pqeff(ig,l,igcm_h2o_ice) + ! ice mass 891 979 & (pdq(ig,l,igcm_h2o_ice) + ! ice mass 892 980 & pdqcloudco2(ig,l,igcm_h2o_ice))*ptimestep, ! ice mass 893 & pq (ig,l,igcm_ccn_mass) + ! ccn mass981 & pqeff(ig,l,igcm_ccn_mass) + ! ccn mass 894 982 & (pdq(ig,l,igcm_ccn_mass) + ! ccn mass 895 983 & pdqcloudco2(ig,l,igcm_ccn_mass))*ptimestep, ! ccn mass 896 & pq (ig,l,igcm_ccn_number) + ! ccn number984 & pqeff(ig,l,igcm_ccn_number) + ! ccn number 897 985 & (pdq(ig,l,igcm_ccn_number) + ! ccn number 898 986 & pdqcloudco2(ig,l,igcm_ccn_number))*ptimestep, ! ccn number … … 901 989 902 990 call updaterdust( 903 & pq (ig,l,igcm_dust_mass) + ! dust mass991 & pqeff(ig,l,igcm_dust_mass) + ! dust mass 904 992 & (pdq(ig,l,igcm_dust_mass) + ! dust mass 905 993 & pdqcloudco2(ig,l,igcm_dust_mass))*ptimestep, ! dust mass 906 & pq (ig,l,igcm_dust_number) + ! dust number994 & pqeff(ig,l,igcm_dust_number) + ! dust number 907 995 & (pdq(ig,l,igcm_dust_number) + ! dust number 908 996 & pdqcloudco2(ig,l,igcm_dust_number))*ptimestep, ! dust number … … 911 999 ENDDO 912 1000 ENDDO 913 914 915 916 ELSE ! no microphys ! not of concern for co2 clouds - listo 917 918 ENDIF ! of IF(microphysco2) 919 920 921 c TO CHEK for relevancy - listo 922 923 c A correction if a lot of subliming CO2 fills the 1st layer FF04/2005 924 c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 925 c Then that should not affect the ice particle radius 926 !do ig=1,ngrid 927 ! if(pdpsrf(ig)*ptimestep.gt.0.9*(pplev(ig,1)-pplev(ig,2)))then 928 ! if(pdpsrf(ig)*ptimestep.gt.0.9*(pplev(ig,1)-pplev(ig,3))) 929 &! riceco2(ig,2)=riceco2(ig,3) 930 ! riceco2(ig,1)=riceco2(ig,2) 931 ! end if 932 !end do 933 1001 tau1mic(:)=0. 1002 do l=1,nlay 1003 do ig=1,ngrid 1004 tau1mic(ig)=tau1mic(ig)+Qext1bins2(ig,l) 1005 enddo 1006 enddo 1007 1008 c$$$ 1009 c$$$ if (riceco2(725,22) .ge. 1.e-10) then 1010 c$$$ 1011 c$$$ write(*,*) 'DIAGJA co2 ',pqeff(725,22,igcm_co2) + 1012 c$$$ & (pdq(725,22,igcm_co2) + 1013 c$$$ & pdqcloudco2(725,22,igcm_co2))*ptimestep 1014 c$$$ write(*,*) 'DIAGJA co2_ice',pqeff(725,22,igcm_co2_ice) + 1015 c$$$ & (pdq(725,22,igcm_co2_ice) + 1016 c$$$ & pdqcloudco2(725,22,igcm_co2_ice))*ptimestep 1017 c$$$ 1018 c$$$ write(*,*) 'DIAGJA riceco2',riceco2(725,22) 1019 c$$$ write(*,*) 'DIAGJA T',zteff(725,22) + 1020 c$$$ & (pdt(725,22) + pdtcloudco2(725,22))*ptimestep 1021 c$$$ write(*,*) 'DIAG pdtcloud',pdtcloudco2(725,22)*ptimestep 1022 c$$$ write(*,*) 'DIAGJA ccnNco2',pqeff(725,22,igcm_ccnco2_number)+ 1023 c$$$ & (pdq(725,22,igcm_ccnco2_number) + 1024 c$$$ & pdqcloudco2(725,22,igcm_ccnco2_number))*ptimestep 1025 c$$$ 1026 c$$$ write(*,*) 'DIAGJA dustN',pqeff(725,22,igcm_dust_number) + 1027 c$$$ & (pdq(725,22,igcm_dust_number) + 1028 c$$$ & pdqcloudco2(725,22,igcm_dust_number))*ptimestep 1029 c$$$ ENDIF 1030 c$$$ 1031 1032 1033 934 1034 c A correction if a lot of subliming CO2 fills the 1st layer FF04/2005 935 1035 c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ … … 958 1058 & (1.+nuiceco2_sed)*(1.+nuiceco2_sed)*(1.+nuiceco2_sed), 959 1059 & rdust(ig,l)) 960 rsedcloudco2(ig,l)=min(rsedcloudco2(ig,l),1.e-5)1060 c rsedcloudco2(ig,l)=min(rsedcloudco2(ig,l),1.e-5) 961 1061 ENDDO 962 1062 ENDDO 963 1063 964 call co2sat(ngrid*nlay, pt,pplay,zqsatco2)965 do ig=1,ngrid966 967 satuco2(ig,l) = (pq(ig,l,igcm_co2) +968 & (pdq(ig,l,igcm_co2) +969 & pdqcloudco2(ig,l,igcm_co2))*ptimestep)*970 & (mmean(ig,l)/44.01)*pplay(ig,l)/zqsatco2(ig,l)971 972 cwrite(*,*) "In CO2 pt,sat ",pt(ig,l),satuco2(ig,l)1064 call co2sat(ngrid*nlay,zteff+(pdt+pdtcloudco2)*ptimestep 1065 & ,pplay,zqsatco2) 1066 do l=1,nlay 1067 do ig=1,ngrid 1068 satuco2(ig,l) = (pqeff(ig,l,igcm_co2) + 1069 & (pdq(ig,l,igcm_co2) + 1070 & pdqcloudco2(ig,l,igcm_co2))*ptimestep)* 1071 & (mmean(ig,l)/44.01)*pplay(ig,l)/zqsatco2(ig,l) 1072 !write(*,*) "In CO2 pt,sat ",pt(ig,l),satuco2(ig,l) 973 1073 enddo 974 1074 enddo 975 call WRITEDIAGFI(ngrid,"satuco2","vap in satu","kg/kg",3, 1075 !Tout ce qui est modifié par la microphysique de CO2 doit être rapporté a cloudfrac 1076 IF (CLFvaryingCO2) THEN 1077 DO l=1,nlay 1078 DO ig=1,ngrid 1079 pdqcloudco2(ig,l,igcm_ccn_mass)= 1080 & pdqcloudco2(ig,l,igcm_ccn_mass)*cloudfrac(ig,l) 1081 pdqcloudco2(ig,l,igcm_ccnco2_mass)= 1082 & pdqcloudco2(ig,l,igcm_ccnco2_mass)*cloudfrac(ig,l) 1083 pdqcloudco2(ig,l,igcm_ccn_number)= 1084 & pdqcloudco2(ig,l,igcm_ccn_number)*cloudfrac(ig,l) 1085 pdqcloudco2(ig,l,igcm_ccnco2_number)= 1086 & pdqcloudco2(ig,l,igcm_ccnco2_number)*cloudfrac(ig,l) 1087 pdqcloudco2(ig,l,igcm_dust_mass)= 1088 & pdqcloudco2(ig,l,igcm_dust_mass)*cloudfrac(ig,l) 1089 pdqcloudco2(ig,l,igcm_dust_number)= 1090 & pdqcloudco2(ig,l,igcm_dust_number)*cloudfrac(ig,l) 1091 pdqcloudco2(ig,l,igcm_h2o_ice)= 1092 & pdqcloudco2(ig,l,igcm_h2o_ice)*cloudfrac(ig,l) 1093 pdqcloudco2(ig,l,igcm_co2_ice)= 1094 & pdqcloudco2(ig,l,igcm_co2_ice)*cloudfrac(ig,l) 1095 pdqcloudco2(ig,l,igcm_co2)= 1096 & pdqcloudco2(ig,l,igcm_co2)*cloudfrac(ig,l) 1097 pdtcloudco2(ig,l)=pdtcloudco2(ig,l)*cloudfrac(ig,l) 1098 ENDDO 1099 ENDDO 1100 ENDIF 1101 1102 call WRITEDIAGFI(ngrid,"satuco2","vap in satu","kg/kg",3, 976 1103 & satuco2) 977 call WRITEdiagfi(ngrid,"riceco2","ice radius","m"1104 call WRITEdiagfi(ngrid,"riceco2","ice radius","m" 978 1105 & ,3,riceco2) 979 ! or output in diagfi.nc (for testphys1d)1106 ! or output in diagfi.nc (for testphys1d) 980 1107 c call WRITEDIAGFI(ngrid,'ps','Surface pressure','Pa',0,ps) 981 1108 c call WRITEDIAGFI(ngrid,'temp','Temperature ', … … 985 1112 & rsedcloudco2) 986 1113 1114 call WRITEdiagfi(ngrid,"tau1mic","co2 ice opacity 1 micron"," ",2, 1115 & tau1mic) 1116 call WRITEdiagfi(ngrid,"cloudfrac","co2 cloud fraction"," ",3, 1117 & cloudfrac) 987 1118 ! used for rad. transfer calculations 988 1119 ! nuice is constant because a lognormal distribution is prescribed -
trunk/LMDZ.MARS/libf/phymars/conf_phys.F
r1711 r1720 36 36 USE ioipsl_getincom, only : getin 37 37 use tracer_mod, only : nuice_sed, ccn_factor, nuiceco2_sed, 38 & nuice_ref,nuiceco2_ref, 39 & meteo_flux_mass, 40 & meteo_flux_number,meteo_alt 38 & nuice_ref,nuiceco2_ref 41 39 use surfdat_h, only: albedo_h2o_ice, inert_h2o_ice, 42 40 & frost_albedo_threshold … … 450 448 451 449 !CO2 clouds scheme? 452 write(*,*) "Compute CO2 clouds ?"450 write(*,*) "Compute CO2 clouds (implies microphysical scheme)?" 453 451 co2clouds=.false. ! default value 454 452 call getin("co2clouds",co2clouds) 455 453 write(*,*) " co2clouds = ",co2clouds 454 !Can water ice particles serve as CCN for CO2clouds 455 write(*,*) "Use water ice as CO2 clouds CCN ?" 456 co2useh2o=.false. ! default value 457 call getin("co2useh2o",co2useh2o) 458 write(*,*) " co2useh2o = ",co2useh2o 459 !Do we allow a supply of meteoritic paricles to serve as CO2 ice CCN? 460 write(*,*) "Supply meteoritic particle for CO2 clouds ?" 461 meteo_flux=.false. !Default value 462 call getin("meteo_flux",meteo_flux) 463 write(*,*) " meteo_flux = ",meteo_flux 464 !Do we allow a sub-grid temperature distribution for the CO2 microphysics 465 write(*,*) "sub-grid temperature distribution for CO2 clouds?" 466 CLFvaryingCO2=.false. !Default value 467 call getin("CLFvaryingCO2",CLFvaryingCO2) 468 write(*,*) " CLFvaryingCO2 = ",CLFvaryingCO2 469 !Amplitude of the sub-grid temperature distribution for the CO2 microphysics 470 write(*,*) "sub-grid temperature amplitude for CO2 clouds?" 471 spantCO2=0 !Default value 472 call getin("spantCO2",spantCO2) 473 write(*,*) " spantCO2 = ",spantCO2 474 456 475 ! thermal inertia feedback 457 476 write(*,*) "Activate the thermal inertia feedback ?" … … 514 533 write(*,*)"the contact parameter is used instead;" 515 534 516 !JAud Meteoritic flux of dust : number, mass and altitude 517 !Jaud dust tracers are incremented by thoses values in co2cloud.F 518 !Jaud Thus only apply if co2clouds=.true. in callphys.def 519 write(*,*) "mass mixing ratio for meteoritic dust?" 520 meteo_flux_mass=0. 521 call getin("meteo_flux_mass",meteo_flux_mass) 522 write(*,*) "number density for meteoritic dust?" 523 meteo_flux_number=0. 524 call getin("meteo_flux_number",meteo_flux_number) 525 write(*,*) "altitude of injection?" 526 meteo_alt=5. 527 call getin("meteo_alt",meteo_alt) 528 WRITE(*,*) "" 529 WRITE(*,*) "meteoritic flux of dust" 530 WRITE(*,*) "dust mass = ",meteo_flux_mass 531 WRITE(*,*) "dust number = ",meteo_flux_number 532 WRITE(*,*) "at altitude = ",meteo_alt 533 534 ! microphys 535 write(*,*)"Microphysical scheme for water-ice clouds?" 536 microphys=.false. ! default value 537 call getin("microphys",microphys) 538 write(*,*)" microphys = ",microphys 539 540 write(*,*)"Microphysical scheme for CO2 clouds?" 541 microphysco2=.false. ! default value 542 call getin("microphysco2",microphysco2) 543 write(*,*)" microphysco2 = ",microphysco2 544 ! supersat 545 write(*,*)"Allow super-saturation of water vapor?" 546 supersat=.true. ! default value 547 call getin("supersat",supersat) 548 write(*,*)"supersat = ",supersat 535 ! microphys 536 write(*,*)"Microphysical scheme for water-ice clouds?" 537 microphys=.false. ! default value 538 call getin("microphys",microphys) 539 write(*,*)" microphys = ",microphys 540 541 ! supersat 542 write(*,*)"Allow super-saturation of water vapor?" 543 supersat=.true. ! default value 544 call getin("supersat",supersat) 545 write(*,*)"supersat = ",supersat 549 546 550 547 ! microphysical parameter contact 551 552 553 554 555 548 write(*,*) "water contact parameter ?" 549 mteta = 0.95 550 call getin("mteta",mteta) 551 write(*,*) "mteta = ", mteta 552 556 553 ! scavenging 557 558 scavenging=.false.! default value559 560 554 write(*,*)"Dust scavenging by H2O/CO2 snowfall ?" 555 scavenging=.false. ! default value 556 call getin("scavenging",scavenging) 557 write(*,*)" scavenging = ",scavenging 561 558 562 559 … … 564 561 ! if scavenging is used, then dustbin should be > 0 565 562 566 563 if ((microphys.and..not.doubleq).or. 567 564 & (microphys.and..not.water)) then 568 569 570 571 572 573 574 575 576 577 578 579 580 if ((scavenging.and..not.microphys.and..not. microphysco2).or.565 print*,'if microphys is used, then doubleq,' 566 print*,'and water must be used!' 567 stop 568 endif 569 if (microphys.and..not.scavenging) then 570 print*,'' 571 print*,'----------------WARNING-----------------' 572 print*,'microphys is used without scavenging !!!' 573 print*,'----------------WARNING-----------------' 574 print*,'' 575 endif 576 577 if ((scavenging.and..not.microphys).or. 581 578 & (scavenging.and.(dustbin.lt.1)))then 582 583 584 585 579 print*,'if scavenging is used, then microphys' 580 print*,'must be used!' 581 stop 582 endif 586 583 587 584 ! Test of incompatibility: -
trunk/LMDZ.MARS/libf/phymars/improvedCO2clouds.F
r1685 r1720 1 1 subroutine improvedCO2clouds(ngrid,nlay,ptimestep, 2 & pplay,p t,pdt,2 & pplay,pzlev,pt,pdt, 3 3 & pq,pdq,pdqcloudco2,pdtcloudco2, 4 4 & nq,tauscaling, … … 15 15 ! & igcm_ccnco2_number 16 16 use conc_mod, only: mmean 17 17 18 implicit none 18 19 19 20 20 c------------------------------------------------------------------ … … 53 53 !#include "dimradmars.h" 54 54 #include "microphys.h" 55 #include "datafile.h" 55 56 !#include "microphysCO2.h" 56 57 !#include "conc.h" … … 62 63 REAL ptimestep ! pas de temps physique (s) 63 64 REAL pplay(ngrid,nlay) ! pression au milieu des couches (Pa) 64 65 REAL pzlev(ngrid,nlay) ! altitude au milieu des couches () 66 65 67 REAL pt(ngrid,nlay) ! temperature at the middle of the 66 68 ! layers (K) … … 74 76 REAL rice(ngrid,nlay) ! Water Ice mass mean radius (m) 75 77 ! used for nucleation of CO2 on ice-coated ccns 78 REAL rccnh2o(ngrid,nlay) ! Water Ice mass mean radius (m) 76 79 77 80 c Outputs: … … 95 98 !external massflowrateCO2 96 99 97 INTEGER ig,l,i 100 INTEGER ig,l,i,flag_pourri 98 101 99 102 REAL zq(ngrid,nlay,nq) ! local value of tracers … … 102 105 REAL zqsat(ngrid,nlay) ! saturation vapor pressure for CO2 103 106 REAL lw !Latent heat of sublimation (J.kg-1) 104 REAL l0,l1,l2,l3,l4107 REAL,save :: l0,l1,l2,l3,l4 105 108 REAL cste 106 109 DOUBLE PRECISION dMice ! mass of condensed ice … … 109 112 DOUBLE PRECISION satu ! Co2 vapor saturation ratio over ice 110 113 DOUBLE PRECISION Mo,No 111 DOUBLE PRECISION Rn, Rm, dev2, n_derf, m_derf112 DOUBLE PRECISION 113 DOUBLE PRECISION 114 DOUBLE PRECISION Rn, Rm, dev2,dev3, n_derf, m_derf 115 DOUBLE PRECISION memdMMccn(ngrid,nlay) 116 DOUBLE PRECISION memdMMh2o(ngrid,nlay) 114 117 DOUBLE PRECISION memdNNccn(ngrid,nlay) 115 118 … … 117 120 DOUBLE PRECISION n_aer(nbinco2_cld) ! number concentration volume-1 of particle/each size bin 118 121 DOUBLE PRECISION m_aer(nbinco2_cld) ! mass mixing ratio of particle/each size bin 122 DOUBLE PRECISION m_aer2(nbinco2_cld) ! mass mixing ratio of particle/each size bin 123 DOUBLE PRECISION m_aer_h2oice2(nbinco2_cld) ! mass mixing ratio of particle/each size bin 119 124 120 125 DOUBLE PRECISION n_aer_h2oice(nbinco2_cld) ! Same - for CO2 nucleation … … 126 131 127 132 DOUBLE PRECISION dN,dM, dNh2o, dMh2o, dNN,dMM,dNNh2o,dMMh2o 133 DOUBLE PRECISION dMh2o_ice,dMh2o_ccn 134 128 135 DOUBLE PRECISION rate(nbinco2_cld) ! nucleation rate 129 136 DOUBLE PRECISION rateh2o(nbinco2_cld) ! nucleation rate … … 142 149 c Parameters of the size discretization 143 150 c used by the microphysical scheme 144 DOUBLE PRECISION, PARAMETER :: rmin_cld = 1.e-1 0! Minimum radius (m)145 DOUBLE PRECISION, PARAMETER :: rmax_cld = 5.e-4! Maximum radius (m)146 DOUBLE PRECISION, PARAMETER :: rbmin_cld =1.e-1 1147 ! Minimum boun dary radius (m)148 DOUBLE PRECISION, PARAMETER :: rbmax_cld = 1.e-3! Maximum boundary radius (m)151 DOUBLE PRECISION, PARAMETER :: rmin_cld = 1.e-11 ! Minimum radius (m) 152 DOUBLE PRECISION, PARAMETER :: rmax_cld = 3.e-6 ! Maximum radius (m) 153 DOUBLE PRECISION, PARAMETER :: rbmin_cld =1.e-12 154 ! Minimum bounary radius (m) 155 DOUBLE PRECISION, PARAMETER :: rbmax_cld = 5.e-6 ! Maximum boundary radius (m) 149 156 DOUBLE PRECISION vrat_cld ! Volume ratio 150 157 DOUBLE PRECISION rb_cldco2(nbinco2_cld+1) ! boundary values of each rad_cldco2 bin (m) … … 157 164 REAL sigma_iceco2 ! Variance of the ice and CCN distributions 158 165 SAVE sigma_iceco2 166 167 168 REAL sigma_ice ! Variance of the ice and CCN distributions 169 SAVE sigma_ice 159 170 160 171 DOUBLE PRECISION Niceco2,Qccnco2,Nccnco2 161 172 173 integer coeffh2o 174 !Variables for the meteoritic flux 175 176 double precision meteo_ccn(ngrid,nlay,100) !100=nbinco2_cld !!! 177 double precision,save :: meteo(130,100) 178 double precision mtemp(100) 179 logical file_ok 180 integer altitudes_meteo(130),nelem,lebon1,lebon2 181 double precision :: ltemp1(130),ltemp2(130) 182 integer ibin,uMeteo,j 162 183 c---------------------------------- 163 184 c TESTS 164 185 165 INTEGER countcells166 167 LOGICAL test_flag ! flag for test/debuging outputs168 SAVE test_flag169 186 170 187 … … 228 245 !! at each timestep and gridpoint 229 246 enddo 230 231 247 c Contact parameter of co2 ice on dst ( m=cos(theta) ) 232 248 c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ … … 240 256 c Variance of the ice and CCN distributions 241 257 sigma_iceco2 = sqrt(log(1.+nuiceco2_sed)) 242 258 sigma_ice = sqrt(log(1.+nuice_sed)) 259 243 260 244 261 c write(*,*) 'Variance of ice & CCN distribs :', sigma_iceco2 … … 249 266 write(*,*) 'nuice for co2 ice sedimentation:', nuiceco2_sed 250 267 write(*,*) 'Volume of a co2 molecule:', vo1co2 251 252 253 254 test_flag = .false. 268 269 coeffh2o=0 270 if (co2useh2o) then 271 write(*,*) 272 write(*,*) "co2useh2o=.true. in callphys.def" 273 write(*,*) "This means water ice particles can " 274 write(*,*) "serve as CCN for CO2 microphysics" 275 coeffh2o=1 276 endif 277 meteo_ccn(:,:,:)=0. 278 279 if (meteo_flux) then 280 write(*,*) 281 write(*,*) "meteo_flux=.true. in callphys.def" 282 write(*,*) "meteoritic dust particles are available" 283 write(*,*) "for co2 ice nucleation! " 284 write(*,*) "Flux given by J. Plane (altitude,size bins)" 285 ! Initialisation of the flux: it is constant and is it saved 286 !We must interpolate the table to the GCM altitudes 287 INQUIRE(FILE=datafile(1:LEN_TRIM(datafile))// 288 & '/Meteo_flux_Plane.dat', EXIST=file_ok) 289 IF (.not. file_ok) THEN 290 write(*,*) 'file Meteo_flux_Plane.dat should be in ' 291 & ,datafile 292 STOP 293 endif 294 !used Variables 295 open(newunit=uMeteo,file=trim(datafile)// 296 & '/Meteo_flux_Plane.dat' 297 & ,FORM='formatted') 298 !13000 records (130 altitudes x 100 bin sizes) 299 do i=1,130 300 do j=1,100 ! les mêmes 100 bins size que la distri nuclea: on touche pas 301 read(uMeteo,'(F12.6)') meteo(i,j) 302 enddo 303 altitudes_meteo(i)=i ! On doit maintenant réinterpoler le tableau(130,100) sur les altitudes du GCM (nlay,100) 304 enddo 305 close(uMeteo) 306 endif !of if meteo_flux 307 308 l0=595594d0 309 l1=903.111d0 310 l2=-11.5959d0 311 l3=0.0528288d0 312 l4=-0.000103183d0 255 313 firstcall=.false. 256 314 315 257 316 END IF 258 259 317 c write(*,*) "max memdNN =",maxval(memdNNccn) 260 318 !============================================================= 261 319 ! 1. Initialisation 262 320 !============================================================= 263 321 !cste = 4*pi*rho_ice*ptimestep !not used for co2 322 flag_pourri=0 264 323 265 324 res_out(:,:) = 0 … … 295 354 !============================================================= 296 355 356 297 357 298 358 dev2 = 1. / ( sqrt(2.) * sigma_iceco2 ) 359 dev3 = 1. / ( sqrt(2.) * sigma_ice ) 299 360 300 361 call co2sat(ngrid*nlay,zt,pplay,zqsat) !zqsat is psat(co2) 301 362 302 countcells = 0 303 304 c Faire rice co2 update en n-1 puis a chaque microdt, mettre a jour riceco2 305 363 364 !============================================================= 365 ! Bonus: additional meteoritic particles for nucleation 366 !============================================================= 367 if (meteo_flux) then 368 !altitude_meteo(130) 369 !pzlev(ngrid,nlay) 370 !meteo(130,100) 371 !resultat: meteo_ccn(ngrid,nlay,100) 372 do l=1,nlay 373 do ig=1,ngrid 374 ltemp1=abs(altitudes_meteo(:)-pzlev(ig,l)) 375 ltemp2=abs(altitudes_meteo(:)-pzlev(ig,l+1)) 376 lebon1=minloc(ltemp1,DIM=1) 377 lebon2=minloc(ltemp2,DIM=1) 378 nelem=lebon2-lebon1+1. 379 mtemp(:)=0d0 !mtemp(100) : valeurs pour les 100bins 380 do ibin=1,100 381 mtemp(ibin)=sum(meteo(lebon1:lebon2,ibin)) 382 enddo 383 meteo_ccn(ig,l,:)=mtemp(:)/nelem 384 enddo 385 enddo 386 endif 387 306 388 c Main loop over the GCM's grid 307 389 DO l=1,nlay 308 390 DO ig=1,ngrid 309 310 311 391 c Get the partial pressure of co2 vapor and its saturation ratio 312 392 pco2 = zq(ig,l,igcm_co2) * (mmean(ig,l)/44.01) * pplay(ig,l) … … 318 398 rho_ice_co2T(ig,l)=1000.*(1.72391-2.53e-4*zt(ig,l) 319 399 & -2.87e-6*zt(ig,l)*zt(ig,l)) 320 vo2co2 = m0co2 / dble(rho_ice_co2T(ig,l)) ! m0co2400 vo2co2 = m0co2 / dble(rho_ice_co2T(ig,l)) 321 401 rho_ice_co2=rho_ice_co2T(ig,l) 322 c call updaterccn(zq(ig,l,igcm_dust_mass), 323 c & zq(ig,l,igcm_dust_number),rdust(ig,l),tauscaling(ig)) 324 c write(*,*) "l, pco2, satu= ",l,pco2,satu 402 325 403 IF ( satu .ge. 1d0 ) THEN ! if there is condensation 326 c write(*,*)327 328 !write(*,*) "Zt, rho=",zt(ig,l),rho_ice_co2329 c Masse_atm=mmean(ig,l)*1.e-3*pplay(ig,l)/rgp/zt(ig,l) !Kg par couche330 331 332 c call updaterccn(zq(ig,l,igcm_dust_mass),333 c & zq(ig,l,igcm_dust_number),rdust(ig,l),tauscaling(ig))334 335 c call updaterccn(zq(ig,l,igcm_dust_mass),336 c & zq(ig,l,igcm_dust_number),rdust(ig,l),tauscaling(ig))337 c write(*,*) "Improved, l,Rdust = ",l,rdust(ig,l)338 404 339 405 rdust(ig,l)= zq(ig,l,igcm_dust_mass) … … 341 407 & / zq(ig,l,igcm_dust_number) 342 408 rdust(ig,l)= rdust(ig,l)**(1./3.) 343 !write(*,*) "Improved2, l,Rdust = ",l,rdust(ig,l)344 rdust(ig,l)=max(1.e- 10,rdust(ig,l))345 rdust(ig,l)=min(5.e-5,rdust(ig,l))409 c write(*,*) "Improved2, l,Rdust = ",l,rdust(ig,l) 410 rdust(ig,l)=max(1.e-9,rdust(ig,l)) 411 c rdust(ig,l)=min(1.e-5,rdust(ig,l)) 346 412 ! write(*,*) "Improved3,Rdust = ",rdust(ig,l) 347 413 … … 349 415 Mo = zq(ig,l,igcm_dust_mass)* tauscaling(ig)+1.e-30 350 416 No = zq(ig,l,igcm_dust_number)* tauscaling(ig)+1.e-30 351 c write(*,*) "Improved dust number, mass = ",352 c & zq(ig,l,igcm_dust_number)* tauscaling(ig),353 c & zq(ig,l,igcm_dust_mass)* tauscaling(ig)354 c write(*,*) "No, Mo = ",No, Mo355 417 Rn = rdust(ig,l) 356 418 Rn = -log(Rn) … … 364 426 n_derf = derf((rb_cldco2(i+1)+Rn) *dev2) 365 427 m_derf = derf((rb_cldco2(i+1)+Rm) *dev2) 366 n_aer(i) = n_aer(i) + 0.5 * No * n_derf 428 n_aer(i) = n_aer(i) + 0.5 * No * n_derf + 429 & meteo_ccn(ig,l,i) !Ajout meteo_ccn 367 430 m_aer(i) = m_aer(i) + 0.5 * Mo * m_derf 368 c write(*,*) "i, rad_cldco2(i) = ",i, rad_cldco2(i), 369 c & n_aer(i) 431 m_aer2(i)=4./3.*pi*rho_dust 432 & *n_aer(i)*rad_cldco2(i)*rad_cldco2(i) 433 & *rad_cldco2(i) 434 c write(*,*) "diff =",rad_cldco2(i),m_aer(i),m_aer2(i) 370 435 enddo 371 372 373 sumcheck = 0 374 do i = 1, nbinco2_cld 375 sumcheck = sumcheck + n_aer(i) 376 enddo 377 sumcheck = abs(sumcheck/No - 1) 378 if ((sumcheck .gt. 1e-5).and. (1./Rn .gt. rmin_cld)) then 379 print*, "WARNING, No sumcheck PROBLEM" 380 print*, "sumcheck, No",sumcheck, No 381 print*, "rdust =",rdust(ig,l) 382 print*, "min radius, Rn, ig, l", rmin_cld, 1./Rn, ig, l 383 print*, "Dust binned distribution", n_aer 384 STOP 385 endif 386 387 sumcheck = 0 388 do i = 1, nbinco2_cld 389 sumcheck = sumcheck + m_aer(i) 390 enddo 391 sumcheck = abs(sumcheck/Mo - 1) 392 if ((sumcheck .gt. 1e-5) .and. (1./Rn .gt. rmin_cld)) 393 & then 394 print*, "WARNING, Mo sumcheck PROBLEM" 395 print*, "sumcheck, Mo",sumcheck, Mo 396 print*, "min radius, Rm, ig, l", rmin_cld, 1./Rm, ig,l 397 print*, "Dust binned distribution", m_aer 398 STOP 399 endif 400 401 call updaterice_micro( 402 & zq(ig,l,igcm_h2o_ice), ! ice mass 403 & zq(ig,l,igcm_ccn_mass), ! ccn mass 404 & zq(ig,l,igcm_ccn_number), ! ccn number 405 & tauscaling(ig),rice(ig,l),rhocloud(ig,l)) 436 c write(*,*) "Bilan =",sum(m_aer),sum(m_aer2) 437 438 c sumcheck = 0 439 c do i = 1, nbinco2_cld 440 c sumcheck = sumcheck + n_aer(i) 441 c enddo 442 c sumcheck = abs(sumcheck/No - 1) 443 c if ((sumcheck .gt. 1e-5).and. (1./Rn .gt. rmin_cld)) then 444 c print*, "WARNING, No sumcheck PROBLEM" 445 c print*, "sumcheck, No, rdust",sumcheck, No,rdust(ig,l) 446 c print*, "min radius, Rn, ig, l", rmin_cld, 1./Rn, ig, l 447 c print*, "Dust binned distribution", n_aer 448 c STOP 449 c endif 450 451 c sumcheck = 0 452 c do i = 1, nbinco2_cld 453 c sumcheck = sumcheck + m_aer(i) 454 c enddo 455 c sumcheck = abs(sumcheck/Mo - 1) 456 c if ((sumcheck .gt. 1e-5) .and. (1./Rn .gt. rmin_cld)) 457 c & then 458 c print*, "WARNING, Mo sumcheck PROBLEM" 459 c print*, "sumcheck, Mo",sumcheck, Mo 460 c print*, "min radius, Rm, ig, l", rmin_cld, 1./Rm, ig,l 461 c print*, "Dust binned distribution", m_aer 462 c STOP 463 c endif 464 m_aer(:)=m_aer2(:) 465 466 467 rccnh2o(ig,l)= zq(ig,l,igcm_ccn_mass) 468 & *0.75/pi/rho_dust 469 & / zq(ig,l,igcm_ccn_number) 470 471 rice(ig,l)=( zq(ig,l,igcm_h2o_ice)*3.0/ 472 & (4.0*rho_ice_co2*zq(ig,l,igcm_ccn_number) 473 & *pi*tauscaling(ig)) +rccnh2o(ig,l)*rccnh2o(ig,l) 474 & *rccnh2o(ig,l) )**(1.0/3.0) 475 rhocloud(ig,l)=( zq(ig,l,igcm_h2o_ice)*rho_ice 476 & +zq(ig,l,igcm_ccn_mass) *tauscaling(ig)*rho_dust) 477 & / (zq(ig,l,igcm_h2o_ice)+zq(ig,l,igcm_ccn_mass) 478 & *tauscaling(ig)) 479 c call updaterice_micro( 480 c & zq(ig,l,igcm_h2o_ice), ! ice mass 481 c & zq(ig,l,igcm_ccn_mass), ! ccn mass 482 c & zq(ig,l,igcm_ccn_number), ! ccn number 483 c & tauscaling(ig),rice(ig,l),rhocloud(ig,l)) 406 484 ! rice radius of CCN + H20 crystal 407 !write(*,*) "Improved1 Rice=",rice(ig,l)485 c write(*,*) "Improved1 Rice=",rice(ig,l) 408 486 rice(ig,l)=max(1.e-10,rice(ig,l)) 409 rice(ig,l)=min(5.e-5,rice(ig,l))487 c rice(ig,l)=min(1.e-5,rice(ig,l)) 410 488 !write(*,*) "Improved2 Rice=",rice(ig,l) 411 489 Mo = zq(ig,l,igcm_h2o_ice) + 412 & zq(ig,l,igcm_ccn_mass)*tauscaling(ig) 413 & + 1.e-30 !Total mass of H20 crystals,CCN included490 & zq(ig,l,igcm_ccn_mass)*tauscaling(ig)+1.e-30!*tauscaling(ig) 491 ! & + 1.e-30 !Total mass of H20 crystals,CCN included 414 492 No = zq(ig,l,igcm_ccn_number)* tauscaling(ig) + 1.e-30 415 493 Rn = rice(ig,l) 416 494 Rn = -log(Rn) 417 Rm = Rn - 3. * sigma_ice co2*sigma_iceco2418 n_derf = erf( (rb_cldco2(1)+Rn) *dev 2)419 m_derf = erf( (rb_cldco2(1)+Rm) *dev 2)495 Rm = Rn - 3. * sigma_ice*sigma_ice 496 n_derf = erf( (rb_cldco2(1)+Rn) *dev3) 497 m_derf = erf( (rb_cldco2(1)+Rm) *dev3) 420 498 do i = 1, nbinco2_cld 421 n_aer_h2oice(i) = -0.5 * No * n_derf !! this ith previously computed422 m_aer_h2oice(i) = -0.5 * Mo * m_derf !! this ith previously computed423 n_derf = derf( (rb_cldco2(i+1)+Rn) *dev 2)424 m_derf = derf( (rb_cldco2(i+1)+Rm) *dev 2)425 n_aer_h2oice(i) = n_aer_h2oice(i) + 0.5 * No * n_derf ! vector not really needed - temp var - listo426 m_aer_h2oice(i) = m_aer_h2oice(i) + 0.5 * Mo * m_derf ! vector not really needed - temp var499 n_aer_h2oice(i) = -0.5 * No * n_derf 500 m_aer_h2oice(i) = -0.5 * Mo * m_derf 501 n_derf = derf( (rb_cldco2(i+1)+Rn) *dev3) 502 m_derf = derf( (rb_cldco2(i+1)+Rm) *dev3) 503 n_aer_h2oice(i) = n_aer_h2oice(i) + 0.5 * No * n_derf 504 m_aer_h2oice(i) = m_aer_h2oice(i) + 0.5 * Mo * m_derf 427 505 rad_h2oice(i) = rad_cldco2(i) 506 m_aer_h2oice2(i)=4./3.*pi*rhocloud(ig,l) 507 & *n_aer_h2oice(i)*rad_h2oice(i)*rad_h2oice(i) 508 & *rad_h2oice(i) 509 510 428 511 c write(*,*) "before nuc, i,rad_h2o(i)= ",i,rad_cldco2(i) 429 512 c & ,m_aer_h2oice(i),n_aer_h2oice(i) 430 513 enddo 431 sumcheck = 0 432 do i = 1, nbinco2_cld 433 sumcheck = sumcheck + n_aer_h2oice(i) 434 enddo 435 sumcheck = abs(sumcheck/No - 1) 436 if ((sumcheck .gt. 1e-5).and. (1./Rn .gt. rmin_cld)) then 437 print*, "WARNING, Noh2o sumcheck PROBLEM" 438 print*, "sumcheck, No",sumcheck, No 439 print*, "rice =",rice(ig,l) 440 print*, "min radius, Rn, ig, l", rmin_cld, 1./Rn, ig, l 441 print*, "Dust binned distribution", n_aer_h2oice 442 STOP 443 endif 444 445 514 c sumcheck = 0 515 c do i = 1, nbinco2_cld 516 c sumcheck = sumcheck + n_aer_h2oice(i) 517 c enddo 518 c sumcheck = abs(sumcheck/No - 1) 519 c if ((sumcheck .gt. 1e-5).and. (1./Rn .gt. rmin_cld)) then 520 c print*, "WARNING, Noh2o sumcheck PROBLEM" 521 c print*, "sumcheck, No, rice",sumcheck, No,rice(ig,l) 522 c print*, "min radius, Rn, ig, l", rmin_cld, 1./Rn, ig, l 523 c print*, "Dust binned distribution", n_aer_h2oice 524 c STOP 525 c endif 526 527 c sumcheck = 0 528 c do i = 1, nbinco2_cld 529 c sumcheck = sumcheck + m_aer_h2oice(i) 530 c enddo 531 c sumcheck = abs(sumcheck/Mo - 1) 532 c if ((sumcheck .gt. 1e-5) .and. (1./Rn .gt. rmin_cld)) 533 c & then 534 c print*, "WARNING, Moh2o sumcheck PROBLEM" 535 c print*, "sumcheck, Mo",sumcheck, Mo 536 c print*, "min radius, Rm, ig, l", rmin_cld, 1./Rm, ig,l 537 c print*, "Dust binned distribution", m_aer_h2oice 538 c STOP 539 c endif 446 540 c Get the rates of nucleation 541 447 542 call nucleaCO2(dble(pco2),zt(ig,l),dble(satu) 448 543 & ,n_aer,rate,n_aer_h2oice 449 544 & ,rad_h2oice,rateh2o,vo2co2) 450 ! regarder rateh20, et mettre = 0 si non nul pour le moment545 m_aer_h2oice(:)=m_aer_h2oice2(:) 451 546 dN = 0. 452 547 dM = 0. … … 455 550 do i = 1, nbinco2_cld 456 551 Proba=1.0-dexp(-1.*ptimestep*rate(i)) 457 Probah2o=1.0-dexp(-1.*ptimestep*rateh2o(i)) 458 552 Probah2o=coeffh2o*(1.0-dexp(-1.*ptimestep*rateh2o(i))) 459 553 dNh2o = dNh2o + n_aer_h2oice(i) * Probah2o 460 554 dMh2o = dMh2o + m_aer_h2oice(i) * Probah2o 461 462 555 dN = dN + n_aer(i) * Proba 463 556 dM = dM + m_aer(i) * Proba 464 c write(*,*) "i, dNi, dN= ",i,n_aer(i)*Proba,dN 557 465 558 enddo 466 467 559 ! dM masse activée (kg) et dN nb particules par kg d'air 468 469 c write(*,*) " nuclea dM = ",dM/tauscaling(ig),470 c & " nuclea dN = ", dN/tauscaling(ig)471 560 472 561 dNN= dN/tauscaling(ig) 473 562 dMM= dM/tauscaling(ig) 563 dNN=min(dNN,zq(ig,l,igcm_dust_number)) 564 dMM=min(dMM,zq(ig,l,igcm_dust_mass)) 565 ! if (dNN .gt. 0) then 566 ! write(*,*) "Nuclea dNN crees=",dNN 567 ! write(*,*) "Nuclea dMM crees=",dMM 568 ! endif 474 569 dNNh2o=dNh2o/tauscaling(ig) 475 dMMh2o=dMh2o/tauscaling(ig) 476 477 dNN=min(dNN,abs(zq(ig,l,igcm_dust_number))) 478 dMM=min(dMM,abs(zq(ig,l,igcm_dust_mass))) 479 c 480 c write(*,*) "Nuclea dNN crees=",dNN 481 dNNh2o=min(dNNh2o,abs(zq(ig,l,igcm_ccn_number))) 482 dMMh2o=min(dMMh2o,abs(zq(ig,l,igcm_h2o_ice)/tauscaling(ig) 483 & +zq(ig,l,igcm_ccn_mass))) !Total mass of H2O crystals available 484 485 c write(*,*) "Nuclea dNNh2o crees=",dNNh2o 486 487 c Update CCNs for CO2 crystals 488 ! WARNING dM dMh2o, interaction nuages eau-co2 -- h20 set to 0 for now 570 dNNh2o=min(dNNh2o,zq(ig,l,igcm_ccn_number)) 571 572 ratioh2o_ccn=1./(zq(ig,l,igcm_h2o_ice) 573 & +zq(ig,l,igcm_ccn_mass)*tauscaling(ig)) 574 575 dMh2o_ice=dMh2o*zq(ig,l,igcm_h2o_ice)*ratioh2o_ccn 576 dMh2o_ccn=dMh2o*zq(ig,l,igcm_ccn_mass)* 577 & tauscaling(ig)*ratioh2o_ccn 578 579 dMh2o_ccn=dMh2o_ccn/tauscaling(ig) 580 dMh2o_ccn=min(dMh2o_ccn,zq(ig,l,igcm_ccn_mass)) 581 dMh2o_ice=min(dMh2o_ice,zq(ig,l,igcm_h2o_ice)) 582 489 583 zq(ig,l,igcm_ccnco2_mass) = 490 584 & zq(ig,l,igcm_ccnco2_mass) + dMM 491 zq(ig,l,igcm_ccnco2_number) = 585 zq(ig,l,igcm_ccnco2_number) = 492 586 & zq(ig,l,igcm_ccnco2_number) + dNN 493 zq(ig,l,igcm_dust_mass) = 494 & zq(ig,l,igcm_dust_mass) - dMM 495 zq(ig,l,igcm_dust_number) = 496 & zq(ig,l,igcm_dust_number) - dNN 587 588 zq(ig,l,igcm_dust_mass)= zq(ig,l,igcm_dust_mass)-dMM 589 zq(ig,l,igcm_dust_number)=zq(ig,l,igcm_dust_number)-dNN 497 590 498 591 c Update CCN for CO2 nucleating on H2O CCN : 499 592 ! Warning: must keep memory of it 500 593 zq(ig,l,igcm_ccnco2_mass) = 501 & zq(ig,l,igcm_ccnco2_mass) + dMMh2o594 & zq(ig,l,igcm_ccnco2_mass) + dMh2o_ice+dMh2o_ccn 502 595 zq(ig,l,igcm_ccnco2_number) = 503 596 & zq(ig,l,igcm_ccnco2_number) + dNNh2o 504 597 505 506 zq(ig,l,igcm_ccn_number) = 507 & zq(ig,l,igcm_ccn_number) - dNNh2o 508 509 ratioh2o_ccn=1./(zq(ig,l,igcm_h2o_ice) 510 & +zq(ig,l,igcm_ccn_mass)*tauscaling(ig)) 511 512 513 memdMMh2o(ig,l)= memdMMh2o(ig,l)+zq(ig,l,igcm_h2o_ice)* 514 & dMMh2o*ratioh2o_ccn 515 memdMMccn(ig,l)= memdMMccn(ig,l)+zq(ig,l,igcm_ccn_mass)* 516 & tauscaling(ig)*dMMh2o*ratioh2o_ccn 598 zq(ig,l,igcm_ccn_number)=zq(ig,l,igcm_ccn_number)-dNNh2o 599 zq(ig,l,igcm_h2o_ice) = zq(ig,l,igcm_h2o_ice)-dMh2o_ice 600 zq(ig,l,igcm_ccn_mass)= zq(ig,l,igcm_ccn_mass)-dMh2o_ccn 601 602 603 memdMMh2o(ig,l)=memdMMh2o(ig,l)+dMh2o_ice 604 memdMMccn(ig,l)=memdMMccn(ig,l)+dMh2o_ccn 517 605 memdNNccn(ig,l)=memdNNccn(ig,l)+dNNh2o 518 c if (dMMh2o .gt. 0) then519 c write(*,*) 'test h2o'520 c write(*,*) "dMMh2o=",dMMh2o521 c write(*,*) "2 =",zq(ig,l,igcm_ccn_mass)*tauscaling(ig)*522 c & dMMh2o*ratioh2o_ccn+zq(ig,l,igcm_h2o_ice)*523 c & dMMh2o*ratioh2o_ccn524 c write(*,*) "3=",zq(ig,l,igcm_ccn_mass)*tauscaling(ig)*525 c & dMMh2o*ratioh2o_ccn526 c write(*,*) "4=",zq(ig,l,igcm_h2o_ice)*527 c & dMMh2o*ratioh2o_ccn528 c write(*,*) "tauscaling=",tauscaling(ig)529 c endif530 zq(ig,l,igcm_h2o_ice) = zq(ig,l,igcm_h2o_ice)*531 & (1.-dMMh2o*ratioh2o_ccn)532 zq(ig,l,igcm_ccn_mass) = zq(ig,l,igcm_ccn_mass)*533 & tauscaling(ig)*(1.-dMMh2o*ratioh2o_ccn)534 535 536 606 ENDIF ! of is satu >1 537 607 !============================================================= … … 542 612 c Indeed, if we are supersaturated and still don't have at least one nuclei, we should better wait 543 613 c to avoid unrealistic value for nuclei radius and so on for cases that remain negligible. 544 c IF ( zq(ig,l,igcm_ccnco2_number)*tauscaling(ig).ge. 1.0) THEN 545 546 IF (zq(ig,l,igcm_ccnco2_number)*tauscaling(ig) .ge. 1.) 547 & THEN 614 IF (zq(ig,l,igcm_ccnco2_number)*tauscaling(ig) .ge. 1)THEN 548 615 ! we trigger crystal growth 549 616 c 617 c write(*,*) "ccn number mass=",zq(ig,l,igcm_ccnco2_number), 618 c & zq(ig,l,igcm_ccnco2_mass) 550 619 551 620 c Niceco2 = zq(ig,l,igcm_co2_ice) … … 561 630 rdust(ig,l)= rdust(ig,l)**(1./3.) 562 631 rdust(ig,l)=max(1.e-10,rdust(ig,l)) 563 !rdust(ig,l)=min(5.e-6,rdust(ig,l))564 632 c rdust(ig,l)=min(5.e-6,rdust(ig,l)) 633 565 634 riceco2(ig,l)=( zq(ig,l,igcm_co2_ice)*3.0/ 566 635 & (4.0*rho_ice_co2*zq(ig,l,igcm_ccnco2_number) … … 581 650 c call updaterice_microco2(Niceco2,Qccnco2,Nccnco2, 582 651 c & tauscaling(ig),riceco2(ig,l),rhocloudco2(ig,l)) 583 !write(*,*) "Riceco2 updatebefore growth = ",riceco2(ig,l)652 c write(*,*) "Riceco2 before growth = ",riceco2(ig,l) 584 653 c write(*,*) "rdust before growth = ",rdust(ig,l) 585 c write(*,*) "co2 before growth=",zq(ig,l,igcm_co2) 654 c write(*,*) "co2ice before growth=",zq(ig,l,igcm_co2_ice) 655 c write(*,*) "co2 before growth=",zq(ig,l,igcm_co2) 586 656 c write(*,*) "pplay before growth=",pplay(ig,l) 587 657 c write(*,*) "zt before growth =",zt(ig,l) 588 !write(*,*) "Satu before growth=",satu658 c write(*,*) "Satu before growth=",satu 589 659 c riceco2(ig,l)=max(riceco2(ig,l),rdust(ig,l)) 590 660 No = zq(ig,l,igcm_ccnco2_number)* tauscaling(ig)+1.e-30 … … 601 671 & satu,riceco2(ig,l),mmean(ig,l),Ic_rice) ! Mass transfer rate (kg/s) for a rice particle 602 672 ! Ic_rice mass flux kg.s-1 <0 si croissance ! 603 drsurdt=-1.0/(4.0*pi*riceco2(ig,l)* 604 & riceco2(ig,l)*rho_ice_co2)*Ic_rice 673 if (isnan(Ic_rice)) then 674 Ic_rice=0. 675 flag_pourri=1 676 endif 677 c drsurdt=-1.0/(4.0*pi*riceco2(ig,l)* 678 c & riceco2(ig,l)*rho_ice_co2)*Ic_rice 605 679 dMice = No * Ic_rice*ptimestep ! Kg par kg d'air, <0 si croissance ! 606 680 c write(*,*) "dMicev0 in improved = " , dMice … … 611 685 dMice =-1.* min(abs(dMice),abs(zq(ig,l,igcm_co2))) 612 686 endif 613 riceco2(ig,l)=riceco2(ig,l)+drsurdt*ptimestep687 c riceco2(ig,l)=riceco2(ig,l)+drsurdt*ptimestep 614 688 c write(*,*) "riceco2+dr/dt = ", riceco2(ig,l) 615 c write(*,*) "dMice in improved = " , dMice 616 617 zq(ig,l,igcm_co2_ice) = zq(ig,l,igcm_co2_ice) 618 & -dMice 689 !write(*,*) "dMice in improved = " , dMice 690 691 zq(ig,l,igcm_co2_ice) = zq(ig,l,igcm_co2_ice)-dMice 619 692 zq(ig,l,igcm_co2) = zq(ig,l,igcm_co2)+dMice 620 c write(*,*) "Improved zq co2 ice = ", zq(ig,l,igcm_co2_ice) 621 ! countcells = countcells + 1 622 623 c riceco2(ig,l)=( zq(ig,l,igcm_co2_ice)*3.0/ 624 c & (4.0*rho_ice_co2*pi*zq(ig,l,igcm_ccnco2_number) 625 c & *tauscaling(ig)) +rdust(ig,l)*rdust(ig,l) 626 c & *rdust(ig,l) )**(1.0/3.0) 627 c write(*,*) "Improved new riceco2 = ",riceco2(ig,l) 628 629 c write(*,*) "new riceco2 improvedupdaterad= ",riceco2(ig,l) 630 631 ! latent heat release 632 633 l0=595594. 634 l1=903.111 635 l2=-11.5959 636 l3=0.0528288 637 l4=-0.000103183 693 694 ! latent heat release >0 if growth i.e. if dMice <0 695 638 696 639 697 lw = l0 + l1 * zt(ig,l) + l2 * zt(ig,l)**2 + 640 698 & l3 * zt(ig,l)**3 + l4 * zt(ig,l)**4 !J.kg-1 641 699 c write(*,*) "CPP= ",cpp ! = 744.5 642 643 pdtcloudco2(ig,l)= dMice*lw/cpp/ptimestep ! kg par couche * J par kg /J par K / s = K par seconde 644 645 c write(*,*) "pdtcloudco2 after growth = ",pdtcloudco2(ig,l) 646 700 !Peut etre un probleme de signe ici! 701 pdtcloudco2(ig,l)= -1.*dMice*lw/cpp/ptimestep ! kg par couche * J par kg /J par K / s = K par seconde 702 703 !write(*,*) "pdtcloudco2 after growth = ",pdtcloudco2(ig,l) 704 705 !write(*,*) "co2 after growth = ",zq(ig,l,igcm_co2) 706 !write(*,*) "co2_ice after growth = ",zq(ig,l,igcm_co2_ice) 647 707 648 708 !deltaT par condens/subli. qui remplace le dT du CO2 du newcondens pré-constantino … … 652 712 !============================================================= 653 713 654 c If all the ice particles sublimate, all the condensation 655 c nuclei are released: 656 657 c !!! with CO2 ice nuclei: dust/H2O nuclei are not released because 658 c they were not subtracted to dust_number 659 c Their counter is just set to "0". 660 c (see end of section 3.) : On peut les enlever à dust 661 662 c interaction ho-co2 ici, dans la mise a jour des traceurs WARNING reflechir 663 !! Niceco2,Qccnco2,Nccnco2 664 665 666 714 c if (abs(dMice) .ge. 1.e-10) then 715 716 c write(*,*) "DIAG zq pdt",(zq(ig,l,igcm_co2_ice)- 717 c & zq0(ig,l,igcm_co2_ice))/ptimestep,pdtcloudco2(ig,l) 718 c endif 719 ENDIF ! of if NCCN > 1 720 667 721 rdust(ig,l)= zq(ig,l,igcm_ccnco2_mass) 668 722 & *0.75/pi/rho_dust 669 723 & / zq(ig,l,igcm_ccnco2_number) 670 724 rdust(ig,l)= rdust(ig,l)**(1./3.) 671 rdust(ig,l)=max(1.e- 10,rdust(ig,l))672 !rdust(ig,l)=min(5.e-6,rdust(ig,l))725 rdust(ig,l)=max(1.e-9,rdust(ig,l)) 726 c rdust(ig,l)=min(5.e-6,rdust(ig,l)) 673 727 674 728 riceco2(ig,l)=( zq(ig,l,igcm_co2_ice)*3.0/ … … 682 736 c call updaterice_microCO2(Niceco2,Qccnco2,Nccnco2, 683 737 c & tauscaling(ig),riceco2(ig,l),rhocloudco2(ig,l)) 684 685 if ((zq(ig,l,igcm_co2_ice).le. 1.e-23) 686 & .or.(riceco2(ig,l) .le. rdust(ig,l))) then 687 738 739 c If there is no more co2 ice, all the ice particles sublimate, all the condensation nuclei are released: 740 741 if ((zq(ig,l,igcm_co2_ice).le. 1.e-25)) then 742 c & .or. flag_pourri .eq. 1 743 c & .or.(riceco2(ig,l) .le. rdust(ig,l)) 744 c & .or. (l .ge.1 .and. l .le. 5) 745 c & .or. (zq(ig,l,igcm_co2_ice) .ge. 0.1) 746 lw = l0 + l1 * zt(ig,l) + l2 * zt(ig,l)**2 + 747 & l3 * zt(ig,l)**3 + l4 * zt(ig,l)**4 !J.kg-1 748 c write(*,*) "CPP= ",cpp ! = 744.5 749 750 pdtcloudco2(ig,l)=pdtcloudco2(ig,l)-1. 751 & *zq(ig,l,igcm_co2_ice)*lw/cpp/ptimestep ! 752 !On sublime tout ! 688 753 c write(*,*) "Riceco2 improved before reset=",riceco2(ig,l) 689 754 c write(*,*) "Niceco2 improved before reset=", … … 705 770 endif 706 771 707 if (zq(ig,l,igcm_ccnco2_mass) .gt. 1.e-30) then772 c if (zq(ig,l,igcm_ccnco2_mass) .gt. 1.e-30) then 708 773 zq(ig,l,igcm_dust_mass) = 709 774 & zq(ig,l,igcm_dust_mass) 710 775 & + zq(ig,l,igcm_ccnco2_mass)- 711 776 & (memdMMh2o(ig,l)+memdMMccn(ig,l)) 712 endif713 if (zq(ig,l,igcm_ccnco2_number) .gt. 1.e-30) then777 c endif 778 c if (zq(ig,l,igcm_ccnco2_number) .gt. 1.e-30) then 714 779 zq(ig,l,igcm_dust_number) = 715 780 & zq(ig,l,igcm_dust_number) 716 781 & + zq(ig,l,igcm_ccnco2_number)-memdNNccn(ig,l) 717 endif782 c endif 718 783 719 if (zq(ig,l,igcm_co2_ice) .gt. 1.e-30) then784 c if (zq(ig,l,igcm_co2_ice) .gt. 1.e-30) then 720 785 zq(ig,l,igcm_co2) = zq(ig,l,igcm_co2) 721 786 & + zq(ig,l,igcm_co2_ice) 722 endif787 c endif 723 788 724 789 zq(ig,l,igcm_ccnco2_mass)=0. … … 729 794 memdMMccn(ig,l)=0. 730 795 riceco2(ig,l)=0. 731 pdtcloudco2(ig,l)=0. 796 732 797 endif 733 734 ENDIF ! of if NCCN > 1 735 ENDDO ! of ig loop 736 ENDDO ! of nlayer loop 737 738 739 ! Get cloud tendencies 798 ENDDO ! of ig loop 799 ENDDO ! of nlayer loop 800 ! Get cloud tendencies 740 801 pdqcloudco2(1:ngrid,1:nlay,igcm_co2) = 741 802 & (zq(1:ngrid,1:nlay,igcm_co2) - 742 803 & zq0(1:ngrid,1:nlay,igcm_co2))/ptimestep 743 744 804 pdqcloudco2(1:ngrid,1:nlay,igcm_co2_ice) = 745 805 & (zq(1:ngrid,1:nlay,igcm_co2_ice) - 746 806 & zq0(1:ngrid,1:nlay,igcm_co2_ice))/ptimestep 747 748 807 pdqcloudco2(1:ngrid,1:nlay,igcm_h2o_ice) = 749 808 & (zq(1:ngrid,1:nlay,igcm_h2o_ice) - 750 809 & zq0(1:ngrid,1:nlay,igcm_h2o_ice))/ptimestep 751 752 810 pdqcloudco2(1:ngrid,1:nlay,igcm_ccn_mass) = 753 811 & (zq(1:ngrid,1:nlay,igcm_ccn_mass) - 754 812 & zq0(1:ngrid,1:nlay,igcm_ccn_mass))/ptimestep 755 756 813 pdqcloudco2(1:ngrid,1:nlay,igcm_ccn_number) = 757 814 & (zq(1:ngrid,1:nlay,igcm_ccn_number) - 758 815 & zq0(1:ngrid,1:nlay,igcm_ccn_number))/ptimestep 759 760 816 pdqcloudco2(1:ngrid,1:nlay,igcm_ccnco2_mass) = 761 817 & (zq(1:ngrid,1:nlay,igcm_ccnco2_mass) - 762 818 & zq0(1:ngrid,1:nlay,igcm_ccnco2_mass))/ptimestep 763 764 819 pdqcloudco2(1:ngrid,1:nlay,igcm_ccnco2_number) = 765 820 & (zq(1:ngrid,1:nlay,igcm_ccnco2_number) - 766 821 & zq0(1:ngrid,1:nlay,igcm_ccnco2_number))/ptimestep 767 768 769 c if (scavenging) then 770 771 pdqcloudco2(1:ngrid,1:nlay,igcm_dust_mass) = 772 & (zq(1:ngrid,1:nlay,igcm_dust_mass) - 773 & zq0(1:ngrid,1:nlay,igcm_dust_mass))/ptimestep 774 775 pdqcloudco2(1:ngrid,1:nlay,igcm_dust_number) = 776 & (zq(1:ngrid,1:nlay,igcm_dust_number) - 777 & zq0(1:ngrid,1:nlay,igcm_dust_number))/ptimestep 778 c endif 779 822 pdqcloudco2(1:ngrid,1:nlay,igcm_dust_mass) = 823 & (zq(1:ngrid,1:nlay,igcm_dust_mass) - 824 & zq0(1:ngrid,1:nlay,igcm_dust_mass))/ptimestep 825 pdqcloudco2(1:ngrid,1:nlay,igcm_dust_number) = 826 & (zq(1:ngrid,1:nlay,igcm_dust_number) - 827 & zq0(1:ngrid,1:nlay,igcm_dust_number))/ptimestep 780 828 return 781 829 end -
trunk/LMDZ.MARS/libf/phymars/initracer.F
r1685 r1720 320 320 count=count+1 321 321 endif 322 if ( microphysco2) then322 if (co2clouds) then 323 323 if (noms(iq).eq."ccnco2_mass") then 324 324 igcm_ccnco2_mass=iq … … 722 722 & "a ho2 tracer !" 723 723 stop 724 724 endif 725 725 if (igcm_h2o2.eq.0) then 726 726 write(*,*) "initracer: error !!" -
trunk/LMDZ.MARS/libf/phymars/massflowrateCO2.F
r1685 r1720 52 52 53 53 54 Tcm = 0!dble(T) ! initialize pourquoi 0 et pas t(i)54 Tcm = 80!dble(T) ! initialize pourquoi 0 et pas t(i) 55 55 56 56 T_inf = 0d0 57 57 T_sup = 800d0 58 58 59 T_dT = 0. 1 ! precision - mettre petit et limiter nb iteration?59 T_dT = 0.01 ! precision - mettre petit et limiter nb iteration? 60 60 61 61 666 CONTINUE … … 114 114 EXTERNAL funcd 115 115 116 PARAMETER (MAXIT= 10000) !Maximum allowed number of iterations. Using a combination of Newton-Raphson and bisection,116 PARAMETER (MAXIT=50000) !Maximum allowed number of iterations. Using a combination of Newton-Raphson and bisection, 117 117 !find the root of a function bracketed between x1 and x2. The root, returned as the function value rtsafe, 118 118 !will be refined until its accuracy is known within !!±xacc. funcd is a user-supplied subroutine which -
trunk/LMDZ.MARS/libf/phymars/microphys.h
r1685 r1720 76 76 ! bachnar 2016 value :0.78 77 77 !old value 0.952 78 REAL, parameter :: mtetaco2 = 0. 7878 REAL, parameter :: mtetaco2 = 0.952 79 79 ! Volume of a co2 molecule (m3) 80 80 DOUBLE PRECISION :: vo1co2 -
trunk/LMDZ.MARS/libf/phymars/nucleaCO2.F
r1685 r1720 65 65 mtetalocal = dble(mtetaco2) !! use mtetalocal for better performance 66 66 mtetalocalh=dble(mteta) 67 68 cccccccccccccccccccccccccccccccccccccccccccccccccc69 ccccccccccc ESSAIS TN MTETA = F (T) cccccccccccccc70 c if (temp .gt. 200) then71 c mtetalocal = mtetalocal72 c else if (temp .lt. 190) then73 c mtetalocal = mtetalocal-0.0574 c else75 c mtetalocal = mtetalocal - (190-temp)*0.00576 c endif77 c----------------exp law, see Trainer 2008, J. Phys. Chem. C 2009, 113, 2036\u2013204078 !mtetalocal = max(mtetalocal - 6005*exp(-0.065*temp),0.1)79 !mtetalocal = max(mtetalocal - 6005*exp(-0.068*temp),0.1)80 !print*, mtetalocal, temp81 cccccccccccccccccccccccccccccccccccccccccccccccccc82 cccccccccccccccccccccccccccccccccccccccccccccccccc83 c IF (firstcall) THEN84 c print*, ' '85 c print*, 'dear user, please keep in mind that'86 c print*, 'contact parameter IS constant'87 !print*, 'contact parameter IS NOT constant:'88 !print*, 'max(mteta - 6005*exp(-0.065*temp),0.1)'89 !print*, 'max(mteta - 6005*exp(-0.068*temp),0.1)'90 c print*, ' '91 c firstcall=.false.92 c END IF93 cccccccccccccccccccccccccccccccccccccccccccccccccc94 cccccccccccccccccccccccccccccccccccccccccccccccccc95 96 c write(*,*) "IN nuc, SAT = ",sat97 c write(*,*) "IN nuc, mtetalocal = ",mtetalocal98 67 99 68 -
trunk/LMDZ.MARS/libf/phymars/physiq_mod.F
r1713 r1720 353 353 REAL nccn(ngrid,nlayer) ! true n ccn (kg/kg) 354 354 REAL qccn(ngrid,nlayer) ! true q ccn (kg/kg) 355 REAL nccnco2(ngrid,nlayer) ! true n ccnco2 (kg/kg) 356 REAL qccnco2(ngrid,nlayer) ! true q ccnco2 (kg/kg) 355 357 356 358 c Test 1d/3d scavenging … … 1240 1242 ! We need to check that we have Nccn & Ndust > 0 1241 1243 ! This is due to single precision rounding problems 1242 if (microphysco2) then1243 1244 1244 1245 pdq(1:ngrid,1:nlayer,igcm_co2) = … … 1254 1255 & pdq(1:ngrid,1:nlayer,igcm_ccnco2_number) + 1255 1256 & zdqcloudco2(1:ngrid,1:nlayer,igcm_ccnco2_number) 1256 c Negative values? 1257 !Update water ice clouds values as well 1258 pdq(1:ngrid,1:nlayer,igcm_h2o_ice) = 1259 & pdq(1:ngrid,1:nlayer,igcm_h2o_ice) + 1260 & zdqcloudco2(1:ngrid,1:nlayer,igcm_h2o_ice) 1261 1262 ! increment dust and ccn masses and numbers 1263 ! We need to check that we have Nccn & Ndust > 0 1264 ! This is due to single precision rounding problems 1265 pdq(1:ngrid,1:nlayer,igcm_ccn_mass) = 1266 & pdq(1:ngrid,1:nlayer,igcm_ccn_mass) + 1267 & zdqcloudco2(1:ngrid,1:nlayer,igcm_ccn_mass) 1268 pdq(1:ngrid,1:nlayer,igcm_ccn_number) = 1269 & pdq(1:ngrid,1:nlayer,igcm_ccn_number) + 1270 & zdqcloudco2(1:ngrid,1:nlayer,igcm_ccn_number) 1271 where (pq(:,:,igcm_ccn_mass) + 1272 & ptimestep*pdq(:,:,igcm_ccn_mass) < 0.) 1273 pdq(:,:,igcm_ccn_mass) = 1274 & - pq(:,:,igcm_ccn_mass)/ptimestep + 1.e-30 1275 pdq(:,:,igcm_ccn_number) = 1276 & - pq(:,:,igcm_ccn_number)/ptimestep + 1.e-30 1277 end where 1278 where (pq(:,:,igcm_ccn_number) + 1279 & ptimestep*pdq(:,:,igcm_ccn_number) < 0.) 1280 pdq(:,:,igcm_ccn_mass) = 1281 & - pq(:,:,igcm_ccn_mass)/ptimestep + 1.e-30 1282 pdq(:,:,igcm_ccn_number) = 1283 & - pq(:,:,igcm_ccn_number)/ptimestep + 1.e-30 1284 end where 1285 c Negative values? 1257 1286 where (pq(:,:,igcm_ccnco2_mass) + 1258 1287 & ptimestep*pdq(:,:,igcm_ccnco2_mass) < 0.) … … 1261 1290 pdq(:,:,igcm_ccnco2_number) = 1262 1291 & - pq(:,:,igcm_ccnco2_number)/ptimestep + 1.e-30 1292 end where 1293 where (pq(:,:,igcm_ccnco2_number) + 1294 & ptimestep*pdq(:,:,igcm_ccnco2_number) < 0.) 1295 pdq(:,:,igcm_ccnco2_mass) = 1296 & - pq(:,:,igcm_ccnco2_mass)/ptimestep + 1.e-30 1297 pdq(:,:,igcm_ccnco2_number) = 1298 & - pq(:,:,igcm_ccnco2_number)/ptimestep + 1.e-30 1263 1299 end where 1264 where (pq(:,:,igcm_ccnco2_number) + 1265 & ptimestep*pdq(:,:,igcm_ccnco2_number) < 0.) 1266 pdq(:,:,igcm_ccnco2_mass) = 1267 & - pq(:,:,igcm_ccnco2_mass)/ptimestep + 1.e-30 1268 pdq(:,:,igcm_ccnco2_number) = 1269 & - pq(:,:,igcm_ccnco2_number)/ptimestep + 1.e-30 1270 end where 1271 endif !(of if micropĥysco2) 1272 1300 1273 1301 ! increment dust tracers tendancies 1274 if (scavenging) then1275 1302 pdq(1:ngrid,1:nlayer,igcm_dust_mass) = 1276 1303 & pdq(1:ngrid,1:nlayer,igcm_dust_mass) + … … 1294 1321 & - pq(:,:,igcm_dust_number)/ptimestep + 1.e-30 1295 1322 end where 1296 endif ! of if scavenging1323 1297 1324 END IF ! of IF (co2clouds) 1298 1325 … … 1821 1848 & zdqssed(ig,igcm_dust_number)*tauscaling(ig) 1822 1849 ndust(ig,:) = 1823 & pq(ig,:,igcm_dust_number)*tauscaling(ig)1850 & zq(ig,:,igcm_dust_number)*tauscaling(ig) 1824 1851 qdust(ig,:) = 1825 & pq(ig,:,igcm_dust_mass)*tauscaling(ig)1852 & zq(ig,:,igcm_dust_mass)*tauscaling(ig) 1826 1853 enddo 1827 1854 if (scavenging) then … … 1832 1859 & zdqssed(ig,igcm_ccn_number)*tauscaling(ig) 1833 1860 nccn(ig,:) = 1834 & pq(ig,:,igcm_ccn_number)*tauscaling(ig)1861 & zq(ig,:,igcm_ccn_number)*tauscaling(ig) 1835 1862 qccn(ig,:) = 1836 & pq(ig,:,igcm_ccn_mass)*tauscaling(ig)1863 & zq(ig,:,igcm_ccn_mass)*tauscaling(ig) 1837 1864 enddo 1838 1865 endif 1839 1866 endif 1867 if (co2clouds) then 1868 do ig=1,ngrid 1869 nccnco2(ig,:) = 1870 & zq(ig,:,igcm_ccnco2_number)*tauscaling(ig) 1871 qccnco2(ig,:) = 1872 & zq(ig,:,igcm_ccnco2_mass)*tauscaling(ig) 1873 enddo 1874 1875 1876 call WRITEDIAGFI(ngrid,'ccnqco2','CCNco2 mass mr', 1877 & 'kg/kg',3,qccnco2) 1878 call WRITEDIAGFI(ngrid,'ccnNco2','CCNco2 number', 1879 & 'part/kg',3,nccnco2) 1880 call WRITEDIAGFI(ngrid,'co2_ice','co2_ice','kg/kg', 1881 & 3,pq(:,:,igcm_co2_ice)) 1882 call WRITEDIAGFI(ngrid,'co2','co2','kg/kg', 1883 & 3,pq(:,:,igcm_co2)) 1884 call WRITEDIAGFI(ngrid,'h2o_ice','h2o_ice','kg/kg', 1885 & 3,pq(:,:,igcm_h2o_ice)) 1886 endif 1887 1840 1888 1841 1889 if (water) then … … 1921 1969 1922 1970 1923 1924 if (co2clouds) then1925 c mtotco2(:)=01926 c icetotco2(:)=01927 c raveco2(:)=01928 c do ig=1,ngrid1929 c do l=1,nlayer1930 c mtotco2(ig) = mtotco2(ig) +1931 c & zq(ig,l,igcm_co2) *1932 c & (zplev(ig,l) - zplev(ig,l+1)) / g1933 c icetotco2(ig) = icetotco2(ig) +1934 c & zq(ig,l,igcm_co2_ice) *1935 c & (zplev(ig,l) - zplev(ig,l+1)) / g1936 1937 c Computing abs optical depth at 825 cm-1 in each ! for now commented for CO2 - listo layer to simulate NEW TES retrieval1938 c Qabsice = min( max(0.4e6*riceco2(ig,l)*1939 c & (1.+nuiceco2_ref)-0.05 ,0.),1.2)1940 c opTESco2(ig,l)= 0.75 * Qabsice *1941 c & zq(ig,l,igcm_co2_ice) *1942 c & (zplev(ig,l) - zplev(ig,l+1)) / g1943 c & / (rho_ice_co2 * riceco2(ig,l)1944 c & * (1.+nuiceco2_ref))1945 c tauTESco2(ig)=tauTESco2(ig)+ opTESco2(ig,l)1946 c enddo1947 c enddo1948 call co2sat(ngrid*nlayer,zt,zplay,zqsatco2)1949 do ig=1,ngrid1950 do l=1,nlayer1951 satuco2(ig,l) = zq(ig,l,igcm_co2)*1952 & (mmean(ig,l)/44.01)*zplay(ig,l)/zqsatco2(ig,l)1953 enddo1954 enddo1955 1956 if (scavenging) then1957 Nccntot(:)= 01958 Mccntot(:)= 01959 raveco2(:)=01960 icetotco2(:)=01961 do ig=1,ngrid1962 do l=1,nlayer1963 icetotco2(ig) = icetotco2(ig) +1964 & zq(ig,l,igcm_co2_ice) *1965 & (zplev(ig,l) - zplev(ig,l+1)) / g1966 Nccntot(ig) = Nccntot(ig) +1967 & zq(ig,l,igcm_ccnco2_number)*tauscaling(ig)1968 & *(zplev(ig,l) - zplev(ig,l+1)) / g1969 Mccntot(ig) = Mccntot(ig) +1970 & zq(ig,l,igcm_ccnco2_mass)*tauscaling(ig)1971 & *(zplev(ig,l) - zplev(ig,l+1)) / g1972 cccc Column integrated effective ice radius1973 cccc is weighted by total ice surface area (BETTER than total ice mass)1974 raveco2(ig) = raveco2(ig) +1975 & tauscaling(ig) *1976 & zq(ig,l,igcm_ccnco2_number) *1977 & (zplev(ig,l) - zplev(ig,l+1)) / g *1978 & riceco2(ig,l) * riceco2(ig,l)*1979 & (1.+nuiceco2_ref)1980 enddo1981 raveco2(ig)=(icetotco2(ig)/rho_ice_co2+Mccntot(ig)/1982 & rho_dust)*0.751983 & /max(pi*raveco2(ig),1.e-30) ! surface weight1984 if (icetotco2(ig)*1e3.lt.0.01) raveco2(ig)=0.1985 enddo1986 else ! of if (scavenging)1987 raveco2(:)=01988 do ig=1,ngrid1989 do l=1,nlayer1990 raveco2(ig) = raveco2(ig) +1991 & zq(ig,l,igcm_co2_ice) *1992 & (zplev(ig,l) - zplev(ig,l+1)) / g *1993 & riceco2(ig,l) * (1.+nuiceco2_ref)1994 enddo1995 raveco2(ig) = max(raveco2(ig) /1996 & max(icetotco2(ig),1.e-30),1.e-30) ! mass weight1997 enddo1998 endif ! of if (scavenging)1999 endif ! of if (co2clouds)2000 1971 endif ! of if (tracer) 2001 1972 … … 2341 2312 & fluxtop_sw_tot) 2342 2313 call WRITEDIAGFI(ngrid,"temp","temperature","K",3,zt) 2314 call WRITEDIAGFI(ngrid,"Time","Time","sols",0,zday) 2315 2343 2316 call WRITEDIAGFI(ngrid,"u","Zonal wind","m.s-1",3,zu) 2344 2317 call WRITEDIAGFI(ngrid,"v","Meridional wind","m.s-1",3,zv) … … 2355 2328 c & 'w.m-2',3,zdtlw) 2356 2329 2357 if (co2clouds) then2358 do iq=1, nq2359 call WRITEDIAGFI(ngrid,trim(noms(iq)),2360 & trim(noms(iq)),'kg/kg',3,zq(:,:,iq))2361 end do2362 endif2363 2330 if (.not.activice) then 2364 2331 CALL WRITEDIAGFI(ngrid,'tauTESap', … … 2915 2882 endif ! of if (scavenging) 2916 2883 2917 if (co2clouds) then2918 if (scavenging) then2919 Nccntot(:)= 02920 Mccntot(:)= 02921 raveco2(:)=02922 icetotco2(:)=02923 do ig=1,ngrid2924 do l=1,nlayer2925 icetotco2(ig) = icetotco2(ig) +2926 & zq(ig,l,igcm_co2_ice) *2927 & (zplev(ig,l) - zplev(ig,l+1)) / g2928 Nccntot(ig) = Nccntot(ig) +2929 & zq(ig,l,igcm_ccnco2_number)*tauscaling(ig)2930 & *(zplev(ig,l) - zplev(ig,l+1)) / g2931 Mccntot(ig) = Mccntot(ig) +2932 & zq(ig,l,igcm_ccnco2_mass)*tauscaling(ig)2933 & *(zplev(ig,l) - zplev(ig,l+1)) / g2934 cccc Column integrated effective ice radius2935 cccc is weighted by total ice surface area (BETTER than total ice mass)2936 raveco2(ig) = raveco2(ig) +2937 & tauscaling(ig) *2938 & zq(ig,l,igcm_ccnco2_number) *2939 & (zplev(ig,l) - zplev(ig,l+1)) / g *2940 & riceco2(ig,l) * riceco2(ig,l)*2941 & (1.+nuiceco2_ref)2942 enddo2943 raveco2(ig)=(icetotco2(ig)/rho_ice_co2+Mccntot(ig)/2944 & rho_dust)*0.752945 & /max(pi*raveco2(ig),1.e-30) ! surface weight2946 if (icetotco2(ig)*1e3.lt.0.01) raveco2(ig)=0.2947 enddo2948 else ! of if (scavenging)2949 raveco2(:)=02950 do ig=1,ngrid2951 do l=1,nlayer2952 raveco2(ig) = raveco2(ig) +2953 & zq(ig,l,igcm_co2_ice) *2954 & (zplev(ig,l) - zplev(ig,l+1)) / g *2955 & riceco2(ig,l) * (1.+nuiceco2_ref)2956 enddo2957 raveco2(ig) = max(raveco2(ig) /2958 & max(icetotco2(ig),1.e-30),1.e-30) ! mass weight2959 enddo2960 endif ! of if (scavenging)2961 2962 call WRITEdiagfi(ngrid,"raveco2","ice eff radius","m",02963 & ,raveco2)2964 2965 endif ! of if (co2clouds)2966 2884 2967 2885 CALL WRITEDIAGFI(ngrid,'reffice', -
trunk/LMDZ.MARS/libf/phymars/tracer_mod.F90
r1660 r1720 25 25 26 26 real,save :: ccn_factor ! ratio of nuclei for water ice particles 27 28 real,save :: meteo_flux_mass ! meteoritic flux of dust - mmr29 real,save :: meteo_flux_number ! meteoritic flux of dust - nd30 integer,save :: meteo_alt31 27 32 28 INTEGER,ALLOCATABLE,SAVE :: nqdust(:) ! to store the indexes of dust tracers (cf aeropacity)
Note: See TracChangeset
for help on using the changeset viewer.
