Changeset 2350


Ignore:
Timestamp:
Jun 9, 2020, 3:27:39 PM (5 years ago)
Author:
aslmd
Message:

CO2 cloud microphysics

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

Legend:

Unmodified
Added
Removed

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

    r2349 r2350  
    1       SUBROUTINE co2sat(naersize,t,p,psat)
    2 c       SUBROUTINE co2sat(naersize,t,p,qsat) JA
    3       IMPLICIT NONE
     1c=======================================================================
     2c     SUBROUTINE co2sat
     3c-----------------------------------------------------------------------
     4c     Aim:
     5c     ----
     6c     Compute saturated steam pressure (from James et al, 1992)
     7c=======================================================================
     8      subroutine co2sat(naersize, t, psat)
     9     
     10      implicit none
     11c-----------------------------------------------------------------------
     12c     VARIABLES
     13c-----------------------------------------------------------------------
     14c     Inputs:
     15c     -------
     16      integer, intent(in) ::
     17     &   naersize ! dimension of tables t and psat
     18     
     19      real, intent(in) ::
     20     &   t(naersize) ! temperature table
    421
     22c     Output:
     23c     -------
     24      real, intent(out) ::
     25     &   psat(naersize) ! Saturated steam pressure (Pa)
     26
     27c     Local:
     28c     ------
     29      integer ::
     30     &   i ! loop on naersize
    531c=======================================================================
    6 c
    7 c
    8 c  now:  straight psat of CO2 (or qsat of CO2 but need of mmean)
    9 c
     32c===== BEGIN
    1033c=======================================================================
    11 
    12 c   declarations:
    13 c   -------------
    14 c   arguments:
    15 c   ----------
    16 
    17 c   INPUT
    18       integer naersize
    19       real t(naersize) , p(naersize)
    20 c   OUTPUT
    21 c      real qsat(naersize) JA
    22       real psat(naersize)
    23 
    24 c   local:
    25 c   ------
    26       INTEGER i
    27       REAL r2,r3,r4 , To, es
    28       SAVE r2,r3,r4
    29       DATA r2,r3,r4/611.14,21.875,7.66/
    30       SAVE To
    31       DATA To/273.16/
    32          
    33       do i=1,naersize
    34 
    35 
    36 c        pression de vapeur saturante (James et al. 1992):
    37 
    38           psat(i)  = 1.382 * 1e12 * exp(-3182.48/t(i)) !; (Pa)
    39 
    40 c         OR:
    41 
    42 c         qsat(i) = psat/p(i)*44.01/mmean ! Need of updated information on mmean
    43 c         qsat(i) = max(qsat(i), 1.e-30)
    44 
    45 
    46       enddo
    47 c      qsat=psat JA
    48          
    49 
     34      do i = 1, naersize
     35        psat(i) = 1.382 * 1e12 * exp(-3182.48/t(i))
     36      end do
     37c=======================================================================
     38c===== END
     39c=======================================================================
    5040      RETURN
    5141      END

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

    r2349 r2350  
    1       SUBROUTINE co2snow (ngrid,nlayer,ptimestep,emisref,condsub
    2      &         ,pplev,pcondicea,pcondices,pfallice,pemisurf)
     1c=======================================================================
     2c     Program for simulate the impact of the CO2 snow fall on
     3c     the surface infrared emission (emissivity)
     4c-----------------------------------------------------------------------
     5c     Algorithme:
     6c     1 - Initialization
     7c     2 - Compute the surface emissivity
     8c-----------------------------------------------------------------------
     9c     Author: F. Forget
     10c-----------------------------------------------------------------------
     11c     Reference paper:
     12c     Francois Forget and James B. Pollack
     13c
     14c     "Thermal infrared observations of the condensing Martian polar
     15c     caps: CO2 ice temperatures and radiative budget"
     16c
     17c     JOURNAL OF GEOPHYSICAL RESEARCH, VOL. 101, NO. E7, PAGES 16,
     18c     865-16,879, JULY 25, 1996
     19c=======================================================================
    320
     21      SUBROUTINE co2snow(ngrid, nlayer, ptimestep, emisref, condsub,
     22     &                   pplev, pcondicea, pcondices, pfallice,
     23     &                   pemisurf)
     24     
    425      use surfdat_h, only: iceradius, dtemisice
    526      use geometry_mod, only: latitude ! grid point latitudes (rad)
    627      use time_phylmdz_mod, only: daysec
     28
    729      IMPLICIT NONE
    8 
    9 c=======================================================================
    10 c     Program for simulate the impact of the CO2 snow fall on
    11 c     the surface infrared emission (emissivity)  and on
    12 c     the airborne dust
    13 c     F.Forget 1996
    14 c=======================================================================
    15 
    16 c Declarations
    17 c ------------
    1830
    1931#include "callkeys.h"
    2032
    21 c     Arguments
    22 c     ---------
    2333
    24       INTEGER,INTENT(IN) :: ngrid ! number of atmospheric columns
    25       INTEGER,INTENT(IN) :: nlayer ! number of atmospheric layers
    26       REAL,INTENT(IN) :: ptimestep ! timestep of the physics (s)
    27       REAL,INTENT(IN) :: emisref(ngrid) ! grd or ice  emissivity without snow
    28       logical,intent(in) :: condsub(ngrid) ! true if there is CO2 condensation
    29                                            ! or sublimation in the column
    30       REAL,INTENT(IN) :: pplev(ngrid,nlayer+1) ! interlayer pressure (Pa)
    31       REAL,INTENT(IN) :: pcondicea(ngrid,nlayer) ! CO2 condensation rate (kg/m2/s)
    32       REAL,INTENT(IN) :: pcondices(ngrid) ! CO2 condensation rate (kg/m2/s)
    33                                           ! on the surface
    34       REAL,INTENT(IN) :: pfallice(ngrid,nlayer+1) ! falling CO2 ice (kg/m2/s)
     34c=======================================================================
     35c     Variables
     36c=======================================================================
     37c     Inputs:
     38c     -------
     39      integer, intent(in) ::
     40     &   ngrid, ! number of atmospheric columns
     41     &   nlayer ! number of atmospheric layers
     42                       
     43      real, intent(in) ::
     44     &   ptimestep,               ! timestep of the physics (s)
     45     &   emisref(ngrid),          ! grd or ice emissivity without snow
     46     &   pplev(ngrid,nlayer+1),   ! interlayer pressure (Pa)
     47     &   pcondicea(ngrid,nlayer), ! CO2 condensation rate in layers
     48                                  !   (kg/m2/s)
     49     &   pcondices(ngrid),        ! CO2 condensation rate on the surface
     50                                  !   (kg/m2/s)       
     51     &   pfallice(ngrid,nlayer+1) ! falling CO2 ice (kg/m2/s)
     52                           
     53      logical, intent(in) ::
     54     &   condsub(ngrid) ! true if there is CO2 condensation or
     55                        ! sublimation in the column
    3556
    36       REAL,INTENT(OUT) :: pemisurf(ngrid) ! surface emissivity
     57c     Output:
     58c     -------
     59      real, intent(out) ::
     60     &   pemisurf(ngrid) ! surface emissivity
    3761
    38 c     local
    39 c     -----
    40       integer ig , l , icap
     62c     local:
     63c     ------
     64      integer ::
     65     &   l,
     66     &   ig,
     67     &   icap
    4168
    42       REAL zdemisurf ,dtemis
    43       REAL sumdaer
     69      real ::
     70     &   sumdaer,
     71     &   zdemisurf
    4472
    45 c     saved
    46 c     -----
    47       REAL,save :: Kscat(2), scaveng
    48       LOGICAL,SAVE :: firstcall=.true.
     73      real, parameter ::
     74     &   alpha = 0.45
    4975
    50 c --------------
    51 c Initialisation
    52 c --------------
     76c     saved:
     77c     ------
     78      real, save ::
     79     &   Kscat(2) ! Kscat: coefficient for decreasing the surface
     80c                 !        emissivity
     81c                 !
     82c                 ! Kscat = (0.001/3.) * alpha / iceradius
     83c                 !          with 0.3 < alpha < 0.6
     84c                 !
     85c                 ! alpha set to 0.45 (coeff from emis = f (tau)) and
     86c                 ! iceradius the mean radius of the scaterring
     87c                 ! particles (200.e-6 < iceradius < 10.e-6)
     88c                 !
     89c                 ! (2) = N and S hemispheres
    5390
    54       ! AS: firstcall OK absolute
     91      logical, save ::
     92     &   firstcall = .true.
     93c=======================================================================
     94c BEGIN
     95c=======================================================================
     96c 1 - Initialization
     97c=======================================================================
    5598      if (firstcall) then
     99        Kscat(1) = (0.001/3.) * alpha / iceradius(1)
     100        Kscat(2) = (0.001/3.) * alpha / iceradius(2)
     101                   
     102        firstcall = .false.
     103      end if
     104c=======================================================================
     105c 2 - Compute the surface emissivity
     106c=======================================================================
     107      do ig = 1, ngrid
     108        if (condsub(ig)) then
     109          if (latitude(ig).lt.0.) then
     110              icap = 2 ! Southern hemisphere
     111          else
     112              icap = 1 ! Northern Hemisphere
     113          end if
    56114
    57 c   Kscat : coefficient for decreasing the surface emissivity
    58 c           =(0.001/3.)*alpha/iceradius ,
    59 c           with 0.3< alpha < 0.6, set to 0.45 (coeff from emis = f (tau))
    60 c           and iceradius the mean radius of the
    61 c           scaterring particles  (200.e-6<iceradius<10.e-6)
    62 
    63             Kscat(1)=(0.001/3.)*0.45/iceradius(1)
    64             Kscat(2)=(0.001/3.)*0.45/iceradius(2)
    65 
    66 c          Scavenging Ratio (dust concentration in the air / in the snow)
    67            scaveng = 100.0
    68            
    69 c          Collision Scavenging coefficient  (m2.kg-1)
    70 c          Csca  = 2.3  ! not used yet !!!!!!!!!!!
    71            firstcall = .false.
    72 
    73       end if
    74 
    75 
    76 c     LOOP on grid points
    77 c     -------------------
    78       do ig=1,ngrid
    79          if (condsub(ig)) then
    80 
    81 c         IF (scavenging) then
    82 c          Airborne Dust
    83 c          -------------
    84 c          sumdaer=0.
    85 c          do l=nlayer, 1, -1
    86 c             pdaerosol(ig,l)= -paerosol(ig,l,1)*
    87 c    &              (1-exp(-scaveng*pcondicea(ig,l)*ptimestep*g/
    88 c    &               (pplev(ig,l)-pplev(ig,l+1))))/ptimestep 
    89 
    90 c    &           - Csca*paerosol(ig,l,1) ! Scavenging by collision
    91 c    &           * 0.5*(pfallice(ig,l)) ! not included
    92 
    93 c             test to avoid releasing to much dust when subliming:
    94 c             if(pdaerosol(ig,l).gt.-sumdaer)pdaerosol(ig,l)=-sumdaer
    95 c             sumdaer=sumdaer + pdaerosol(ig,l)
    96 
    97 c            if (-pdaerosol(ig,l)*ptimestep.gt.paerosol(ig,l,1)) then
    98 c                write(*,*) 'ds co2snow: aerosol < 0.0 !!!'
    99 c                write(*,*) 'ig =' , ig
    100 c            end if
    101 c          end do
    102 c         END IF
    103 
    104 c          Surface emissivity
    105 c          ------------------
    106 c   dtemis: Time scale for increasing the ice emissivity
    107 
    108            IF(latitude(ig).LT. 0.) THEN
    109               icap=2 ! Southern hemisphere
    110            ELSE
    111               icap=1 ! Northern Hemisphere
    112            ENDIF
    113 
    114            zdemisurf =
    115      &    (emisref(ig)-pemisurf(ig))/(dtemisice(icap)*daysec)
    116 c Using directly the diferential equation:
    117 c    &       -Kscat(icap)*emisref(ig)*
    118 c    &      (pemisurf(ig)/emisref(ig))**4 *pfallice(ig,1)
    119 c Using an integrated form for numerical safety instead
    120      & +(emisref(ig)* ((pemisurf(ig)/emisref(ig))**(-3)+3.*Kscat(icap)*
    121      & pfallice(ig,1)*ptimestep)**(-1/3.) -pemisurf(ig))/ptimestep   
    122 
    123 
    124            pemisurf(ig) = pemisurf(ig) + zdemisurf*ptimestep
    125 
    126            if (pemisurf(ig).lt.0.1) then
    127                  write(*,*) 'ds co2snow: emis < 0.1 !!!'
    128                  write(*,*) 'ig =' , ig
    129                  write(*,*)'pemisurf(ig)',pemisurf(ig)
    130                  write(*,*) 'zdemisurf*ptimestep',zdemisurf*ptimestep
    131            end if
    132          else
     115          ! compute zdemisurf using an integrated form for numerical
     116          ! safety instead
     117          zdemisurf =
     118     &     (emisref(ig) - pemisurf(ig)) / (dtemisice(icap) * daysec)
     119     &     +
     120     &     ( emisref(ig) *
     121     &              ( (pemisurf(ig)/emisref(ig))**(-3) +
     122     &                3.*Kscat(icap)*pfallice(ig,1)*ptimestep )**(-1/3.)
     123     &       - pemisurf(ig)
     124     &      ) / ptimestep
     125         
     126          pemisurf(ig) = pemisurf(ig) + zdemisurf * ptimestep
     127          if (pemisurf(ig).lt.0.1) then
     128            write(*,*)'ds co2snow: emis < 0.1 !!!'
     129            write(*,*)'ig =' , ig
     130            write(*,*)'pemisurf(ig)', pemisurf(ig)
     131            write(*,*)'zdemisurf*ptimestep', zdemisurf*ptimestep
     132          end if
     133        else ! if condsub(ig) is false
    133134           pemisurf(ig) = emisref(ig)
    134          end if
    135       end do
    136 
     135        end if
     136      end do ! ngrid
     137c=======================================================================
     138c END
     139c=======================================================================
    137140      return
    138141      end

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

    r2349 r2350  
    1818      use watersat_mod, only: watersat
    1919      use co2condens_mod, only: co2condens
     20      use co2condens_mod4micro, only: co2condens4micro
    2021      use co2cloud_mod, only: co2cloud, mem_Mccn_co2, mem_Mh2o_co2,
    2122     &                        mem_Nccn_co2
    … …  
    251252c     Variables used by the CO2 clouds microphysical scheme:
    252253      DOUBLE PRECISION riceco2(ngrid,nlayer)   ! co2 ice geometric mean radius (m)
    253       real rsedcloudco2(ngrid,nlayer) !CO2 Cloud sedimentation radius
    254       real rhocloudco2(ngrid,nlayer) !co2 Cloud density (kg.m-3)
    255254      real zdqssed_co2(ngrid)  ! CO2 flux at the surface  (kg.m-2.s-1)
    256255      real zcondicea_co2microp(ngrid,nlayer)
    … …  
    303302
    304303c     Tendancies due to various processes:
    305       REAL dqsurf(ngrid,nq)
     304      REAL dqsurf(ngrid,nq)  ! tendency for tracers on surface (Kg/m2/s)
    306305      REAL zdtlw(ngrid,nlayer)     ! (K/s)
    307306      REAL zdtsw(ngrid,nlayer)     ! (K/s)
    … …  
    384383      REAL mtot(ngrid)          ! Total mass of water vapor (kg/m2)
    385384      REAL mstormdtot(ngrid)    ! Total mass of stormdust tracer (kg/m2)
    386       REAL mdusttot(ngrid)    ! Total mass of dust tracer (kg/m2)
     385      REAL mdusttot(ngrid)      ! Total mass of dust tracer (kg/m2)
    387386      REAL icetot(ngrid)        ! Total mass of water ice (kg/m2)
    388       REAL mtotco2(ngrid)       ! Total mass of co2 vapor (kg/m2)
    389       REAL icetotco2(ngrid)     ! Total mass of co2 ice (kg/m2)
     387      REAL mtotco2       ! Total mass of co2, including ice at the surface (kg/m2)
     388      REAL vaptotco2     ! Total mass of co2 vapor (kg/m2)
     389      REAL icetotco2     ! Total mass of co2 ice (kg/m2)
    390390      REAL Nccntot(ngrid)       ! Total number of ccn (nbr/m2)
     391      REAL NccnCO2tot(ngrid)    ! Total number of ccnCO2 (nbr/m2)
    391392      REAL Mccntot(ngrid)       ! Total mass of ccn (kg/m2)
    392393      REAL rave(ngrid)          ! Mean water ice effective radius (m)
    … …  
    403404      REAL DoH_ice(ngrid,nlayer) !D/H ratio
    404405      REAL DoH_surf(ngrid) !D/H ratio surface
    405                                
     406
    406407      REAL dqdustsurf(ngrid) ! surface q dust flux (kg/m2/s)
    407408      REAL dndustsurf(ngrid) ! surface n dust flux (number/m2/s)
    … …  
    483484      real tf_clf, ntf_clf ! tf: fraction of clouds, ntf: fraction without clouds
    484485      real rave2(ngrid), totrave2(ngrid) ! Mean water ice mean radius (m)
     486C test de conservation de la masse de CO2
     487      REAL co2totA
     488      REAL co2totB
    485489
    486490c entrainment by slope winds above sb-grid scale topography
    … …  
    496500
    497501c=======================================================================
     502c---------- begin added by CM
     503      pdq(:,:,:) = 0.
     504c---------end added by CM
     505
    498506
    499507c 1. Initialisation:
    500508c -----------------
    501 
    502509c  1.1   Initialisation only at first call
    503510c  ---------------------------------------
    … …  
    556563           ! starting without startfi.nc and with callsoil
    557564           ! is not yet possible as soildepth default is not defined
    558            if (callsoil) then                     
     565           if (callsoil) then
    559566              ! default mlayer distribution, following a power law:
    560567              !  mlayer(k)=lay1*alpha**(k-1/2)
    … …  
    725732      dsotop(:,:)=0.
    726733      dwatercap(:)=0
    727      
     734
    728735#ifdef DUSTSTORM
    729736      pq_tmp(:,:,:)=0
    … …  
    810817      enddo
    811818
     819      ! calculates the amount of co2 at the beginning of physics
     820      co2totA = 0.
     821      do ig=1,ngrid
     822        do l=1,nlayer
     823           co2totA = co2totA + (zplev(ig,l)-zplev(ig,l+1))/g*
     824     &           (pq(ig,l,igcm_co2)+pq(ig,l,igcm_co2_ice)
     825     &      +(pdq(ig,l,igcm_co2)+pdq(ig,l,igcm_co2_ice))*ptimestep)
     826        end do
     827        co2totA = co2totA + co2ice(ig)
     828      end do
    812829c-----------------------------------------------------------------------
    813830c    2. Compute radiative tendencies :
    … …  
    852869     &                      pq(1:ngrid,1:nlayer,igcm_o)*
    853870     &                      mmean(1:ngrid,1:nlayer)/mmol(igcm_o)
    854                  
     871
    855872                 CALL nlte_tcool(ngrid,nlayer,zplay*9.869e-6,
    856873     $                pt,zzlay,co2vmr_gcm, n2vmr_gcm, covmr_gcm,
    … …  
    902919c          ------------------
    903920           ! callradite for the part with clouds
    904            clearsky=.false. ! part with clouds for both cases CLFvarying true/false 
     921           clearsky=.false. ! part with clouds for both cases CLFvarying true/false
    905922           CALL callradite(icount,ngrid,nlayer,nq,zday,zls,pq,albedo,
    906923     &     emis,mu0,zplev,zplay,pt,tsurf,fract,dist_sol,igout,
    … …  
    917934               ! if
    918935               ! CLFvarying ...) ?? AP ??
    919                clearsky=.true. ! 
     936               clearsky=.true. !
    920937               CALL callradite(icount,ngrid,nlayer,nq,zday,zls,pq,
    921938     &              albedo,emis,mu0,zplev,zplay,pt,tsurf,fract,
    … …  
    11191136     &                      pdqrds,wspeed,dsodust,dsords,dsotop,
    11201137     &                      tauref)
    1121                    
     1138
    11221139c      update the tendencies of both dust after vertical transport
    11231140         DO l=1,nlayer
    … …  
    11711188           hsummit(:)=14000.
    11721189         endif
    1173          clearatm=.true. ! stormdust is not accounted in the extra heating on top of the mountains 
     1190         clearatm=.true. ! stormdust is not accounted in the extra heating on top of the mountains
    11741191         nohmons=.false.
    1175          pdqtop(:,:,:)=0. 
     1192         pdqtop(:,:,:)=0.
    11761193         CALL topmons(ngrid,nlayer,nq,ptime,ptimestep,
    11771194     &                pq,pdq,pt,pdt,zplev,zplay,zzlev,
    … …  
    12231240        CALL calldrag_noro(ngrid,nlayer,ptimestep,
    12241241     &                 zplay,zplev,pt,pu,pv,zdtgw,zdugw,zdvgw)
    1225  
     1242
    12261243        DO l=1,nlayer
    12271244          DO ig=1,ngrid
    … …  
    13091326            DO l=1,nlayer
    13101327              DO ig=1,ngrid
    1311                  pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqdif(ig,l,iq) 
     1328                 pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqdif(ig,l,iq)
    13121329              ENDDO
    13131330            ENDDO
    … …  
    13641381
    13651382      if(calltherm .and. .not.turb_resolved) then
    1366  
     1383
    13671384        call calltherm_interface(ngrid,nlayer,nq,
    13681385     $ tracer,igcm_co2,
    … …  
    13721389     $ pdu_th,pdv_th,pdt_th,pdq_th,lmax_th,zmax_th,
    13731390     $ dtke_th,zdhdif,hfmax_th,wstar,sensibFlux)
    1374      
     1391
    13751392         DO l=1,nlayer
    13761393           DO ig=1,ngrid
    … …  
    13811398           ENDDO
    13821399        ENDDO
    1383  
     1400
    13841401        DO ig=1,ngrid
    13851402          q2(ig,nlayer+1)=q2(ig,nlayer+1)+dtke_th(ig,nlayer+1)*ptimestep
    1386         ENDDO     
     1403        ENDDO
    13871404
    13881405        if (tracer) then
    … …  
    14931510     &                nq,tau,tauscaling,rdust,rice,nuice,
    14941511     &                rsedcloud,rhocloud,totcloudfrac)
    1495      
     1512
    14961513c Temperature variation due to latent heat release
    14971514           if (activice) then
    1498                pdt(1:ngrid,1:nlayer) = 
    1499      &          pdt(1:ngrid,1:nlayer) + 
     1515               pdt(1:ngrid,1:nlayer) =
     1516     &          pdt(1:ngrid,1:nlayer) +
    15001517     &          zdtcloud(1:ngrid,1:nlayer)
    15011518           endif
    … …  
    15111528           if (hdo) then
    15121529! increment HDO vapour and ice atmospheric tracers tendencies
    1513            pdq(1:ngrid,1:nlayer,igcm_hdo_vap) = 
    1514      &       pdq(1:ngrid,1:nlayer,igcm_hdo_vap) + 
     1530           pdq(1:ngrid,1:nlayer,igcm_hdo_vap) =
     1531     &       pdq(1:ngrid,1:nlayer,igcm_hdo_vap) +
    15151532     &       zdqcloud(1:ngrid,1:nlayer,igcm_hdo_vap)
    1516            pdq(1:ngrid,1:nlayer,igcm_hdo_ice) = 
    1517      &       pdq(1:ngrid,1:nlayer,igcm_hdo_ice) + 
     1533           pdq(1:ngrid,1:nlayer,igcm_hdo_ice) =
     1534     &       pdq(1:ngrid,1:nlayer,igcm_hdo_ice) +
    15181535     &       zdqcloud(1:ngrid,1:nlayer,igcm_hdo_ice)
    15191536           endif !hdo
    … …  
    15881605c                                nucleation
    15891606c               imicroco2=50     micro-timestep is 1/50 of physical timestep
    1590          
    1591          IF (co2clouds ) THEN
    1592 
    1593 
    1594             call co2cloud(ngrid,nlayer,ptimestep,
     1607        zdqssed_co2(:) = 0.
     1608
     1609         IF (co2clouds) THEN
     1610           call co2cloud(ngrid,nlayer,ptimestep,
    15951611     &           zplev,zplay,pdpsrf,zzlay,pt,pdt,
    15961612     &           pq,pdq,zdqcloudco2,zdtcloudco2,
    15971613     &           nq,tau,tauscaling,rdust,rice,riceco2,nuice,
    1598      &           rsedcloudco2,rhocloudco2,
    15991614     &           rsedcloud,rhocloud,zzlev,zdqssed_co2,
    1600      &           pdu,pu,zcondicea_co2microp)
     1615     &           pdu,pu,zcondicea_co2microp, co2ice)
    16011616
    16021617
    … …  
    16121627! We need to check that we have Nccn & Ndust > 0
    16131628! This is due to single precision rounding problems
    1614              
    16151629! increment dust tracers tendancies
    1616                pdq(1:ngrid,1:nlayer,igcm_dust_mass) =
    1617      &              pdq(1:ngrid,1:nlayer,igcm_dust_mass) +
    1618      &              zdqcloudco2(1:ngrid,1:nlayer,igcm_dust_mass)         
    1619                pdq(1:ngrid,1:nlayer,igcm_dust_number) =
    1620      &              pdq(1:ngrid,1:nlayer,igcm_dust_number) +
    1621      &              zdqcloudco2(1:ngrid,1:nlayer,igcm_dust_number)
    1622                pdq(1:ngrid,1:nlayer,igcm_co2) =
    1623      &              pdq(1:ngrid,1:nlayer,igcm_co2) +
    1624      &              zdqcloudco2(1:ngrid,1:nlayer,igcm_co2)
    1625                pdq(1:ngrid,1:nlayer,igcm_co2_ice) =
    1626      &              pdq(1:ngrid,1:nlayer,igcm_co2_ice) +
    1627      &              zdqcloudco2(1:ngrid,1:nlayer,igcm_co2_ice)
    1628                pdq(1:ngrid,1:nlayer,igcm_ccnco2_mass) =
    1629      &              pdq(1:ngrid,1:nlayer,igcm_ccnco2_mass) +
    1630      &              zdqcloudco2(1:ngrid,1:nlayer,igcm_ccnco2_mass)
    1631                pdq(1:ngrid,1:nlayer,igcm_ccnco2_number) =
    1632      &              pdq(1:ngrid,1:nlayer,igcm_ccnco2_number) +
    1633      &              zdqcloudco2(1:ngrid,1:nlayer,igcm_ccnco2_number)
     1630               pdq(:,:,igcm_dust_mass) = pdq(:,:,igcm_dust_mass)
     1631     &                                 + zdqcloudco2(:,:,igcm_dust_mass)
     1632
     1633               pdq(:,:,igcm_dust_number) = pdq(:,:,igcm_dust_number)
     1634     &                               + zdqcloudco2(:,:,igcm_dust_number)
     1635
     1636               pdq(:,:,igcm_co2) = pdq(:,:,igcm_co2)
     1637     &                             + zdqcloudco2(:,:,igcm_co2)
     1638
     1639               pdq(:,:,igcm_co2_ice) = pdq(:,:,igcm_co2_ice)
     1640     &                                 + zdqcloudco2(:,:,igcm_co2_ice)
     1641
     1642               pdq(:,:,igcm_ccnco2_mass) = pdq(:,:,igcm_ccnco2_mass)
     1643     &                               + zdqcloudco2(:,:,igcm_ccnco2_mass)
     1644
     1645               pdq(:,:,igcm_ccnco2_number) = pdq(:,:,igcm_ccnco2_number)
     1646     &                             + zdqcloudco2(:,:,igcm_ccnco2_number)
     1647
    16341648!Update water ice clouds values as well
    16351649             if (co2useh2o) then
    … …  
    16431657     &               pdq(1:ngrid,1:nlayer,igcm_ccn_number) +
    16441658     &               zdqcloudco2(1:ngrid,1:nlayer,igcm_ccn_number)
    1645                 where (pq(:,:,igcm_ccn_mass) +
     1659
     1660c     Negative values?
     1661                where (pq(:,:,igcm_ccn_mass) +
    16461662     &                 ptimestep*pdq(:,:,igcm_ccn_mass) < 0.)
    16471663                  pdq(:,:,igcm_ccn_mass) =
    … …  
    16501666     &               - pq(:,:,igcm_ccn_number)/ptimestep + 1.e-30
    16511667                end where
    1652                 where (pq(:,:,igcm_ccn_number) +
     1668c     Negative values?
     1669                where (pq(:,:,igcm_ccn_number) +
    16531670     &                 ptimestep*pdq(:,:,igcm_ccn_number) < 0.)
    16541671                  pdq(:,:,igcm_ccn_mass) =
    … …  
    16891706     &           - pq(:,:,igcm_dust_number)/ptimestep + 1.e-30
    16901707             end where
    1691      
    1692       END IF                    ! of IF (co2clouds)
     1708      END IF ! of IF (co2clouds)
    16931709
    16941710c   9b. Aerosol particles
    … …  
    17341750     &                tau,tauscaling)
    17351751c Flux at the surface of co2 ice computed in co2cloud microtimestep
    1736            IF (co2clouds) THEN
    1737               zdqssed(1:ngrid,igcm_co2_ice)=zdqssed_co2(1:ngrid)
    1738            ENDIF
    1739 
    17401752           IF (rdstorm) THEN
    17411753c             Storm dust cannot sediment to the surface
    … …  
    18671879          DO ig=1,ngrid
    18681880            dtrad(ig,l)=dtrad(ig,l)+zdteuv(ig,l)
    1869             pdt(ig,l)=pdt(ig,l)+zdtconduc(ig,l)
    1870      &                         +zdteuv(ig,l)
     1881            pdt(ig,l)=pdt(ig,l)+zdtconduc(ig,l)+zdteuv(ig,l)
    18711882            pdv(ig,l)=pdv(ig,l)+zdvmolvis(ig,l)
    18721883            pdu(ig,l)=pdu(ig,l)+zdumolvis(ig,l)
    … …  
    18831894c     (should be the last atmospherical physical process to be computed)
    18841895c   -------------------------------------------
    1885 
    18861896      IF (tituscap) THEN
    18871897         !!! get the actual co2 seasonal cap from Titus observations
    1888          CALL geticecover( ngrid, 180.*zls/pi,
     1898         CALL geticecover(ngrid, 180.*zls/pi,
    18891899     .                  180.*longitude/pi, 180.*latitude/pi, co2ice )
    18901900         co2ice = co2ice * 10000.
    18911901      ENDIF
    18921902     
    1893      
    1894       pdpsrf(:) = 0
    18951903
    18961904      IF (callcond) THEN
    1897          CALL co2condens(ngrid,nlayer,nq,ptimestep,
     1905         zdtc(:,:) = 0.
     1906         zdtsurfc(:) = 0.
     1907         zduc(:,:) = 0.
     1908         zdvc(:,:) = 0.
     1909         zdqc(:,:,:) = 0.
     1910
     1911        if (co2clouds) then
     1912          CALL co2condens4micro(ngrid,nlayer,nq,ptimestep,
     1913     $                  capcal,zplay,zplev,tsurf,pt,
     1914     $                  pphi,pdt,pdu,pdv,zdtsurf,pu,pv,pq,pdq,
     1915     $                  co2ice,albedo,emis,
     1916     $                  zdtc,zdtsurfc,pdpsrf,zduc,zdvc,zdqc,
     1917     $                  fluxsurf_sw,zls,
     1918     $                  zdqssed_co2,zcondicea_co2microp,
     1919     &                  zdtcloudco2)
     1920        else
     1921          CALL co2condens(ngrid,nlayer,nq,ptimestep,
    18981922     $              capcal,zplay,zplev,tsurf,pt,
    18991923     $              pphi,pdt,pdu,pdv,zdtsurf,pu,pv,pq,pdq,
    … …  
    19031927     $              zdqssed_co2,zcondicea_co2microp,
    19041928     &              zdtcloudco2,zdqsc)
    1905 
    1906          DO l=1,nlayer
     1929           DO iq=1, nq
     1930           DO ig=1,ngrid
     1931              dqsurf(ig,iq)=dqsurf(ig,iq)+zdqsc(ig,iq)
     1932           ENDDO  ! (ig)
     1933           ENDDO    ! (iq)
     1934        end if
     1935        DO l=1,nlayer
    19071936           DO ig=1,ngrid
    19081937             pdt(ig,l)=pdt(ig,l)+zdtc(ig,l)
    … …  
    19101939             pdu(ig,l)=pdu(ig,l)+zduc(ig,l)
    19111940           ENDDO
    1912          ENDDO
    1913          DO ig=1,ngrid
    1914             zdtsurf(ig) = zdtsurf(ig) + zdtsurfc(ig)
    1915          ENDDO
    1916 
    1917          IF (tracer) THEN
     1941        ENDDO
     1942        DO ig=1,ngrid
     1943           zdtsurf(ig) = zdtsurf(ig) + zdtsurfc(ig)
     1944        ENDDO
     1945
     1946        IF (tracer) THEN
    19181947           DO iq=1, nq
    19191948            DO l=1,nlayer
    19201949              DO ig=1,ngrid
    1921                 pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqc(ig,l,iq) 
     1950                pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqc(ig,l,iq)
    19221951              ENDDO
    19231952            ENDDO
    … …  
    19371966          ps(ig) = zplev(ig,1) + pdpsrf(ig)*ptimestep
    19381967        ENDDO
    1939        
    19401968        ! update pressure levels
    19411969        DO l=1,nlayer
    … …  
    19461974        ENDDO
    19471975        zplev(:,nlayer+1) = 0.
    1948        
    19491976        ! update layers altitude
    19501977        DO l=2,nlayer
    … …  
    19561983        ENDDO
    19571984#endif
    1958      
    19591985      ENDIF  ! of IF (callcond)
    19601986
    … …  
    19711997        ENDDO    ! (iq)
    19721998      ENDIF
    1973      
    19741999c-----------------------------------------------------------------------
    19752000c   12. Surface  and sub-surface soil temperature
    … …  
    20742099
    20752100      DO ig=1,ngrid
    2076          watercap(ig)=watercap(ig)+ptimestep*dwatercap(ig) 
     2101         watercap(ig)=watercap(ig)+ptimestep*dwatercap(ig)
    20772102      ENDDO
    20782103
    … …  
    21032128      ENDDO
    21042129
    2105       ! Density 
     2130      ! Density
    21062131      DO l=1,nlayer
    21072132         DO ig=1,ngrid
    … …  
    21182143       ENDDO
    21192144
     2145      ! Total mass of co2
     2146      mtotco2 = 0.
     2147      icetotco2 = 0.
     2148      vaptotco2 = 0.
     2149      do ig=1,ngrid
     2150        do l=1,nlayer
     2151          vaptotco2 = vaptotco2 + zq(ig,l,igcm_co2) *
     2152     &                    (zplev(ig,l) - zplev(ig,l+1)) / g
     2153          icetotco2 = icetotco2 + zq(ig,l,igcm_co2_ice) *
     2154     &                    (zplev(ig,l) - zplev(ig,l+1)) / g
     2155        end do
     2156        mtotco2 = mtotco2 + co2ice(ig)
     2157      end do
     2158      mtotco2 = mtotco2 + vaptotco2 + icetotco2
    21202159
    21212160c    Compute surface stress : (NB: z0 is a common in surfdat.h)
    … …  
    22502289     .          icount,' date=',ztime_fin
    22512290           
    2252            
     2291
    22532292          call physdem1("restartfi.nc",nsoilmx,ngrid,nlayer,nq,
    22542293     .                ptimestep,ztime_fin,
    … …  
    22672306         if (tracer) then
    22682307           ! Density-scaled opacities
    2269               do ig=1,ngrid 
    2270                 dsodust(ig,:) = 
     2308              do ig=1,ngrid
     2309                dsodust(ig,:) =
    22712310     &                dsodust(ig,:)*tauscaling(ig)
    2272                 dsords(ig,:) = 
     2311                dsords(ig,:) =
    22732312     &                dsords(ig,:)*tauscaling(ig)
    22742313                dsotop(ig,:) =
    22752314     &                dsotop(ig,:)*tauscaling(ig)
    22762315              enddo
    2277          
     2316
    22782317           if(doubleq) then
    22792318              do ig=1,ngrid
    … …  
    23332372     &                zq(ig,:,igcm_ccnco2_mass)*tauscaling(ig)
    23342373              enddo
    2335 c D. BARDET compute integrated CO2 vapor and ice content
    2336               mtotco2(:)=0
    2337               icetotco2(:)=0
    2338               do ig=1,ngrid
    2339                do l=1,nlayer
    2340                  mtotco2(ig) = mtotco2(ig) +
    2341      &                      zq(ig,l,igcm_co2) *
    2342      &                      (zplev(ig,l) - zplev(ig,l+1)) / g
    2343                  icetotco2(ig) = icetotco2(ig) +
    2344      &                        zq(ig,l,igcm_co2_ice) *
    2345      &                        (zplev(ig,l) - zplev(ig,l+1)) / g
    2346                enddo
    2347               enddo
    23482374           endif ! of if (co2clouds)
    23492375                             
    … …  
    23592385             ENDIF !hdo
    23602386
    2361              do ig=1,ngrid 
     2387             do ig=1,ngrid
    23622388               do l=1,nlayer
    23632389                 mtot(ig) = mtot(ig) +
    … …  
    23682394     &                        (zplev(ig,l) - zplev(ig,l+1)) / g
    23692395                 IF (hdo) then
    2370                  mtotD(ig) = mtotD(ig) + 
    2371      &                      zq(ig,l,igcm_hdo_vap) * 
     2396                 mtotD(ig) = mtotD(ig) +
     2397     &                      zq(ig,l,igcm_hdo_vap) *
    23722398     &                      (zplev(ig,l) - zplev(ig,l+1)) / g
    2373                  icetotD(ig) = icetotD(ig) + 
    2374      &                        zq(ig,l,igcm_hdo_ice) * 
     2399                 icetotD(ig) = icetotD(ig) +
     2400     &                        zq(ig,l,igcm_hdo_ice) *
    23752401     &                        (zplev(ig,l) - zplev(ig,l+1)) / g
    23762402                 ENDIF !hdo
    … …  
    24442470
    24452471           endif ! of if (water)
    2446 
    2447 
    24482472        endif                   ! of if (tracer)
    2449 
    24502473#ifndef MESOSCALE
    24512474c        -----------------------------------------------------------------
    … …  
    28262849
    28272850#else
    2828      !!! this is to ensure correct initialisation of mesoscale model
    2829         call WRITEDIAGFI(ngrid,"tsurf","Surface temperature","K",2,
    2830      &                  tsurf)
    2831         call WRITEDIAGFI(ngrid,"ps","surface pressure","Pa",2,ps)
    2832         call WRITEDIAGFI(ngrid,"co2ice","co2 ice thickness","kg.m-2",2,
    2833      &                  co2ice)
    2834         call WRITEDIAGFI(ngrid,"temp","temperature","K",3,zt)
    2835         call WRITEDIAGFI(ngrid,"u","Zonal wind","m.s-1",3,zu)
    2836         call WRITEDIAGFI(ngrid,"v","Meridional wind","m.s-1",3,zv)
    2837         call WRITEDIAGFI(ngrid,"emis","Surface emissivity","w.m-1",2,
    2838      &                  emis)
    2839         call WRITEDIAGFI(ngrid,"tsoil","Soil temperature",
    2840      &                       "K",3,tsoil)
    2841         call WRITEDIAGFI(ngrid,"inertiedat","Soil inertia",
    2842      &                       "K",3,inertiedat)
     2851      !!! this is to ensure correct initialisation of mesoscale model
     2852      call WRITEDIAGFI(ngrid,"tsurf","Surface temperature","K",2,
     2853     &                 tsurf)
     2854      call WRITEDIAGFI(ngrid,"ps","surface pressure","Pa",2,ps)
     2855      call WRITEDIAGFI(ngrid,"co2ice","co2 ice thickness","kg.m-2",2,
     2856     &                 co2ice)
     2857      call WRITEDIAGFI(ngrid,"temp","temperature","K",3,zt)
     2858      call WRITEDIAGFI(ngrid,"u","Zonal wind","m.s-1",3,zu)
     2859      call WRITEDIAGFI(ngrid,"v","Meridional wind","m.s-1",3,zv)
     2860      call WRITEDIAGFI(ngrid,"emis","Surface emissivity","w.m-1",2,
     2861     &                 emis)
     2862      call WRITEDIAGFI(ngrid,"tsoil","Soil temperature",
     2863     &                 "K",3,tsoil)
     2864      call WRITEDIAGFI(ngrid,"inertiedat","Soil inertia",
     2865     &                 "K",3,inertiedat)
    28432866#endif
    2844 
    28452867
    28462868c        ----------------------------------------------------------
    28472869c        Outputs of the CO2 cycle
    28482870c        ----------------------------------------------------------
    2849            
    2850          if (tracer.and.(igcm_co2.ne.0)) then
    2851 !          call WRITEDIAGFI(ngrid,"co2_l1","co2 mix. ratio in 1st layer",
    2852 !    &                     "kg/kg",2,zq(1,1,igcm_co2))
    2853           call WRITEDIAGFI(ngrid,"co2","co2 mass mixing ratio",
    2854      &                     "kg/kg",3,zq(1,1,igcm_co2))
    2855           if (co2clouds) then
    2856             CALL WRITEDIAGFI(ngrid,'mtotco2',
    2857      &                       'total mass of CO2 vapor',
    2858      &                       'kg/m2',2,mtotco2)
    2859             CALL WRITEDIAGFI(ngrid,'zdtcloudco2',
    2860      &                       'temperature variation of CO2 latent heat',
    2861      &                       'K/s',3,zdtcloudco2)
    2862 
    2863             CALL WRITEDIAGFI(ngrid,'icetotco2',
    2864      &                       'total mass of CO2 ice',
    2865      &                       'kg/m2',2,icetotco2)
     2871
     2872      if (tracer.and.(igcm_co2.ne.0)) then
     2873!       call WRITEDIAGFI(ngrid,"co2_l1","co2 mix. ratio in 1st layer",
     2874!    &                   "kg/kg",2,zq(1,1,igcm_co2))
     2875        call WRITEDIAGFI(ngrid,"co2","co2 mass mixing ratio",
     2876     &                   "kg/kg",3,zq(:,:,igcm_co2))
     2877
     2878        if (co2clouds) then
     2879          call WRITEDIAGFI(ngrid,'zdtcloudco2',
     2880     &                     'temperature variation of CO2 latent heat',
     2881     &                     'K/s',3,zdtcloudco2)
    28662882 
    2867             call WRITEDIAGFI(ngrid,'ccnqco2','CCNco2 mass mr',
    2868      &                       'kg/kg',3,qccnco2)
    2869             call WRITEDIAGFI(ngrid,'ccnNco2','CCNco2 number',
    2870      &                       'part/kg',3,nccnco2)
    2871             call WRITEDIAGFI(ngrid,'co2_ice','co2_ice','kg/kg',
    2872      &                       3,zq(:,:,igcm_co2_ice))
    2873             call WRITEDIAGFI(ngrid,'precip_co2_ice',
    2874      &                       'surface deposition of co2 ice',
    2875      &                       'kg.m-2.s-1',2,
    2876      &                       zdqssed(1:ngrid,igcm_co2_ice))
    2877            endif ! of if (co2clouds)
    2878          endif ! of if (tracer.and.(igcm_co2.ne.0))
    2879         ! Output He tracer, if there is one
    2880         if (tracer.and.(igcm_he.ne.0)) then
    2881           call WRITEDIAGFI(ngrid,"he","helium mass mixing ratio",
    2882      &                     "kg/kg",3,zq(1,1,igcm_he))
    2883           vmr=zq(1:ngrid,1:nlayer,igcm_he)
    2884      &       *mmean(1:ngrid,1:nlayer)/mmol(igcm_he)
    2885           call WRITEDIAGFI(ngrid,'vmr_he','helium vol. mixing ratio',
    2886      &                       'mol/mol',3,vmr)
    2887         endif
     2883          call WRITEDIAGFI(ngrid,'ccnqco2','CCNco2 mass mr',
     2884     &                     'kg/kg',3,qccnco2)
     2885
     2886          call WRITEDIAGFI(ngrid,'ccnNco2','CCNco2 number',
     2887     &                     'part/kg',3,nccnco2)
     2888
     2889          call WRITEDIAGFI(ngrid,'co2_ice','co2_ice','kg/kg',
     2890     &                     3,zq(:,:,igcm_co2_ice))
     2891
     2892          call WRITEDIAGFI(ngrid,'precip_co2_ice',
     2893     &                     'surface deposition of co2 ice',
     2894     &                     'kg.m-2.s-1',2,
     2895     &                     zdqssed(1:ngrid,igcm_co2_ice))
     2896        end if ! of if (co2clouds)
     2897      end if ! of if (tracer.and.(igcm_co2.ne.0))
     2898
     2899      ! Output He tracer, if there is one
     2900      if (tracer.and.(igcm_he.ne.0)) then
     2901        call WRITEDIAGFI(ngrid,"he","helium mass mixing ratio",
     2902     &                   "kg/kg",3,zq(1,1,igcm_he))
     2903        vmr = zq(1:ngrid,1:nlayer,igcm_he)
     2904     &        * mmean(1:ngrid,1:nlayer)/mmol(igcm_he)
     2905        call WRITEDIAGFI(ngrid,'vmr_he','helium vol. mixing ratio',
     2906     &                   'mol/mol',3,vmr)
     2907      end if
    28882908
    28892909c        ----------------------------------------------------------
    28902910c        Outputs of the water cycle
    28912911c        ----------------------------------------------------------
    2892          if (tracer) then
    2893            if (water) then
    2894 
     2912      if (tracer) then
     2913        if (water) then
    28952914#ifdef MESOINI
    28962915            !!!! waterice = q01, voir readmeteo.F90
    2897             call WRITEDIAGFI(ngrid,'q01',noms(igcm_h2o_ice),
    2898      &                      'kg/kg',3,
    2899      &                       zq(1:ngrid,1:nlayer,igcm_h2o_ice))
     2916          call WRITEDIAGFI(ngrid,'q01',noms(igcm_h2o_ice),
     2917     &                     'kg/kg',3,
     2918     &                     zq(1:ngrid,1:nlayer,igcm_h2o_ice))
    29002919            !!!! watervapor = q02, voir readmeteo.F90
    2901             call WRITEDIAGFI(ngrid,'q02',noms(igcm_h2o_vap),
    2902      &                      'kg/kg',3,
    2903      &                       zq(1:ngrid,1:nlayer,igcm_h2o_vap))
     2920          call WRITEDIAGFI(ngrid,'q02',noms(igcm_h2o_vap),
     2921     &                     'kg/kg',3,
     2922     &                     zq(1:ngrid,1:nlayer,igcm_h2o_vap))
    29042923            !!!! surface waterice qsurf02 (voir readmeteo)
    2905             call WRITEDIAGFI(ngrid,'qsurf02','surface tracer',
    2906      &                      'kg.m-2',2,
    2907      &                       qsurf(1:ngrid,igcm_h2o_ice))
     2924          call WRITEDIAGFI(ngrid,'qsurf02','surface tracer',
     2925     &                     'kg.m-2',2,
     2926     &                     qsurf(1:ngrid,igcm_h2o_ice))
    29082927#endif
    2909             CALL WRITEDIAGFI(ngrid,'mtot',
    2910      &                       'total mass of water vapor',
    2911      &                       'kg/m2',2,mtot)
    2912             CALL WRITEDIAGFI(ngrid,'icetot',
    2913      &                       'total mass of water ice',
    2914      &                       'kg/m2',2,icetot)
    2915             vmr=zq(1:ngrid,1:nlayer,igcm_h2o_ice)
    2916      &      *mmean(1:ngrid,1:nlayer)/mmol(igcm_h2o_ice)
    2917             CALL WRITEDIAGFI(ngrid,'vmr_h2oice','h2o ice vmr',
    2918      &                       'mol/mol',3,vmr)
    2919             vmr=zq(1:ngrid,1:nlayer,igcm_h2o_vap)
    2920      &      *mmean(1:ngrid,1:nlayer)/mmol(igcm_h2o_vap)
    2921             CALL WRITEDIAGFI(ngrid,'vmr_h2ovap','h2o vap vmr',
    2922      &                       'mol/mol',3,vmr)
    2923             CALL WRITEDIAGFI(ngrid,'reffice',
    2924      &                       'Mean reff',
    2925      &                       'm',2,rave)
    2926 
    2927             call WRITEDIAGFI(ngrid,'h2o_ice','h2o_ice','kg/kg',
    2928      &             3,zq(:,:,igcm_h2o_ice))
    2929             call WRITEDIAGFI(ngrid,'h2o_vap','h2o_vap','kg/kg',
    2930      &             3,zq(:,:,igcm_h2o_vap))
     2928          call WRITEDIAGFI(ngrid,'mtot',
     2929     &                     'total mass of water vapor',
     2930     &                     'kg/m2',2,mtot)
     2931          call WRITEDIAGFI(ngrid,'icetot',
     2932     &                     'total mass of water ice',
     2933     &                     'kg/m2',2,icetot)
     2934          vmr = zq(1:ngrid,1:nlayer,igcm_h2o_ice)
     2935     &          * mmean(1:ngrid,1:nlayer)/mmol(igcm_h2o_ice)
     2936          call WRITEDIAGFI(ngrid,'vmr_h2oice','h2o ice vmr',
     2937     &                     'mol/mol',3,vmr)
     2938          vmr = zq(1:ngrid,1:nlayer,igcm_h2o_vap)
     2939     &          * mmean(1:ngrid,1:nlayer)/mmol(igcm_h2o_vap)
     2940          call WRITEDIAGFI(ngrid,'vmr_h2ovap','h2o vap vmr',
     2941     &                     'mol/mol',3,vmr)
     2942          call WRITEDIAGFI(ngrid,'reffice',
     2943     &                     'Mean reff',
     2944     &                     'm',2,rave)
     2945          call WRITEDIAGFI(ngrid,'h2o_ice','h2o_ice','kg/kg',
     2946     &                     3,zq(:,:,igcm_h2o_ice))
     2947          call WRITEDIAGFI(ngrid,'h2o_vap','h2o_vap','kg/kg',
     2948     &                     3,zq(:,:,igcm_h2o_vap))
    29312949
    29322950            if (hdo) then
    … …  
    29542972            do l=1,nlayer
    29552973                do ig=1,ngrid
    2956                 if (zq(ig,l,igcm_h2o_vap).gt.qperemin) then 
     2974                if (zq(ig,l,igcm_h2o_vap).gt.qperemin) then
    29572975                    DoH_vap(ig,l) = ( zq(ig,l,igcm_hdo_vap)/
    29582976     &              zq(ig,l,igcm_h2o_vap) )*1./(2.*155.76e-6)
    … …  
    29762994            CALL WRITEDIAGFI(ngrid,'DoH_vap',
    29772995     &                       'D/H ratio in vapor',
    2978      &                       ' ',3,DoH_vap) 
     2996     &                       ' ',3,DoH_vap)
    29792997            CALL WRITEDIAGFI(ngrid,'DoH_ice',
    29802998     &                       'D/H ratio in ice',
    … …  
    30333051              endif
    30343052!A. Pottier
    3035               if (CLFvarying) then !AP14 nebulosity
    3036                  call WRITEDIAGFI(ngrid,'totcloudfrac',
    3037      &                'Total cloud fraction',
     3053          if (CLFvarying) then !AP14 nebulosity
     3054            call WRITEDIAGFI(ngrid,'totcloudfrac',
     3055     &                       'Total cloud fraction',
    30383056     &                       ' ',2,totcloudfrac)
    3039               endif !clf varying
    3040 
    3041            endif !(water)
    3042 
    3043 
    3044            if (water.and..not.photochem) then
    3045              iq=nq
    3046 c            write(str2(1:2),'(i2.2)') iq
    3047 c            call WRITEDIAGFI(ngrid,'dqs'//str2,'dqscloud',
    3048 c    &                       'kg.m-2',2,zdqscloud(1,iq))
    3049 c            call WRITEDIAGFI(ngrid,'dqch'//str2,'var chim',
    3050 c    &                       'kg/kg',3,zdqchim(1,1,iq))
    3051 c            call WRITEDIAGFI(ngrid,'dqd'//str2,'var dif',
    3052 c    &                       'kg/kg',3,zdqdif(1,1,iq))
    3053 c            call WRITEDIAGFI(ngrid,'dqa'//str2,'var adj',
    3054 c    &                       'kg/kg',3,zdqadj(1,1,iq))
    3055 c            call WRITEDIAGFI(ngrid,'dqc'//str2,'var c',
    3056 c    &                       'kg/kg',3,zdqc(1,1,iq))
    3057            endif  !(water.and..not.photochem)
    3058          endif
     3057          end if !clf varying
     3058        end if !(water)
     3059
     3060        if (water.and..not.photochem) then
     3061          iq = nq
     3062c          write(str2(1:2),'(i2.2)') iq
     3063c          call WRITEDIAGFI(ngrid,'dqs'//str2,'dqscloud',
     3064c    &                      'kg.m-2',2,zdqscloud(1,iq))
     3065c          call WRITEDIAGFI(ngrid,'dqch'//str2,'var chim',
     3066c    &                      'kg/kg',3,zdqchim(1,1,iq))
     3067c          call WRITEDIAGFI(ngrid,'dqd'//str2,'var dif',
     3068c    &                      'kg/kg',3,zdqdif(1,1,iq))
     3069c          call WRITEDIAGFI(ngrid,'dqa'//str2,'var adj',
     3070c    &                      'kg/kg',3,zdqadj(1,1,iq))
     3071c          call WRITEDIAGFI(ngrid,'dqc'//str2,'var c',
     3072c    &                      'kg/kg',3,zdqc(1,1,iq))
     3073        end if  !(water.and..not.photochem)
     3074      end if !tracer
    30593075
    30603076c        ----------------------------------------------------------
    … …  
    30623078c        ----------------------------------------------------------
    30633079
    3064         call WRITEDIAGFI(ngrid,'tauref',
    3065      &                    'Dust ref opt depth','NU',2,tauref)
    3066 
    3067          if (tracer.and.(dustbin.ne.0)) then
    3068 
    3069           call WRITEDIAGFI(ngrid,'tau','taudust','SI',2,tau(1,1))
     3080      call WRITEDIAGFI(ngrid,'tauref',
     3081     &                 'Dust ref opt depth','NU',2,tauref)
     3082
     3083      if (tracer.and.(dustbin.ne.0)) then
     3084
     3085        call WRITEDIAGFI(ngrid,'tau','taudust','SI',2,tau(1,1))
    30703086
    30713087#ifndef MESOINI
    … …  
    30973113             call WRITEDIAGFI(ngrid,'dustN','Dust number',
    30983114     &                        'part/kg',3,ndust)
    3099              
     3115
    31003116             select case (trim(dustiropacity))
    31013117              case ("tes")
    … …  
    31703186             call WRITEDIAGFI(ngrid,'totaldustq','total dust mass',
    31713187     &                        'kg/kg',3,qdusttotal)
    3172              
     3188
    31733189             select case (trim(dustiropacity))
    31743190              case ("tes")
    … …  
    31953211             end select
    31963212           endif ! (slpwind)
    3197            
     3213
    31983214           if (scavenging) then
    31993215             call WRITEDIAGFI(ngrid,'ccnq','CCN mass mr',
    … …  
    32613277c        ----------------------------------------------------------
    32623278
    3263          if(callthermos) then 
     3279         if(callthermos) then
    32643280
    32653281            call WRITEDIAGFI(ngrid,"q15um","15 um cooling","K/s",
    … …  
    32833299         endif  !(callthermos)
    32843300
     3301            call WRITEDIAGFI(ngrid,"q15um","15 um cooling","K/s",
     3302     $           3,zdtnlte)
     3303            call WRITEDIAGFI(ngrid,"quv","UV heating","K/s",
     3304     $           3,zdteuv)
     3305            call WRITEDIAGFI(ngrid,"cond","Thermal conduction","K/s",
     3306     $           3,zdtconduc)
     3307            call WRITEDIAGFI(ngrid,"qnir","NIR heating","K/s",
     3308     $           3,zdtnirco2)
     3309
    32853310c        ----------------------------------------------------------
    32863311c        ----------------------------------------------------------
    … …  
    32923317c        Outputs of thermals
    32933318c        ----------------------------------------------------------
    3294          if (calltherm) then
    3295 
     3319      if (calltherm) then
    32963320!        call WRITEDIAGFI(ngrid,'dtke',
    3297 !     &              'tendance tke thermiques','m**2/s**2',
    3298 !     &                         3,dtke_th)
     3321!     &                   'tendance tke thermiques','m**2/s**2',
     3322!     &                   3,dtke_th)
    32993323!        call WRITEDIAGFI(ngrid,'d_u_ajs',
    3300 !     &              'tendance u thermiques','m/s',
    3301 !     &                         3,pdu_th*ptimestep)
     3324!     &                   'tendance u thermiques','m/s',
     3325!     &                   3,pdu_th*ptimestep)
    33023326!        call WRITEDIAGFI(ngrid,'d_v_ajs',
    3303 !     &              'tendance v thermiques','m/s',
    3304 !     &                         3,pdv_th*ptimestep)
     3327!     &                   'tendance v thermiques','m/s',
     3328!     &                   3,pdv_th*ptimestep)
    33053329!        if (tracer) then
    3306 !        if (nq .eq. 2) then
    3307 !        call WRITEDIAGFI(ngrid,'deltaq_th',
    3308 !     &              'delta q thermiques','kg/kg',
    3309 !     &                         3,ptimestep*pdq_th(:,:,2))
    3310 !        endif
    3311 !        endif
     3330!          if (nq .eq. 2) then
     3331!            call WRITEDIAGFI(ngrid,'deltaq_th',
     3332!     &                       'delta q thermiques','kg/kg',
     3333!     &                       3,ptimestep*pdq_th(:,:,2))
     3334!          end if
     3335!        end if
    33123336
    33133337        call WRITEDIAGFI(ngrid,'zmax_th',
    3314      &              'hauteur du thermique','m',
    3315      &                         2,zmax_th)
     3338     &                   'hauteur du thermique','m',
     3339     &                    2,zmax_th)
    33163340        call WRITEDIAGFI(ngrid,'hfmax_th',
    3317      &              'maximum TH heat flux','K.m/s',
    3318      &                         2,hfmax_th)
     3341     &                   'maximum TH heat flux','K.m/s',
     3342     &                   2,hfmax_th)
    33193343        call WRITEDIAGFI(ngrid,'wstar',
    3320      &              'maximum TH vertical velocity','m/s',
    3321      &                         2,wstar)
    3322 
    3323          endif
     3344     &                   'maximum TH vertical velocity','m/s',
     3345     &                   2,wstar)
     3346      end if
    33243347
    33253348c        ----------------------------------------------------------
    … …  
    33593382c
    33603383c        CALL writeg1d(ngrid,1,fluxsurf_lw,'Fs_ir','W.m-2')
    3361 c         CALL writeg1d(ngrid,1,tsurf,'tsurf','K')
    3362 c         CALL writeg1d(ngrid,1,ps,'ps','Pa')
    3363        
    3364 c         CALL writeg1d(ngrid,nlayer,zt,'T','K')
     3384c        CALL writeg1d(ngrid,1,tsurf,'tsurf','K')
     3385c        CALL writeg1d(ngrid,1,ps,'ps','Pa')
     3386c        CALL writeg1d(ngrid,nlayer,zt,'T','K')
    33653387c        CALL writeg1d(ngrid,nlayer,pu,'u','m.s-1')
    33663388c        CALL writeg1d(ngrid,nlayer,pv,'v','m.s-1')
    … …  
    33833405     &                         0,wstar)
    33843406
    3385          co2col(:)=0.
    3386          if (tracer) then
    3387          do l=1,nlayer
    3388            do ig=1,ngrid
    3389              co2col(ig)=co2col(ig)+
    3390      &                  zq(ig,l,1)*(zplev(ig,l)-zplev(ig,l+1))/g
    3391          enddo
    3392          enddo
    3393 
    3394          end if
    3395          call WRITEDIAGFI(ngrid,'co2col','integrated co2 mass'          &
    3396      &                                      ,'kg/m-2',0,co2col)
    3397          endif ! of if (calltherm)
    3398 
    3399          call WRITEDIAGFI(ngrid,'w','vertical velocity'                 &
     3407         end if ! of if (calltherm)
     3408
     3409         call WRITEDIAGFI(ngrid,'w','vertical velocity'
    34003410     &                              ,'m/s',1,pw)
    34013411         call WRITEDIAGFI(ngrid,"q2","q2","kg.m-3",1,q2)
    … …  
    34083418         call WRITEDIAGFI(ngrid,"pplev","Pressure","Pa",1,zplev)
    34093419         call WRITEDIAGFI(ngrid,"rho","rho","kg.m-3",1,rho)
    3410          call WRITEDIAGFI(ngrid,"dtrad","rad. heat. rate",              &
     3420         call WRITEDIAGFI(ngrid,"dtrad","rad. heat. rate",
    34113421     &              "K.s-1",1,dtrad)
    3412         call WRITEDIAGFI(ngrid,'sw_htrt','sw heat. rate',
     3422         call WRITEDIAGFI(ngrid,'sw_htrt','sw heat. rate',
    34133423     &                   'w.m-2',1,zdtsw)
    3414         call WRITEDIAGFI(ngrid,'lw_htrt','lw heat. rate',
     3424         call WRITEDIAGFI(ngrid,'lw_htrt','lw heat. rate',
    34153425     &                   'w.m-2',1,zdtlw)
    34163426         call WRITEDIAGFI(ngrid,"co2ice","co2 ice thickness"
    … …  
    34183428
    34193429         if (igcm_co2.ne.0) then
    3420           call co2sat(ngrid*nlayer,zt,zplay,zqsatco2)
     3430          call co2sat(ngrid*nlayer,zt,zqsatco2)
    34213431           do ig=1,ngrid
    34223432            do l=1,nlayer
    … …  
    34373447c         call WRITEDIAGFI(ngrid,"satuco2","vap in satu","kg/kg",1,
    34383448c     &        satuco2)
    3439 c         call WRITEdiagfi(ngrid,"riceco2","ice radius","m"
     3449c         call WRITEDIAGFI(ngrid,"riceco2","ice radius","m"
    34403450c     &        ,1,riceco2)
    34413451! or output in diagfi.nc (for testphys1d)
    … …  
    35563566     &                 * (zplev(1,l) - zplev(1,l+1)) / g
    35573567               if (hdo) THEN
    3558                  mtotD = mtotD +  zq(1,l,igcm_hdo_vap) 
     3568                 mtotD = mtotD +  zq(1,l,igcm_hdo_vap)
    35593569     &                 * (zplev(1,l) - zplev(1,l+1)) / g
    3560                  icetotD = icetotD +  zq(1,l,igcm_hdo_ice) 
     3570                 icetotD = icetotD +  zq(1,l,igcm_hdo_ice)
    35613571     &                 * (zplev(1,l) - zplev(1,l+1)) / g
    35623572               ENDIF !hdo
    … …  
    35663576                   hdotot = hdotot+mtotD+icetotD
    35673577               ENDIF ! hdo
    3568                
     3578
    35693579
    35703580             CALL WRITEDIAGFI(ngrid,'h2otot',
    … …  
    36103620            CALL WRITEDIAGFI(ngrid,'DoH_vap',
    36113621     &                       'D/H ratio in vapor',
    3612      &                       ' ',1,DoH_vap) 
     3622     &                       ' ',1,DoH_vap)
    36133623            CALL WRITEDIAGFI(ngrid,'DoH_ice',
    36143624     &                       'D/H ratio in ice',
    … …  
    37053715
    37063716        ENDIF ! hdo
    3707      
     3717
    37083718        call WRITEDIAGFI(ngrid,"rice","ice radius","m",1,
    37093719     &                    rice)
    … …  
    37153725              endif !clfvarying
    37163726
    3717          ENDIF ! of IF (water)
     3727        ENDIF ! of IF (water)
     3728
     3729
     3730      ! co2clouds
     3731      if (co2clouds) then
     3732        call WRITEDIAGFI(ngrid,'mtotco2',
     3733     &                   'total mass of CO2, including ice surf',
     3734     &                   'kg/m2',0,mtotco2)
     3735
     3736        call WRITEDIAGFI(ngrid,'vaptotco2',
     3737     &                   'total mass of CO2 vapor',
     3738     &                   'kg/m2',0,vaptotco2)
     3739
     3740        call WRITEDIAGFI(ngrid,'zdtcloudco2',
     3741     &                   'temperature variation of CO2 latent heat',
     3742     &                   'K/s',3,zdtcloudco2)
     3743
     3744        call WRITEDIAGFI(ngrid,'icetotco2',
     3745     &                   'total mass of CO2 ice',
     3746     &                   'kg/m2',0,icetotco2)
     3747      end if
    37183748         
    37193749ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    37203750ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    3721 
    37223751
    37233752         zlocal(1)=-log(zplay(1,1)/zplev(1,1))* Rnew(1,1)*zt(1,1)/g
    … …  
    37363765
    37373766      END IF       ! if(ngrid.ne.1)
     3767
     3768      co2totB = 0. ! added by C.M.
     3769      do ig=1,ngrid
     3770        do l=1,nlayer
     3771          co2totB = co2totB + (zplev(ig,l)-zplev(ig,l+1))/g*
     3772     &           (pq(ig,l,igcm_co2)+pq(ig,l,igcm_co2_ice)
     3773     &      +(pdq(ig,l,igcm_co2)+pdq(ig,l,igcm_co2_ice))*ptimestep)
     3774        enddo
     3775        co2totB = co2totB + co2ice(ig)
     3776      enddo
     3777
     3778      call WRITEDIAGFI(ngrid,'co2conservation',
     3779     &                   'Total CO2 mass conservation in physic',
     3780     &                   '%',0,(co2totA-co2totB)/co2totA)
     3781      call test_vmr(pq(:,:,igcm_co2), pdq(:,:,igcm_co2), ptimestep,
     3782     &               ngrid, nlayer, 'end physiq_mod')
    37383783
    37393784! XIOS outputs

  • trunk/LMDZ.MARS/libf/phymars/tcondco2.F90

    r2349 r2350  
    88! Condensation temperature for co2 ice; based on
    99! the saturation in co2sat.F JA17
    10 !--------------------------------------------------i
     10!--------------------------------------------------
    1111
    1212integer, intent(in) :: ngrid,nlay
    … …  
    1414double precision, intent(out), dimension(ngrid,nlay):: tcond ! CO2 condensation temperature     (atm)
    1515double precision:: A,B,pco2
     16real :: qco2
    1617integer:: ig,l
    1718
    1819A=dlog(1.382d12)
    1920B=-3182.48
    20 
     21qco2=0.
    2122DO l=1,nlay
    2223   DO ig=1,ngrid
     24!      qco2 = min(1., q(ig,l))!added by CM not sure, to be verified before
     25                             !  cleaning
     26!      pco2 = qco2 * (mmean(ig,l)/44.01) * p(ig,l)
    2327      pco2 = q(ig,l) * (mmean(ig,l)/44.01) * p(ig,l)
    2428      tcond(ig,l)=B/(dlog(pco2)-A)
Note: See TracChangeset for help on using the changeset viewer.