Changeset 1820


Ignore:
Timestamp:
Nov 15, 2017, 8:32:34 AM (7 years ago)
Author:
emillour
Message:

Mars GCM:

  • Bug fix in co2cloud.F (in the computation of the saturation index) along with some code tidying.

EM

Location:
trunk/LMDZ.MARS
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.MARS/README

    r1818 r1820  
    24842484- 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.
    24852485 
     2486== 15/11/2017 == EM
     2487- Bug fix in co2cloud.F (in the computation of the saturation index) along
     2488  with some code tidying.

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

    r1818 r1820  
    66     &                rsedcloud,rhocloud,pzlev,pdqs_sedco2,
    77     &                pdu,pu)
    8 ! to use  'getin'
     8      USE ioipsl_getincom, only: getin
    99      use dimradmars_mod, only: naerkind
    10       USE comcstfi_h
    11       USE ioipsl_getincom
    12       USE updaterad
     10      USE comcstfi_h, only: pi, g, cpp
     11      USE updaterad, only: updaterice_microco2, updaterice_micro,
     12     &                     updaterdust
    1313      use conc_mod, only: mmean,rnew
    1414      use tracer_mod, only: nqmx, igcm_co2, igcm_co2_ice,
    … …  
    2121      IMPLICIT NONE
    2222
    23 #include "datafile.h"
    24 #include "callkeys.h"
    25 #include "microphys.h"
     23      include "datafile.h"
     24      include "callkeys.h"
     25      include "microphys.h"
    2626
    2727c=======================================================================
    … …  
    6868
    6969c-----------------------------------------------------------------------
    70 c   declarations:
     70c   arguments:
    7171c   -------------
    7272
    73 c   Inputs:
    74 c   ------
    75 
    76       INTEGER ngrid,nlay
    77       INTEGER nq                 ! nombre de traceurs
    78       REAL ptimestep            ! pas de temps physique (s)
    79       REAL pplev(ngrid,nlay+1)   ! pression aux inter-couches (Pa)
    80       REAL pplay(ngrid,nlay)     ! pression au milieu des couches (Pa)
    81       REAL pdpsrf(ngrid)         ! tendence surf pressure
    82       REAL pzlay(ngrid,nlay)     ! altitude at the middle of the layers
    83       REAL pt(ngrid,nlay)        ! temperature at the middle of the layers (K)
    84       REAL pdt(ngrid,nlay)       ! tendence temperature des autres param.
    85       real,intent(in) :: pzlev(ngrid,nlay+1) ! altitude at the boundaries of the layers
    86       real pq(ngrid,nlay,nq)     ! traceur (kg/kg)
    87       real pdq(ngrid,nlay,nq)    ! tendance avant condensation  (kg/kg.s-1)
    88 
    89       real rice(ngrid,nlay)    ! Water Ice mass mean radius (m)
     73      INTEGER, INTENT(IN) :: ngrid,nlay
     74      REAL, INTENT(IN) :: ptimestep            ! pas de temps physique (s)
     75      REAL, INTENT(IN) :: pplev(ngrid,nlay+1)   ! Inter-layer pressures (Pa)
     76      REAL, INTENT(IN) :: pplay(ngrid,nlay)     ! mid-layer pressures (Pa)
     77      REAL, INTENT(IN) :: pdpsrf(ngrid)         ! tendency on surface pressure
     78      REAL, INTENT(IN) :: pzlay(ngrid,nlay)     ! altitude at the middle of the layers
     79      REAL, INTENT(IN) :: pt(ngrid,nlay)        ! temperature at the middle of the layers (K)
     80      REAL, INTENT(IN) :: pdt(ngrid,nlay)       ! tendency on temperature from other parametrizations
     81      real, INTENT(IN) :: pq(ngrid,nlay,nq)     ! tracers (kg/kg)
     82      real, INTENT(IN) :: pdq(ngrid,nlay,nq)    ! tendencies before condensation  (kg/kg.s-1)
     83      real, intent(out) :: pdqcloudco2(ngrid,nlay,nq) ! tendency due to CO2 condensation (kg/kg.s-1)
     84      real, intent(out) :: pdtcloudco2(ngrid,nlay)    ! tendency on temperature due to latent heat
     85      INTEGER, INTENT(IN) :: nq                 ! number of tracers
     86      REAL, INTENT(IN) :: tau(ngrid,naerkind) ! Column dust optical depth at each point
     87      REAL, INTENT(IN) :: tauscaling(ngrid)   ! Convertion factor for dust amount
     88      REAL, INTENT(OUT) :: rdust(ngrid,nlay)   ! Dust geometric mean radius (m)
     89      real, intent(OUT) :: rice(ngrid,nlay)    ! Water Ice mass mean radius (m)
    9090                                ! used for nucleation of CO2 on ice-coated ccns
    91       DOUBLE PRECISION rho_ice_co2T(ngrid,nlay) !T-dependant CO2 ice density
    92       DOUBLE PRECISION :: myT   ! temperature scalar for co2 density computation
    93       REAL tau(ngrid,naerkind) ! Column dust optical depth at each point
    94       REAL tauscaling(ngrid)   ! Convertion factor for dust amount
    95       real rdust(ngrid,nlay)   ! Dust geometric mean radius (m)
    96 
    97 c   Outputs:
    98 c   -------
    99 
    100       real pdqcloudco2(ngrid,nlay,nq) ! tendence de la condensation H2O(kg/kg.s-1)
    101       REAL pdtcloudco2(ngrid,nlay)    ! tendence temperature due
    102                                    ! a la chaleur latente
    103       DOUBLE PRECISION riceco2(ngrid,nlay)    ! Ice mass mean radius (m)
     91      DOUBLE PRECISION, INTENT(out) :: riceco2(ngrid,nlay)    ! Ice mass mean radius (m)
    10492                               ! (r_c in montmessin_2004)
    105       REAL nuice(ngrid,nlay)   ! Estimated effective variance
     93      REAL, INTENT(in) :: nuice(ngrid,nlay)   ! Estimated effective variance
    10694                               !   of the size distribution
    107       real rsedcloudco2(ngrid,nlay) ! Cloud sedimentation radius
    108       real rhocloudco2(ngrid,nlay)  ! Cloud density (kg.m-3)
    109       real rhocloudco2t(ngrid,nlay)  ! Cloud density (kg.m-3)
    110       real pdqs_sedco2(ngrid) ! CO2 flux at the surface
     95      real, intent(out) :: rsedcloudco2(ngrid,nlay) ! Cloud sedimentation radius
     96      real, intent(out) :: rhocloudco2(ngrid,nlay)  ! Cloud density (kg.m-3)
     97      real, intent(out) :: rsedcloud(ngrid,nlay) ! Water Cloud sedimentation radius
     98      real, intent(out) :: rhocloud(ngrid,nlay)  ! Water Cloud density (kg.m-3)
     99      real, intent(in) :: pzlev(ngrid,nlay+1) ! altitude at the boundaries of the layers
     100      real, intent(out) :: pdqs_sedco2(ngrid) ! CO2 flux at the surface
     101      REAL, INTENT(IN) :: pdu(ngrid,nlay),pu(ngrid,nlay) !Zonal Wind: zu=pu+pdu*ptimestep
    111102
    112103c   local:
    113104c   ------
    114       !water
    115       real rsedcloud(ngrid,nlay) ! Cloud sedimentation radius
    116       real rhocloud(ngrid,nlay)  ! Cloud density (kg.m-3)
    117       ! for ice radius computation
    118105       
    119106      ! for time loop
    120107      INTEGER microstep  ! time subsampling step variable
    121       INTEGER imicroco2     ! time subsampling for coupled water microphysics & sedimentation
    122       SAVE imicroco2
    123       REAL microtimestep ! integration timestep for coupled water microphysics & sedimentation
    124       SAVE microtimestep
     108      INTEGER, SAVE :: imicroco2     ! time subsampling for coupled water microphysics & sedimentation
     109      REAL, SAVE :: microtimestep ! integration timestep for coupled water microphysics & sedimentation
    125110     
    126111      ! tendency given by clouds (inside the micro loop)
    … …  
    135120      REAL satuco2(ngrid,nlay)  ! co2 satu ratio for output
    136121      REAL zqsatco2(ngrid,nlay) ! saturation co2
     122
     123      DOUBLE PRECISION rho_ice_co2T(ngrid,nlay) !T-dependant CO2 ice density
     124      DOUBLE PRECISION :: myT   ! temperature scalar for co2 density computation
    137125
    138126      INTEGER iq,ig,l,i
    … …  
    154142!     What we need for Qext reading and tau computation : size distribution
    155143      DOUBLE PRECISION vrat_cld ! Volume ratio
    156       DOUBLE PRECISION rb_cldco2(nbinco2_cld+1) ! boundary values of each rad_cldco2 bin (m)
    157       SAVE rb_cldco2
     144      DOUBLE PRECISION, SAVE :: rb_cldco2(nbinco2_cld+1) ! boundary values of each rad_cldco2 bin (m)
    158145      DOUBLE PRECISION, PARAMETER :: rmin_cld = 1.e-9 ! Minimum radius (m)
    159146      DOUBLE PRECISION, PARAMETER :: rmax_cld = 5.e-6 ! Maximum radius (m)
    … …  
    162149      DOUBLE PRECISION dr_cld(nbinco2_cld) ! width of each rad_cldco2 bin (m)
    163150      DOUBLE PRECISION vol_cld(nbinco2_cld) ! particle volume for each bin (m3)
    164       REAL sigma_iceco2 ! Variance of the ice and CCN distributions
     151      REAL, SAVE :: sigma_iceco2 ! Variance of the ice and CCN distributions
    165152      logical :: file_ok !Qext file reading
    166153      double precision :: radv(10000),Qextv1mic(10000)
    167       double precision :: Qext1bins(100),Qtemp
     154      double precision, save :: Qext1bins(100)
     155      double precision :: Qtemp
    168156      double precision :: ltemp1(10000),ltemp2(10000)
    169157      integer :: nelem,lebon1,lebon2,uQext
    170       save Qext1bins
    171158      DOUBLE PRECISION n_aer(nbinco2_cld),Rn,No,n_derf,dev2
    172159      DOUBLE PRECISION Qext1bins2(ngrid,nlay)   
    … …  
    177164      REAL zt(ngrid,nlay)       ! local value of temperature
    178165      REAL :: zq(ngrid, nlay,nq)
     166
     167      real :: rhocloudco2t(ngrid,nlay)  ! Cloud density (kg.m-3)
    179168
    180169      DOUBLE PRECISION :: tcond(ngrid,nlay) !CO2 condensation temperature
    … …  
    187176      REAL ::  mincloud ! min cloud frac
    188177      DOUBLE PRECISION:: rho,zu,NN,gradT !For Saturation Index computation
    189       REAL :: pdu(ngrid,nlay),pu(ngrid,nlay) !Wind field zu=pu+pdu*ptimestep
    190178      DOUBLE PRECISION :: SatIndex(ngrid,nlay),SatIndexmap(ngrid)
    191179     
    … …  
    378366
    379367         if (satindexco2) then !logical in callphys.def
    380          DO l=12,26
    381             ! layers 12 --> 26 ~ 12->85 km
    382             DO ig=1,ngrid
    383          ! Calcul de N^2 static stability
    384             gradT=(zt(ig,l+1)-zt(ig,l))/(pzlev(ig,l+1)-pzlev(ig,l))
    385             NN=sqrt(g/zt(iq,l)*(max(gradT,-g/cpp)+g/cpp))
    386             !calcul of wind field
    387             zu=pu(ig,l)  + pdu(ig,l)*ptimestep
    388 !     calcul of background density
    389             rho=pplay(ig,l)/(rnew(ig,l)*zt(ig,l))
    390             !saturation index
    391             SatIndex(ig,l)=sqrt(7.5e-7*150.e3/(2.*pi)*NN/(rho*zu*zu*zu))
    392          enddo
    393          enddo
     368           DO l=12,26
     369             ! layers 12 --> 26 ~ 12->85 km
     370             DO ig=1,ngrid
     371             ! compute N^2 static stability
     372             gradT=(zt(ig,l+1)-zt(ig,l))/(pzlev(ig,l+1)-pzlev(ig,l))
     373             NN=sqrt(g/zt(iq,l)*(max(gradT,-g/cpp)+g/cpp))
     374             ! compute absolute value of zonal wind field
     375             zu=abs(pu(ig,l)  + pdu(ig,l)*ptimestep)
     376             ! compute background density
     377             rho=pplay(ig,l)/(rnew(ig,l)*zt(ig,l))
     378             !saturation index
     379             SatIndex(ig,l)=sqrt(7.5e-7*150.e3/(2.*pi)*NN/
     380     &                                               (rho*zu*zu*zu))
     381             ENDDO
     382           ENDDO
    394383            !Then compute Satindex map
    395384            ! layers 12 --> 26 ~ 12->85 km
    396          DO ig=1,ngrid
    397             SatIndexmap(ig)=maxval(SatIndex(ig,12:26))
    398          enddo
    399 
    400          call WRITEDIAGFI(ngrid,"SatIndexmap","SatIndexmap","km",2,
    401      &        SatIndexmap)
     385           DO ig=1,ngrid
     386             SatIndexmap(ig)=maxval(SatIndex(ig,12:26))
     387           ENDDO
     388
     389           call WRITEDIAGFI(ngrid,"SatIndexmap","SatIndexmap","km",2,
     390     &         SatIndexmap)
    402391         else
    403             DO ig=1,ngrid
    404                SatIndexmap(ig)=0.05 !maxval(SatIndex(ig,12:26))
    405             enddo
    406          endif
     392           do ig=1,ngrid
     393             SatIndexmap(ig)=0.05 !maxval(SatIndex(ig,12:26))
     394           enddo
     395         endif ! of if (satindexco2)
    407396
    408397!Modulate the DeltaT by GW propagation index :
    … …  
    414403       
    415404         CALL tcondco2(ngrid,nlay,pplay,zqvap,tcond)
    416 ! A tester:            CALL tcondco2(ngrid,nlay,pplay,zqvap,tcond)
    417405         zdelt=spant 
    418            DO ig=1,ngrid
     406         DO ig=1,ngrid
    419407             
    420             if (SatIndexmap(ig) .le. 0.1) THEN
    421               DO l=1,nlay-1
     408           IF (SatIndexmap(ig) .le. 0.1) THEN
     409             DO l=1,nlay-1
    422410         
    423411               IF (tcond(ig,l) .ge. (zt(ig,l)+zdelt)
    424      &             .or. tcond(ig,l) .le. 0 ) THEN !Toute la fraction est saturée
     412     &             .or. tcond(ig,l) .le. 0 ) THEN !The entire fraction is saturated
    425413                  zteff(ig,l)=zt(ig,l)
    426414                  cloudfrac(ig,l)=1.
    427                ELSE IF (tcond(ig,l) .le. (zt(ig,l)-zdelt)) THEN !Rien n'est saturé
     415               ELSE IF (tcond(ig,l) .le. (zt(ig,l)-zdelt)) THEN ! No saturation at all
    428416                  zteff(ig,l)=zt(ig,l)-zdelt
    429417                  cloudfrac(ig,l)=mincloud
    … …  
    431419                  cloudfrac(ig,l)=(tcond(ig,l)-zt(ig,l)+zdelt)/
    432420     &                 (2.0*zdelt)
    433                   zteff(ig,l)=(tcond(ig,l)+zt(ig,l)-zdelt)/2. !Temperature moyenne de la fraction nuageuse
     421                  zteff(ig,l)=(tcond(ig,l)+zt(ig,l)-zdelt)/2. !Mean temperature of the cloud fraction
    434422               END IF           !ig if (tcond(ig,l) ...
    435423               zteff(ig,l)=zteff(ig,l)-pdt(ig,l)*ptimestep
    … …  
    439427                  cloudfrac(ig,l)=1.
    440428               END IF
    441             ENDDO
    442          ELSE
     429             ENDDO
     430           ELSE
    443431!SatIndex not favorable for GW : leave pt untouched
    444             zteff(ig,l)=pt(ig,l)
    445             cloudfrac(ig,l)=mincloud
    446          END IF                 ! of if(SatIndexmap...
    447       ENDDO
     432             zteff(ig,l)=pt(ig,l)
     433             cloudfrac(ig,l)=mincloud
     434           ENDIF                 ! of if(SatIndexmap...
     435         ENDDO ! of DO ig=1,ngrid
    448436!     Totalcloud frac of the column missing here
    449437c-----------------------
    … …  
    466454c------ Temperature tendency subpdt
    467455        ! If imicro=1 subpdt is the same as pdt
    468          DO l=1,nlay
    469             DO ig=1,ngrid
     456        DO l=1,nlay
     457          DO ig=1,ngrid
    470458               subpdt(ig,l) = subpdt(ig,l)
    471459     &              + pdt(ig,l) ! At each micro timestep we add pdt in order to have a stepped entry
    … …  
    500488     &              subpdq(ig,l,igcm_ccn_number)
    501489     &              + pdq(ig,l,igcm_ccn_number)
    502             ENDDO
    503          ENDDO
     490          ENDDO
     491        ENDDO
    504492c- Effective tracers quantities in the cloud fraction
    505          IF (CLFvaryingCO2) THEN     
     493        IF (CLFvaryingCO2) THEN     
    506494            pqeff(:,:,:)=pq(:,:,:) ! prevent from buggs (A. Pottier)
    507495            pqeff(:,:,igcm_ccnco2_mass) =pq(:,:,igcm_ccnco2_mass)/
    … …  
    511499            pqeff(:,:,igcm_co2_ice)= pq(:,:,igcm_co2_ice)/
    512500     &           cloudfrac(:,:)
    513          ELSE
     501        ELSE
    514502            pqeff(:,:,:)=pq(:,:,:)
    515          END IF 
     503        END IF 
    516504       
    517505c------------------------------------------------------
    … …  
    519507c   call to sedimentation routine, update tendancies
    520508c------------------------------------------------------
    521        DO l=1, nlay
     509        DO l=1, nlay
    522510          DO ig=1,ngrid             
    523511             tempo_traceur_t(ig,l)=zteff(ig,l)+subpdt(ig,l)
    … …  
    548536     &            riceco2(ig,l))
    549537          ENDDO
    550        ENDDO
     538        ENDDO
    551539!     Gravitational sedimentation       
    552        sav_trac(:,:,igcm_co2_ice)=tempo_traceurs(:,:,igcm_co2_ice)
    553        sav_trac(:,:,igcm_ccnco2_mass)=
     540        sav_trac(:,:,igcm_co2_ice)=tempo_traceurs(:,:,igcm_co2_ice)
     541        sav_trac(:,:,igcm_ccnco2_mass)=
    554542     &      tempo_traceurs(:,:,igcm_ccnco2_mass)
    555        sav_trac(:,:,igcm_ccnco2_number)=
     543        sav_trac(:,:,igcm_ccnco2_number)=
    556544     &      tempo_traceurs(:,:,igcm_ccnco2_number)
    557545       !We save actualized tracer values to compute sedimentation tendancies
    558        call newsedim(ngrid,nlay,ngrid*nlay,ngrid*nlay,
     546        call newsedim(ngrid,nlay,ngrid*nlay,ngrid*nlay,
    559547     &     microtimestep,pplev,masse,epaisseur,tempo_traceur_t,
    560548     &     rsedcloudco2,rhocloudco2t,
    561549     &     tempo_traceurs(:,:,igcm_co2_ice),wq,beta) !  3 traceurs
    562550!     sedim at the surface of co2 ice : keep track of it for physiq_mod
    563       do ig=1,ngrid
    564          pdqs_sedco2(ig)=pdqs_sedco2(ig)+ wq(ig,1)/microtimestep
    565       end do
    566       call newsedim(ngrid,nlay,ngrid*nlay,ngrid*nlay,
     551        do ig=1,ngrid
     552          pdqs_sedco2(ig)=pdqs_sedco2(ig)+ wq(ig,1)/microtimestep
     553        end do
     554        call newsedim(ngrid,nlay,ngrid*nlay,ngrid*nlay,
    567555     &     microtimestep,pplev,masse,epaisseur,tempo_traceur_t,
    568556     &     rsedcloudco2,rhocloudco2t,
    569557     &     tempo_traceurs(:,:,igcm_ccnco2_mass),wq,beta)
    570       call newsedim(ngrid,nlay,ngrid*nlay,ngrid*nlay,
     558        call newsedim(ngrid,nlay,ngrid*nlay,ngrid*nlay,
    571559     &     microtimestep,pplev,masse,epaisseur,tempo_traceur_t,
    572560     &     rsedcloudco2,rhocloudco2t,
    573561     &     tempo_traceurs(:,:,igcm_ccnco2_number),wq,beta)
    574       DO l = 1, nlay            !Compute tendencies
    575          DO ig=1,ngrid
     562        DO l = 1, nlay            !Compute tendencies
     563          DO ig=1,ngrid
    576564            pdqsed(ig,l,igcm_ccnco2_mass)=
    577565     &           (tempo_traceurs(ig,l,igcm_ccnco2_mass)-
    … …  
    583571     &           (tempo_traceurs(ig,l,igcm_co2_ice)-
    584572     &           sav_trac(ig,l,igcm_co2_ice))/microtimestep
    585          ENDDO
    586       ENDDO
     573          ENDDO
     574        ENDDO
    587575      !update subtimestep tendencies with sedimentation input
    588576        DO l=1,nlay
    … …  
    598586     &           +pdqsed(ig,l,igcm_co2_ice)
    599587         ENDDO
    600       ENDDO   
     588        ENDDO   
    601589c------------------------------------------------------
    602590c      2.  Main call to the cloud schemes:
    603591c------------------------------------------------------
    604       CALL improvedCO2clouds(ngrid,nlay,microtimestep,
     592        CALL improvedCO2clouds(ngrid,nlay,microtimestep,
    605593     &     pplay,pplev,zteff,subpdt,
    606594     &     pqeff,subpdq,subpdqcloudco2,subpdtcloudco2,
    … …  
    609597c      3.  Updating tendencies after cloud scheme:
    610598c-----------------------------------------------------
    611           DO l=1,nlay
    612             DO ig=1,ngrid
     599        DO l=1,nlay
     600          DO ig=1,ngrid
    613601               subpdt(ig,l) =
    614602     &              subpdt(ig,l) + subpdtcloudco2(ig,l)
    … …  
    644632     &              subpdq(ig,l,igcm_ccn_number)
    645633     &              + subpdqcloudco2(ig,l,igcm_ccn_number)
    646             ENDDO
    647          ENDDO
     634          ENDDO
     635        ENDDO
    648636      ENDDO                     ! of DO microstep=1,imicro
    649637     
    … …  
    656644      enddo
    657645c------ Temperature tendency pdtcloud
    658        DO l=1,nlay
    659          DO ig=1,ngrid
     646      DO l=1,nlay
     647        DO ig=1,ngrid
    660648             pdtcloudco2(ig,l) =
    661649     &         subpdt(ig,l)/real(imicroco2)-pdt(ig,l)
    662           ENDDO
    663        ENDDO
     650        ENDDO
     651      ENDDO
    664652c------ Tracers tendencies pdqcloud
    665        DO l=1,nlay
    666           DO ig=1,ngrid       
     653      DO l=1,nlay
     654        DO ig=1,ngrid       
    667655             pdqcloudco2(ig,l,igcm_co2_ice) =
    668656     &            subpdq(ig,l,igcm_co2_ice)/real(imicroco2)
    … …  
    674662     &            subpdq(ig,l,igcm_h2o_ice)/real(imicroco2)
    675663     &            - pdq(ig,l,igcm_h2o_ice)
    676           ENDDO
    677        ENDDO
    678        DO l=1,nlay
    679           DO ig=1,ngrid
     664        ENDDO
     665      ENDDO
     666      DO l=1,nlay
     667        DO ig=1,ngrid
    680668             pdqcloudco2(ig,l,igcm_ccnco2_mass) =
    681669     &            subpdq(ig,l,igcm_ccnco2_mass)/real(imicroco2)
    … …  
    690678     &            subpdq(ig,l,igcm_ccn_number)/real(imicroco2)
    691679     &            - pdq(ig,l,igcm_ccn_number)
    692           ENDDO
    693        ENDDO
    694        DO l=1,nlay
    695           DO ig=1,ngrid
     680        ENDDO
     681      ENDDO
     682      DO l=1,nlay
     683        DO ig=1,ngrid
    696684             pdqcloudco2(ig,l,igcm_dust_mass) =
    697685     &            subpdq(ig,l,igcm_dust_mass)/real(imicroco2)
    … …  
    700688     &            subpdq(ig,l,igcm_dust_number)/real(imicroco2)
    701689     &            - pdq(ig,l,igcm_dust_number)
    702           ENDDO
    703        ENDDO
     690        ENDDO
     691      ENDDO
    704692c-------Due to stepped entry, other processes tendencies can add up to negative values
    705693c-------Therefore, enforce positive values and conserve mass
    706        DO l=1,nlay
    707           DO ig=1,ngrid
     694      DO l=1,nlay
     695        DO ig=1,ngrid
    708696             IF ((pqeff(ig,l,igcm_ccnco2_number) +
    709697     &            ptimestep* (pdq(ig,l,igcm_ccnco2_number) +
    … …  
    725713     &               -pdqcloudco2(ig,l,igcm_ccnco2_mass)
    726714             ENDIF
    727           ENDDO
    728        ENDDO
    729        DO l=1,nlay
    730           DO ig=1,ngrid
     715        ENDDO
     716      ENDDO
     717      DO l=1,nlay
     718        DO ig=1,ngrid
    731719             IF ( (pqeff(ig,l,igcm_dust_number) +
    732720     &            ptimestep* (pdq(ig,l,igcm_dust_number) +
    … …  
    747735     &               -pdqcloudco2(ig,l,igcm_dust_mass)
    748736             ENDIF
    749           ENDDO
    750        ENDDO
     737        ENDDO
     738      ENDDO
    751739        !pq+ptime*(pdq+pdqc)=1 ! pdqc=1-pq/ptime-pdq     
    752        DO l=1,nlay
    753           DO ig=1,ngrid
     740      DO l=1,nlay
     741        DO ig=1,ngrid
    754742             IF (pqeff(ig,l,igcm_co2_ice) + ptimestep*
    755743     &       (pdq(ig,l,igcm_co2_ice) + pdqcloudco2(ig,l,igcm_co2_ice))
    … …  
    766754           pdqcloudco2(ig,l,igcm_co2_ice)= -pdqcloudco2(ig,l,igcm_co2)
    767755          ENDIF
    768          ENDDO
    769         ENDDO
     756        ENDDO
     757      ENDDO
    770758
    771759c Update clouds parameters values in the cloud fraction (for output)
    772         DO l=1, nlay
    773            DO ig=1,ngrid
     760      DO l=1, nlay
     761        DO ig=1,ngrid
    774762
    775763              Niceco2=pqeff(ig,l,igcm_co2_ice) +                   
    … …  
    801789     &             Nccnco2*tauscaling(ig) .le. 1.) )THEN
    802790                 riceco2(ig,l)=0.
     791                 Qext1bins2(ig,l)=0.
     792              else
     793c     Compute opacities
     794                No=Nccnco2*tauscaling(ig)
     795                Rn=-dlog(riceco2(ig,l))
     796                n_derf = derf( (rb_cldco2(1)+Rn) *dev2)
     797                Qext1bins2(ig,l)=0.
     798                do i = 1, nbinco2_cld
     799                 n_aer(i) = -0.5 * No * n_derf !! this ith previously computed
     800                 n_derf = derf((rb_cldco2(i+1)+Rn) *dev2)
     801                 n_aer(i) = n_aer(i) + 0.5 * No * n_derf
     802                 Qext1bins2(ig,l)=Qext1bins2(ig,l)+Qext1bins(i)*n_aer(i)
     803                enddo
    803804              endif
    804 c     Compute opacities
    805            No=Nccnco2*tauscaling(ig)
    806            Rn=-dlog(riceco2(ig,l))
    807            n_derf = derf( (rb_cldco2(1)+Rn) *dev2)
    808            Qext1bins2(ig,l)=0.
    809            do i = 1, nbinco2_cld
    810               n_aer(i) = -0.5 * No * n_derf !! this ith previously computed
    811               n_derf = derf((rb_cldco2(i+1)+Rn) *dev2)
    812               n_aer(i) = n_aer(i) + 0.5 * No * n_derf
    813               Qext1bins2(ig,l)=Qext1bins2(ig,l)+Qext1bins(i)*n_aer(i)
    814            enddo
    815805     
    816806      !update rice water
    817         call updaterice_micro(
     807          call updaterice_micro(
    818808     &    pqeff(ig,l,igcm_h2o_ice) +                    ! ice mass
    819809     &   (pdq(ig,l,igcm_h2o_ice) +                   ! ice mass
    … …  
    827817     &    tauscaling(ig),rice(ig,l),rhocloud(ig,l))
    828818         
    829         call updaterdust(
     819          call updaterdust(
    830820     &    pqeff(ig,l,igcm_dust_mass) +                   ! dust mass
    831821     &   (pdq(ig,l,igcm_dust_mass) +                  ! dust mass
    … …  
    836826     &    rdust(ig,l))
    837827
    838          ENDDO
    839        ENDDO
     828        ENDDO
     829      ENDDO
    840830     
    841831c     A correction if a lot of subliming CO2 fills the 1st layer FF04/2005
    … …  
    843833c     Then that should not affect the ice particle radius
    844834       
    845        do ig=1,ngrid
    846           if(pdpsrf(ig)*ptimestep.gt.0.9*(pplev(ig,1)-pplev(ig,2)))then
     835      do ig=1,ngrid
     836        if(pdpsrf(ig)*ptimestep.gt.0.9*(pplev(ig,1)-pplev(ig,2)))then
    847837             if(pdpsrf(ig)*ptimestep.gt.0.9*(pplev(ig,1)-pplev(ig,3)))
    848838     &       riceco2(ig,2)=riceco2(ig,3)
    849839             riceco2(ig,1)=riceco2(ig,2)
    850         end if
     840        endif
    851841      end do
    852842       
    … …  
    860850      ENDDO
    861851       
    862        DO l=1,nlay
     852      DO l=1,nlay
    863853         DO ig=1,ngrid
    864854           rsedcloudco2(ig,l)=max(riceco2(ig,l)*
    … …  
    867857c          rsedcloudco2(ig,l)=min(rsedcloudco2(ig,l),1.e-5)
    868858         ENDDO
    869        ENDDO
     859      ENDDO
    870860       
    871        call co2sat(ngrid*nlay,zteff+(pdt+pdtcloudco2)*ptimestep
     861      call co2sat(ngrid*nlay,zteff+(pdt+pdtcloudco2)*ptimestep
    872862     &      ,pplay,zqsatco2)
    873        do l=1,nlay
    874           do ig=1,ngrid
     863      do l=1,nlay
     864        do ig=1,ngrid
    875865             satuco2(ig,l) = (pqeff(ig,l,igcm_co2)  +
    876866     &            (pdq(ig,l,igcm_co2) +
    877867     &            pdqcloudco2(ig,l,igcm_co2))*ptimestep)*
    878868     &            (mmean(ig,l)/44.01)*pplay(ig,l)/zqsatco2(ig,l)
    879             enddo
    880          enddo
    881 !Tout ce qui est modifié par la microphysique de CO2 doit être rapporté a cloudfrac
    882          IF (CLFvaryingCO2) THEN
     869        enddo
     870      enddo
     871!Everything modified by CO2 microphysics must be wrt cloudfrac
     872      IF (CLFvaryingCO2) THEN
    883873            DO l=1,nlay
    884874               DO ig=1,ngrid
    … …  
    905895               ENDDO
    906896            ENDDO   
    907          ENDIF
    908 !l'opacité de la case ig est la somme sur l de Qext1bins2: Est-ce-vrai?
    909          tau1mic(:)=0.
    910          do l=1,nlay
    911             do ig=1,ngrid
    912                tau1mic(ig)=tau1mic(ig)+Qext1bins2(ig,l)
    913             enddo
    914          enddo
     897      ENDIF
     898! opacity in mesh ig is the sum over l of Qext1bins2: Is this true ?
     899      tau1mic(:)=0.
     900      do l=1,nlay
     901        do ig=1,ngrid
     902          tau1mic(ig)=tau1mic(ig)+Qext1bins2(ig,l)
     903        enddo
     904      enddo
    915905!Outputs:
    916          call WRITEDIAGFI(ngrid,"SatIndex","SatIndex"," ",3,
     906      call WRITEDIAGFI(ngrid,"SatIndex","SatIndex"," ",3,
    917907     &        SatIndex)
    918          call WRITEDIAGFI(ngrid,"satuco2","vap in satu","kg/kg",3,
     908      call WRITEDIAGFI(ngrid,"satuco2","vap in satu","kg/kg",3,
    919909     &        satuco2)
    920          call WRITEdiagfi(ngrid,"riceco2","ice radius","m"
     910      call WRITEdiagfi(ngrid,"riceco2","ice radius","m"
    921911     &        ,3,riceco2)         
    922          call WRITEdiagfi(ngrid,"cloudfrac","co2 cloud fraction"
     912      call WRITEdiagfi(ngrid,"cloudfrac","co2 cloud fraction"
    923913     &        ," ",3,cloudfrac)
    924          call WRITEdiagfi(ngrid,"rsedcloudco2","rsed co2"
     914      call WRITEdiagfi(ngrid,"rsedcloudco2","rsed co2"
    925915     &        ,"m",3,rsedcloudco2)
    926          call WRITEdiagfi(ngrid,"Tau3D1mic"," co2 ice opacities"
     916      call WRITEdiagfi(ngrid,"Tau3D1mic"," co2 ice opacities"
    927917     &        ," ",3,Qext1bins2)
    928          call WRITEdiagfi(ngrid,"tau1mic","co2 ice opacity 1 micron"
     918      call WRITEdiagfi(ngrid,"tau1mic","co2 ice opacity 1 micron"
    929919     &        ," ",2,tau1mic)
    930920         
    931          END
     921      END
Note: See TracChangeset for help on using the changeset viewer.