Ignore:
Timestamp:
Jun 9, 2020, 1:56:52 PM (5 years ago)
Author:
aslmd
Message:

Cancelled commits 2341 and 2342

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

Legend:

Unmodified
Added
Removed

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

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

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

    r2341 r2343  
    1 c=======================================================================
    2 c     Program for simulate the impact of the CO2 snow fall on
    3 c     the surface infrared emission (emissivity)
    4 c-----------------------------------------------------------------------
    5 c     Algorithme:
    6 c     1 - Initialization
    7 c     2 - Compute the surface emissivity
    8 c-----------------------------------------------------------------------
    9 c     Author: F. Forget
    10 c-----------------------------------------------------------------------
    11 c     Reference paper:
    12 c     Francois Forget and James B. Pollack
    13 c
    14 c     "Thermal infrared observations of the condensing Martian polar
    15 c     caps: CO2 ice temperatures and radiative budget"
    16 c
    17 c     JOURNAL OF GEOPHYSICAL RESEARCH, VOL. 101, NO. E7, PAGES 16,
    18 c     865-16,879, JULY 25, 1996
    19 c=======================================================================
     1      SUBROUTINE co2snow (ngrid,nlayer,ptimestep,emisref,condsub
     2     &         ,pplev,pcondicea,pcondices,pfallice,pemisurf)
    203
    21       SUBROUTINE co2snow(ngrid, nlayer, ptimestep, emisref, condsub,
    22      &                   pplev, pcondicea, pcondices, pfallice,
    23      &                   pemisurf)
    24      
    254      use surfdat_h, only: iceradius, dtemisice
    265      use geometry_mod, only: latitude ! grid point latitudes (rad)
    276      use time_phylmdz_mod, only: daysec
     7      IMPLICIT NONE
    288
    29       IMPLICIT NONE
     9c=======================================================================
     10c     Program for simulate the impact of the CO2 snow fall on
     11c     the surface infrared emission (emissivity)  and on
     12c     the airborne dust
     13c     F.Forget 1996
     14c=======================================================================
     15
     16c Declarations
     17c ------------
    3018
    3119#include "callkeys.h"
    3220
     21c     Arguments
     22c     ---------
    3323
    34 c=======================================================================
    35 c     Variables
    36 c=======================================================================
    37 c     Inputs:
    38 c     -------
    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
     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)
    5635
    57 c     Output:
    58 c     -------
    59       real, intent(out) ::
    60      &   pemisurf(ngrid) ! surface emissivity
     36      REAL,INTENT(OUT) :: pemisurf(ngrid) ! surface emissivity
    6137
    62 c     local:
    63 c     ------
    64       integer ::
    65      &   l,
    66      &   ig,
    67      &   icap
     38c     local
     39c     -----
     40      integer ig , l , icap
    6841
    69       real ::
    70      &   sumdaer,
    71      &   zdemisurf
     42      REAL zdemisurf ,dtemis
     43      REAL sumdaer
    7244
    73       real, parameter ::
    74      &   alpha = 0.45
     45c     saved
     46c     -----
     47      REAL,save :: Kscat(2), scaveng
     48      LOGICAL,SAVE :: firstcall=.true.
    7549
    76 c     saved:
    77 c     ------
    78       real, save ::
    79      &   Kscat(2) ! Kscat: coefficient for decreasing the surface
    80 c                 !        emissivity
    81 c                 !
    82 c                 ! Kscat = (0.001/3.) * alpha / iceradius
    83 c                 !          with 0.3 < alpha < 0.6
    84 c                 !
    85 c                 ! alpha set to 0.45 (coeff from emis = f (tau)) and
    86 c                 ! iceradius the mean radius of the scaterring
    87 c                 ! particles (200.e-6 < iceradius < 10.e-6)
    88 c                 !
    89 c                 ! (2) = N and S hemispheres
     50c --------------
     51c Initialisation
     52c --------------
    9053
    91       logical, save ::
    92      &   firstcall = .true.
    93 c=======================================================================
    94 c BEGIN
    95 c=======================================================================
    96 c 1 - Initialization
    97 c=======================================================================
     54      ! AS: firstcall OK absolute
    9855      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.
     56
     57c   Kscat : coefficient for decreasing the surface emissivity
     58c           =(0.001/3.)*alpha/iceradius ,
     59c           with 0.3< alpha < 0.6, set to 0.45 (coeff from emis = f (tau))
     60c           and iceradius the mean radius of the
     61c           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
     66c          Scavenging Ratio (dust concentration in the air / in the snow)
     67           scaveng = 100.0
     68           
     69c          Collision Scavenging coefficient  (m2.kg-1)
     70c          Csca  = 2.3  ! not used yet !!!!!!!!!!!
     71           firstcall = .false.
     72
    10373      end if
    104 c=======================================================================
    105 c 2 - Compute the surface emissivity
    106 c=======================================================================
    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
    11474
    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
     75
     76c     LOOP on grid points
     77c     -------------------
     78      do ig=1,ngrid
     79         if (condsub(ig)) then
     80
     81c         IF (scavenging) then
     82c          Airborne Dust
     83c          -------------
     84c          sumdaer=0.
     85c          do l=nlayer, 1, -1
     86c             pdaerosol(ig,l)= -paerosol(ig,l,1)*
     87c    &              (1-exp(-scaveng*pcondicea(ig,l)*ptimestep*g/
     88c    &               (pplev(ig,l)-pplev(ig,l+1))))/ptimestep 
     89
     90c    &           - Csca*paerosol(ig,l,1) ! Scavenging by collision
     91c    &           * 0.5*(pfallice(ig,l)) ! not included
     92
     93c             test to avoid releasing to much dust when subliming:
     94c             if(pdaerosol(ig,l).gt.-sumdaer)pdaerosol(ig,l)=-sumdaer
     95c             sumdaer=sumdaer + pdaerosol(ig,l)
     96
     97c            if (-pdaerosol(ig,l)*ptimestep.gt.paerosol(ig,l,1)) then
     98c                write(*,*) 'ds co2snow: aerosol < 0.0 !!!'
     99c                write(*,*) 'ig =' , ig
     100c            end if
     101c          end do
     102c         END IF
     103
     104c          Surface emissivity
     105c          ------------------
     106c   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)
     116c Using directly the diferential equation:
     117c    &       -Kscat(icap)*emisref(ig)*
     118c    &      (pemisurf(ig)/emisref(ig))**4 *pfallice(ig,1)
     119c 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
    134133           pemisurf(ig) = emisref(ig)
    135         end if
    136       end do ! ngrid
    137 c=======================================================================
    138 c END
    139 c=======================================================================
     134         end if
     135      end do
     136
    140137      return
    141138      end

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

    r2341 r2343  
    1818      use watersat_mod, only: watersat
    1919      use co2condens_mod, only: co2condens
    20       use co2condens_mod4micro, only: co2condens4micro
    2120      use co2cloud_mod, only: co2cloud, mem_Mccn_co2, mem_Mh2o_co2,
    2221     &                        mem_Nccn_co2
    … …  
    252251c     Variables used by the CO2 clouds microphysical scheme:
    253252      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)
    254255      real zdqssed_co2(ngrid)  ! CO2 flux at the surface  (kg.m-2.s-1)
    255256      real zcondicea_co2microp(ngrid,nlayer)
    … …  
    302303
    303304c     Tendancies due to various processes:
    304       REAL dqsurf(ngrid,nq)  ! tendency for tracers on surface (Kg/m2/s)
     305      REAL dqsurf(ngrid,nq)
    305306      REAL zdtlw(ngrid,nlayer)     ! (K/s)
    306307      REAL zdtsw(ngrid,nlayer)     ! (K/s)
    … …  
    383384      REAL mtot(ngrid)          ! Total mass of water vapor (kg/m2)
    384385      REAL mstormdtot(ngrid)    ! Total mass of stormdust tracer (kg/m2)
    385       REAL mdusttot(ngrid)      ! Total mass of dust tracer (kg/m2)
     386      REAL mdusttot(ngrid)    ! Total mass of dust tracer (kg/m2)
    386387      REAL icetot(ngrid)        ! Total mass of water 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)
     388      REAL mtotco2(ngrid)       ! Total mass of co2 vapor (kg/m2)
     389      REAL icetotco2(ngrid)     ! 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)
    392391      REAL Mccntot(ngrid)       ! Total mass of ccn (kg/m2)
    393392      REAL rave(ngrid)          ! Mean water ice effective radius (m)
    … …  
    404403      REAL DoH_ice(ngrid,nlayer) !D/H ratio
    405404      REAL DoH_surf(ngrid) !D/H ratio surface
    406 
     405                               
    407406      REAL dqdustsurf(ngrid) ! surface q dust flux (kg/m2/s)
    408407      REAL dndustsurf(ngrid) ! surface n dust flux (number/m2/s)
    … …  
    484483      real tf_clf, ntf_clf ! tf: fraction of clouds, ntf: fraction without clouds
    485484      real rave2(ngrid), totrave2(ngrid) ! Mean water ice mean radius (m)
    486 C test de conservation de la masse de CO2
    487       REAL co2totA
    488       REAL co2totB
    489485
    490486c entrainment by slope winds above sb-grid scale topography
    … …  
    500496
    501497c=======================================================================
    502 c---------- begin added by CM
    503       pdq(:,:,:) = 0.
    504 c---------end added by CM
    505 
    506498
    507499c 1. Initialisation:
    508500c -----------------
     501
    509502c  1.1   Initialisation only at first call
    510503c  ---------------------------------------
    … …  
    563556           ! starting without startfi.nc and with callsoil
    564557           ! is not yet possible as soildepth default is not defined
    565            if (callsoil) then
     558           if (callsoil) then                     
    566559              ! default mlayer distribution, following a power law:
    567560              !  mlayer(k)=lay1*alpha**(k-1/2)
    … …  
    732725      dsotop(:,:)=0.
    733726      dwatercap(:)=0
    734 
     727     
    735728#ifdef DUSTSTORM
    736729      pq_tmp(:,:,:)=0
    … …  
    817810      enddo
    818811
    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
    829812c-----------------------------------------------------------------------
    830813c    2. Compute radiative tendencies :
    … …  
    869852     &                      pq(1:ngrid,1:nlayer,igcm_o)*
    870853     &                      mmean(1:ngrid,1:nlayer)/mmol(igcm_o)
    871 
     854                 
    872855                 CALL nlte_tcool(ngrid,nlayer,zplay*9.869e-6,
    873856     $                pt,zzlay,co2vmr_gcm, n2vmr_gcm, covmr_gcm,
    … …  
    919902c          ------------------
    920903           ! callradite for the part with clouds
    921            clearsky=.false. ! part with clouds for both cases CLFvarying true/false
     904           clearsky=.false. ! part with clouds for both cases CLFvarying true/false 
    922905           CALL callradite(icount,ngrid,nlayer,nq,zday,zls,pq,albedo,
    923906     &     emis,mu0,zplev,zplay,pt,tsurf,fract,dist_sol,igout,
    … …  
    934917               ! if
    935918               ! CLFvarying ...) ?? AP ??
    936                clearsky=.true. !
     919               clearsky=.true. ! 
    937920               CALL callradite(icount,ngrid,nlayer,nq,zday,zls,pq,
    938921     &              albedo,emis,mu0,zplev,zplay,pt,tsurf,fract,
    … …  
    11361119     &                      pdqrds,wspeed,dsodust,dsords,dsotop,
    11371120     &                      tauref)
    1138 
     1121                   
    11391122c      update the tendencies of both dust after vertical transport
    11401123         DO l=1,nlayer
    … …  
    11881171           hsummit(:)=14000.
    11891172         endif
    1190          clearatm=.true. ! stormdust is not accounted in the extra heating on top of the mountains
     1173         clearatm=.true. ! stormdust is not accounted in the extra heating on top of the mountains 
    11911174         nohmons=.false.
    1192          pdqtop(:,:,:)=0.
     1175         pdqtop(:,:,:)=0. 
    11931176         CALL topmons(ngrid,nlayer,nq,ptime,ptimestep,
    11941177     &                pq,pdq,pt,pdt,zplev,zplay,zzlev,
    … …  
    12401223        CALL calldrag_noro(ngrid,nlayer,ptimestep,
    12411224     &                 zplay,zplev,pt,pu,pv,zdtgw,zdugw,zdvgw)
    1242 
     1225 
    12431226        DO l=1,nlayer
    12441227          DO ig=1,ngrid
    … …  
    13261309            DO l=1,nlayer
    13271310              DO ig=1,ngrid
    1328                  pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqdif(ig,l,iq)
     1311                 pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqdif(ig,l,iq) 
    13291312              ENDDO
    13301313            ENDDO
    … …  
    13811364
    13821365      if(calltherm .and. .not.turb_resolved) then
    1383 
     1366 
    13841367        call calltherm_interface(ngrid,nlayer,nq,
    13851368     $ tracer,igcm_co2,
    … …  
    13891372     $ pdu_th,pdv_th,pdt_th,pdq_th,lmax_th,zmax_th,
    13901373     $ dtke_th,zdhdif,hfmax_th,wstar,sensibFlux)
    1391 
     1374     
    13921375         DO l=1,nlayer
    13931376           DO ig=1,ngrid
    … …  
    13981381           ENDDO
    13991382        ENDDO
    1400 
     1383 
    14011384        DO ig=1,ngrid
    14021385          q2(ig,nlayer+1)=q2(ig,nlayer+1)+dtke_th(ig,nlayer+1)*ptimestep
    1403         ENDDO
     1386        ENDDO     
    14041387
    14051388        if (tracer) then
    … …  
    15101493     &                nq,tau,tauscaling,rdust,rice,nuice,
    15111494     &                rsedcloud,rhocloud,totcloudfrac)
    1512 
     1495     
    15131496c Temperature variation due to latent heat release
    15141497           if (activice) then
    1515                pdt(1:ngrid,1:nlayer) =
    1516      &          pdt(1:ngrid,1:nlayer) +
     1498               pdt(1:ngrid,1:nlayer) = 
     1499     &          pdt(1:ngrid,1:nlayer) + 
    15171500     &          zdtcloud(1:ngrid,1:nlayer)
    15181501           endif
    … …  
    15281511           if (hdo) then
    15291512! increment HDO vapour and ice atmospheric tracers tendencies
    1530            pdq(1:ngrid,1:nlayer,igcm_hdo_vap) =
    1531      &       pdq(1:ngrid,1:nlayer,igcm_hdo_vap) +
     1513           pdq(1:ngrid,1:nlayer,igcm_hdo_vap) = 
     1514     &       pdq(1:ngrid,1:nlayer,igcm_hdo_vap) + 
    15321515     &       zdqcloud(1:ngrid,1:nlayer,igcm_hdo_vap)
    1533            pdq(1:ngrid,1:nlayer,igcm_hdo_ice) =
    1534      &       pdq(1:ngrid,1:nlayer,igcm_hdo_ice) +
     1516           pdq(1:ngrid,1:nlayer,igcm_hdo_ice) = 
     1517     &       pdq(1:ngrid,1:nlayer,igcm_hdo_ice) + 
    15351518     &       zdqcloud(1:ngrid,1:nlayer,igcm_hdo_ice)
    15361519           endif !hdo
    … …  
    16051588c                                nucleation
    16061589c               imicroco2=50     micro-timestep is 1/50 of physical timestep
    1607         zdqssed_co2(:) = 0.
    1608 
    1609          IF (co2clouds) THEN
    1610            call co2cloud(ngrid,nlayer,ptimestep,
     1590         
     1591         IF (co2clouds ) THEN
     1592
     1593
     1594            call co2cloud(ngrid,nlayer,ptimestep,
    16111595     &           zplev,zplay,pdpsrf,zzlay,pt,pdt,
    16121596     &           pq,pdq,zdqcloudco2,zdtcloudco2,
    16131597     &           nq,tau,tauscaling,rdust,rice,riceco2,nuice,
     1598     &           rsedcloudco2,rhocloudco2,
    16141599     &           rsedcloud,rhocloud,zzlev,zdqssed_co2,
    1615      &           pdu,pu,zcondicea_co2microp, co2ice)
     1600     &           pdu,pu,zcondicea_co2microp)
    16161601
    16171602
    … …  
    16271612! We need to check that we have Nccn & Ndust > 0
    16281613! This is due to single precision rounding problems
     1614             
    16291615! increment dust tracers tendancies
    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 
     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)
    16481634!Update water ice clouds values as well
    16491635             if (co2useh2o) then
    … …  
    16571643     &               pdq(1:ngrid,1:nlayer,igcm_ccn_number) +
    16581644     &               zdqcloudco2(1:ngrid,1:nlayer,igcm_ccn_number)
    1659 
    1660 c     Negative values?
    1661                 where (pq(:,:,igcm_ccn_mass) +
     1645                where (pq(:,:,igcm_ccn_mass) +
    16621646     &                 ptimestep*pdq(:,:,igcm_ccn_mass) < 0.)
    16631647                  pdq(:,:,igcm_ccn_mass) =
    … …  
    16661650     &               - pq(:,:,igcm_ccn_number)/ptimestep + 1.e-30
    16671651                end where
    1668 c     Negative values?
    1669                 where (pq(:,:,igcm_ccn_number) +
     1652                where (pq(:,:,igcm_ccn_number) +
    16701653     &                 ptimestep*pdq(:,:,igcm_ccn_number) < 0.)
    16711654                  pdq(:,:,igcm_ccn_mass) =
    … …  
    17061689     &           - pq(:,:,igcm_dust_number)/ptimestep + 1.e-30
    17071690             end where
    1708       END IF ! of IF (co2clouds)
     1691     
     1692      END IF                    ! of IF (co2clouds)
    17091693
    17101694c   9b. Aerosol particles
    … …  
    17501734     &                tau,tauscaling)
    17511735c 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
    17521740           IF (rdstorm) THEN
    17531741c             Storm dust cannot sediment to the surface
    … …  
    18791867          DO ig=1,ngrid
    18801868            dtrad(ig,l)=dtrad(ig,l)+zdteuv(ig,l)
    1881             pdt(ig,l)=pdt(ig,l)+zdtconduc(ig,l)+zdteuv(ig,l)
     1869            pdt(ig,l)=pdt(ig,l)+zdtconduc(ig,l)
     1870     &                         +zdteuv(ig,l)
    18821871            pdv(ig,l)=pdv(ig,l)+zdvmolvis(ig,l)
    18831872            pdu(ig,l)=pdu(ig,l)+zdumolvis(ig,l)
    … …  
    18941883c     (should be the last atmospherical physical process to be computed)
    18951884c   -------------------------------------------
     1885
    18961886      IF (tituscap) THEN
    18971887         !!! get the actual co2 seasonal cap from Titus observations
    1898          CALL geticecover(ngrid, 180.*zls/pi,
     1888         CALL geticecover( ngrid, 180.*zls/pi,
    18991889     .                  180.*longitude/pi, 180.*latitude/pi, co2ice )
    19001890         co2ice = co2ice * 10000.
    19011891      ENDIF
    19021892     
     1893     
     1894      pdpsrf(:) = 0
    19031895
    19041896      IF (callcond) THEN
    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,
     1897         CALL co2condens(ngrid,nlayer,nq,ptimestep,
    19221898     $              capcal,zplay,zplev,tsurf,pt,
    19231899     $              pphi,pdt,pdu,pdv,zdtsurf,pu,pv,pq,pdq,
    … …  
    19271903     $              zdqssed_co2,zcondicea_co2microp,
    19281904     &              zdtcloudco2,zdqsc)
    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
     1905
     1906         DO l=1,nlayer
    19361907           DO ig=1,ngrid
    19371908             pdt(ig,l)=pdt(ig,l)+zdtc(ig,l)
    … …  
    19391910             pdu(ig,l)=pdu(ig,l)+zduc(ig,l)
    19401911           ENDDO
    1941         ENDDO
    1942         DO ig=1,ngrid
    1943            zdtsurf(ig) = zdtsurf(ig) + zdtsurfc(ig)
    1944         ENDDO
    1945 
    1946         IF (tracer) THEN
     1912         ENDDO
     1913         DO ig=1,ngrid
     1914            zdtsurf(ig) = zdtsurf(ig) + zdtsurfc(ig)
     1915         ENDDO
     1916
     1917         IF (tracer) THEN
    19471918           DO iq=1, nq
    19481919            DO l=1,nlayer
    19491920              DO ig=1,ngrid
    1950                 pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqc(ig,l,iq)
     1921                pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqc(ig,l,iq) 
    19511922              ENDDO
    19521923            ENDDO
    … …  
    19661937          ps(ig) = zplev(ig,1) + pdpsrf(ig)*ptimestep
    19671938        ENDDO
     1939       
    19681940        ! update pressure levels
    19691941        DO l=1,nlayer
    … …  
    19741946        ENDDO
    19751947        zplev(:,nlayer+1) = 0.
     1948       
    19761949        ! update layers altitude
    19771950        DO l=2,nlayer
    … …  
    19831956        ENDDO
    19841957#endif
     1958     
    19851959      ENDIF  ! of IF (callcond)
    19861960
    … …  
    19971971        ENDDO    ! (iq)
    19981972      ENDIF
     1973     
    19991974c-----------------------------------------------------------------------
    20001975c   12. Surface  and sub-surface soil temperature
    … …  
    20992074
    21002075      DO ig=1,ngrid
    2101          watercap(ig)=watercap(ig)+ptimestep*dwatercap(ig)
     2076         watercap(ig)=watercap(ig)+ptimestep*dwatercap(ig) 
    21022077      ENDDO
    21032078
    … …  
    21282103      ENDDO
    21292104
    2130       ! Density
     2105      ! Density 
    21312106      DO l=1,nlayer
    21322107         DO ig=1,ngrid
    … …  
    21432118       ENDDO
    21442119
    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
    21592120
    21602121c    Compute surface stress : (NB: z0 is a common in surfdat.h)
    … …  
    22892250     .          icount,' date=',ztime_fin
    22902251           
    2291 
     2252           
    22922253          call physdem1("restartfi.nc",nsoilmx,ngrid,nlayer,nq,
    22932254     .                ptimestep,ztime_fin,
    … …  
    23062267         if (tracer) then
    23072268           ! Density-scaled opacities
    2308               do ig=1,ngrid
    2309                 dsodust(ig,:) =
     2269              do ig=1,ngrid 
     2270                dsodust(ig,:) = 
    23102271     &                dsodust(ig,:)*tauscaling(ig)
    2311                 dsords(ig,:) =
     2272                dsords(ig,:) = 
    23122273     &                dsords(ig,:)*tauscaling(ig)
    23132274                dsotop(ig,:) =
    23142275     &                dsotop(ig,:)*tauscaling(ig)
    23152276              enddo
    2316 
     2277         
    23172278           if(doubleq) then
    23182279              do ig=1,ngrid
    … …  
    23722333     &                zq(ig,:,igcm_ccnco2_mass)*tauscaling(ig)
    23732334              enddo
     2335c 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
    23742348           endif ! of if (co2clouds)
    23752349                             
    … …  
    23852359             ENDIF !hdo
    23862360
    2387              do ig=1,ngrid
     2361             do ig=1,ngrid 
    23882362               do l=1,nlayer
    23892363                 mtot(ig) = mtot(ig) +
    … …  
    23942368     &                        (zplev(ig,l) - zplev(ig,l+1)) / g
    23952369                 IF (hdo) then
    2396                  mtotD(ig) = mtotD(ig) +
    2397      &                      zq(ig,l,igcm_hdo_vap) *
     2370                 mtotD(ig) = mtotD(ig) + 
     2371     &                      zq(ig,l,igcm_hdo_vap) * 
    23982372     &                      (zplev(ig,l) - zplev(ig,l+1)) / g
    2399                  icetotD(ig) = icetotD(ig) +
    2400      &                        zq(ig,l,igcm_hdo_ice) *
     2373                 icetotD(ig) = icetotD(ig) + 
     2374     &                        zq(ig,l,igcm_hdo_ice) * 
    24012375     &                        (zplev(ig,l) - zplev(ig,l+1)) / g
    24022376                 ENDIF !hdo
    … …  
    24702444
    24712445           endif ! of if (water)
     2446
     2447
    24722448        endif                   ! of if (tracer)
     2449
    24732450#ifndef MESOSCALE
    24742451c        -----------------------------------------------------------------
    … …  
    28492826
    28502827#else
    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)
     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)
    28662843#endif
     2844
    28672845
    28682846c        ----------------------------------------------------------
    28692847c        Outputs of the CO2 cycle
    28702848c        ----------------------------------------------------------
    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)
     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)
    28822866 
    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
     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
    29082888
    29092889c        ----------------------------------------------------------
    29102890c        Outputs of the water cycle
    29112891c        ----------------------------------------------------------
    2912       if (tracer) then
    2913         if (water) then
     2892         if (tracer) then
     2893           if (water) then
     2894
    29142895#ifdef MESOINI
    29152896            !!!! waterice = q01, voir readmeteo.F90
    2916           call WRITEDIAGFI(ngrid,'q01',noms(igcm_h2o_ice),
    2917      &                     'kg/kg',3,
    2918      &                     zq(1:ngrid,1:nlayer,igcm_h2o_ice))
     2897            call WRITEDIAGFI(ngrid,'q01',noms(igcm_h2o_ice),
     2898     &                      'kg/kg',3,
     2899     &                       zq(1:ngrid,1:nlayer,igcm_h2o_ice))
    29192900            !!!! watervapor = q02, voir readmeteo.F90
    2920           call WRITEDIAGFI(ngrid,'q02',noms(igcm_h2o_vap),
    2921      &                     'kg/kg',3,
    2922      &                     zq(1:ngrid,1:nlayer,igcm_h2o_vap))
     2901            call WRITEDIAGFI(ngrid,'q02',noms(igcm_h2o_vap),
     2902     &                      'kg/kg',3,
     2903     &                       zq(1:ngrid,1:nlayer,igcm_h2o_vap))
    29232904            !!!! surface waterice qsurf02 (voir readmeteo)
    2924           call WRITEDIAGFI(ngrid,'qsurf02','surface tracer',
    2925      &                     'kg.m-2',2,
    2926      &                     qsurf(1:ngrid,igcm_h2o_ice))
     2905            call WRITEDIAGFI(ngrid,'qsurf02','surface tracer',
     2906     &                      'kg.m-2',2,
     2907     &                       qsurf(1:ngrid,igcm_h2o_ice))
    29272908#endif
    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))
     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))
    29492931
    29502932            if (hdo) then
    … …  
    29722954            do l=1,nlayer
    29732955                do ig=1,ngrid
    2974                 if (zq(ig,l,igcm_h2o_vap).gt.qperemin) then
     2956                if (zq(ig,l,igcm_h2o_vap).gt.qperemin) then 
    29752957                    DoH_vap(ig,l) = ( zq(ig,l,igcm_hdo_vap)/
    29762958     &              zq(ig,l,igcm_h2o_vap) )*1./(2.*155.76e-6)
    … …  
    29942976            CALL WRITEDIAGFI(ngrid,'DoH_vap',
    29952977     &                       'D/H ratio in vapor',
    2996      &                       ' ',3,DoH_vap)
     2978     &                       ' ',3,DoH_vap) 
    29972979            CALL WRITEDIAGFI(ngrid,'DoH_ice',
    29982980     &                       'D/H ratio in ice',
    … …  
    30513033              endif
    30523034!A. Pottier
    3053           if (CLFvarying) then !AP14 nebulosity
    3054             call WRITEDIAGFI(ngrid,'totcloudfrac',
    3055      &                       'Total cloud fraction',
     3035              if (CLFvarying) then !AP14 nebulosity
     3036                 call WRITEDIAGFI(ngrid,'totcloudfrac',
     3037     &                'Total cloud fraction',
    30563038     &                       ' ',2,totcloudfrac)
    3057           end if !clf varying
    3058         end if !(water)
    3059 
    3060         if (water.and..not.photochem) then
    3061           iq = nq
    3062 c          write(str2(1:2),'(i2.2)') iq
    3063 c          call WRITEDIAGFI(ngrid,'dqs'//str2,'dqscloud',
    3064 c    &                      'kg.m-2',2,zdqscloud(1,iq))
    3065 c          call WRITEDIAGFI(ngrid,'dqch'//str2,'var chim',
    3066 c    &                      'kg/kg',3,zdqchim(1,1,iq))
    3067 c          call WRITEDIAGFI(ngrid,'dqd'//str2,'var dif',
    3068 c    &                      'kg/kg',3,zdqdif(1,1,iq))
    3069 c          call WRITEDIAGFI(ngrid,'dqa'//str2,'var adj',
    3070 c    &                      'kg/kg',3,zdqadj(1,1,iq))
    3071 c          call WRITEDIAGFI(ngrid,'dqc'//str2,'var c',
    3072 c    &                      'kg/kg',3,zdqc(1,1,iq))
    3073         end if  !(water.and..not.photochem)
    3074       end if !tracer
     3039              endif !clf varying
     3040
     3041           endif !(water)
     3042
     3043
     3044           if (water.and..not.photochem) then
     3045             iq=nq
     3046c            write(str2(1:2),'(i2.2)') iq
     3047c            call WRITEDIAGFI(ngrid,'dqs'//str2,'dqscloud',
     3048c    &                       'kg.m-2',2,zdqscloud(1,iq))
     3049c            call WRITEDIAGFI(ngrid,'dqch'//str2,'var chim',
     3050c    &                       'kg/kg',3,zdqchim(1,1,iq))
     3051c            call WRITEDIAGFI(ngrid,'dqd'//str2,'var dif',
     3052c    &                       'kg/kg',3,zdqdif(1,1,iq))
     3053c            call WRITEDIAGFI(ngrid,'dqa'//str2,'var adj',
     3054c    &                       'kg/kg',3,zdqadj(1,1,iq))
     3055c            call WRITEDIAGFI(ngrid,'dqc'//str2,'var c',
     3056c    &                       'kg/kg',3,zdqc(1,1,iq))
     3057           endif  !(water.and..not.photochem)
     3058         endif
    30753059
    30763060c        ----------------------------------------------------------
    … …  
    30783062c        ----------------------------------------------------------
    30793063
    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))
     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))
    30863070
    30873071#ifndef MESOINI
    … …  
    31133097             call WRITEDIAGFI(ngrid,'dustN','Dust number',
    31143098     &                        'part/kg',3,ndust)
    3115 
     3099             
    31163100             select case (trim(dustiropacity))
    31173101              case ("tes")
    … …  
    31863170             call WRITEDIAGFI(ngrid,'totaldustq','total dust mass',
    31873171     &                        'kg/kg',3,qdusttotal)
    3188 
     3172             
    31893173             select case (trim(dustiropacity))
    31903174              case ("tes")
    … …  
    32113195             end select
    32123196           endif ! (slpwind)
    3213 
     3197           
    32143198           if (scavenging) then
    32153199             call WRITEDIAGFI(ngrid,'ccnq','CCN mass mr',
    … …  
    32773261c        ----------------------------------------------------------
    32783262
    3279          if(callthermos) then
     3263         if(callthermos) then 
    32803264
    32813265            call WRITEDIAGFI(ngrid,"q15um","15 um cooling","K/s",
    … …  
    32993283         endif  !(callthermos)
    33003284
    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 
    33103285c        ----------------------------------------------------------
    33113286c        ----------------------------------------------------------
    … …  
    33173292c        Outputs of thermals
    33183293c        ----------------------------------------------------------
    3319       if (calltherm) then
     3294         if (calltherm) then
     3295
    33203296!        call WRITEDIAGFI(ngrid,'dtke',
    3321 !     &                   'tendance tke thermiques','m**2/s**2',
    3322 !     &                   3,dtke_th)
     3297!     &              'tendance tke thermiques','m**2/s**2',
     3298!     &                         3,dtke_th)
    33233299!        call WRITEDIAGFI(ngrid,'d_u_ajs',
    3324 !     &                   'tendance u thermiques','m/s',
    3325 !     &                   3,pdu_th*ptimestep)
     3300!     &              'tendance u thermiques','m/s',
     3301!     &                         3,pdu_th*ptimestep)
    33263302!        call WRITEDIAGFI(ngrid,'d_v_ajs',
    3327 !     &                   'tendance v thermiques','m/s',
    3328 !     &                   3,pdv_th*ptimestep)
     3303!     &              'tendance v thermiques','m/s',
     3304!     &                         3,pdv_th*ptimestep)
    33293305!        if (tracer) then
    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
     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
    33363312
    33373313        call WRITEDIAGFI(ngrid,'zmax_th',
    3338      &                   'hauteur du thermique','m',
    3339      &                    2,zmax_th)
     3314     &              'hauteur du thermique','m',
     3315     &                         2,zmax_th)
    33403316        call WRITEDIAGFI(ngrid,'hfmax_th',
    3341      &                   'maximum TH heat flux','K.m/s',
    3342      &                   2,hfmax_th)
     3317     &              'maximum TH heat flux','K.m/s',
     3318     &                         2,hfmax_th)
    33433319        call WRITEDIAGFI(ngrid,'wstar',
    3344      &                   'maximum TH vertical velocity','m/s',
    3345      &                   2,wstar)
    3346       end if
     3320     &              'maximum TH vertical velocity','m/s',
     3321     &                         2,wstar)
     3322
     3323         endif
    33473324
    33483325c        ----------------------------------------------------------
    … …  
    33823359c
    33833360c        CALL writeg1d(ngrid,1,fluxsurf_lw,'Fs_ir','W.m-2')
    3384 c        CALL writeg1d(ngrid,1,tsurf,'tsurf','K')
    3385 c        CALL writeg1d(ngrid,1,ps,'ps','Pa')
    3386 c        CALL writeg1d(ngrid,nlayer,zt,'T','K')
     3361c         CALL writeg1d(ngrid,1,tsurf,'tsurf','K')
     3362c         CALL writeg1d(ngrid,1,ps,'ps','Pa')
     3363       
     3364c         CALL writeg1d(ngrid,nlayer,zt,'T','K')
    33873365c        CALL writeg1d(ngrid,nlayer,pu,'u','m.s-1')
    33883366c        CALL writeg1d(ngrid,nlayer,pv,'v','m.s-1')
    … …  
    34053383     &                         0,wstar)
    34063384
    3407          end if ! of if (calltherm)
    3408 
    3409          call WRITEDIAGFI(ngrid,'w','vertical velocity'
     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'                 &
    34103400     &                              ,'m/s',1,pw)
    34113401         call WRITEDIAGFI(ngrid,"q2","q2","kg.m-3",1,q2)
    … …  
    34183408         call WRITEDIAGFI(ngrid,"pplev","Pressure","Pa",1,zplev)
    34193409         call WRITEDIAGFI(ngrid,"rho","rho","kg.m-3",1,rho)
    3420          call WRITEDIAGFI(ngrid,"dtrad","rad. heat. rate",
     3410         call WRITEDIAGFI(ngrid,"dtrad","rad. heat. rate",              &
    34213411     &              "K.s-1",1,dtrad)
    3422          call WRITEDIAGFI(ngrid,'sw_htrt','sw heat. rate',
     3412        call WRITEDIAGFI(ngrid,'sw_htrt','sw heat. rate',
    34233413     &                   'w.m-2',1,zdtsw)
    3424          call WRITEDIAGFI(ngrid,'lw_htrt','lw heat. rate',
     3414        call WRITEDIAGFI(ngrid,'lw_htrt','lw heat. rate',
    34253415     &                   'w.m-2',1,zdtlw)
    34263416         call WRITEDIAGFI(ngrid,"co2ice","co2 ice thickness"
    … …  
    34283418
    34293419         if (igcm_co2.ne.0) then
    3430           call co2sat(ngrid*nlayer,zt,zqsatco2)
     3420          call co2sat(ngrid*nlayer,zt,zplay,zqsatco2)
    34313421           do ig=1,ngrid
    34323422            do l=1,nlayer
    … …  
    34473437c         call WRITEDIAGFI(ngrid,"satuco2","vap in satu","kg/kg",1,
    34483438c     &        satuco2)
    3449 c         call WRITEDIAGFI(ngrid,"riceco2","ice radius","m"
     3439c         call WRITEdiagfi(ngrid,"riceco2","ice radius","m"
    34503440c     &        ,1,riceco2)
    34513441! or output in diagfi.nc (for testphys1d)
    … …  
    35663556     &                 * (zplev(1,l) - zplev(1,l+1)) / g
    35673557               if (hdo) THEN
    3568                  mtotD = mtotD +  zq(1,l,igcm_hdo_vap)
     3558                 mtotD = mtotD +  zq(1,l,igcm_hdo_vap) 
    35693559     &                 * (zplev(1,l) - zplev(1,l+1)) / g
    3570                  icetotD = icetotD +  zq(1,l,igcm_hdo_ice)
     3560                 icetotD = icetotD +  zq(1,l,igcm_hdo_ice) 
    35713561     &                 * (zplev(1,l) - zplev(1,l+1)) / g
    35723562               ENDIF !hdo
    … …  
    35763566                   hdotot = hdotot+mtotD+icetotD
    35773567               ENDIF ! hdo
    3578 
     3568               
    35793569
    35803570             CALL WRITEDIAGFI(ngrid,'h2otot',
    … …  
    36203610            CALL WRITEDIAGFI(ngrid,'DoH_vap',
    36213611     &                       'D/H ratio in vapor',
    3622      &                       ' ',1,DoH_vap)
     3612     &                       ' ',1,DoH_vap) 
    36233613            CALL WRITEDIAGFI(ngrid,'DoH_ice',
    36243614     &                       'D/H ratio in ice',
    … …  
    37153705
    37163706        ENDIF ! hdo
    3717 
     3707     
    37183708        call WRITEDIAGFI(ngrid,"rice","ice radius","m",1,
    37193709     &                    rice)
    … …  
    37253715              endif !clfvarying
    37263716
    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
     3717         ENDIF ! of IF (water)
    37483718         
    37493719ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    37503720ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
     3721
    37513722
    37523723         zlocal(1)=-log(zplay(1,1)/zplev(1,1))* Rnew(1,1)*zt(1,1)/g
    … …  
    37653736
    37663737      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')
    37833738
    37843739! XIOS outputs

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

    r2342 r2343  
    88! Condensation temperature for co2 ice; based on
    99! the saturation in co2sat.F JA17
    10 !--------------------------------------------------
     10!--------------------------------------------------i
    1111
    1212integer, intent(in) :: ngrid,nlay
    … …  
    1414double precision, intent(out), dimension(ngrid,nlay):: tcond ! CO2 condensation temperature     (atm)
    1515double precision:: A,B,pco2
    16 real :: qco2
    1716integer:: ig,l
    1817
    1918A=dlog(1.382d12)
    2019B=-3182.48
    21 qco2=0.
     20
    2221DO l=1,nlay
    2322   DO ig=1,ngrid
Note: See TracChangeset for help on using the changeset viewer.