Ignore:
Timestamp:
Apr 17, 2018, 3:58:30 PM (7 years ago)
Author:
emillour
Message:

Mars GCM:
CO2 code update:

  • nucleaCO2.F: code cleanup and use co2useh2o flag to handle cases where

h2o ccns have to be tracked and accounted for.

  • co2cloud.F & improvedco2clouds.F: code cleanup and use co2useh2o flag to

control updates on water variables if necessary.

  • physiq.F : cleanup on outputs & compute mtotco2 and icetotco2

DB

Location:
trunk/LMDZ.MARS/libf/phymars
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.MARS/libf/phymars/co2cloud.F

    r1918 r1921  
    2020      USE newsedim_mod, ONLY: newsedim
    2121      USE datafile_mod, ONLY: datadir
     22
    2223      IMPLICIT NONE
    2324
    … …  
    308309        firstcall=.false.
    309310      ENDIF                     ! of IF (firstcall)
    310    
    311 c-----Initialization
     311c ===========================================================================   
     312c     Initialization
     313c ===========================================================================
    312314      dev2 = 1. / ( sqrt(2.) * sigma_iceco2 )
    313315      beta=0.85
    … …  
    335337      enddo
    336338
    337 c-------------------------------------------------------------------
     339c ==========================================================================
    338340c   0.  Representation of sub-grid water ice clouds
    339 c-------------------------------------------------------------------
     341c ==========================================================================
    340342      IF (CLFvaryingCO2) THEN
    341343
    … …  
    345347         co2cloudfrac(:,:)=mincloud
    346348         
    347 c-----Tendencies
     349c Tendencies
    348350         DO l=1,nlay
    349351            DO ig=1,ngrid
    … …  
    439441         ENDDO ! of DO ig=1,ngrid
    440442!     Totalcloud frac of the column missing here
    441 c-----------------------
    442 c-----No sub-grid cloud representation (CLFvarying=false)
     443c
     444c No sub-grid cloud representation (CLFvarying=false)
    443445      ELSE
    444446         DO l=1,nlay
    … …  
    448450         END DO
    449451      END IF                    ! end if (CLFvaryingco2)
    450 c------------------------------------------------------------------
     452c =============================================================================
    451453c microtimestep timeloop for microphysics:
    452454c 0.Stepped entry for tendancies
    … …  
    454456c 2.Call co2clouds microphysics
    455457c 3.Update tendancies
    456 c------------------------------------------------------------------
     458c =============================================================================
    457459      DO microstep=1,imicroco2
    458 c------ Temperature tendency subpdt
     460c Temperature tendency subpdt
    459461        ! If imicro=1 subpdt is the same as pdt
    460462        DO l=1,nlay
    … …  
    482484     &              sum_subpdq(ig,l,igcm_co2)
    483485     &              + pdq(ig,l,igcm_co2)
    484 
     486c D.BARDET :
     487               if (co2useh2o) then
    485488               sum_subpdq(ig,l,igcm_h2o_ice) =
    486489     &              sum_subpdq(ig,l,igcm_h2o_ice)
    487490     &              + pdq(ig,l,igcm_h2o_ice)
     491
    488492               sum_subpdq(ig,l,igcm_ccn_mass) =
    489493     &              sum_subpdq(ig,l,igcm_ccn_mass)
    … …  
    492496     &              sum_subpdq(ig,l,igcm_ccn_number)
    493497     &              + pdq(ig,l,igcm_ccn_number)
     498                endif
    494499          ENDDO
    495500        ENDDO
    496 c- Effective tracers quantities in the cloud fraction
     501c Effective tracers quantities in the cloud fraction
    497502        IF (CLFvaryingCO2) THEN     
    498503            pqeff(:,:,:)=pq(:,:,:) ! prevent from buggs (A. Pottier)
    … …  
    507512        END IF 
    508513       
    509 c------------------------------------------------------
     514c ========================================================================
    510515c 1.SEDIMENTATION : update tracers, compute parameters,
    511516c   call to sedimentation routine, update tendancies
    512 c------------------------------------------------------
     517c ========================================================================
    513518        IF (sedimentation) THEN
    514519       
    … …  
    543548          ENDDO
    544549        ENDDO
    545 !     Gravitational sedimentation       
     550! Gravitational sedimentation       
    546551        zqsed0(:,:,igcm_co2_ice)=zqsed(:,:,igcm_co2_ice)
    547552        zqsed0(:,:,igcm_ccnco2_mass)=zqsed(:,:,igcm_ccnco2_mass)
    548553        zqsed0(:,:,igcm_ccnco2_number)=zqsed(:,:,igcm_ccnco2_number)
    549        !We save actualized tracer values to compute sedimentation tendancies
     554! We save actualized tracer values to compute sedimentation tendancies
    550555        call newsedim(ngrid,nlay,ngrid*nlay,ngrid*nlay,
    551556     &     microtimestep,pplev,masse,epaisseur,ztsed,
    552557     &     rsedcloudco2,rhocloudco2t,
    553558     &     zqsed(:,:,igcm_co2_ice),wq,beta) !  3 traceurs
    554 !     sedim at the surface of co2 ice : keep track of it for physiq_mod
     559
     560! sedim at the surface of co2 ice : keep track of it for physiq_mod
    555561        do ig=1,ngrid
    556562          sum_subpdqs_sedco2(ig)=
    557563     &         sum_subpdqs_sedco2(ig)+ wq(ig,1)/microtimestep
    558564        end do
     565
    559566        call newsedim(ngrid,nlay,ngrid*nlay,ngrid*nlay,
    560567     &     microtimestep,pplev,masse,epaisseur,ztsed,
    561568     &     rsedcloudco2,rhocloudco2t,
    562569     &     zqsed(:,:,igcm_ccnco2_mass),wq,beta)
     570
    563571        call newsedim(ngrid,nlay,ngrid*nlay,ngrid*nlay,
    564572     &     microtimestep,pplev,masse,epaisseur,ztsed,
    565573     &     rsedcloudco2,rhocloudco2t,
    566574     &     zqsed(:,:,igcm_ccnco2_number),wq,beta)
     575
    567576        DO l = 1, nlay            !Compute tendencies
    568577          DO ig=1,ngrid
    … …  
    578587          ENDDO
    579588        ENDDO
    580       !update subtimestep tendencies with sedimentation input
     589! update subtimestep tendencies with sedimentation input
    581590        DO l=1,nlay
    582591         DO ig=1,ngrid
    … …  
    595604        END IF !(end if sedimentation)
    596605       
    597 c------------------------------------------------------
     606c ==============================================================================
    598607c      2.  Main call to the cloud schemes:
    599 c------------------------------------------------------
     608c ==============================================================================
    600609        CALL improvedCO2clouds(ngrid,nlay,microtimestep,
    601610     &     pplay,pplev,pteff,sum_subpdt,
    602611     &     pqeff,sum_subpdq,subpdqcloudco2,subpdtcloudco2,
    603612     &     nq,tauscaling,memdMMccn,memdMMh2o,memdNNccn)
    604 c-----------------------------------------------------
     613c ==============================================================================
    605614c      3.  Updating tendencies after cloud scheme:
    606 c-----------------------------------------------------
     615c ==============================================================================
    607616        DO l=1,nlay
    608617          DO ig=1,ngrid
    … …  
    630639     &              sum_subpdq(ig,l,igcm_co2)
    631640     &              + subpdqcloudco2(ig,l,igcm_co2)
    632 
     641c D.BARDET :
     642               if (co2useh2o) then
    633643               sum_subpdq(ig,l,igcm_h2o_ice) =
    634644     &              sum_subpdq(ig,l,igcm_h2o_ice)
    635645     &              + subpdqcloudco2(ig,l,igcm_h2o_ice)
     646
    636647               sum_subpdq(ig,l,igcm_ccn_mass) =
    637648     &              sum_subpdq(ig,l,igcm_ccn_mass)
    … …  
    640651     &              sum_subpdq(ig,l,igcm_ccn_number)
    641652     &              + subpdqcloudco2(ig,l,igcm_ccn_number)
     653                endif
    642654          ENDDO
    643655        ENDDO
    … …  
    651663         pdqs_sedco2(ig)=sum_subpdqs_sedco2(ig)/real(imicroco2)
    652664      enddo
    653 c------ Temperature tendency pdtcloud
     665c Temperature tendency pdtcloud
    654666      DO l=1,nlay
    655667        DO ig=1,ngrid
    … …  
    658670        ENDDO
    659671      ENDDO
    660 c------ Tracers tendencies pdqcloud
     672c Tracers tendencies pdqcloud
    661673      DO l=1,nlay
    662674        DO ig=1,ngrid       
    … …  
    667679     &            sum_subpdq(ig,l,igcm_co2)/real(imicroco2)
    668680     &            - pdq(ig,l,igcm_co2)
     681c D.BARDET :
     682             if (co2useh2o) then
    669683             pdqcloudco2(ig,l,igcm_h2o_ice) =
    670684     &            sum_subpdq(ig,l,igcm_h2o_ice)/real(imicroco2)
    671685     &            - pdq(ig,l,igcm_h2o_ice)
    672         ENDDO
    673       ENDDO
    674       DO l=1,nlay
    675         DO ig=1,ngrid
     686
     687             pdqcloudco2(ig,l,igcm_ccn_mass) =
     688     &            sum_subpdq(ig,l,igcm_ccn_mass)/real(imicroco2)
     689     &            - pdq(ig,l,igcm_ccn_mass)
     690
     691             pdqcloudco2(ig,l,igcm_ccn_number) =
     692     &            sum_subpdq(ig,l,igcm_ccn_number)/real(imicroco2)
     693     &            - pdq(ig,l,igcm_ccn_number)
     694             endif
     695       
    676696             pdqcloudco2(ig,l,igcm_ccnco2_mass) =
    677697     &            sum_subpdq(ig,l,igcm_ccnco2_mass)/real(imicroco2)
    678698     &            - pdq(ig,l,igcm_ccnco2_mass)
     699       
    679700             pdqcloudco2(ig,l,igcm_ccnco2_number) =
    680701     &            sum_subpdq(ig,l,igcm_ccnco2_number)/real(imicroco2)
    681702     &            - pdq(ig,l,igcm_ccnco2_number)
    682              pdqcloudco2(ig,l,igcm_ccn_mass) =
    683      &            sum_subpdq(ig,l,igcm_ccn_mass)/real(imicroco2)
    684      &            - pdq(ig,l,igcm_ccn_mass)
    685              pdqcloudco2(ig,l,igcm_ccn_number) =
    686      &            sum_subpdq(ig,l,igcm_ccn_number)/real(imicroco2)
    687      &            - pdq(ig,l,igcm_ccn_number)
    688         ENDDO
    689       ENDDO
    690       DO l=1,nlay
    691         DO ig=1,ngrid
     703
    692704             pdqcloudco2(ig,l,igcm_dust_mass) =
    693705     &            sum_subpdq(ig,l,igcm_dust_mass)/real(imicroco2)
    694706     &            - pdq(ig,l,igcm_dust_mass)
     707
    695708             pdqcloudco2(ig,l,igcm_dust_number) =
    696709     &            sum_subpdq(ig,l,igcm_dust_number)/real(imicroco2)
    … …  
    698711        ENDDO
    699712      ENDDO
    700 c-------Due to stepped entry, other processes tendencies can add up to negative values
    701 c-------Therefore, enforce positive values and conserve mass
     713c Due to stepped entry, other processes tendencies can add up to negative values
     714c Therefore, enforce positive values and conserve mass
    702715      DO l=1,nlay
    703716        DO ig=1,ngrid
    … …  
    713726     &               - pqeff(ig,l,igcm_ccnco2_number)/ptimestep
    714727     &               - pdq(ig,l,igcm_ccnco2_number)+1.
     728               
    715729                pdqcloudco2(ig,l,igcm_dust_number) = 
    716730     &               -pdqcloudco2(ig,l,igcm_ccnco2_number)
     731
    717732                pdqcloudco2(ig,l,igcm_ccnco2_mass) =
    718733     &               - pqeff(ig,l,igcm_ccnco2_mass)/ptimestep
    719734     &               - pdq(ig,l,igcm_ccnco2_mass)+1.e-20
     735
    720736                pdqcloudco2(ig,l,igcm_dust_mass) =
    721737     &               -pdqcloudco2(ig,l,igcm_ccnco2_mass)
    … …  
    735751     &               - pqeff(ig,l,igcm_dust_number)/ptimestep
    736752     &               - pdq(ig,l,igcm_dust_number)+1.
     753
    737754                pdqcloudco2(ig,l,igcm_ccnco2_number) = 
    738755     &               -pdqcloudco2(ig,l,igcm_dust_number)
     756
    739757                pdqcloudco2(ig,l,igcm_dust_mass) =
    740758     &               - pqeff(ig,l,igcm_dust_mass)/ptimestep
    741759     &               - pdq(ig,l,igcm_dust_mass) +1.e-20
     760
    742761                pdqcloudco2(ig,l,igcm_ccnco2_mass) =
    743762     &               -pdqcloudco2(ig,l,igcm_dust_mass)
    … …  
    745764        ENDDO
    746765      ENDDO
    747         !pq+ptime*(pdq+pdqc)=1 ! pdqc=1-pq/ptime-pdq     
     766! pq+ptime*(pdq+pdqc)=1 ! pdqc=1-pq/ptime-pdq     
    748767      DO l=1,nlay
    749768        DO ig=1,ngrid
    … …  
    812831              endif
    813832     
    814       !update rice water
    815           call updaterice_micro(
    816      &    pqeff(ig,l,igcm_h2o_ice) +                    ! ice mass
    817      &   (pdq(ig,l,igcm_h2o_ice) +                   ! ice mass
    818      &    pdqcloudco2(ig,l,igcm_h2o_ice))*ptimestep,    ! ice mass
    819      &    pqeff(ig,l,igcm_ccn_mass) +                   ! ccn mass
    820      &   (pdq(ig,l,igcm_ccn_mass) +                  ! ccn mass
    821      &    pdqcloudco2(ig,l,igcm_ccn_mass))*ptimestep,   ! ccn mass
    822      &    pqeff(ig,l,igcm_ccn_number) +                 ! ccn number
    823      &   (pdq(ig,l,igcm_ccn_number) +                ! ccn number
    824      &    pdqcloudco2(ig,l,igcm_ccn_number))*ptimestep, ! ccn number
    825      &    tauscaling(ig),rice(ig,l),rhocloud(ig,l))
    826          
     833c D.BARDET : update rice water only if co2 use h2o ice as CCN
     834          if (co2useh2o) then
     835             call updaterice_micro(
     836     &       pqeff(ig,l,igcm_h2o_ice) +                    ! ice mass
     837     &       (pdq(ig,l,igcm_h2o_ice) +                     ! ice mass
     838     &       pdqcloudco2(ig,l,igcm_h2o_ice))*ptimestep,    ! ice mass
     839     &       pqeff(ig,l,igcm_ccn_mass) +                   ! ccn mass
     840     &       (pdq(ig,l,igcm_ccn_mass) +                    ! ccn mass
     841     &       pdqcloudco2(ig,l,igcm_ccn_mass))*ptimestep,   ! ccn mass
     842     &       pqeff(ig,l,igcm_ccn_number) +                 ! ccn number
     843     &       (pdq(ig,l,igcm_ccn_number) +                  ! ccn number
     844     &       pdqcloudco2(ig,l,igcm_ccn_number))*ptimestep, ! ccn number
     845     &       tauscaling(ig),rice(ig,l),rhocloud(ig,l))
     846          endif
     847
    827848          call updaterdust(
    828849     &    pqeff(ig,l,igcm_dust_mass) +                   ! dust mass
    829      &   (pdq(ig,l,igcm_dust_mass) +                  ! dust mass
     850     &   (pdq(ig,l,igcm_dust_mass) +                     ! dust mass
    830851     &    pdqcloudco2(ig,l,igcm_dust_mass))*ptimestep,   ! dust mass
    831852     &    pqeff(ig,l,igcm_dust_number) +                 ! dust number
    832      &   (pdq(ig,l,igcm_dust_number) +                ! dust number
     853     &   (pdq(ig,l,igcm_dust_number) +                   ! dust number
    833854     &    pdqcloudco2(ig,l,igcm_dust_number))*ptimestep, ! dust number
    834855     &    rdust(ig,l))
    … …  
    881902        DO l=1,nlay
    882903          DO ig=1,ngrid
     904
     905            pdqcloudco2(ig,l,igcm_ccnco2_mass)=
     906     &        pdqcloudco2(ig,l,igcm_ccnco2_mass)*co2cloudfrac(ig,l)
     907
     908            pdqcloudco2(ig,l,igcm_ccnco2_number)=
     909     &        pdqcloudco2(ig,l,igcm_ccnco2_number)*co2cloudfrac(ig,l)
     910
     911            pdqcloudco2(ig,l,igcm_dust_mass)=
     912     &        pdqcloudco2(ig,l,igcm_dust_mass)*co2cloudfrac(ig,l)
     913
     914            pdqcloudco2(ig,l,igcm_dust_number)=
     915     &        pdqcloudco2(ig,l,igcm_dust_number)*co2cloudfrac(ig,l)
     916c D.BARDET
     917            if (co2useh2o) then
     918            pdqcloudco2(ig,l,igcm_h2o_ice)=
     919     &        pdqcloudco2(ig,l,igcm_h2o_ice)*co2cloudfrac(ig,l)
     920
    883921            pdqcloudco2(ig,l,igcm_ccn_mass)=
    884922     &        pdqcloudco2(ig,l,igcm_ccn_mass)*co2cloudfrac(ig,l)
    885             pdqcloudco2(ig,l,igcm_ccnco2_mass)=
    886      &        pdqcloudco2(ig,l,igcm_ccnco2_mass)*co2cloudfrac(ig,l)
     923
    887924            pdqcloudco2(ig,l,igcm_ccn_number)=
    888      &        pdqcloudco2(ig,l,igcm_ccn_number)*co2cloudfrac(ig,l)
    889             pdqcloudco2(ig,l,igcm_ccnco2_number)=
    890      &        pdqcloudco2(ig,l,igcm_ccnco2_number)*co2cloudfrac(ig,l)
    891             pdqcloudco2(ig,l,igcm_dust_mass)=
    892      &        pdqcloudco2(ig,l,igcm_dust_mass)*co2cloudfrac(ig,l)
    893             pdqcloudco2(ig,l,igcm_dust_number)=
    894      &        pdqcloudco2(ig,l,igcm_dust_number)*co2cloudfrac(ig,l)
    895             pdqcloudco2(ig,l,igcm_h2o_ice)=
    896      &        pdqcloudco2(ig,l,igcm_h2o_ice)*co2cloudfrac(ig,l)
     925     &        pdqcloudco2(ig,l,igcm_ccn_number)*co2cloudfrac(ig,l)           
     926            endif
     927
    897928            pdqcloudco2(ig,l,igcm_co2_ice)=
    898929     &        pdqcloudco2(ig,l,igcm_co2_ice)*co2cloudfrac(ig,l)
     930
    899931            pdqcloudco2(ig,l,igcm_co2)=
    900932     &        pdqcloudco2(ig,l,igcm_co2)*co2cloudfrac(ig,l)
     933
    901934            pdtcloudco2(ig,l)=pdtcloudco2(ig,l)*co2cloudfrac(ig,l)
     935
    902936            Qext1bins2(ig,l)=Qext1bins2(ig,l)*co2cloudfrac(ig,l)
    903937          ENDDO
    … …  
    926960      call WRITEdiagfi(ngrid,"tau1mic","co2 ice opacity 1 micron"
    927961     &        ," ",2,tau1mic)
    928          
     962      call WRITEDIAGFI(ngrid,"memdNNccn","Nombre de CCN de glace d eau"
     963     &        ,"kg/kg ",3,memdNNccn)
     964      call WRITEDIAGFI(ngrid,"memdMMccn","Masse de CCN de glace d eau"
     965     &        ,"kg/kg ",3,memdMMccn)
     966      call WRITEDIAGFI(ngrid,"memdMMh2o","Masse de CCN de glace d eau"
     967     &        ,"kg/kg ",3,memdMMh2o)         
    929968      END

  • trunk/LMDZ.MARS/libf/phymars/improvedCO2clouds.F

    r1918 r1921  
    8383      LOGICAL,SAVE :: firstcall=.true.
    8484      REAL*8   derf ! Error function
    85       INTEGER ig,l,i,flag_pourri
    86      
     85      INTEGER ig,l,i
     86
    8787      real masse (ngrid,nlay) ! Layer mass (kg.m-2)
    8888      REAL rice(ngrid,nlay)    ! Water Ice mass mean radius (m)
    … …  
    102102      DOUBLE PRECISION pco2,psat  ! Co2 vapor partial pressure (Pa)
    103103      DOUBLE PRECISION satu ! Co2 vapor saturation ratio over ice
    104       DOUBLE PRECISION Mo,No
     104      DOUBLE PRECISION Mo,No,No_dust,Mo_dust
    105105      DOUBLE PRECISION  Rn, Rm, dev2,dev3, n_derf, m_derf
    106106      DOUBLE PRECISION memdMMccn(ngrid,nlay)
    … …  
    111111      DOUBLE PRECISION n_aer(nbinco2_cld) ! number concentration volume-1 of particle/each size bin
    112112      DOUBLE PRECISION m_aer(nbinco2_cld) ! mass mixing ratio of particle/each size bin
    113       DOUBLE PRECISION m_aer2(nbinco2_cld) ! mass mixing ratio of particle/each size bin
    114113      DOUBLE PRECISION m_aer_h2oice2(nbinco2_cld) ! mass mixing ratio of particle/each size bin
    115114
    … …  
    138137      DOUBLE PRECISION ratioh2o_ccn
    139138      DOUBLE PRECISION vo2co2
    140 c     Parameters of the size discretization
    141 c       used by the microphysical scheme
    142       DOUBLE PRECISION, PARAMETER :: rmin_cld = 1.e-9 ! Minimum radius (m)
    143       DOUBLE PRECISION, PARAMETER :: rmax_cld = 5.e-6 ! Maximum radius (m)
    144       DOUBLE PRECISION, PARAMETER :: rbmin_cld =1.e-10! Minimum bounary radius (m)
    145       DOUBLE PRECISION, PARAMETER :: rbmax_cld = 2.e-4 ! Maximum boundary radius (m)
    146       DOUBLE PRECISION vrat_cld ! Volume ratio
    147       DOUBLE PRECISION rb_cldco2(nbinco2_cld+1) ! boundary values of each rad_cldco2 bin (m)
     139
     140c     Parameters of the size discretization used by the microphysical scheme
     141      DOUBLE PRECISION, PARAMETER :: rmin_cld = 1.e-9   ! Minimum radius (m)
     142      DOUBLE PRECISION, PARAMETER :: rmax_cld = 5.e-6   ! Maximum radius (m)
     143      DOUBLE PRECISION, PARAMETER :: rbmin_cld =1.e-10  ! Minimum bounary radius (m)
     144      DOUBLE PRECISION, PARAMETER :: rbmax_cld = 2.e-4  ! Maximum boundary radius (m)
     145      DOUBLE PRECISION vrat_cld                         ! Volume ratio
     146      DOUBLE PRECISION rb_cldco2(nbinco2_cld+1)         ! boundary values of each rad_cldco2 bin (m)
    148147      SAVE rb_cldco2
    149       DOUBLE PRECISION dr_cld(nbinco2_cld) ! width of each rad_cldco2 bin (m)
    150       DOUBLE PRECISION vol_cld(nbinco2_cld)  ! particle volume for each bin (m3)
     148      DOUBLE PRECISION dr_cld(nbinco2_cld)              ! width of each rad_cldco2 bin (m)
     149      DOUBLE PRECISION vol_cld(nbinco2_cld)             ! particle volume for each bin (m3)
    151150
    152151      DOUBLE PRECISION Proba,Masse_atm,drsurdt,reff,Probah2o
    153       REAL sigma_iceco2 ! Variance of the co2 ice and CCN distributions
     152      REAL sigma_iceco2   ! Variance of the co2 ice and CCN distributions
    154153      SAVE sigma_iceco2
    155       REAL sigma_ice ! Variance of the h2o ice and CCN distributions
     154      REAL sigma_ice      ! Variance of the h2o ice and CCN distributions
    156155      SAVE sigma_ice
    157156      DOUBLE PRECISION Niceco2,Qccnco2,Nccnco2
    158       integer,save :: coeffh2o  !will be =0 is co2useh2o=.false.
     157
    159158!     Variables for the meteoritic flux:
    160159      integer,parameter :: nbin_meteor=100
    … …  
    217216c        mteta  = 0.952
    218217c        mtetaco2 = 0.952
    219         write(*,*) 'co2_param contact parameter:', mtetaco2
     218c        write(*,*) 'co2_param contact parameter:', mtetaco2
     219
    220220c       Volume of a co2 molecule (m3)
    221         vo1 = m0co2 / dble(rho_ice_co2) ! m0co2 et non mco2
    222 c       below, vo1 will be recomputed using real rho ice co2 (T-dependant)
    223         vo1co2=vo1 ! AJOUT JA
     221        vo1co2 = m0co2 / dble(rho_ice_co2) ! m0co2 et non mco2
     222
    224223c       Variance of the ice and CCN distributions
    225224        sigma_iceco2 = sqrt(log(1.+nuiceco2_sed))
    … …  
    231230        write(*,*) 'Volume of a co2 molecule:', vo1co2
    232231       
    233         coeffh2o=0
    234232        if (co2useh2o) then
    235233           write(*,*)
    … …  
    237235           write(*,*) "This means water ice particles can "
    238236           write(*,*) "serve as CCN for CO2 microphysics"
    239            coeffh2o=1
    240237        endif
     238
    241239        if (meteo_flux) then
    242240           write(*,*)
    … …  
    255253           endif
    256254!used Variables
    257 !           open(newunit=uMeteor,file=trim(datadir)//
    258255           open(unit=uMeteor,file=trim(datadir)//
    259256     &          '/Meteo_flux_Plane.dat'
    … …  
    298295! 1. Initialisation
    299296!=============================================================
    300       flag_pourri=0
     297
    301298      meteor_ccn(:,:,:)=0.
    302299      rice(:,:) = 1.e-8
    … …  
    373370           IF ( satu .ge. 1 ) THEN ! if there is condensation
    374371
    375 c     call updaterccn(zq(ig,l,igcm_dust_mass),
    376 c     &             zq(ig,l,igcm_dust_number),rdust(ig,l),tauscaling(ig))
    377 
    378               !We do rdust computation "by hand" because we don't want to
    379               ! change the mininumum rdust value in updaterad... It would
    380               ! mess up various part of the GCM !
     372c We do rdust computation "by hand" because we don't want to
     373c change the mininumum rdust value in updaterad... It would
     374c mess up various part of the GCM !
    381375             
    382376              rdust(ig,l)= zq(ig,l,igcm_dust_mass)
    383      &             *0.75/pi/rho_dust
    384      &             / zq(ig,l,igcm_dust_number)
     377     &             *0.75/(pi*rho_dust
     378     &             * zq(ig,l,igcm_dust_number))
    385379              rdust(ig,l)= rdust(ig,l)**(1./3.)
    386380              if (zq(ig,l,igcm_dust_mass)*tauscaling(ig) .le. 1.e-20
    … …  
    391385              rdust(ig,l)=min(5.e-4,rdust(ig,l))
    392386     
    393 c       Expand the dust moments into a binned distribution
     387c Expand the dust moments into a binned distribution
     388             
     389              n_aer(:)=0.
     390              m_aer(:)=0.
     391
    394392              Mo = zq(ig,l,igcm_dust_mass)* tauscaling(ig)+1.e-30
    395393              No = zq(ig,l,igcm_dust_number)* tauscaling(ig)+1.e-30
     394
     395              No_dust=No
     396              Mo_dust=Mo
     397
    396398              Rn = rdust(ig,l)
    397399              Rn = -dlog(Rn)
    … …  
    399401              n_derf = derf( (rb_cldco2(1)+Rn) *dev2)
    400402              m_derf = derf( (rb_cldco2(1)+Rm) *dev2)
    401      
     403
    402404              do i = 1, nbinco2_cld
    403405                 n_aer(i) = -0.5 * No * n_derf !! this ith previously computed
    … …  
    405407                 n_derf = derf((rb_cldco2(i+1)+Rn) *dev2)
    406408                 m_derf = derf((rb_cldco2(i+1)+Rm) *dev2)
    407                  n_aer(i) = n_aer(i) + 0.5 * No * n_derf +
    408      &                meteor_ccn(ig,l,i) !Ajout meteor_ccn particles
     409                 n_aer(i) = n_aer(i) + 0.5 * No * n_derf
    409410                 m_aer(i) = m_aer(i) + 0.5 * Mo * m_derf
    410                  m_aer2(i)=4./3.*pi*rho_dust
     411              enddo
     412
     413c Ajout meteor_ccn particles aux particules de poussière background
     414              if (meteo_flux) then
     415                 do i = 1, nbinco2_cld
     416                    n_aer(i) = n_aer(i) + meteor_ccn(ig,l,i)
     417                    m_aer(i) = m_aer(i) + 4./3.*pi*rho_dust
    411418     &                *meteor_ccn(ig,l,i)*rad_cldco2(i)*rad_cldco2(i)
    412419     &                *rad_cldco2(i)
    413               enddo
    414               if (meteo_flux) m_aer(:)=m_aer(:)+m_aer2(:)
     420                 enddo
     421              endif
    415422             
    416               !Same but with h2o particles as CCN
    417               call updaterice_micro(zq(ig,l,igcm_h2o_ice),
    418      &             zq(ig,l,igcm_ccn_mass),zq(ig,l,igcm_ccn_number),
    419      &             tauscaling(ig),rice(ig,l),rhocloud(ig,l))
    420               Mo = zq(ig,l,igcm_h2o_ice) +
    421      &             zq(ig,l,igcm_ccn_mass)*tauscaling(ig)+1.e-30!*tauscaling(ig)
    422    !  &             + 1.e-30 !Total mass of H20 crystals,CCN included
    423               No = zq(ig,l,igcm_ccn_number)* tauscaling(ig) + 1.e-30
    424               Rn = rice(ig,l)
    425               Rn = -dlog(Rn)
    426               Rm = Rn - 3. * sigma_ice*sigma_ice 
    427               n_derf = derf( (rb_cldco2(1)+Rn) *dev3)
    428               m_derf = derf( (rb_cldco2(1)+Rm) *dev3)
    429               do i = 1, nbinco2_cld
    430                  n_aer_h2oice(i) = -0.5 * No * n_derf
    431                  m_aer_h2oice(i) = -0.5 * Mo * m_derf
    432                  n_derf = derf( (rb_cldco2(i+1)+Rn) *dev3)
    433                  m_derf = derf( (rb_cldco2(i+1)+Rm) *dev3)
    434                  n_aer_h2oice(i) = n_aer_h2oice(i) + 0.5 * No * n_derf
    435                  m_aer_h2oice(i) = m_aer_h2oice(i) + 0.5 * Mo * m_derf
    436                  rad_h2oice(i) = rad_cldco2(i)
    437               enddo
    438               !Call to nucleation routine
     423c Same but with h2o particles as CCN only if co2useh2o=.true.
     424             
     425              n_aer_h2oice(:)=0.
     426              m_aer_h2oice(:)=0.
     427
     428              if (co2useh2o) then
     429                call updaterice_micro(zq(ig,l,igcm_h2o_ice),
     430     &               zq(ig,l,igcm_ccn_mass),zq(ig,l,igcm_ccn_number),
     431     &                 tauscaling(ig),rice(ig,l),rhocloud(ig,l))
     432                Mo = zq(ig,l,igcm_h2o_ice) +
     433     &               zq(ig,l,igcm_ccn_mass)*tauscaling(ig)+1.e-30
     434                     ! Total mass of H20 crystals,CCN included
     435                No = zq(ig,l,igcm_ccn_number)* tauscaling(ig) + 1.e-30
     436                Rn = rice(ig,l)
     437                Rn = -dlog(Rn)
     438                Rm = Rn - 3. * sigma_ice*sigma_ice 
     439                n_derf = derf( (rb_cldco2(1)+Rn) *dev3)
     440                m_derf = derf( (rb_cldco2(1)+Rm) *dev3)
     441                do i = 1, nbinco2_cld
     442                  n_aer_h2oice(i) = -0.5 * No * n_derf
     443                  m_aer_h2oice(i) = -0.5 * Mo * m_derf
     444                  n_derf = derf( (rb_cldco2(i+1)+Rn) *dev3)
     445                  m_derf = derf( (rb_cldco2(i+1)+Rm) *dev3)
     446                  n_aer_h2oice(i) = n_aer_h2oice(i) + 0.5 * No * n_derf
     447                  m_aer_h2oice(i) = m_aer_h2oice(i) + 0.5 * Mo * m_derf
     448                  rad_h2oice(i) = rad_cldco2(i)
     449                enddo
     450              endif
     451             
     452
     453! Call to nucleation routine
    439454              call nucleaCO2(dble(pco2),zt(ig,l),dble(satu)
    440455     &             ,n_aer,rate,n_aer_h2oice
    … …  
    445460              dMh2o = 0.
    446461              do i = 1, nbinco2_cld
    447                  Proba=1.0-dexp(-1.*microtimestep*rate(i))
    448                  Probah2o=coeffh2o*
    449      &            (1.0-dexp(-1.*microtimestep*rateh2o(i))) !if co2useh2o=.false., this is =0
    450                  dNh2o    = dNh2o + n_aer_h2oice(i) * Probah2o
    451                  dMh2o    = dMh2o + m_aer_h2oice(i) * Probah2o
     462                 Proba    =1.0-dexp(-1.*microtimestep*rate(i))
    452463                 dN       = dN + n_aer(i) * Proba
    453464                 dM       = dM + m_aer(i) * Proba             
    454465              enddo
    455         ! dM  masse activée (kg) et dN nb particules par  kg d'air
    456         ! Now increment CCN tracers and update dust tracers
     466              if (co2useh2o) then
     467                 do i = 1, nbinco2_cld
     468                    Probah2o = 1.0-dexp(-1.*microtimestep*rateh2o(i))
     469                    dNh2o    = dNh2o + n_aer_h2oice(i) * Probah2o
     470                    dMh2o    = dMh2o + m_aer_h2oice(i) * Probah2o
     471                 enddo
     472              endif
     473
     474! dM  masse activée (kg) et dN nb particules par  kg d'air
     475! Now increment CCN tracers and update dust tracers
    457476              dNN= dN/tauscaling(ig)
    458477              dMM= dM/tauscaling(ig)
    … …  
    466485              zq(ig,l,igcm_dust_number)=zq(ig,l,igcm_dust_number)-dNN
    467486
    468 c     Update CCN for CO2 nucleating on H2O CCN :
     487
     488c Update CCN for CO2 nucleating on H2O CCN :
    469489              ! Warning: must keep memory of it
    470490              if (co2useh2o) then
    … …  
    489509                memdMMccn(ig,l)=memdMMccn(ig,l)+dMh2o_ccn
    490510                memdNNccn(ig,l)=memdNNccn(ig,l)+dNNh2o
    491              endif !of if co2useh2o
    492            ENDIF                ! of is satu >1
    493 !=============================================================
    494 ! 4. Ice growth: scheme for radius evolution
     511             endif ! of if co2useh2o
     512           ENDIF   ! of is satu >1
     513
     514!=============================================================
     515! 5. Ice growth: scheme for radius evolution
    495516!=============================================================
    496517
    … …  
    505526
    506527              Ic_rice=0.
    507               flag_pourri=0
    508528              lw = l0 + l1 * zt(ig,l) + l2 * zt(ig,l)**2 +
    509529     &             l3 * zt(ig,l)**3 + l4 * zt(ig,l)**4 !J.kg-1
    … …  
    512532              !specific heat of atm cpp = 744.5 J.kg-1.K-1
    513533
    514 ccccccc call scheme of microphys. mass growth for CO2
     534c call scheme of microphys. mass growth for CO2
    515535              call massflowrateCO2(pplay(ig,l),zt(ig,l),
    516536     &             satu,riceco2(ig,l),mmean(ig,l),Ic_rice)
    517 c     Ic_rice Mass transfer rate (kg/s) for a rice particle >0 si croissance !
     537c Ic_rice Mass transfer rate (kg/s) for a rice particle >0 si croissance !
    518538             
    519539              if (isnan(Ic_rice) .or. Ic_rice .eq. 0.) then
    520540                 Ic_rice=0.
    521                  flag_pourri=1
    522541                 subpdtcloudco2(ig,l)=-sum_subpdt(ig,l)
    523542                 dMice=0
    … …  
    531550              dMice = max(dMice,max(-facteurmax,-zq(ig,l,igcm_co2_ice)))
    532551              dMice = min(dMice,min(facteurmax,zq(ig,l,igcm_co2)))
    533 !facteurmax maximum quantity of CO2 that can sublime/condense according to available thermal energy
     552! facteurmax maximum quantity of CO2 that can sublime/condense according to available thermal energy
    534553! latent heat release       >0 if growth i.e. if dMice >0
    535554              subpdtcloudco2(ig,l)=dMice*lw/cpp/microtimestep
    … …  
    541560
    542561!=============================================================
    543 ! 5. Dust cores releasing if no more co2 ice :
     562! 6. Dust cores releasing if no more co2 ice :
    544563!=============================================================
    545564
    546565              if (zq(ig,l,igcm_co2_ice).le. 1.e-25)THEN
    547  !On sublime tout
     566! On sublime tout
    548567                 if (co2useh2o) then
    549                  if (memdMMccn(ig,l) .gt. 0) then
     568                   if (memdMMccn(ig,l) .gt. 0) then
    550569                    zq(ig,l,igcm_ccn_mass)=zq(ig,l,igcm_ccn_mass)
    551570     &                   +memdMMccn(ig,l)
    552                  endif
    553                  if (memdMMh2o(ig,l) .gt. 0) then
     571                   endif
     572                   if (memdMMh2o(ig,l) .gt. 0) then
    554573                    zq(ig,l,igcm_h2o_ice)=zq(ig,l,igcm_h2o_ice)
    555574     &                   +memdMMh2o(ig,l)
    556                  endif
     575                   endif
    557576                 
    558                  if (memdNNccn(ig,l) .gt. 0) then
     577                   if (memdNNccn(ig,l) .gt. 0) then
    559578                    zq(ig,l,igcm_ccn_number)=zq(ig,l,igcm_ccn_number)
    560579     &                   +memdNNccn(ig,l)
     580                   endif
    561581                 endif
    562               endif
    563582                    zq(ig,l,igcm_dust_mass) =
    564583     &                   zq(ig,l,igcm_dust_mass)
    … …  
    582601              endif !of if co2_ice <1e-25
    583602
    584               ENDIF                ! of if NCCN > 1
     603              ENDIF            ! of if NCCN > 1
    585604          ENDDO                ! of ig loop
    586         ENDDO                   ! of nlayer loop 
    587 
    588 !=============================================================
    589 ! 6. END: get cloud tendencies
     605        ENDDO                  ! of nlayer loop 
     606
     607!=============================================================
     608! 7. END: get cloud tendencies
    590609!=============================================================
    591610
    … …  
    597616     &       (zq(1:ngrid,1:nlay,igcm_co2_ice) -
    598617     &       zq0(1:ngrid,1:nlay,igcm_co2_ice))/microtimestep
    599         subpdqcloudco2(1:ngrid,1:nlay,igcm_h2o_ice) =
    600      &       (zq(1:ngrid,1:nlay,igcm_h2o_ice) -
    601      &       zq0(1:ngrid,1:nlay,igcm_h2o_ice))/microtimestep
    602         subpdqcloudco2(1:ngrid,1:nlay,igcm_ccn_mass) =
    603      &       (zq(1:ngrid,1:nlay,igcm_ccn_mass) -
    604      &       zq0(1:ngrid,1:nlay,igcm_ccn_mass))/microtimestep
    605         subpdqcloudco2(1:ngrid,1:nlay,igcm_ccn_number) =
    606      &       (zq(1:ngrid,1:nlay,igcm_ccn_number) -
    607      &       zq0(1:ngrid,1:nlay,igcm_ccn_number))/microtimestep
     618       
     619        if (co2useh2o) then
     620           subpdqcloudco2(1:ngrid,1:nlay,igcm_h2o_ice) =
     621     &         (zq(1:ngrid,1:nlay,igcm_h2o_ice) -
     622     &         zq0(1:ngrid,1:nlay,igcm_h2o_ice))/microtimestep
     623           subpdqcloudco2(1:ngrid,1:nlay,igcm_ccn_mass) =
     624     &         (zq(1:ngrid,1:nlay,igcm_ccn_mass) -
     625     &         zq0(1:ngrid,1:nlay,igcm_ccn_mass))/microtimestep
     626           subpdqcloudco2(1:ngrid,1:nlay,igcm_ccn_number) =
     627     &         (zq(1:ngrid,1:nlay,igcm_ccn_number) -
     628     &         zq0(1:ngrid,1:nlay,igcm_ccn_number))/microtimestep
     629        endif
     630
    608631        subpdqcloudco2(1:ngrid,1:nlay,igcm_ccnco2_mass) =
    609632     &       (zq(1:ngrid,1:nlay,igcm_ccnco2_mass) -
    … …  
    619642     &       zq0(1:ngrid,1:nlay,igcm_dust_number))/microtimestep
    620643
     644
    621645        end
    622646     

  • trunk/LMDZ.MARS/libf/phymars/nucleaCO2.F

    r1818 r1921  
    2626 ! nucrate_h2o en sortie aussi :
    2727!nucleation sur dust et h2o separement ici
    28 #include "microphys.h"
     28
     29      include "microphys.h"
     30      include "callkeys.h"
    2931
    3032c     Inputs
    … …  
    5355
    5456
    55       if (sat .gt. 1.) then    ! minimum condition to activate nucleation
     57      IF (sat .gt. 1.) THEN    ! minimum condition to activate nucleation
    5658
    5759        nco2   = pco2 / kbz / temp
    … …  
    6365
    6466c       Loop over size bins
    65         do 200 i=1,nbinco2_cld
     67        do i=1,nbinco2_cld
    6668c            write(*,*) "IN NUCLEA, i, RAD_CLDCO2(i) = ",i, rad_cldco2(i),
    6769c     &          n_ccn(i)
    … …  
    7072c           no dust, no need to compute nucleation!
    7173            nucrate(i)=0.
    72             goto 210
    73           endif
     74c            goto 210
     75c          endif
     76          else
     77            if (rad_cldco2(i).gt.3000.*rstar) then
     78              zefshapeco2 = fshapeco2simple
     79            else
     80             zefshapeco2 = fshapeco2(mtetalocal,rad_cldco2(i)/rstar)
     81            endif
    7482
    75           if (rad_cldco2(i).gt.3000.*rstar) then
    76             zefshapeco2 = fshapeco2simple
    77           else
    78             zefshapeco2 = fshapeco2(mtetalocal,rad_cldco2(i)/rstar)
    79           endif
     83            fistar = (4./3.*pi) * sigco2 * (rstar * rstar) *
     84     &             zefshapeco2
     85            deltaf = (2.*desorpco2-surfdifco2-fistar)/
     86     &             (kbz*temp)
     87            deltaf = min( max(deltaf, -100.d0), 100.d0)
    8088
    81           fistar = (4./3.*pi) * sigco2 * (rstar * rstar) *
    82      &             zefshapeco2
    83           deltaf = (2.*desorpco2-surfdifco2-fistar)/
    84      &             (kbz*temp)
    85           deltaf = min( max(deltaf, -100.d0), 100.d0)
    86 
    87           if (deltaf.eq.-100.) then
    88             nucrate(i) = 0.
    89           else
    90             nucrate(i)= dble(sqrt ( fistar /
     89            if (deltaf.eq.-100.) then
     90                nucrate(i) = 0.
     91            else
     92                nucrate(i)= dble(sqrt ( fistar /
    9193     &               (3.*pi*kbz*temp*(gstar*gstar)) )
    9294     &                  * kbz * temp * rstar
    … …  
    98100
    99101           
    100           endif
     102            endif
     103          endif ! if n_ccn(i) .lt. 1e-10
    101104
    102 210     continue
     105          if (co2useh2o) then
    103106
    104           if ( n_ccn_h2oice(i) .lt. 1e-10 ) then
    105 c           no dust, no need to compute nucleation!
    106             nucrate_h2oice(i)=0.
    107             goto 200
    108           endif
     107            if ( n_ccn_h2oice(i) .lt. 1e-10 ) then
     108c           no H2O ice, no need to compute nucleation!
     109               nucrate_h2oice(i)=0.
     110            else   
     111               if (rad_h2oice(i).gt.3000.*rstar) then
     112                 zefshapeco2 = (2.+mtetalocalh)*(1.-mtetalocalh)*
     113     &                (1.-mtetalocalh) / 4.
     114               else  ! same m for dust/h2o ice
     115                 zefshapeco2 = fshapeco2(mtetalocalh,
     116     &                            (rad_h2oice(i)/rstar))
     117               endif
    109118
    110           if (rad_h2oice(i).gt.3000.*rstar) then
    111             zefshapeco2 = (2.+mtetalocalh)*(1.-mtetalocalh)*
    112      &                (1.-mtetalocalh) / 4.
    113           else
    114             zefshapeco2 = fshapeco2(mtetalocalh,rad_h2oice(i)/rstar) ! same m for dust/h2o ice
    115           endif
     119               fistar = (4./3.*pi) * sigco2 * (rstar * rstar) *
     120     &             zefshapeco2
     121              deltaf = (2.*desorpco2-surfdifco2-fistar)/
     122     &             (kbz*temp)
     123              deltaf = min( max(deltaf, -100.d0), 100.d0)
    116124
    117           fistar = (4./3.*pi) * sigco2 * (rstar * rstar) *
    118      &             zefshapeco2
    119           deltaf = (2.*desorpco2-surfdifco2-fistar)/
    120      &             (kbz*temp)
    121           deltaf = min( max(deltaf, -100.d0), 100.d0)
    122 
    123           if (deltaf.eq.-100.) then
    124             nucrate_h2oice(i) = 0.
    125           else
    126             nucrate_h2oice(i)= dble(sqrt ( fistar /
     125              if (deltaf.eq.-100.) then
     126                  nucrate_h2oice(i) = 0.
     127              else
     128                  nucrate_h2oice(i)= dble(sqrt ( fistar /
    127129     &               (3.*pi*kbz*temp*(gstar*gstar)) )
    128130     &                  * kbz * temp * rstar
    … …  
    132134     &                  / ( zefshapeco2 * nusco2 * m0co2 )
    133135     &                  * dexp (deltaf))
     136              endif
     137            endif
    134138          endif
    135          
    136          
     139        enddo
    137140
    138 200     continue
    139 
    140       else
     141      ELSE ! parcelle d'air non saturée
    141142
    142143        do i=1,nbinco2_cld
    … …  
    145146        enddo
    146147
    147       endif
     148      ENDIF ! if (sat .gt. 1.)
    148149
    149       return
    150150      end
    151151
    … …  
    175175          fshapeco2 = 0.5*fshapeco2
    176176
    177       return
    178177      end

  • trunk/LMDZ.MARS/libf/phymars/physiq_mod.F

    r1912 r1921  
    335335      REAL icetot(ngrid)        ! Total mass of water ice (kg/m2)
    336336      REAL mtotco2(ngrid)       ! Total mass of co2 vapor (kg/m2)
    337       REAL icetotco2(ngrid)        ! Total mass of co2 ice (kg/m2)
     337      REAL icetotco2(ngrid)     ! Total mass of co2 ice (kg/m2)
    338338      REAL Nccntot(ngrid)       ! Total number of ccn (nbr/m2)
    339339      REAL Mccntot(ngrid)       ! Total mass of ccn (kg/m2)
    … …  
    12841284     &               zdqcloudco2(1:ngrid,1:nlayer,igcm_ccn_number)
    12851285                where (pq(:,:,igcm_ccn_mass) +
    1286      &               ptimestep*pdq(:,:,igcm_ccn_mass) < 0.)
    1287                 pdq(:,:,igcm_ccn_mass) =
     1286     &                 ptimestep*pdq(:,:,igcm_ccn_mass) < 0.)
     1287                  pdq(:,:,igcm_ccn_mass) =
    12881288     &               - pq(:,:,igcm_ccn_mass)/ptimestep + 1.e-30
    1289                 pdq(:,:,igcm_ccn_number) =
     1289                  pdq(:,:,igcm_ccn_number) =
    12901290     &               - pq(:,:,igcm_ccn_number)/ptimestep + 1.e-30
    1291              end where
    1292              where (pq(:,:,igcm_ccn_number) +
    1293      &            ptimestep*pdq(:,:,igcm_ccn_number) < 0.)
    1294              pdq(:,:,igcm_ccn_mass) =
    1295      &            - pq(:,:,igcm_ccn_mass)/ptimestep + 1.e-30
    1296              pdq(:,:,igcm_ccn_number) =
    1297      &            - pq(:,:,igcm_ccn_number)/ptimestep + 1.e-30
    1298           end where
    1299        endif
     1291                end where
     1292                where (pq(:,:,igcm_ccn_number) +
     1293     &                 ptimestep*pdq(:,:,igcm_ccn_number) < 0.)
     1294                  pdq(:,:,igcm_ccn_mass) =
     1295     &              - pq(:,:,igcm_ccn_mass)/ptimestep + 1.e-30
     1296                  pdq(:,:,igcm_ccn_number) =
     1297     &              - pq(:,:,igcm_ccn_number)/ptimestep + 1.e-30
     1298                end where
     1299             endif ! of if (co2useh2o)
    13001300c     Negative values?
    13011301             where (pq(:,:,igcm_ccnco2_mass) +
    1302      &            ptimestep*pdq(:,:,igcm_ccnco2_mass) < 0.)
    1303              pdq(:,:,igcm_ccnco2_mass) =
     1302     &              ptimestep*pdq(:,:,igcm_ccnco2_mass) < 0.)
     1303                pdq(:,:,igcm_ccnco2_mass) =
    13041304     &            - pq(:,:,igcm_ccnco2_mass)/ptimestep + 1.e-30
    1305              pdq(:,:,igcm_ccnco2_number) =
     1305                pdq(:,:,igcm_ccnco2_number) =
    13061306     &            - pq(:,:,igcm_ccnco2_number)/ptimestep + 1.e-30
    13071307             end where
    13081308             where (pq(:,:,igcm_ccnco2_number) +
    1309      &            ptimestep*pdq(:,:,igcm_ccnco2_number) < 0.)
    1310              pdq(:,:,igcm_ccnco2_mass) =
     1309     &              ptimestep*pdq(:,:,igcm_ccnco2_number) < 0.)
     1310                pdq(:,:,igcm_ccnco2_mass) =
    13111311     &            - pq(:,:,igcm_ccnco2_mass)/ptimestep + 1.e-30
    1312              pdq(:,:,igcm_ccnco2_number) =
     1312                pdq(:,:,igcm_ccnco2_number) =
    13131313     &            - pq(:,:,igcm_ccnco2_number)/ptimestep + 1.e-30
    1314           end where
     1314             end where
    13151315       
    13161316c     Negative values?
    1317           where (pq(:,:,igcm_dust_mass) +
    1318      &         ptimestep*pdq(:,:,igcm_dust_mass) < 0.)
    1319           pdq(:,:,igcm_dust_mass) =
    1320      &         - pq(:,:,igcm_dust_mass)/ptimestep + 1.e-30
    1321           pdq(:,:,igcm_dust_number) =
    1322      &         - pq(:,:,igcm_dust_number)/ptimestep + 1.e-30
    1323        end where
    1324        where (pq(:,:,igcm_dust_number) +
    1325      &      ptimestep*pdq(:,:,igcm_dust_number) < 0.)
    1326        pdq(:,:,igcm_dust_mass) =
    1327      &      - pq(:,:,igcm_dust_mass)/ptimestep + 1.e-30
    1328        pdq(:,:,igcm_dust_number) =
    1329      &      - pq(:,:,igcm_dust_number)/ptimestep + 1.e-30
    1330       end where
     1317             where (pq(:,:,igcm_dust_mass) +
     1318     &              ptimestep*pdq(:,:,igcm_dust_mass) < 0.)
     1319                pdq(:,:,igcm_dust_mass) =
     1320     &            - pq(:,:,igcm_dust_mass)/ptimestep + 1.e-30
     1321                pdq(:,:,igcm_dust_number) =
     1322     &            - pq(:,:,igcm_dust_number)/ptimestep + 1.e-30
     1323             end where
     1324             where (pq(:,:,igcm_dust_number) +
     1325     &              ptimestep*pdq(:,:,igcm_dust_number) < 0.)
     1326                pdq(:,:,igcm_dust_mass) =
     1327     &            - pq(:,:,igcm_dust_mass)/ptimestep + 1.e-30
     1328                pdq(:,:,igcm_dust_number) =
     1329     &            - pq(:,:,igcm_dust_number)/ptimestep + 1.e-30
     1330             end where
    13311331     
    13321332      END IF                    ! of IF (co2clouds)
    … …  
    13671367           zdqssed(1:ngrid,1:nq)=0
    13681368
    1369 cSedimentation for co2 clouds tracers are inside co2cloud microtimestep
    1370 cZdqssed isn't
     1369c Sedimentation for co2 clouds tracers are inside co2cloud microtimestep
     1370c Zdqssed isn't
    13711371           call callsedim(ngrid,nlayer, ptimestep,
    13721372     &                zplev,zzlev, zzlay, pt, pdt, rdust, rice,
    … …  
    13741374     &                pq, pdq, zdqsed, zdqssed,nq,
    13751375     &                tau,tauscaling)
    1376    !Flux at the surface of co2 ice computed in co2cloud microtimestep
     1376c Flux at the surface of co2 ice computed in co2cloud microtimestep
    13771377           IF (co2clouds) THEN
    13781378              zdqssed(1:ngrid,igcm_co2_ice)=zdqssed_co2(1:ngrid)
    … …  
    18741874                enddo
    18751875              endif
    1876            endif
     1876           endif ! of (doubleq)
    18771877           if (co2clouds) then
    18781878              do ig=1,ngrid
    … …  
    18821882     &                zq(ig,:,igcm_ccnco2_mass)*tauscaling(ig)
    18831883              enddo
     1884c D. BARDET compute integrated CO2 vapor and ice content
     1885              mtotco2(:)=0
     1886              icetotco2(:)=0
     1887              do ig=1,ngrid
     1888               do l=1,nlayer
     1889                 mtotco2(ig) = mtotco2(ig) +
     1890     &                      zq(ig,l,igcm_co2) *
     1891     &                      (zplev(ig,l) - zplev(ig,l+1)) / g
     1892                 icetotco2(ig) = icetotco2(ig) +
     1893     &                        zq(ig,l,igcm_co2_ice) *
     1894     &                        (zplev(ig,l) - zplev(ig,l+1)) / g
     1895               enddo
     1896              enddo
    18841897             
    1885              
    1886               call WRITEDIAGFI(ngrid,'ccnqco2','CCNco2 mass mr',
    1887      &             'kg/kg',3,qccnco2)
    1888               call WRITEDIAGFI(ngrid,'ccnNco2','CCNco2 number',
    1889      &             'part/kg',3,nccnco2)
    1890               call WRITEDIAGFI(ngrid,'co2_ice','co2_ice','kg/kg',
    1891      &             3,pq(:,:,igcm_co2_ice))
    1892               call WRITEDIAGFI(ngrid,'co2','co2','kg/kg',
    1893      &             3,pq(:,:,igcm_co2))
    1894                 call WRITEDIAGFI(ngrid,'h2o_ice','h2o_ice','kg/kg',
    1895      &             3,pq(:,:,igcm_h2o_ice))
    1896            endif
     1898           endif ! of if (co2clouds)
    18971899           
    18981900                   
    … …  
    23732375          call WRITEDIAGFI(ngrid,"co2","co2 mass mixing ratio",
    23742376     &                     "kg/kg",3,zq(1,1,igcm_co2))
    2375        
    2376          ! Compute co2 column
    2377          co2col(:)=0
    2378          do l=1,nlayer
    2379            do ig=1,ngrid
    2380              co2col(ig)=co2col(ig)+
    2381      &                  zq(ig,l,igcm_co2)*(zplev(ig,l)-zplev(ig,l+1))/g
    2382            enddo
    2383          enddo
    2384          call WRITEDIAGFI(ngrid,"co2col","CO2 column","kg.m-2",2,
    2385      &                  co2col)
     2377          if (co2clouds) then
     2378            CALL WRITEDIAGFI(ngrid,'mtotco2',
     2379     &                       'total mass of CO2 vapor',
     2380     &                       'kg/m2',2,mtotco2)
     2381            CALL WRITEDIAGFI(ngrid,'zdtcloudco2',
     2382     &                       'temperature variation of CO2 latent heat',
     2383     &                       'K/s',3,zdtcloudco2)
     2384
     2385            CALL WRITEDIAGFI(ngrid,'icetotco2',
     2386     &                       'total mass of CO2 ice',
     2387     &                       'kg/m2',2,icetotco2)
     2388
     2389            call WRITEDIAGFI(ngrid,'ccnqco2','CCNco2 mass mr',
     2390     &                       'kg/kg',3,qccnco2)
     2391            call WRITEDIAGFI(ngrid,'ccnNco2','CCNco2 number',
     2392     &                       'part/kg',3,nccnco2)
     2393            call WRITEDIAGFI(ngrid,'co2_ice','co2_ice','kg/kg',
     2394     &                       3,zq(:,:,igcm_co2_ice))
     2395          endif ! of if (co2clouds)
    23862396         endif ! of if (tracer.and.(igcm_co2.ne.0))
    23872397
    … …  
    24352445     &                       'Mean reff',
    24362446     &                       'm',2,rave)
     2447            call WRITEDIAGFI(ngrid,'h2o_ice','h2o_ice','kg/kg',
     2448     &             3,zq(:,:,igcm_h2o_ice))
     2449            call WRITEDIAGFI(ngrid,'h2o_vap','h2o_vap','kg/kg',
     2450     &             3,zq(:,:,igcm_h2o_vap))
     2451
    24372452!A. Pottier
    24382453             CALL WRITEDIAGFI(ngrid,'rmoym',
Note: See TracChangeset for help on using the changeset viewer.