Index: /trunk/LMDZ.MARS/README =================================================================== --- /trunk/LMDZ.MARS/README (revision 1819) +++ /trunk/LMDZ.MARS/README (revision 1820) @@ -2484,2 +2484,5 @@ - All that you need to launch a GCM run with co2 clouds is described in co2clouds.F comments. Ehouarn Millour and Anni Määttänen have the files. +== 15/11/2017 == EM +- Bug fix in co2cloud.F (in the computation of the saturation index) along + with some code tidying. Index: /trunk/LMDZ.MARS/libf/phymars/co2cloud.F =================================================================== --- /trunk/LMDZ.MARS/libf/phymars/co2cloud.F (revision 1819) +++ /trunk/LMDZ.MARS/libf/phymars/co2cloud.F (revision 1820) @@ -6,9 +6,9 @@ & rsedcloud,rhocloud,pzlev,pdqs_sedco2, & pdu,pu) -! to use 'getin' + USE ioipsl_getincom, only: getin use dimradmars_mod, only: naerkind - USE comcstfi_h - USE ioipsl_getincom - USE updaterad + USE comcstfi_h, only: pi, g, cpp + USE updaterad, only: updaterice_microco2, updaterice_micro, + & updaterdust use conc_mod, only: mmean,rnew use tracer_mod, only: nqmx, igcm_co2, igcm_co2_ice, @@ -21,7 +21,7 @@ IMPLICIT NONE -#include "datafile.h" -#include "callkeys.h" -#include "microphys.h" + include "datafile.h" + include "callkeys.h" + include "microphys.h" c======================================================================= @@ -68,59 +68,44 @@ c----------------------------------------------------------------------- -c declarations: +c arguments: c ------------- -c Inputs: -c ------ - - INTEGER ngrid,nlay - INTEGER nq ! nombre de traceurs - REAL ptimestep ! pas de temps physique (s) - REAL pplev(ngrid,nlay+1) ! pression aux inter-couches (Pa) - REAL pplay(ngrid,nlay) ! pression au milieu des couches (Pa) - REAL pdpsrf(ngrid) ! tendence surf pressure - REAL pzlay(ngrid,nlay) ! altitude at the middle of the layers - REAL pt(ngrid,nlay) ! temperature at the middle of the layers (K) - REAL pdt(ngrid,nlay) ! tendence temperature des autres param. - real,intent(in) :: pzlev(ngrid,nlay+1) ! altitude at the boundaries of the layers - real pq(ngrid,nlay,nq) ! traceur (kg/kg) - real pdq(ngrid,nlay,nq) ! tendance avant condensation (kg/kg.s-1) - - real rice(ngrid,nlay) ! Water Ice mass mean radius (m) + INTEGER, INTENT(IN) :: ngrid,nlay + REAL, INTENT(IN) :: ptimestep ! pas de temps physique (s) + REAL, INTENT(IN) :: pplev(ngrid,nlay+1) ! Inter-layer pressures (Pa) + REAL, INTENT(IN) :: pplay(ngrid,nlay) ! mid-layer pressures (Pa) + REAL, INTENT(IN) :: pdpsrf(ngrid) ! tendency on surface pressure + REAL, INTENT(IN) :: pzlay(ngrid,nlay) ! altitude at the middle of the layers + REAL, INTENT(IN) :: pt(ngrid,nlay) ! temperature at the middle of the layers (K) + REAL, INTENT(IN) :: pdt(ngrid,nlay) ! tendency on temperature from other parametrizations + real, INTENT(IN) :: pq(ngrid,nlay,nq) ! tracers (kg/kg) + real, INTENT(IN) :: pdq(ngrid,nlay,nq) ! tendencies before condensation (kg/kg.s-1) + real, intent(out) :: pdqcloudco2(ngrid,nlay,nq) ! tendency due to CO2 condensation (kg/kg.s-1) + real, intent(out) :: pdtcloudco2(ngrid,nlay) ! tendency on temperature due to latent heat + INTEGER, INTENT(IN) :: nq ! number of tracers + REAL, INTENT(IN) :: tau(ngrid,naerkind) ! Column dust optical depth at each point + REAL, INTENT(IN) :: tauscaling(ngrid) ! Convertion factor for dust amount + REAL, INTENT(OUT) :: rdust(ngrid,nlay) ! Dust geometric mean radius (m) + real, intent(OUT) :: rice(ngrid,nlay) ! Water Ice mass mean radius (m) ! used for nucleation of CO2 on ice-coated ccns - DOUBLE PRECISION rho_ice_co2T(ngrid,nlay) !T-dependant CO2 ice density - DOUBLE PRECISION :: myT ! temperature scalar for co2 density computation - REAL tau(ngrid,naerkind) ! Column dust optical depth at each point - REAL tauscaling(ngrid) ! Convertion factor for dust amount - real rdust(ngrid,nlay) ! Dust geometric mean radius (m) - -c Outputs: -c ------- - - real pdqcloudco2(ngrid,nlay,nq) ! tendence de la condensation H2O(kg/kg.s-1) - REAL pdtcloudco2(ngrid,nlay) ! tendence temperature due - ! a la chaleur latente - DOUBLE PRECISION riceco2(ngrid,nlay) ! Ice mass mean radius (m) + DOUBLE PRECISION, INTENT(out) :: riceco2(ngrid,nlay) ! Ice mass mean radius (m) ! (r_c in montmessin_2004) - REAL nuice(ngrid,nlay) ! Estimated effective variance + REAL, INTENT(in) :: nuice(ngrid,nlay) ! Estimated effective variance ! of the size distribution - real rsedcloudco2(ngrid,nlay) ! Cloud sedimentation radius - real rhocloudco2(ngrid,nlay) ! Cloud density (kg.m-3) - real rhocloudco2t(ngrid,nlay) ! Cloud density (kg.m-3) - real pdqs_sedco2(ngrid) ! CO2 flux at the surface + real, intent(out) :: rsedcloudco2(ngrid,nlay) ! Cloud sedimentation radius + real, intent(out) :: rhocloudco2(ngrid,nlay) ! Cloud density (kg.m-3) + real, intent(out) :: rsedcloud(ngrid,nlay) ! Water Cloud sedimentation radius + real, intent(out) :: rhocloud(ngrid,nlay) ! Water Cloud density (kg.m-3) + real, intent(in) :: pzlev(ngrid,nlay+1) ! altitude at the boundaries of the layers + real, intent(out) :: pdqs_sedco2(ngrid) ! CO2 flux at the surface + REAL, INTENT(IN) :: pdu(ngrid,nlay),pu(ngrid,nlay) !Zonal Wind: zu=pu+pdu*ptimestep c local: c ------ - !water - real rsedcloud(ngrid,nlay) ! Cloud sedimentation radius - real rhocloud(ngrid,nlay) ! Cloud density (kg.m-3) - ! for ice radius computation ! for time loop INTEGER microstep ! time subsampling step variable - INTEGER imicroco2 ! time subsampling for coupled water microphysics & sedimentation - SAVE imicroco2 - REAL microtimestep ! integration timestep for coupled water microphysics & sedimentation - SAVE microtimestep + INTEGER, SAVE :: imicroco2 ! time subsampling for coupled water microphysics & sedimentation + REAL, SAVE :: microtimestep ! integration timestep for coupled water microphysics & sedimentation ! tendency given by clouds (inside the micro loop) @@ -135,4 +120,7 @@ REAL satuco2(ngrid,nlay) ! co2 satu ratio for output REAL zqsatco2(ngrid,nlay) ! saturation co2 + + DOUBLE PRECISION rho_ice_co2T(ngrid,nlay) !T-dependant CO2 ice density + DOUBLE PRECISION :: myT ! temperature scalar for co2 density computation INTEGER iq,ig,l,i @@ -154,6 +142,5 @@ ! What we need for Qext reading and tau computation : size distribution DOUBLE PRECISION vrat_cld ! Volume ratio - DOUBLE PRECISION rb_cldco2(nbinco2_cld+1) ! boundary values of each rad_cldco2 bin (m) - SAVE rb_cldco2 + DOUBLE PRECISION, SAVE :: rb_cldco2(nbinco2_cld+1) ! boundary values of each rad_cldco2 bin (m) DOUBLE PRECISION, PARAMETER :: rmin_cld = 1.e-9 ! Minimum radius (m) DOUBLE PRECISION, PARAMETER :: rmax_cld = 5.e-6 ! Maximum radius (m) @@ -162,11 +149,11 @@ DOUBLE PRECISION dr_cld(nbinco2_cld) ! width of each rad_cldco2 bin (m) DOUBLE PRECISION vol_cld(nbinco2_cld) ! particle volume for each bin (m3) - REAL sigma_iceco2 ! Variance of the ice and CCN distributions + REAL, SAVE :: sigma_iceco2 ! Variance of the ice and CCN distributions logical :: file_ok !Qext file reading double precision :: radv(10000),Qextv1mic(10000) - double precision :: Qext1bins(100),Qtemp + double precision, save :: Qext1bins(100) + double precision :: Qtemp double precision :: ltemp1(10000),ltemp2(10000) integer :: nelem,lebon1,lebon2,uQext - save Qext1bins DOUBLE PRECISION n_aer(nbinco2_cld),Rn,No,n_derf,dev2 DOUBLE PRECISION Qext1bins2(ngrid,nlay) @@ -177,4 +164,6 @@ REAL zt(ngrid,nlay) ! local value of temperature REAL :: zq(ngrid, nlay,nq) + + real :: rhocloudco2t(ngrid,nlay) ! Cloud density (kg.m-3) DOUBLE PRECISION :: tcond(ngrid,nlay) !CO2 condensation temperature @@ -187,5 +176,4 @@ REAL :: mincloud ! min cloud frac DOUBLE PRECISION:: rho,zu,NN,gradT !For Saturation Index computation - REAL :: pdu(ngrid,nlay),pu(ngrid,nlay) !Wind field zu=pu+pdu*ptimestep DOUBLE PRECISION :: SatIndex(ngrid,nlay),SatIndexmap(ngrid) @@ -378,31 +366,32 @@ if (satindexco2) then !logical in callphys.def - DO l=12,26 - ! layers 12 --> 26 ~ 12->85 km - DO ig=1,ngrid - ! Calcul de N^2 static stability - gradT=(zt(ig,l+1)-zt(ig,l))/(pzlev(ig,l+1)-pzlev(ig,l)) - NN=sqrt(g/zt(iq,l)*(max(gradT,-g/cpp)+g/cpp)) - !calcul of wind field - zu=pu(ig,l) + pdu(ig,l)*ptimestep -! calcul of background density - rho=pplay(ig,l)/(rnew(ig,l)*zt(ig,l)) - !saturation index - SatIndex(ig,l)=sqrt(7.5e-7*150.e3/(2.*pi)*NN/(rho*zu*zu*zu)) - enddo - enddo + DO l=12,26 + ! layers 12 --> 26 ~ 12->85 km + DO ig=1,ngrid + ! compute N^2 static stability + gradT=(zt(ig,l+1)-zt(ig,l))/(pzlev(ig,l+1)-pzlev(ig,l)) + NN=sqrt(g/zt(iq,l)*(max(gradT,-g/cpp)+g/cpp)) + ! compute absolute value of zonal wind field + zu=abs(pu(ig,l) + pdu(ig,l)*ptimestep) + ! compute background density + rho=pplay(ig,l)/(rnew(ig,l)*zt(ig,l)) + !saturation index + SatIndex(ig,l)=sqrt(7.5e-7*150.e3/(2.*pi)*NN/ + & (rho*zu*zu*zu)) + ENDDO + ENDDO !Then compute Satindex map ! layers 12 --> 26 ~ 12->85 km - DO ig=1,ngrid - SatIndexmap(ig)=maxval(SatIndex(ig,12:26)) - enddo - - call WRITEDIAGFI(ngrid,"SatIndexmap","SatIndexmap","km",2, - & SatIndexmap) + DO ig=1,ngrid + SatIndexmap(ig)=maxval(SatIndex(ig,12:26)) + ENDDO + + call WRITEDIAGFI(ngrid,"SatIndexmap","SatIndexmap","km",2, + & SatIndexmap) else - DO ig=1,ngrid - SatIndexmap(ig)=0.05 !maxval(SatIndex(ig,12:26)) - enddo - endif + do ig=1,ngrid + SatIndexmap(ig)=0.05 !maxval(SatIndex(ig,12:26)) + enddo + endif ! of if (satindexco2) !Modulate the DeltaT by GW propagation index : @@ -414,16 +403,15 @@ CALL tcondco2(ngrid,nlay,pplay,zqvap,tcond) -! A tester: CALL tcondco2(ngrid,nlay,pplay,zqvap,tcond) zdelt=spant - DO ig=1,ngrid + DO ig=1,ngrid - if (SatIndexmap(ig) .le. 0.1) THEN - DO l=1,nlay-1 + IF (SatIndexmap(ig) .le. 0.1) THEN + DO l=1,nlay-1 IF (tcond(ig,l) .ge. (zt(ig,l)+zdelt) - & .or. tcond(ig,l) .le. 0 ) THEN !Toute la fraction est saturée + & .or. tcond(ig,l) .le. 0 ) THEN !The entire fraction is saturated zteff(ig,l)=zt(ig,l) cloudfrac(ig,l)=1. - ELSE IF (tcond(ig,l) .le. (zt(ig,l)-zdelt)) THEN !Rien n'est saturé + ELSE IF (tcond(ig,l) .le. (zt(ig,l)-zdelt)) THEN ! No saturation at all zteff(ig,l)=zt(ig,l)-zdelt cloudfrac(ig,l)=mincloud @@ -431,5 +419,5 @@ cloudfrac(ig,l)=(tcond(ig,l)-zt(ig,l)+zdelt)/ & (2.0*zdelt) - zteff(ig,l)=(tcond(ig,l)+zt(ig,l)-zdelt)/2. !Temperature moyenne de la fraction nuageuse + zteff(ig,l)=(tcond(ig,l)+zt(ig,l)-zdelt)/2. !Mean temperature of the cloud fraction END IF !ig if (tcond(ig,l) ... zteff(ig,l)=zteff(ig,l)-pdt(ig,l)*ptimestep @@ -439,11 +427,11 @@ cloudfrac(ig,l)=1. END IF - ENDDO - ELSE + ENDDO + ELSE !SatIndex not favorable for GW : leave pt untouched - zteff(ig,l)=pt(ig,l) - cloudfrac(ig,l)=mincloud - END IF ! of if(SatIndexmap... - ENDDO + zteff(ig,l)=pt(ig,l) + cloudfrac(ig,l)=mincloud + ENDIF ! of if(SatIndexmap... + ENDDO ! of DO ig=1,ngrid ! Totalcloud frac of the column missing here c----------------------- @@ -466,6 +454,6 @@ c------ Temperature tendency subpdt ! If imicro=1 subpdt is the same as pdt - DO l=1,nlay - DO ig=1,ngrid + DO l=1,nlay + DO ig=1,ngrid subpdt(ig,l) = subpdt(ig,l) & + pdt(ig,l) ! At each micro timestep we add pdt in order to have a stepped entry @@ -500,8 +488,8 @@ & subpdq(ig,l,igcm_ccn_number) & + pdq(ig,l,igcm_ccn_number) - ENDDO - ENDDO + ENDDO + ENDDO c- Effective tracers quantities in the cloud fraction - IF (CLFvaryingCO2) THEN + IF (CLFvaryingCO2) THEN pqeff(:,:,:)=pq(:,:,:) ! prevent from buggs (A. Pottier) pqeff(:,:,igcm_ccnco2_mass) =pq(:,:,igcm_ccnco2_mass)/ @@ -511,7 +499,7 @@ pqeff(:,:,igcm_co2_ice)= pq(:,:,igcm_co2_ice)/ & cloudfrac(:,:) - ELSE + ELSE pqeff(:,:,:)=pq(:,:,:) - END IF + END IF c------------------------------------------------------ @@ -519,5 +507,5 @@ c call to sedimentation routine, update tendancies c------------------------------------------------------ - DO l=1, nlay + DO l=1, nlay DO ig=1,ngrid tempo_traceur_t(ig,l)=zteff(ig,l)+subpdt(ig,l) @@ -548,30 +536,30 @@ & riceco2(ig,l)) ENDDO - ENDDO + ENDDO ! Gravitational sedimentation - sav_trac(:,:,igcm_co2_ice)=tempo_traceurs(:,:,igcm_co2_ice) - sav_trac(:,:,igcm_ccnco2_mass)= + sav_trac(:,:,igcm_co2_ice)=tempo_traceurs(:,:,igcm_co2_ice) + sav_trac(:,:,igcm_ccnco2_mass)= & tempo_traceurs(:,:,igcm_ccnco2_mass) - sav_trac(:,:,igcm_ccnco2_number)= + sav_trac(:,:,igcm_ccnco2_number)= & tempo_traceurs(:,:,igcm_ccnco2_number) !We save actualized tracer values to compute sedimentation tendancies - call newsedim(ngrid,nlay,ngrid*nlay,ngrid*nlay, + call newsedim(ngrid,nlay,ngrid*nlay,ngrid*nlay, & microtimestep,pplev,masse,epaisseur,tempo_traceur_t, & rsedcloudco2,rhocloudco2t, & tempo_traceurs(:,:,igcm_co2_ice),wq,beta) ! 3 traceurs ! sedim at the surface of co2 ice : keep track of it for physiq_mod - do ig=1,ngrid - pdqs_sedco2(ig)=pdqs_sedco2(ig)+ wq(ig,1)/microtimestep - end do - call newsedim(ngrid,nlay,ngrid*nlay,ngrid*nlay, + do ig=1,ngrid + pdqs_sedco2(ig)=pdqs_sedco2(ig)+ wq(ig,1)/microtimestep + end do + call newsedim(ngrid,nlay,ngrid*nlay,ngrid*nlay, & microtimestep,pplev,masse,epaisseur,tempo_traceur_t, & rsedcloudco2,rhocloudco2t, & tempo_traceurs(:,:,igcm_ccnco2_mass),wq,beta) - call newsedim(ngrid,nlay,ngrid*nlay,ngrid*nlay, + call newsedim(ngrid,nlay,ngrid*nlay,ngrid*nlay, & microtimestep,pplev,masse,epaisseur,tempo_traceur_t, & rsedcloudco2,rhocloudco2t, & tempo_traceurs(:,:,igcm_ccnco2_number),wq,beta) - DO l = 1, nlay !Compute tendencies - DO ig=1,ngrid + DO l = 1, nlay !Compute tendencies + DO ig=1,ngrid pdqsed(ig,l,igcm_ccnco2_mass)= & (tempo_traceurs(ig,l,igcm_ccnco2_mass)- @@ -583,6 +571,6 @@ & (tempo_traceurs(ig,l,igcm_co2_ice)- & sav_trac(ig,l,igcm_co2_ice))/microtimestep - ENDDO - ENDDO + ENDDO + ENDDO !update subtimestep tendencies with sedimentation input DO l=1,nlay @@ -598,9 +586,9 @@ & +pdqsed(ig,l,igcm_co2_ice) ENDDO - ENDDO + ENDDO c------------------------------------------------------ c 2. Main call to the cloud schemes: c------------------------------------------------------ - CALL improvedCO2clouds(ngrid,nlay,microtimestep, + CALL improvedCO2clouds(ngrid,nlay,microtimestep, & pplay,pplev,zteff,subpdt, & pqeff,subpdq,subpdqcloudco2,subpdtcloudco2, @@ -609,6 +597,6 @@ c 3. Updating tendencies after cloud scheme: c----------------------------------------------------- - DO l=1,nlay - DO ig=1,ngrid + DO l=1,nlay + DO ig=1,ngrid subpdt(ig,l) = & subpdt(ig,l) + subpdtcloudco2(ig,l) @@ -644,6 +632,6 @@ & subpdq(ig,l,igcm_ccn_number) & + subpdqcloudco2(ig,l,igcm_ccn_number) - ENDDO - ENDDO + ENDDO + ENDDO ENDDO ! of DO microstep=1,imicro @@ -656,13 +644,13 @@ enddo c------ Temperature tendency pdtcloud - DO l=1,nlay - DO ig=1,ngrid + DO l=1,nlay + DO ig=1,ngrid pdtcloudco2(ig,l) = & subpdt(ig,l)/real(imicroco2)-pdt(ig,l) - ENDDO - ENDDO + ENDDO + ENDDO c------ Tracers tendencies pdqcloud - DO l=1,nlay - DO ig=1,ngrid + DO l=1,nlay + DO ig=1,ngrid pdqcloudco2(ig,l,igcm_co2_ice) = & subpdq(ig,l,igcm_co2_ice)/real(imicroco2) @@ -674,8 +662,8 @@ & subpdq(ig,l,igcm_h2o_ice)/real(imicroco2) & - pdq(ig,l,igcm_h2o_ice) - ENDDO - ENDDO - DO l=1,nlay - DO ig=1,ngrid + ENDDO + ENDDO + DO l=1,nlay + DO ig=1,ngrid pdqcloudco2(ig,l,igcm_ccnco2_mass) = & subpdq(ig,l,igcm_ccnco2_mass)/real(imicroco2) @@ -690,8 +678,8 @@ & subpdq(ig,l,igcm_ccn_number)/real(imicroco2) & - pdq(ig,l,igcm_ccn_number) - ENDDO - ENDDO - DO l=1,nlay - DO ig=1,ngrid + ENDDO + ENDDO + DO l=1,nlay + DO ig=1,ngrid pdqcloudco2(ig,l,igcm_dust_mass) = & subpdq(ig,l,igcm_dust_mass)/real(imicroco2) @@ -700,10 +688,10 @@ & subpdq(ig,l,igcm_dust_number)/real(imicroco2) & - pdq(ig,l,igcm_dust_number) - ENDDO - ENDDO + ENDDO + ENDDO c-------Due to stepped entry, other processes tendencies can add up to negative values c-------Therefore, enforce positive values and conserve mass - DO l=1,nlay - DO ig=1,ngrid + DO l=1,nlay + DO ig=1,ngrid IF ((pqeff(ig,l,igcm_ccnco2_number) + & ptimestep* (pdq(ig,l,igcm_ccnco2_number) + @@ -725,8 +713,8 @@ & -pdqcloudco2(ig,l,igcm_ccnco2_mass) ENDIF - ENDDO - ENDDO - DO l=1,nlay - DO ig=1,ngrid + ENDDO + ENDDO + DO l=1,nlay + DO ig=1,ngrid IF ( (pqeff(ig,l,igcm_dust_number) + & ptimestep* (pdq(ig,l,igcm_dust_number) + @@ -747,9 +735,9 @@ & -pdqcloudco2(ig,l,igcm_dust_mass) ENDIF - ENDDO - ENDDO + ENDDO + ENDDO !pq+ptime*(pdq+pdqc)=1 ! pdqc=1-pq/ptime-pdq - DO l=1,nlay - DO ig=1,ngrid + DO l=1,nlay + DO ig=1,ngrid IF (pqeff(ig,l,igcm_co2_ice) + ptimestep* & (pdq(ig,l,igcm_co2_ice) + pdqcloudco2(ig,l,igcm_co2_ice)) @@ -766,10 +754,10 @@ pdqcloudco2(ig,l,igcm_co2_ice)= -pdqcloudco2(ig,l,igcm_co2) ENDIF - ENDDO - ENDDO + ENDDO + ENDDO c Update clouds parameters values in the cloud fraction (for output) - DO l=1, nlay - DO ig=1,ngrid + DO l=1, nlay + DO ig=1,ngrid Niceco2=pqeff(ig,l,igcm_co2_ice) + @@ -801,19 +789,21 @@ & Nccnco2*tauscaling(ig) .le. 1.) )THEN riceco2(ig,l)=0. + Qext1bins2(ig,l)=0. + else +c Compute opacities + No=Nccnco2*tauscaling(ig) + Rn=-dlog(riceco2(ig,l)) + n_derf = derf( (rb_cldco2(1)+Rn) *dev2) + Qext1bins2(ig,l)=0. + do i = 1, nbinco2_cld + n_aer(i) = -0.5 * No * n_derf !! this ith previously computed + n_derf = derf((rb_cldco2(i+1)+Rn) *dev2) + n_aer(i) = n_aer(i) + 0.5 * No * n_derf + Qext1bins2(ig,l)=Qext1bins2(ig,l)+Qext1bins(i)*n_aer(i) + enddo endif -c Compute opacities - No=Nccnco2*tauscaling(ig) - Rn=-dlog(riceco2(ig,l)) - n_derf = derf( (rb_cldco2(1)+Rn) *dev2) - Qext1bins2(ig,l)=0. - do i = 1, nbinco2_cld - n_aer(i) = -0.5 * No * n_derf !! this ith previously computed - n_derf = derf((rb_cldco2(i+1)+Rn) *dev2) - n_aer(i) = n_aer(i) + 0.5 * No * n_derf - Qext1bins2(ig,l)=Qext1bins2(ig,l)+Qext1bins(i)*n_aer(i) - enddo !update rice water - call updaterice_micro( + call updaterice_micro( & pqeff(ig,l,igcm_h2o_ice) + ! ice mass & (pdq(ig,l,igcm_h2o_ice) + ! ice mass @@ -827,5 +817,5 @@ & tauscaling(ig),rice(ig,l),rhocloud(ig,l)) - call updaterdust( + call updaterdust( & pqeff(ig,l,igcm_dust_mass) + ! dust mass & (pdq(ig,l,igcm_dust_mass) + ! dust mass @@ -836,6 +826,6 @@ & rdust(ig,l)) - ENDDO - ENDDO + ENDDO + ENDDO c A correction if a lot of subliming CO2 fills the 1st layer FF04/2005 @@ -843,10 +833,10 @@ c Then that should not affect the ice particle radius - do ig=1,ngrid - if(pdpsrf(ig)*ptimestep.gt.0.9*(pplev(ig,1)-pplev(ig,2)))then + do ig=1,ngrid + if(pdpsrf(ig)*ptimestep.gt.0.9*(pplev(ig,1)-pplev(ig,2)))then if(pdpsrf(ig)*ptimestep.gt.0.9*(pplev(ig,1)-pplev(ig,3))) & riceco2(ig,2)=riceco2(ig,3) riceco2(ig,1)=riceco2(ig,2) - end if + endif end do @@ -860,5 +850,5 @@ ENDDO - DO l=1,nlay + DO l=1,nlay DO ig=1,ngrid rsedcloudco2(ig,l)=max(riceco2(ig,l)* @@ -867,18 +857,18 @@ c rsedcloudco2(ig,l)=min(rsedcloudco2(ig,l),1.e-5) ENDDO - ENDDO + ENDDO - call co2sat(ngrid*nlay,zteff+(pdt+pdtcloudco2)*ptimestep + call co2sat(ngrid*nlay,zteff+(pdt+pdtcloudco2)*ptimestep & ,pplay,zqsatco2) - do l=1,nlay - do ig=1,ngrid + do l=1,nlay + do ig=1,ngrid satuco2(ig,l) = (pqeff(ig,l,igcm_co2) + & (pdq(ig,l,igcm_co2) + & pdqcloudco2(ig,l,igcm_co2))*ptimestep)* & (mmean(ig,l)/44.01)*pplay(ig,l)/zqsatco2(ig,l) - enddo - enddo -!Tout ce qui est modifié par la microphysique de CO2 doit être rapporté a cloudfrac - IF (CLFvaryingCO2) THEN + enddo + enddo +!Everything modified by CO2 microphysics must be wrt cloudfrac + IF (CLFvaryingCO2) THEN DO l=1,nlay DO ig=1,ngrid @@ -905,27 +895,27 @@ ENDDO ENDDO - ENDIF -!l'opacité de la case ig est la somme sur l de Qext1bins2: Est-ce-vrai? - tau1mic(:)=0. - do l=1,nlay - do ig=1,ngrid - tau1mic(ig)=tau1mic(ig)+Qext1bins2(ig,l) - enddo - enddo + ENDIF +! opacity in mesh ig is the sum over l of Qext1bins2: Is this true ? + tau1mic(:)=0. + do l=1,nlay + do ig=1,ngrid + tau1mic(ig)=tau1mic(ig)+Qext1bins2(ig,l) + enddo + enddo !Outputs: - call WRITEDIAGFI(ngrid,"SatIndex","SatIndex"," ",3, + call WRITEDIAGFI(ngrid,"SatIndex","SatIndex"," ",3, & SatIndex) - call WRITEDIAGFI(ngrid,"satuco2","vap in satu","kg/kg",3, + call WRITEDIAGFI(ngrid,"satuco2","vap in satu","kg/kg",3, & satuco2) - call WRITEdiagfi(ngrid,"riceco2","ice radius","m" + call WRITEdiagfi(ngrid,"riceco2","ice radius","m" & ,3,riceco2) - call WRITEdiagfi(ngrid,"cloudfrac","co2 cloud fraction" + call WRITEdiagfi(ngrid,"cloudfrac","co2 cloud fraction" & ," ",3,cloudfrac) - call WRITEdiagfi(ngrid,"rsedcloudco2","rsed co2" + call WRITEdiagfi(ngrid,"rsedcloudco2","rsed co2" & ,"m",3,rsedcloudco2) - call WRITEdiagfi(ngrid,"Tau3D1mic"," co2 ice opacities" + call WRITEdiagfi(ngrid,"Tau3D1mic"," co2 ice opacities" & ," ",3,Qext1bins2) - call WRITEdiagfi(ngrid,"tau1mic","co2 ice opacity 1 micron" + call WRITEdiagfi(ngrid,"tau1mic","co2 ice opacity 1 micron" & ," ",2,tau1mic) - END + END