Changeset 3193


Ignore:
Timestamp:
Jan 30, 2024, 6:36:45 PM (11 months ago)
Author:
afalco
Message:

Pluto PCM:
Added .def examples for Pluto-generic.
AF

Location:
trunk/LMDZ.PLUTO
Files:
7 added
9 edited

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.PLUTO/libf/phypluto/aerosol_mod.F90

    r3184 r3193  
    3838  use radinc_h, only: naerkind
    3939  use tracer_h, only: n_rgcs, nqtot, is_rgcs
    40   use callkeys_mod, only: aeron2, dusttau, aeroh2so4, &
     40  use callkeys_mod, only: aeron2, dusttau, &
    4141                          aeroback2lay, aeronh3, nlayaero, aeronlay, &
    42                           aeroaurora, aerogeneric, &
    43                           aerovenus1, aerovenus2, aerovenus2p, &
    44                           aerovenus3, aerovenusUV
     42                          aeroaurora, aerogeneric
    4543
    4644  IMPLICIT NONE
    … …  
    8987  if (is_master) write(*,*) '--- Dust aerosol = ', iaero_dust
    9088
    91   if (aeroh2so4) then
    92      ia=ia+1
    93      iaero_h2so4=ia
    94   endif
    95   if (is_master) write(*,*) '--- H2SO4 aerosol = ', iaero_h2so4
    96      
     89     
    9790  if (aeroback2lay) then
    9891     ia=ia+1

  • trunk/LMDZ.PLUTO/libf/phypluto/callcorrk.F90

    r3184 r3193  
    3737                              diagdtau,kastprof,strictboundcorrk,specOLR, &
    3838                              tplanckmin,tplanckmax,global1d, &
    39                               generic_condensation, aerovenus
     39                              generic_condensation
    4040      use optcv_mod, only: optcv
    4141      use optci_mod, only: optci
    … …  
    808808      plevrad(1) = 0.
    809809!      plevrad(2) = 0.   !! JL18 enabling this line puts the radiative top at p=0 which was the idea before, but does not seem to perform best after all.
    810       if (aerovenus) then
    811 !!  GG19 modified below after SL routines
    812         plevrad(2) = 0.
    813       endif
    814 
     810     
    815811      tlevrad(1) = tlevrad(2)
    816812      tlevrad(2*nlayer+1)=tsurf(ig)
    817813     
    818814      pmid(1) = pplay(ig,nlayer)/scalep   
    819       if (aerovenus) then
    820 !! GG19 modified below after SL routines
    821         pmid(1) = max(pgasmin,0.0001*plevrad(3))
    822       endif
    823815      pmid(2) =  pmid(1)
    824816

  • trunk/LMDZ.PLUTO/libf/phypluto/callkeys_mod.F90

    r3184 r3193  
    66      logical,save :: calladj,calltherm,n2cond,callsoil
    77!$OMP THREADPRIVATE(calladj,calltherm,n2cond,callsoil)
    8       logical,save :: season,diurnal,tlocked,rings_shadow,lwrite
    9 !$OMP THREADPRIVATE(season,diurnal,tlocked,rings_shadow,lwrite)
     8      logical,save :: season,diurnal,lwrite
     9!$OMP THREADPRIVATE(season,diurnal,lwrite)
    1010      logical,save :: callgasvis,continuum,graybody
    1111!$OMP THREADPRIVATE(callgasvis,continuum,graybody)
    … …  
    3737      logical,save :: sedimentation
    3838      logical,save :: generic_condensation
    39       logical,save :: generic_rain
    4039!$OMP THREADPRIVATE(varactive,varfixed,sedimentation,generic_condensation,generic_rain)
    41       logical,save :: water ,watercond, waterrain, moistadjustment
    42 !$OMP THREADPRIVATE(water, watercond, waterrain, moistadjustment)
    4340      logical,save :: aeron2, aeroh2o, aeroh2so4, aeroback2lay
    4441!$OMP THREADPRIVATE(aeron2, aeroh2o, aeroh2so4, aeroback2lay)
    4542      logical,save :: aeronh3, aeronlay, aeroaurora
    4643!$OMP THREADPRIVATE(aeronh3,aeronlay,aeroaurora)
    47 
    48       logical,save :: aerovenus ! master flag for "Venus-like" aerosol additions
    49 !$OMP THREADPRIVATE(aerovenus)
    50       ! detailed sub-options when with "Venus-like" aerosol additions
    51       logical,save :: aerovenus1, aerovenus2, aerovenus2p, aerovenus3, aerovenusUV
    52 !$OMP THREADPRIVATE(aerovenus1, aerovenus2, aerovenus2p, aerovenus3, aerovenusUV)
    5344
    5445      logical,save :: aerofixn2
    … …  
    6354      logical,save :: jonline
    6455      logical,save :: depos
    65       logical,save :: haze
    6656!$OMP THREADPRIVATE(photoheat,jonline,depos)
    67       logical,save :: calllott_nonoro
    68 !$OMP THREADPRIVATE(calllott_nonoro)
     57
     58!! Pluto-specific variables
     59      logical,save :: haze_radproffix,haze_proffix
     60!$OMP THREADPRIVATE(haze_radproffix,haze_proffix)     
     61      logical,save :: haze,fasthaze,changeti,changetid,aerohaze         
     62!$OMP THREADPRIVATE(haze,fasthaze,changeti,changetid,aerohaze)     
     63      logical,save :: fast,metcloud,monoxcloud,glaflow,triton,paleo,nlte,strobel
     64!$OMP THREADPRIVATE(fast,metcloud,monoxcloud,glaflow,triton,paleo,nlte,strobel     
     65      logical,save :: kbo,nbsub,cooling,mode_n2,thres_n2ice,thres_coice
     66!$OMP THREADPRIVATE(kbo,nbsub,cooling,mode_n2,thres_n2ice,thres_coice)     
     67      logical,save :: deltab,nb_monomer                                 
     68!$OMP THREADPRIVATE(deltab,nb_monomer)     
     69      logical,save :: metlateq,tholateq,metls1,metls2                   
     70!$OMP THREADPRIVATE(metlateq,tholateq,metls1,metls2)     
     71      logical,save :: mode_ch4,mode_tholins                             
     72!$OMP THREADPRIVATE(mode_ch4,mode_tholins)     
     73      logical,save :: source_haze,ch4lag,latlag,vmrlag,kfix,hazeconservch4     
     74!$OMP THREADPRIVATE(source_haze,ch4lag,latlag,vmrlag,kfix,hazeconservch4)     
     75      logical,save :: fracsource,latsource,lonsource                   
     76!$OMP THREADPRIVATE(fracsource,latsource,lonsource)     
     77      logical,save :: spelon1,spelon2,spelat1,spelat2,specalb           
     78!$OMP THREADPRIVATE(spelon1,spelon2,spelat1,spelat2,specalb)     
     79      logical,save :: assymflux,mode_hs,tsurfmax,albmin_ch4
     80!$OMP THREADPRIVATE(assymflux,mode_hs,tsurfmax,albmin_ch4)     
     81      logical,save :: feedback_met,thres_ch4ice,fdch4_latn,fdch4_lats   
     82!$OMP THREADPRIVATE(feedback_met,thres_ch4ice,fdch4_latn,fdch4_lats)     
     83      logical,save :: globmean1d,kmix_proffix,rad_haze         
     84!$OMP THREADPRIVATE(globmean1d,kmix_proffix,rad_haze)     
     85      logical,save :: fdch4_lone,fdch4_lonw,fdch4_maxalb,fdch4_depalb,fdch4_finalb
     86!$OMP THREADPRIVATE(fdch4_lone,fdch4_lonw,fdch4_maxalb,fdch4_depalb,fdch4_finalb     
     87      logical,save :: tholatn,tholats,tholone,tholonw                   
     88!$OMP THREADPRIVATE(tholatn,tholats,tholone,tholonw)     
     89      logical,save :: fdch4_ampl,fdch4_maxice,condmetsurf,condcosurf,vertdiff   
     90!$OMP THREADPRIVATE(fdch4_ampl,fdch4_maxice,condmetsurf,condcosurf,vertdiff)     
     91      logical,save :: convergeps,conservn2,patchflux,condensn2,no_n2frost
     92!$OMP THREADPRIVATE(convergeps,conservn2,patchflux,condensn2,no_n2frost     
     93
     94
    6995      logical,save :: global1d
    7096      real,save    :: szangle

  • trunk/LMDZ.PLUTO/libf/phypluto/callsedim.F

    r3184 r3193  
    88      USE tracer_h, only : igcm_n2_ice,radius,rho_q
    99      use comcstfi_mod, only: g
    10       use callkeys_mod, only : water
    1110
    1211      IMPLICIT NONE
    … …  
    7574        firstcall=.false.
    7675        ! add some tests on presence of required tracers/aerosols:
    77     !     if (water.and.(igcm_h2o_ice.eq.0)) then
    78     !         write(*,*) "callsedim error: water=.true.",
    79     !  &                 " but igcm_h2o_ice=0"
    80     !       stop
    81     !     endif
    8276        ! allocate "naerkind" size local arrays (which are also
    8377        ! "saved" so that this is done only once in for all even if
    … …  
    120114! Sedimentation
    121115!======================================================================
    122 ! Water         
    123 !        if (water.and.(iaero_h2o.ne.0).and.(iq.eq.igcm_h2o_ice)) then
    124 !             ! compute radii for h2o_ice
    125 !              call h2o_reffrad(ngrid,nlay,zqi(1,1,igcm_h2o_ice),zt,
    126 !      &                reffrad(1,1,iaero_h2o),nueffrad(1,1,iaero_h2o))
    127 !             ! call sedimentation for h2o_ice
    128 !              call newsedim(ngrid,nlay,ngrid*nlay,ptimestep,
    129 !      &            pplev,masse,epaisseur,zt,reffrad(1,1,iaero_h2o),
    130 !      &            rho_q(iq),zqi(1,1,igcm_h2o_ice),wq,iq)
     116! Water      !AF24: deleted
    131117
    132118! ! General Case

  • trunk/LMDZ.PLUTO/libf/phypluto/condense_n2.F90

    r3187 r3193  
    1       subroutine condense_n2(ngrid,nlayer,nq,ptimestep, &
    2           pcapcal,pplay,pplev,ptsrf,pt,                  &
    3           pdt,pdtsrf,pq,pdq,                            &
    4           pqsurf,pdqsurfc,albedo,pemisurf,              &
    5           albedo_bareground,albedo_n2_ice_SPECTV,      &
    6           pdtc,pdtsrfc,pdpsrfc,pdqc)
    7 
    8       use radinc_h, only : L_NSPECTV
    9       use gases_h, only: gfrac, igas_n2
    10       use radii_mod, only : n2_reffrad
    11       use aerosol_mod, only : iaero_n2
    12       USE surfdat_h, only: emisice, emissiv
    13       USE geometry_mod, only: latitude ! in radians
    14       USE tracer_h, only: noms, rho_n2
    15       use comcstfi_mod, only: g, r, cpp
    16 
    17       implicit none
     1subroutine condens_n2(klon,nlayer,nq,ptimestep, &
     2      pcapcal,pplay,pplev,ptsrf,pt, &
     3      pphi,pdt,pdu,pdv,pdtsrf,pu,pv,pq,pdq, &
     4      picen2,psolaralb,pemisurf, &
     5      pdtc,pdtsrfc,pdpsrf,pduc,pdvc, &
     6      pdqc,pdicen2)
     7
     8  use radinc_h, only : naerkind
     9  use comgeomfi_h
     10  use comcstfi_mod, only: g, r, cpp
     11  USE surfdat_h, only: phisfi
     12  USE tracer_h, only: noms, igcm_n2
     13  USE callkeys_mod, only: fast,ch4lag,latlag,lw_n2,nbsub,no_n2frost,tsurfmax,kmixmin,source_haze,vmrlag
     14  USE dimphy, only : klon,klev
     15
     16
     17  implicit none
    1818
    1919!==================================================================
    2020!     Purpose
    2121!     -------
    22 !     Condense and/or sublime N2 ice on the ground and in the atmosphere, and sediment the ice.
     22!     Condense and/or sublime N2 ice on the ground and in the
     23!     atmosphere, and sediment the ice.
    2324! 
    2425!     Inputs
    25 !     ------
    26 !     ngrid                 Number of vertical columns.
    27 !     nlayer                Number of vertical layers.
    28 !     nq                    Number of tracers.
    29 !     ptimestep             Duration of the physical timestep (s).
    30 !     pplay(ngrid,nlayer)   Pressure layers (Pa).
    31 !     pplev(ngrid,nlayer+1) Pressure levels (Pa).
    32 !     pt(ngrid,nlayer)      Atmospheric Temperatures (K).
    33 !     ptsrf(ngrid)          Surface temperatures (K).
    34 !     pq(ngrid,nlayer,nq)   Atmospheric tracers mixing ratios (kg/kg of air).
    35 !     pqsurf(ngrid,nq)      Surface tracers (kg/m2).
     26!     ------
     27!     klon                 Number of vertical columns
     28!     nlayer                Number of layers
     29!     pplay(klon,nlayer)   Pressure layers
     30!     pplev(klon,nlayer+1) Pressure levels
     31!     pt(klon,nlayer)      Temperature (in K)
     32!     ptsrf(klon)          Surface temperature
    3633!     
    37 !     pdt(ngrid,nlayer)     Time derivative before condensation/sublimation of pt.
    38 !     pdtsrf(ngrid)         Time derivative before condensation/sublimation of ptsrf.
    39 !     pdq(ngrid,nlayer,nq)  Time derivative before condensation/sublimation of
    40 !
    41 !     albedo_bareground(ngrid)           Albedo of the bare ground.
    42 !     albedo_n2_ice_SPECTV(L_NSPECTV)   Spectral albedo of N2 ice.
     34!     pdt(klon,nlayermx)   Time derivative before condensation/sublimation of pt
     35!     pdtsrf(klon)         Time derivative before condensation/sublimation of ptsrf
     36!     picen2(klon)         n2 ice at the surface (kg/m2)
    4337!     
    4438!     Outputs
    4539!     -------
    46 !     pdpsrfc(ngrid)          \       Contribution of condensation/sublimation
    47 !     pdtc(ngrid,nlayer)       \            to the time derivatives of
    48 !     pdtsrfc(ngrid)           /     Surface Pressure, Atmospheric Temperatures,
    49 !     pdqsurfc(ngrid)         /         Surface Temperatures, Surface Tracers,
    50 !     pdqc(ngrid,nlayer,nq)  /                and Atmospheric Tracers.*
    51 !
    52 !     pemisurf(ngrid)              Emissivity of the surface.
     40!     pdpsrf(klon)         \  Contribution of condensation/sublimation
     41!     pdtc(klon,nlayermx)  /  to the time derivatives of Ps, pt, and ptsrf
     42!     pdtsrfc(klon)       /
     43!     pdicen2(klon)         Tendancy n2 ice at the surface (kg/m2)
    5344!     
    5445!     Both
    5546!     ----
    56 !     albedo(ngrid,L_NSPECTV)      Spectral albedo of the surface.
     47!     psolaralb(klon)      Albedo at the surface
     48!     pemisurf(klon)       Emissivity of the surface
    5749!
    5850!     Authors
    5951!     -------
    60 !     Francois Forget (1996)
     52!     Francois Forget (1996,2013)
    6153!     Converted to Fortran 90 and slightly modified by R. Wordsworth (2009)
    62 !     Includes simplifed nucleation by J. Leconte (2011)
     54!     
    6355!     
    6456!==================================================================
    6557
    66 !--------------------------
    67 !        Arguments
    68 !--------------------------
    69 
    70 
    71       INTEGER,INTENT(IN) :: ngrid
    72       INTEGER,INTENT(IN) :: nlayer
    73       INTEGER,INTENT(IN) :: nq
    74       REAL,INTENT(IN) :: ptimestep
    75       REAL,INTENT(IN) :: pcapcal(ngrid)
    76       REAL,INTENT(IN) :: pplay(ngrid,nlayer)
    77       REAL,INTENT(IN) :: pplev(ngrid,nlayer+1)
    78       REAL,INTENT(IN) :: ptsrf(ngrid)
    79       REAL,INTENT(IN) :: pt(ngrid,nlayer)
    80       REAL,INTENT(IN) :: pdt(ngrid,nlayer)
    81       REAL,INTENT(IN) :: pdtsrf(ngrid)
    82       REAL,INTENT(IN) :: pq(ngrid,nlayer,nq)
    83       REAL,INTENT(IN) :: pqsurf(ngrid,nq)
    84       REAL,INTENT(IN) :: pdq(ngrid,nlayer,nq)
    85       REAL,INTENT(IN) :: albedo_bareground(ngrid)
    86       REAL,INTENT(IN) :: albedo_n2_ice_SPECTV(L_NSPECTV)
    87       REAL,INTENT(INOUT) :: albedo(ngrid,L_NSPECTV)
    88       REAL,INTENT(OUT) :: pemisurf(ngrid)
    89       REAL,INTENT(OUT) :: pdtc(ngrid,nlayer)
    90       REAL,INTENT(OUT) :: pdtsrfc(ngrid)
    91       REAL,INTENT(OUT) :: pdpsrfc(ngrid)
    92       REAL,INTENT(OUT) :: pdqc(ngrid,nlayer,nq)
    93       REAL,INTENT(OUT) :: pdqsurfc(ngrid)
    94 
    95 !------------------------------
    96 !       Local variables
    97 !------------------------------
    98 
    99       INTEGER l,ig,icap,ilay,iq,nw,igas,it
    100 
    101       REAL reffrad(ngrid,nlayer) ! Radius (m) of the N2 ice particles.
    102       REAL*8 zt(ngrid,nlayer)    ! Updated Atmospheric Temperatures (K).
    103       REAL ztsrf(ngrid)          ! Updated Surface Temperatures (K).
    104       REAL zq(ngrid,nlayer,nq)   ! Updated Atmospheric tracers mixing ratios (kg/kg of air).
    105       REAL picen2(ngrid)        ! Updated Surface Tracer (kg/m2).
    106       REAL ztcond (ngrid,nlayer) ! Atmospheric Temperatures of condensation of N2.
    107       REAL ztnuc (ngrid,nlayer)  ! Atmospheric Nucleation Temperatures.
    108       REAL ztcondsol(ngrid)      ! Temperatures of condensation of N2 at the surface.
    109       REAL zcondices(ngrid)      ! Condensation rate on the ground (kg/m2/s).
    110       REAL zfallice(ngrid)       ! Flux of ice falling on the surface (kg/m2/s).
    111       REAL Mfallice(ngrid)       ! Total amount of ice fallen to the ground during the timestep (kg/m2).
    112       REAL wq(ngrid,nlayer+1)    ! Total amount of ice fallen to the ground during the timestep (kg/m2).
    113       REAL subptimestep          ! Duration of the subtimestep (s) for the sedimentation.
    114       Integer Ntime              ! Number of subtimesteps.
    115       REAL masse (ngrid,nlayer)  ! Mass of atmospheric layers (kg/m2)
    116       REAL w(ngrid,nlayer,nq)    !
    117       REAL vstokes,reff          !
    118       REAL ppn2                 !
    119 
    120 
    121 !------------------------------------------
    122 !         Saved local variables
    123 !------------------------------------------
    124 
    125 
    126       REAL,SAVE :: latcond=2.5e5
    127       REAL,SAVE :: ccond
    128       REAL,SAVE,ALLOCATABLE,DIMENSION(:) :: emisref
    129 !$OMP THREADPRIVATE(latcond,ccond,emisref)
    130       LOGICAL,SAVE :: firstcall=.true.
    131 !$OMP THREADPRIVATE(firstcall)
    132       INTEGER,SAVE :: i_n2ice=0      ! n2 ice
    133 !$OMP THREADPRIVATE(i_n2ice)
    134       CHARACTER(LEN=20) :: tracername ! to temporarily store text
    135 
    136 
    137 !------------------------------------------------
    138 !       Initialization at the first call
    139 !------------------------------------------------
    140 
    141 
    142       IF (firstcall) THEN
    143 
    144 !         ALLOCATE(emisref(ngrid))
    145 !         ! Find N2 ice tracer.
    146 !         do iq=1,nq
    147 !            tracername=noms(iq)
    148 !            if (tracername.eq."n2_ice") then
    149 !               i_n2ice=iq
    150 !            endif
    151 !         enddo
    152          
    153 !         write(*,*) "condense_n2: i_n2ice=",i_n2ice       
    154 
    155 !         if((i_n2ice.lt.1))then
    156 !            print*,'In condens_cloud but no N2 ice tracer, exiting.'
    157 !            print*,'Still need generalisation to arbitrary species!'
    158 !            stop
    159 !         endif
    160 
    161          ccond=cpp/(g*latcond)
    162          print*,'In condens_cloud: ccond=',ccond,' latcond=',latcond
    163 
    164          ! Prepare special treatment if gas is not pure N2
    165          ! if (addn2) then
    166          !    m_n2   = 44.01E-3 ! N2 molecular mass (kg/mol)   
    167          !    m_non2 = 28.02E-3 ! N2 molecular mass (kg/mol) 
    168          ! Compute A and B coefficient use to compute
    169          ! mean molecular mass Mair defined by
    170          ! 1/Mair = q(in2)/m_n2 + (1-q(in2))/m_non2
    171          ! 1/Mair = A*q(in2) + B
    172          ! A = (1/m_n2 - 1/m_non2)
    173          ! B = 1/m_non2
    174          ! endif
    175 
    176          ! Minimum N2 mixing ratio below which mixing occurs with layer above : qn2min =0.75 
    177 
    178            firstcall=.false.
    179       ENDIF
    180 
    181 
    182 !------------------------------------------------
    183 !        Tendencies initially set to 0
    184 !------------------------------------------------
    185 
    186 
    187       pdqc(1:ngrid,1:nlayer,1:nq)     = 0.
    188       pdtc(1:ngrid,1:nlayer)          = 0.
    189       zq(1:ngrid,1:nlayer,1:nq)       = 0.
    190       zt(1:ngrid,1:nlayer)            = 0.
    191       Mfallice(1:ngrid)               = 0.     
    192       zfallice(1:ngrid)               = 0.
    193       zcondices(1:ngrid)              = 0.
    194       pdtsrfc(1:ngrid)                = 0.
    195       pdpsrfc(1:ngrid)                = 0.
    196       pdqsurfc(1:ngrid)               = 0.
    197 
    198 
    199 !----------------------------------
     58!-----------------------------------------------------------------------
     59!     Arguments
     60
     61  INTEGER klon, nlayer, nq
     62
     63  REAL ptimestep
     64  REAL pplay(klon,nlayer),pplev(klon,nlayer+1)
     65  REAL pcapcal(klon)
     66  REAL pt(klon,nlayer)
     67  REAL ptsrf(klon),flu1(klon),flu2(klon),flu3(klon)
     68  REAL pphi(klon,nlayer)
     69  REAL pdt(klon,nlayer),pdtsrf(klon),pdtc(klon,nlayer)
     70  REAL pdtsrfc(klon),pdpsrf(klon)
     71  REAL picen2(klon),psolaralb(klon),pemisurf(klon)
     72 
     73
     74  REAL pu(klon,nlayer) ,  pv(klon,nlayer)
     75  REAL pdu(klon,nlayer) , pdv(klon,nlayer)
     76  REAL pduc(klon,nlayer) , pdvc(klon,nlayer)
     77  REAL pq(ngridmx,nlayer,nq),pdq(klon,nlayer,nq)
     78  REAL pdqc(klon,nlayer,nq)
     79
     80!-----------------------------------------------------------------------
     81!     Local variables
     82
     83  INTEGER l,ig,ilay,it,iq,ich4_gas
     84
     85  REAL*8 zt(ngridmx,nlayermx)
     86  REAL tcond_n2
     87  REAL pcond_n2
     88  REAL glob_average2d         ! function
     89  REAL zqn2(ngridmx,nlayermx) ! N2 MMR used to compute Tcond/zqn2
     90  REAL ztcond (ngridmx,nlayermx)
     91  REAL ztcondsol(ngridmx),zfallheat
     92  REAL pdicen2(ngridmx)
     93  REAL zcondicea(ngridmx,nlayermx), zcondices(ngridmx)
     94  REAL zfallice(ngridmx,nlayermx+1)
     95  REAL zmflux(nlayermx+1)
     96  REAL zu(nlayermx),zv(nlayermx)
     97  REAL zq(nlayermx,nqmx),zq1(nlayermx)
     98  REAL ztsrf(ngridmx)
     99  REAL ztc(nlayermx), ztm(nlayermx+1)
     100  REAL zum(nlayermx+1) , zvm(nlayermx+1)
     101  REAL zqm(nlayermx+1,nqmx),zqm1(nlayermx+1)
     102  LOGICAL condsub(ngridmx)
     103  REAL subptimestep
     104  Integer Ntime
     105  real masse (nlayermx),w(nlayermx+1)
     106  real wq(ngridmx,nlayermx+1)
     107  real vstokes,reff
     108  real dWtotsn2
     109  real condnconsn2(ngridmx)
     110  real nconsMAXn2
     111!     Special diagnostic variables
     112  real tconda1(ngridmx,nlayermx)
     113  real tconda2(ngridmx,nlayermx)
     114  real zdtsig (ngridmx,nlayermx)
     115  real zdtlatent (ngridmx,nlayermx)
     116  real zdt (ngridmx,nlayermx)
     117  REAL albediceF(ngridmx)
     118  SAVE albediceF
     119  INTEGER nsubtimestep,itsub    !number of subtimestep when calling vl1d
     120  REAL subtimestep        !ptimestep/nsubtimestep
     121  REAL dtmax             
     122
     123  REAL zplevhist(ngridmx)
     124  REAL zplevnew(ngridmx)
     125  REAL zplev(ngridmx)
     126  REAL zpicen2(ngridmx)
     127  REAL ztsrfhist(ngridmx)
     128  REAL zdtsrf(ngridmx)
     129  REAL globzplevnew
     130
     131  REAL vmrn2(ngridmx)
     132  SAVE vmrn2
     133  REAL stephan   
     134  DATA stephan/5.67e-08/  ! Stephan Boltzman constant
     135  SAVE stephan
     136!-----------------------------------------------------------------------
     137!     Saved local variables
     138
     139!      REAL latcond
     140  REAL ccond
     141  REAL cpice  ! for atm condensation
     142  SAVE cpice
     143!      SAVE latcond,ccond
     144  SAVE ccond
     145
     146  LOGICAL firstcall
     147  SAVE firstcall
     148  REAL SSUM
     149  EXTERNAL SSUM
     150
     151!      DATA latcond /2.5e5/
     152!      DATA latcond /1.98e5/
     153  DATA cpice /1300./
     154  DATA firstcall/.true./
     155
     156  INTEGER,SAVE :: i_n2ice=0       ! n2 ice
     157  CHARACTER(LEN=20) :: tracername ! to temporarily store text
     158  logical olkin  ! option to prevent N2 ice effect in the south
     159  DATA olkin/.false./
     160  save olkin
     161
     162  CHARACTER(len=10) :: tname
     163
     164!-----------------------------------------------------------------------
     165
     166!     Initialisation
     167  IF (firstcall) THEN
     168     ccond=cpp/(g*lw_n2)
     169     print*,'In condens_n2cloud: ccond=',ccond,' latcond=',lw_n2
     170
     171     ! calculate global mean surface pressure for the fast mode
     172     IF (fast) THEN
     173        DO ig=1,ngridmx
     174           kp(ig) = exp(-phisfi(ig)/(r*38.))
     175        ENDDO
     176        p00=glob_average2d(kp) ! mean pres at ref level
     177     ENDIF
     178
     179     vmrn2(:) = 1.
     180     IF (ch4lag) then
     181        DO ig=1,ngridmx
     182           if (lati(ig)*180./pi.ge.latlag) then
     183              vmrn2(ig) = vmrlag
     184           endif
     185        ENDDO
     186     ENDIF   
     187     IF (no_n2frost) then
     188        DO ig=1,ngridmx
     189           if (picen2(ig).eq.0.) then
     190              vmrn2(ig) = 1.e-15
     191           endif
     192        ENDDO
     193     ENDIF   
     194     firstcall=.false.
     195  ENDIF
     196
     197!-----------------------------------------------------------------------
     198!     Calculation of n2 condensation / sublimation
     199
     200!     Variables used:
     201!     picen2(klon)            : amount of n2 ice on the ground     (kg/m2)
     202!     zcondicea(klon,nlayermx): condensation rate in layer l       (kg/m2/s)
     203!     zcondices(klon)         : condensation rate on the ground    (kg/m2/s)
     204!     zfallice(klon,nlayermx) : amount of ice falling from layer l (kg/m2/s)
     205!     zdtlatent(klon,nlayermx): dT/dt due to phase changes         (K/s)
     206
     207!     Tendencies initially set to 0
     208  zcondices(1:klon) = 0.
     209  pdtsrfc(1:klon) = 0.
     210  pdpsrf(1:klon) = 0.
     211  ztsrfhist(1:klon) = 0.
     212  condsub(1:klon) = .false.
     213  pdicen2(1:klon) = 0.
     214  zfallheat=0
     215  pdqc(1:klon,1:nlayer,1:nq)=0
     216  pdtc(1:klon,1:nlayer)=0
     217  pduc(1:klon,1:nlayer)=0
     218  pdvc(1:klon,1:nlayer)=0
     219  zfallice(1:klon,1:nlayer+1)=0
     220  zcondicea(1:klon,1:nlayer)=0
     221  zdtlatent(1:klon,1:nlayer)=0
     222  zt(1:klon,1:nlayer)=0.
     223
     224!-----------------------------------------------------------------------
    200225!     Atmospheric condensation
    201 !----------------------------------
    202 
    203 
    204 !     Compute N2 Volume mixing ratio
    205 !     -------------------------------
    206 !      if (addn2) then
    207 !         DO l=1,nlayer
    208 !            DO ig=1,ngrid
    209 !              qn2=pq(ig,l,in2)+pdq(ig,l,in2)*ptimestep
    210 !              Mean air molecular mass = 1/(q(in2)/m_n2 + (1-q(in2))/m_non2)
    211 !              mmean=1/(A*qn2 +B)
    212 !              vmr_n2(ig,l) = qn2*mmean/m_n2
    213 !            ENDDO
    214 !         ENDDO
    215 !      else
    216 !         DO l=1,nlayer
    217 !            DO ig=1,ngrid
    218 !              vmr_n2(ig,l)=0.5
    219 !            ENDDO
    220 !         ENDDO
    221 !      end if
    222 
    223 
    224       ! Forecast the atmospheric frost temperature 'ztcond' and nucleation temperature 'ztnuc'.
    225 !      DO l=1,nlayer
    226 !         DO ig=1,ngrid
    227 !            ppn2=gfrac(igas_N2)*pplay(ig,l)
    228 !            call get_tcond_n2(ppn2,ztcond(ig,l))
    229 !            call get_tnuc_n2(ppn2,ztnuc(ig,l))
    230 !         ENDDO
    231 !      ENDDO
     226
     227!     Condensation / sublimation in the atmosphere
     228!     --------------------------------------------
     229!      (calcul of zcondicea , zfallice and pdtc)
     230
     231  zt(1:klon,1:nlayer)=pt(1:klon,1:nlayer)+ pdt(1:klon,1:nlayer)*ptimestep
     232  if (igcm_n2.ne.0) then
     233     zqn2(1:klon,1:nlayer) = 1. ! & temporaire
     234!        zqn2(1:klon,1:nlayer)= pq(1:klon,1:nlayer,igcm_n2) + pdq(1:klon,1:nlayer,igcm_n2)*ptimestep
     235  else
     236     zqn2(1:klon,1:nlayer) = 1.
     237  end if
     238 
     239  if (.not.fast) then
     240!     Forecast the atmospheric frost temperature 'ztcond' with function tcond_n2
     241    DO l=1,nlayer
     242      DO ig=1,klon
     243          ztcond (ig,l) = tcond_n2(pplay(ig,l),zqn2(ig,l))
     244      ENDDO
     245    ENDDO
     246
     247    DO l=nlayer,1,-1
     248      DO ig=1,klon
     249       pdtc(ig,l)=0.  ! final tendancy temperature set to 0
     250
     251       IF((zt(ig,l).LT.ztcond(ig,l)).or.(zfallice(ig,l+1).gt.0))THEN
     252           condsub(ig)=.true.  !condensation at level l
     253           IF (zfallice(ig,l+1).gt.0) then
     254             zfallheat=zfallice(ig,l+1)*& 
     255            (pphi(ig,l+1)-pphi(ig,l) +&
     256           cpice*(ztcond(ig,l+1)-ztcond(ig,l)))/lw_n2
     257           ELSE
     258               zfallheat=0.
     259           ENDIF
     260           zdtlatent(ig,l)=(ztcond(ig,l) - zt(ig,l))/ptimestep
     261           zcondicea(ig,l)=(pplev(ig,l)-pplev(ig,l+1))&
     262                         *ccond*zdtlatent(ig,l)- zfallheat
     263!              Case when the ice from above sublimes entirely
     264!              """""""""""""""""""""""""""""""""""""""""""""""
     265           IF ((zfallice(ig,l+1).lt.-zcondicea(ig,l)) &
     266                .AND. (zfallice(ig,l+1).gt.0)) THEN
     267
     268              zdtlatent(ig,l)=(-zfallice(ig,l+1)+zfallheat)/&     
     269               (ccond*(pplev(ig,l)-pplev(ig,l+1)))
     270              zcondicea(ig,l)= -zfallice(ig,l+1)
     271           END IF
     272
     273           zfallice(ig,l) = zcondicea(ig,l)+zfallice(ig,l+1)
    232274     
    233       ! Initialize zq and zt at the beginning of the sub-timestep loop and qsurf.
    234 !      DO ig=1,ngrid
    235 !         picen2(ig)=pqsurf(ig,i_n2ice)
    236 !         DO l=1,nlayer
    237 !            zt(ig,l)=pt(ig,l)
    238 !            zq(ig,l,i_n2ice)=pq(ig,l,i_n2ice)
    239 !            IF( zq(ig,l,i_n2ice).lt.-1.e-6 ) THEN
    240 !               print*,'Uh-oh, zq = ',zq(ig,l,i_n2ice),'at ig,l=',ig,l
    241 !               if(l.eq.1)then
    242 !                  print*,'Perhaps the atmosphere is collapsing on surface...?'
    243 !               endif
    244 !            END IF
    245 !         ENDDO
    246 !      ENDDO
    247 
    248       ! Calculate the mass of each atmospheric layer (kg.m-2)
    249 !      do  ilay=1,nlayer
    250 !         DO ig=1,ngrid
    251 !            masse(ig,ilay)=(pplev(ig,ilay) - pplev(ig,ilay+1)) /g
    252 !         end do
    253 !      end do
     275        END IF
     276     
     277      ENDDO
     278    ENDDO
     279  endif  ! not fast
     280
     281!-----------------------------------------------------------------------
     282!     Condensation/sublimation on the ground
     283!     (calculation of zcondices and pdtsrfc)
     284
     285!     Adding subtimesteps : in fast version, pressures too low lead to
     286!     instabilities.
     287  IF (fast) THEN
     288     IF (pplev(1,1).gt.0.3) THEN
     289         nsubtimestep= 1 
     290     ELSE
     291         nsubtimestep= nbsub !max(nint(ptimestep/dtmax),1)
     292     ENDIF
     293  ELSE
     294     nsubtimestep= 1 ! if more then kp and p00 have to be calculated
     295                     ! for nofast mode
     296  ENDIF
     297  subtimestep=ptimestep/float(nsubtimestep)
     298
     299  do itsub=1,nsubtimestep
     300     ! first loop : getting zplev, ztsurf
     301     IF (itsub.eq.1) then
     302       DO ig=1,klon
     303          zplev(ig)=pplev(ig,1)
     304          ztsrfhist(ig)=ptsrf(ig) + pdtsrf(ig)*ptimestep
     305          ztsrf(ig)=ptsrf(ig) + pdtsrf(ig)*subtimestep    !!
     306          zpicen2(ig)=picen2(ig)
     307       ENDDO
     308     ELSE
     309     ! additional loop :
     310     ! 1)  pressure updated
     311     ! 2)  direct redistribution of pressure over the globe
     312     ! 3)  modification pressure for unstable cases
     313     ! 4)  pressure update to conserve mass
     314     ! 5)  temperature updated with radiative tendancies
     315       DO ig=1,klon
     316          zplevhist(ig)=zplev(ig)
     317          zplevnew(ig)=zplev(ig)+pdpsrf(ig)*subtimestep   ! 1)
     318          !IF (zplevnew(ig).le.0.0001) then
     319          !   zplevnew(ig)=0.0001*kp(ig)/p00
     320          !ENDIF
     321       ENDDO
     322       ! intermediaire de calcul: valeur moyenne de zplevnew (called twice in the code)
     323       globzplevnew=glob_average2d(zplevnew)
     324       DO ig=1,klon
     325         zplev(ig)=kp(ig)*globzplevnew/p00  ! 2)
     326       ENDDO
     327       DO ig=1,klon     ! 3) unstable case condensation and sublimation: zplev=zplevhist
     328          IF (((pdpsrf(ig).lt.0.).and.(tcond_n2(zplev(ig),zqn2(ig,1)).le.ztsrf(ig))).or.  &
     329          ((pdpsrf(ig).gt.0.).and.(tcond_n2(zplev(ig),zqn2(ig,1)).ge.ztsrf(ig)))) then
     330             zplev(ig)=zplevhist(ig)
     331          ENDIF
     332          zplevhist(ig)=zplev(ig)
     333       ENDDO
     334       zplev=zplev*globzplevnew/glob_average2d(zplevhist)   ! 4)
     335       DO ig=1,klon    ! 5)
     336         zdtsrf(ig)=pdtsrf(ig) + (stephan/pcapcal(ig))*(ptsrf(ig)**4-ztsrf(ig)**4)
     337         ztsrf(ig)=ztsrf(ig)+pdtsrfc(ig)*subtimestep+zdtsrf(ig)*subtimestep
     338         zpicen2(ig)=zpicen2(ig)+pdicen2(ig)*subtimestep
     339       ENDDO
     340     ENDIF  ! (itsub=1)
     341   
     342   DO ig=1,klon
     343     ! forecast of frost temperature ztcondsol
     344     !ztcondsol(ig) = tcond_n2(zplev(ig),zqn2(ig,1))
     345     ztcondsol(ig) = tcond_n2(zplev(ig),vmrn2(ig))
     346
     347!     Loop over where we have condensation / sublimation
     348     IF ((ztsrf(ig) .LT. ztcondsol(ig)) .OR.           &    ! ground cond
     349                ((ztsrf(ig) .GT. ztcondsol(ig)) .AND.  &    ! ground sublim
     350               (zpicen2(ig) .GT. 0.))) THEN
     351        condsub(ig) = .true.    ! condensation or sublimation
     352
     353!     Condensation or partial sublimation of N2 ice
     354        if (ztsrf(ig) .LT. ztcondsol(ig)) then  ! condensation
     355!             Include a correction to account for the cooling of air near
     356!             the surface before condensing:
     357         zcondices(ig)=pcapcal(ig)*(ztcondsol(ig)-ztsrf(ig)) &
     358         /((lw_n2+cpp*(zt(ig,1)-ztcondsol(ig)))*subtimestep)
     359        else                                    ! sublimation
     360          zcondices(ig)=pcapcal(ig)*(ztcondsol(ig)-ztsrf(ig)) &
     361          /(lw_n2*subtimestep)
     362        end if
     363
     364        pdtsrfc(ig) = (ztcondsol(ig) - ztsrf(ig))/subtimestep
     365
     366!     partial sublimation of N2 ice
     367!     If the entire N_2 ice layer sublimes
     368!     (including what has just condensed in the atmosphere)
     369        IF((zpicen2(ig)/subtimestep).LE. &
     370            -zcondices(ig))THEN
     371           zcondices(ig) = -zpicen2(ig)/subtimestep
     372           pdtsrfc(ig)=(lw_n2/pcapcal(ig))*       &
     373               (zcondices(ig))
     374        END IF
     375
     376!     Changing N2 ice amount and pressure
     377
     378        pdicen2(ig) = zcondices(ig)
     379        pdpsrf(ig)   = -pdicen2(ig)*g
     380    !    pdpsrf(ig)   = 0. ! OPTION to check impact N2 sub/cond
     381        IF (zplev(ig)+pdpsrf(ig)*subtimestep.le.0.0000001) then
     382            pdpsrf(ig)=(0.0000001*kp(ig)/p00-zplev(ig))/subtimestep
     383            pdicen2(ig)=-pdpsrf(ig)/g
     384        ENDIF
     385
     386     ELSE  ! no condsub
     387        pdpsrf(ig)=0.
     388        pdicen2(ig)=0.
     389        pdtsrfc(ig)=0.
     390     ENDIF
     391   ENDDO ! ig
     392  enddo ! subtimestep
     393
     394!     Updating pressure, temperature and ice reservoir
     395  DO ig=1,klon
     396     pdpsrf(ig)=(zplev(ig)+pdpsrf(ig)*subtimestep-pplev(ig,1))/ptimestep
     397     ! Two options here : 1 ok, 2 is wrong
     398     pdicen2(ig)=(zpicen2(ig)+pdicen2(ig)*subtimestep-picen2(ig))/ptimestep
     399     !pdicen2(ig)=-pdpsrf(ig)/g
     400
     401     pdtsrfc(ig)=((ztsrf(ig)+pdtsrfc(ig)*subtimestep)-(ztsrfhist(ig)))/ptimestep
     402
     403! security
     404     if (picen2(ig) + pdicen2(ig)*ptimestep.lt.0.) then
     405       write(*,*) 'WARNING in condense_n2:'
     406       write(*,*) picen2(ig),pdicen2(ig)*ptimestep
     407       pdicen2(ig)= -picen2(ig)/ptimestep
     408       pdpsrf(ig)=-pdicen2(ig)*g
     409     endif
     410
     411     if(.not.picen2(ig).ge.0.) THEN
     412!           if(picen2(ig) + pdicen2(ig)*ptimestep.le.-1.e-8) then
     413          print*, 'WARNING NEG RESERVOIR in condense_n2: picen2(',ig,')=', picen2(ig) + pdicen2(ig)*ptimestep 
     414!              pdicen2(ig)= -picen2(ig)/ptimestep
     415!           else
     416          picen2(ig)=0.0
     417!           endif
     418     endif
     419  ENDDO
     420
     421! ***************************************************************
     422!  Correction to account for redistribution between sigma or hybrid
     423!  layers when changing surface pressure (and warming/cooling
     424!  of the n2 currently changing phase).
     425! *************************************************************
     426  if (.not.fast) then
     427   DO ig=1,klon
     428    if (condsub(ig)) then
     429
     430!  Mass fluxes through the sigma levels (kg.m-2.s-1)  (>0 when up)
     431!  """"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
     432        zmflux(1) = -zcondices(ig)
     433        DO l=1,nlayer
     434         zmflux(l+1) = zmflux(l) -zcondicea(ig,l)  &
     435         + (bp(l)-bp(l+1))*(zfallice(ig,1)-zmflux(1))     
     436! zmflux set to 0 if very low to avoid: top layer is disappearing in v1ld 
     437      if (abs(zmflux(l+1)).lt.1E-13.OR.bp(l+1).eq.0.) zmflux(l+1)=0.
     438        END DO
     439
     440! Mass of each layer
     441! ------------------
     442       DO l=1,nlayer
     443         masse(l)=(pplev(ig,l) - pplev(ig,l+1))/g
     444       END DO
     445
     446
     447!  Corresponding fluxes for T,U,V,Q
     448!  """"""""""""""""""""""""""""""""
     449!           averaging operator for TRANSPORT 
     450!           """"""""""""""""""""""""""""""""
     451
     452!           Subtimestep loop to perform the redistribution gently and simultaneously with
     453!           the other tendencies 
     454!           Estimation of subtimestep (using only  the first layer, the most critical)
     455        dtmax=ptimestep
     456        if (zmflux(1).gt.1.e-20) then
     457            dtmax=min(dtmax,masse(1)*zqn2(ig,1)/abs(zmflux(1)))
     458        endif
     459        nsubtimestep= max(nint(ptimestep/dtmax),nint(2.))
     460        subtimestep=ptimestep/float(nsubtimestep)
     461
     462!          New flux for each subtimestep
     463       do l=1,nlayer+1
     464            w(l)=-zmflux(l)*subtimestep
     465       enddo
     466!          initializing variables that will vary during subtimestep:
     467       do l=1,nlayer
     468           ztc(l)  =pt(ig,l)   
     469           zu(l)   =pu(ig,l)
     470           zv(l)   =pv(ig,l) 
     471           do iq=1,nqmx
     472              zq(l,iq) = pq(ig,l,iq)
     473           enddo
     474       end do
     475
     476!          loop over nsubtimestep
     477!          """"""""""""""""""""""
     478       do itsub=1,nsubtimestep
     479!             Progressively adding tendancies from other processes.
     480          do l=1,nlayer
     481             ztc(l)  =ztc(l)   +(pdt(ig,l) + zdtlatent(ig,l))*subtimestep
     482             zu(l)   =zu(l)   +pdu( ig,l) * subtimestep
     483             zv(l)   =zv(l)   +pdv( ig,l) * subtimestep
     484             do iq=1,nqmx
     485                zq(l,iq) = zq(l,iq) + pdq(ig,l,iq)* subtimestep
     486             enddo
     487           end do
     488
     489!             Change of mass in each layer
     490          do l=1,nlayer
     491             masse(l)=masse(l)+pdpsrf(ig)*subtimestep*(pplev(ig,l) - pplev(ig,l+1))&
     492                                             /(g*pplev(ig,1))
     493          end do
     494
     495           ztm(2:nlayermx+1)=0.
     496           zum(2:nlayermx+1)=0.
     497           zvm(2:nlayermx+1)=0.
     498           zqm1(1:nlayermx+1)=0.
     499
     500!             Van Leer scheme:
     501          call vl1d(nlayer,ztc,2.,masse,w,ztm)
     502          call vl1d(nlayer,zu ,2.,masse,w,zum)
     503          call vl1d(nlayer,zv ,2.,masse,w,zvm)
     504         do iq=1,nqmx
     505           do l=1,nlayer
     506            zq1(l)=zq(l,iq)
     507           enddo
     508           zqm1(1)=zqm(1,iq)
     509           call vl1d(nlayer,zq1,2.,masse,w,zqm1)
     510           do l=2,nlayer
     511            zqm(l,iq)=zqm1(l)
     512           enddo
     513          enddo
     514
     515!             Correction to prevent negative value for qn2
     516          if (igcm_n2.ne.0) then
     517            zqm(1,igcm_n2)=1.
     518            do l=1,nlayer-1
     519              if (w(l)*zqm(l,igcm_n2).gt.zq(l,igcm_n2)*masse(l)) then
     520                 zqm(l+1,igcm_n2)=max(zqm(l+1,igcm_n2),    &
     521                 (zqm(l,igcm_n2)*w(l) -zq(l,igcm_n2)*masse(l))/w(l+1) )
     522              else
     523                exit
     524              endif
     525            end do
     526           end if
     527
     528!             Value transfert at the surface interface when condensation sublimation:
     529
     530          if (zmflux(1).lt.0) then
     531!               Surface condensation
     532            zum(1)= zu(1)
     533            zvm(1)= zv(1)
     534            ztm(1) = ztc(1)
     535          else
     536!               Surface sublimation:
     537            ztm(1) = ztsrf(ig) + pdtsrfc(ig)*ptimestep
     538            zum(1) = 0 
     539            zvm(1) = 0 
     540          end if
     541          do iq=1,nqmx
     542             zqm(1,iq)=0. ! most tracer do not condense !
     543          enddo
     544!             Special case if the tracer is n2 gas
     545          if (igcm_n2.ne.0) zqm(1,igcm_n2)=1.
     546
     547          !!! Source haze: 0.02 pourcent when n2 sublimes
     548          IF (source_haze) THEN
     549           IF (pdicen2(ig).lt.0) THEN
     550            DO iq=1,nq
     551             tname=noms(iq)
     552             if (tname(1:4).eq."haze") then
     553                   !zqm(1,iq)=0.02
     554                   !zqm(1,iq)=-pdicen2(ig)*0.02
     555                   zqm(1,iq)=-pdicen2(ig)*ptimestep*0.02
     556                   !zqm(10,iq)=-pdicen2(ig)*ptimestep*100.
     557                   !zqm(1,iq)=-pdicen2(ig)*1000000.
     558
     559             endif
     560            ENDDO
     561           ENDIF     
     562          ENDIF     
     563          ztm(nlayer+1)= ztc(nlayer) ! should not be used, but...
     564          zum(nlayer+1)= zu(nlayer)  ! should not be used, but...
     565          zvm(nlayer+1)= zv(nlayer)  ! should not be used, but...
     566          do iq=1,nqmx
     567           zqm(nlayer+1,iq)= zq(nlayer,iq)
     568          enddo
     569
     570!             Tendencies on T, U, V, Q
     571!             """""""""""""""""""""""
     572          DO l=1,nlayer
     573
     574!               Tendencies on T
     575            zdtsig(ig,l) = (1/masse(l)) *   &
     576           ( zmflux(l)*(ztm(l) - ztc(l))    &
     577           - zmflux(l+1)*(ztm(l+1) - ztc(l)) &
     578           + zcondicea(ig,l)*(ztcond(ig,l)-ztc(l))  )
     579
     580!               Tendencies on U
     581              pduc(ig,l)   = (1/masse(l)) *   &
     582           ( zmflux(l)*(zum(l) - zu(l))&
     583           - zmflux(l+1)*(zum(l+1) - zu(l)) )
     584
     585!               Tendencies on V
     586              pdvc(ig,l)   = (1/masse(l)) *   &
     587           ( zmflux(l)*(zvm(l) - zv(l))   &
     588           - zmflux(l+1)*(zvm(l+1) - zv(l)) )
     589
     590          END DO
     591
     592!             Tendencies on Q
     593          do iq=1,nqmx
     594            if (iq.eq.igcm_n2) then
     595!                 SPECIAL Case when the tracer IS N2 :
     596              DO l=1,nlayer
     597              pdqc(ig,l,iq)= (1/masse(l)) *        &
     598               ( zmflux(l)*(zqm(l,iq) - zq(l,iq))    &
     599               - zmflux(l+1)*(zqm(l+1,iq) - zq(l,iq))&
     600               + zcondicea(ig,l)*(zq(l,iq)-1.) )
     601              END DO
     602            else
     603              DO l=1,nlayer
     604                pdqc(ig,l,iq)= (1/masse(l)) *         &
     605               ( zmflux(l)*(zqm(l,iq) - zq(l,iq))     &
     606               - zmflux(l+1)*(zqm(l+1,iq) - zq(l,iq)) &
     607               + zcondicea(ig,l)*zq(l,iq) )
     608              END DO
     609            end if
     610          enddo
     611!             Update variables at the end of each subtimestep.
     612          do l=1,nlayer
     613            ztc(l)  =ztc(l)  + zdtsig(ig,l)  *subtimestep
     614            zu(l)   =zu(l)   + pduc(ig,l)  *subtimestep
     615            zv(l)   =zv(l)   + pdvc(ig,l)  *subtimestep
     616            do iq=1,nqmx
     617              zq(l,iq) = zq(l,iq) + pdqc(ig,l,iq)*subtimestep
     618            enddo
     619          end do
     620       enddo   ! loop on nsubtimestep
     621!          Recomputing Total tendencies
     622       do l=1,nlayer
     623         pdtc(ig,l)   = (ztc(l) - zt(ig,l) )/ptimestep
     624         pduc(ig,l)   = (zu(l) - (pu(ig,l) + pdu(ig,l)*ptimestep))/ptimestep
     625         pdvc(ig,l)   = (zv(l) - (pv(ig,l) + pdv(ig,l)*ptimestep))/ptimestep
     626         do iq=1,nqmx
     627           pdqc(ig,l,iq) = (zq(l,iq) - (pq(ig,l,iq) + pdq(ig,l,iq)*ptimestep))/ptimestep
     628
     629
     630!           Correction temporaire
     631            if (iq.eq.igcm_n2) then
     632             if((pq(ig,l,iq) +(pdqc(ig,l,iq)+ pdq(ig,l,iq))*ptimestep) &
     633                     .lt.0.01) then ! if n2 < 1 %  !
     634               pdqc(ig,l,iq)=(0.01-pq(ig,l,iq))/ptimestep-pdq(ig,l,iq) 
     635             end if
     636            end if
     637
     638         enddo
     639       end do
     640! *******************************TEMPORAIRE ******************
     641       if (klon.eq.1) then
     642              write(*,*) 'nsubtimestep=' ,nsubtimestep
     643           write(*,*) 'masse avant' , (pplev(ig,1) - pplev(ig,2))/g
     644              write(*,*) 'masse apres' , masse(1)
     645              write(*,*) 'zmflux*DT, l=1' ,  zmflux(1)*ptimestep
     646              write(*,*) 'zmflux*DT, l=2' ,  zmflux(2)*ptimestep
     647          write(*,*) 'pq, l=1,2,3' ,  pq(1,1,1), pq(1,2,1),pq(1,3,1)
     648          write(*,*) 'zq, l=1,2,3' ,  zq(1,1), zq(2,1),zq(3,1)
     649              write(*,*) 'dq*Dt l=1' ,  pdq(1,1,1)*ptimestep
     650              write(*,*) 'dqcond*Dt l=1' ,  pdqc(1,1,1)*ptimestep
     651       end if
     652
     653! ***********************************************************
     654    end if ! if (condsub)
     655   END DO  ! loop on ig
     656  endif ! not fast
     657
     658! ************ Option Olkin to prevent N2 effect in the south********
     659112   continue
     660  if (olkin) then
     661  DO ig=1,klon
     662     if (lati(ig).lt.0) then
     663       pdtsrfc(ig) = max(0.,pdtsrfc(ig))
     664       pdpsrf(ig) = 0.
     665       pdicen2(ig)  = 0.
     666       do l=1,nlayer
     667         pdtc(ig,l)  =  max(0.,zdtlatent(ig,l))
     668         pduc(ig,l)  = 0.
     669         pdvc(ig,l)  = 0.
     670         do iq=1,nqmx
     671           pdqc(ig,l,iq) = 0
     672         enddo
     673       end do
     674     end if
     675  END DO
     676  end if
     677! *******************************************************************
     678
     679! ***************************************************************
     680! Ecriture des diagnostiques
     681! ***************************************************************
     682
     683!     DO l=1,nlayer
     684!        DO ig=1,klon
     685!         Taux de cond en kg.m-2.pa-1.s-1
     686!          tconda1(ig,l)=zcondicea(ig,l)/(pplev(ig,l)-pplev(ig,l+1))
     687!          Taux de cond en kg.m-3.s-1
     688!          tconda2(ig,l)=tconda1(ig,l)*pplay(ig,l)*g/(r*pt(ig,l))
     689!        END DO
     690!     END DO
     691!     call WRITEDIAGFI(ngridmx,'tconda1',              &
     692!     'Taux de condensation N2 atmospherique /Pa',    &
     693!      'kg.m-2.Pa-1.s-1',3,tconda1)
     694!     call WRITEDIAGFI(ngridmx,'tconda2',              &
     695!     'Taux de condensation N2 atmospherique /m',     &
     696!      'kg.m-3.s-1',3,tconda2)
     697
     698
     699  return
     700  end subroutine condens_n2
     701
     702!-------------------------------------------------------------------------
     703
     704  real function tcond_n2(p,vmr)
     705!   Calculates the condensation temperature for N2 at pressure P and vmr
     706  implicit none
     707  real, intent(in):: p,vmr
     708
     709!       tcond_n2 = (1.)/(0.026315-0.0011877*log(.7143*p*vmr))
     710  IF (p.ge.0.529995) then
     711!     tcond Fray and Schmitt for N2 phase beta (T>35.6 K) FIT TB
     712  ! tcond_n2 = (1.)/(1./63.1470-296.925/(2.5e5*0.98)*log(1./(0.125570*1.e5)*p*vmr))
     713    tcond_n2 = (1.)/(0.01583606505-1.211938776e-3*log(7.963685594e-5*p*vmr))
     714  ELSE
     715!     tcond Fray and Schmitt for N2 phase alpha(T<35.6 K) FIT BT
     716  ! tcond_n2 = (1.)/(1./35.6-296.925/(2.5e5*1.09)*log(1./(0.508059)*p*vmr))
     717    tcond_n2 = (1.)/(1./35.6-1.089633028e-3*log(1.968275338*p*vmr))
     718  ENDIF
     719  return
     720  end function tcond_n2
     721
     722!-------------------------------------------------------------------------
     723
     724  real function pcond_n2(t,vmr)
     725!   Calculates the condensation pressure for N2 at temperature T and vmr
     726  implicit none
     727  real, intent(in):: t,vmr
     728
     729!       tcond_n2 = (1.)/(0.026315-0.0011877*log(.7143*p*vmr))
     730  IF (t.ge.35.6) then
     731!     tcond Fray and Schmitt for N2 phase beta (T>35.6 K) FIT TB
     732  ! pcond_n2 = 0.125570*1.e5/vmr*exp((2.5e5*0.98)/296.925*(1./63.1470-1./t))
     733    pcond_n2 = 0.125570e5/vmr*exp(825.1241896*(1./63.147-1./t))
     734  ELSE
     735!     tcond Fray and Schmitt for N2 phase alpha(T<35.6 K) FIT TB
     736  ! pcond_n2 = 0.508059/vmr*exp((2.5e5*1.09)/296.925*(1./35.6-1./t))
     737    pcond_n2 = 0.508059/vmr*exp(917.7401701*(1./35.6-1./t))
     738  ENDIF
     739  return
     740  end function pcond_n2
     741
     742!-------------------------------------------------------------------------
     743
     744  real function glob_average2d(var)
     745!   Calculates the global average of variable var     
     746  use comgeomfi_h
     747  implicit none
     748
     749!      INTEGER klon
     750  REAL var(ngridmx)
     751  INTEGER ig
     752
     753  glob_average2d = 0.
     754  DO ig=2,ngridmx-1
     755       glob_average2d = glob_average2d + var(ig)*area(ig)
     756  END DO
     757  glob_average2d = glob_average2d + var(1)*area(1)*iim
     758  glob_average2d = glob_average2d + var(ngridmx)*area(ngridmx)*iim
     759  glob_average2d = glob_average2d/(totarea+(area(1)+area(ngridmx))*(iim-1))
     760
     761  end function glob_average2d
     762
     763! *****************************************************************
     764
     765  subroutine vl1d(nlayer,q,pente_max,masse,w,qm)
     766!   
     767!     Operateur de moyenne inter-couche pour calcul de transport type
     768!     Van-Leer " pseudo amont " dans la verticale
     769!    q,w sont des arguments d'entree  pour le s-pg ....
     770!    masse : masse de la couche Dp/g
     771!    w : masse d'atm ``transferee'' a chaque pas de temps (kg.m-2)
     772!    pente_max = 2 conseillee
     773!   --------------------------------------------------------------------
     774  IMPLICIT NONE
     775
     776!   Arguments:
     777!   ----------
     778  integer nlayer
     779  real masse(nlayer),pente_max
     780  REAL q(nlayer),qm(nlayer+1)
     781  REAL w(nlayer+1)
    254782!
    255 
    256 !-----------------------------------------------------------
    257 !     START CONDENSATION/SEDIMENTATION SUB-TIME LOOP
    258 !-----------------------------------------------------------
    259 
    260 
    261 !      Ntime =  20  ! number of sub-timestep
    262 !      subptimestep = ptimestep/float(Ntime)           
    263 
    264       ! Add the tendencies from other physical processes at each subtimstep.
    265 !      DO it=1,Ntime
    266 !         DO l=1,nlayer
    267 !            DO ig=1,ngrid
    268 !               zt(ig,l)   = zt(ig,l)   + pdt(ig,l)   * subptimestep
    269 !               zq(ig,l,i_n2ice) = zq(ig,l,i_n2ice) + pdq(ig,l,i_n2ice) * subptimestep
    270 !            END DO
    271 !         END DO
    272 
    273          ! Gravitational sedimentation starts.
    274            
    275          ! Sedimentation computed from radius computed from q in module radii_mod.
    276 !        call n2_reffrad(ngrid,nlayer,nq,zq,reffrad)
    277 !       
    278 !         DO  ilay=1,nlayer
    279 !            DO ig=1,ngrid
    280 
    281 !               reff = reffrad(ig,ilay)
    282 
    283 !              call stokes                      &
    284 !                   (pplev(ig,ilay),pt(ig,ilay), &
    285 !                    reff,vstokes,rho_n2)
    286 
    287 !               !w(ig,ilay,i_n2ice) = 0.0
    288 !               w(ig,ilay,i_n2ice) = vstokes *  subptimestep * &
    289 !                   pplev(ig,ilay)/(r*pt(ig,ilay))
    290 
    291 !            END DO
    292 !         END DO
    293 
    294          ! Computing q after sedimentation
    295 !         call vlz_fi(ngrid,nlayer,zq(1,1,i_n2ice),2.,masse,w(1,1,i_n2ice),wq)
    296 
    297 
    298          ! Progressively accumulating the flux to the ground.
    299          ! Mfallice is the total amount of ice fallen to the ground.
    300 !         DO ig=1,ngrid
    301 !            Mfallice(ig) =  Mfallice(ig) + wq(ig,i_n2ice)
    302 !         END DO
    303 
    304 !----------------------------------------------------------
    305 !       Condensation / sublimation in the atmosphere
    306 !----------------------------------------------------------
    307 !     (MODIFICATIONS FOR EARLY MARS: falling heat neglected, condensation of N2 into tracer i_n2ice)
    308 
    309 
    310 !         DO l=nlayer , 1, -1
    311 !            DO ig=1,ngrid
    312 !               pdtc(ig,l)=0.
    313 
    314                ! ztcond-> ztnuc in test beneath to nucleate only when super saturation occurs(JL 2011)
    315 !               IF ((zt(ig,l).LT.ztnuc(ig,l)).or.(zq(ig,l,i_n2ice).gt.1.E-10)) THEN
    316 !                  pdtc(ig,l)   = (ztcond(ig,l) - zt(ig,l))/subptimestep
    317 !                  pdqc(ig,l,i_n2ice) = pdtc(ig,l)*ccond*g
    318 
    319                   ! Case when the ice from above sublimes entirely
    320 !                  IF ((zq(ig,l,i_n2ice).lt.-pdqc(ig,l,i_n2ice)*subptimestep) &
    321 !                      .AND. (zq(ig,l,i_n2ice).gt.0)) THEN
    322 
    323  !                    pdqc(ig,l,i_n2ice) = -zq(ig,l,i_n2ice)/subptimestep
    324  !                    pdtc(ig,l)   =-zq(ig,l,i_n2ice)/(ccond*g*subptimestep)
    325 
    326 !                  END IF
    327 
    328                   ! Temperature and q after condensation
    329 !                  zt(ig,l)   = zt(ig,l)   + pdtc(ig,l)   * subptimestep
    330 !                  zq(ig,l,i_n2ice) = zq(ig,l,i_n2ice) + pdqc(ig,l,i_n2ice) * subptimestep
    331 !               END IF
    332 
    333 !            ENDDO
    334 !         ENDDO
    335          
    336 !      ENDDO! end of subtimestep loop.
    337 
    338       ! Computing global tendencies after the subtimestep.
    339  !     DO l=1,nlayer
    340  !        DO ig=1,ngrid
    341  !           pdtc(ig,l) = &
    342  !             (zt(ig,l) - (pt(ig,l) + pdt(ig,l)*ptimestep))/ptimestep
    343  !           pdqc(ig,l,i_n2ice) = &
    344  !             (zq(ig,l,i_n2ice)-(pq(ig,l,i_n2ice)+pdq(ig,l,i_n2ice)*ptimestep))/ptimestep
    345  !        END DO
    346  !     END DO
    347  !     DO ig=1,ngrid
    348 !         zfallice(ig) = Mfallice(ig)/ptimestep
    349 !      END DO
    350 
    351 
    352 !-----------------------------------------------------------------------
    353 !              Condensation/sublimation on the ground
    354 !-----------------------------------------------------------------------
    355 
    356 
    357       ! Forecast of ground temperature ztsrf and frost temperature ztcondsol.
    358       DO ig=1,ngrid
    359          picen2(ig)=pqsurf(ig,i_n2ice)
    360          ppn2=gfrac(igas_N2)*pplay(ig,1)
    361          call get_tcond_n2(ppn2,ztcondsol(ig))
    362          
    363          ztsrf(ig) = ptsrf(ig)
    364 
    365          if((ztsrf(ig).le.ztcondsol(ig)+2.0).and.(ngrid.eq.1))then
    366             print*,'N2 is condensing on the surface in 1D. This atmosphere is doomed.'
    367             print*,'T_surf = ',ztsrf,'K'
    368             print*,'T_cond = ',ztcondsol,'K'
    369             open(116,file='surf_vals.out')
    370             write(116,*) 0.0, pplev(1,1), 0.0, 0.0
    371             close(116)
    372             call abort
     783!      Local
     784!   ---------
     785!
     786  INTEGER l
     787!
     788  real dzq(nlayer),dzqw(nlayer),adzqw(nlayer),dzqmax
     789  real sigw, Mtot, MQtot
     790  integer m
     791
     792
     793!    On oriente tout dans le sens de la pression
     794!     W > 0 WHEN DOWN !!!!!!!!!!!!!
     795
     796  do l=2,nlayer
     797        dzqw(l)=q(l-1)-q(l)
     798        adzqw(l)=abs(dzqw(l))
     799  enddo
     800
     801  do l=2,nlayer-1
     802        if(dzqw(l)*dzqw(l+1).gt.0.) then
     803            dzq(l)=0.5*(dzqw(l)+dzqw(l+1))
     804        else
     805            dzq(l)=0.
     806        endif
     807        dzqmax=pente_max*min(adzqw(l),adzqw(l+1))
     808        dzq(l)=sign(min(abs(dzq(l)),dzqmax),dzq(l))
     809  enddo
     810
     811     dzq(1)=0.
     812     dzq(nlayer)=0.
     813
     814   do l = 1,nlayer-1
     815
     816!         Regular scheme (transfered mass < layer mass)
     817!         ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     818      if(w(l+1).gt.0. .and. w(l+1).le.masse(l+1)) then
     819         sigw=w(l+1)/masse(l+1)
     820         qm(l+1)=(q(l+1)+0.5*(1.-sigw)*dzq(l+1))
     821      else if(w(l+1).le.0. .and. -w(l+1).le.masse(l)) then
     822         sigw=w(l+1)/masse(l)
     823         qm(l+1)=(q(l)-0.5*(1.+sigw)*dzq(l))
     824
     825!         Extended scheme (transfered mass > layer mass)
     826!         ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     827      else if(w(l+1).gt.0.) then
     828         m=l+1
     829         Mtot = masse(m)
     830         MQtot = masse(m)*q(m)
     831         !if (m.lt.nlayer) then ! because some compilers will have problems
     832         !                      ! evaluating masse(nlayer+1)
     833         do while ((m.lt.nlayer).and.(w(l+1).gt.(Mtot+masse(m+1))))
     834            m=m+1
     835            Mtot = Mtot + masse(m)
     836            MQtot = MQtot + masse(m)*q(m)
     837         !   if (m.eq.nlayer) exit
     838         end do
     839         !endif
     840         if (m.lt.nlayer) then
     841            sigw=(w(l+1)-Mtot)/masse(m+1)
     842            qm(l+1)= (1/w(l+1))*(MQtot + (w(l+1)-Mtot)* &
     843           (q(m+1)+0.5*(1.-sigw)*dzq(m+1)) )
     844         else
     845            w(l+1) = Mtot
     846            qm(l+1) = Mqtot / Mtot
     847            write(*,*) 'top layer is disapearing !'
     848            stop
     849         end if
     850      else      ! if(w(l+1).lt.0)
     851         m = l-1
     852         Mtot = masse(m+1)
     853         MQtot = masse(m+1)*q(m+1)
     854         if (m.gt.0) then ! because some compilers will have problems
     855                          ! evaluating masse(0)
     856          do while ((m.gt.0).and.(-w(l+1).gt.(Mtot+masse(m))))
     857            m=m-1
     858            Mtot = Mtot + masse(m+1)
     859            MQtot = MQtot + masse(m+1)*q(m+1)
     860            if (m.eq.0) exit
     861          end do
    373862         endif
    374 
    375          ztsrf(ig) = ptsrf(ig) + pdtsrf(ig)*ptimestep
    376 
    377       ENDDO
    378      
    379       DO ig=1,ngrid
    380      
    381          IF(ig.GT.ngrid/2+1) THEN
    382             icap=2
    383          ELSE
    384             icap=1
    385          ENDIF
    386 
    387          ! Loop over where we have condensation / sublimation
    388          IF ((ztsrf(ig) .LT. ztcondsol(ig)) .OR.           &        ! ground condensation
    389                     (zfallice(ig).NE.0.) .OR.              &        ! falling snow
    390                     ((ztsrf(ig) .GT. ztcondsol(ig)) .AND.  &        ! ground sublimation
    391                     ((picen2(ig)+zfallice(ig)*ptimestep) .NE. 0.))) THEN
    392 
    393 
    394             ! Condensation or partial sublimation of N2 ice
    395             zcondices(ig)=pcapcal(ig)*(ztcondsol(ig)-ztsrf(ig)) &
    396                 /(latcond*ptimestep)         
    397             pdtsrfc(ig) = (ztcondsol(ig) - ztsrf(ig))/ptimestep
    398 
    399             ! If the entire CO_2 ice layer sublimes
    400             ! (including what has just condensed in the atmosphere)
    401             IF((picen2(ig)/ptimestep+zfallice(ig)).LE. &
    402                 -zcondices(ig))THEN
    403                zcondices(ig) = -picen2(ig)/ptimestep - zfallice(ig)
    404                pdtsrfc(ig)=(latcond/pcapcal(ig))*       &
    405                    (zcondices(ig))
    406             END IF
    407 
    408             ! Changing N2 ice amount and pressure
    409             picen2(ig)  =  picen2(ig)   + pdqsurfc(ig)*ptimestep
    410             pdqsurfc(ig) =  zcondices(ig) + zfallice(ig)
    411             pdpsrfc(ig)  = -pdqsurfc(ig)*g
    412 
    413             IF(ABS(pdpsrfc(ig)*ptimestep).GT.pplev(ig,1)) THEN
    414                PRINT*,'STOP in condens in condense_n2'
    415                PRINT*,'condensing more than total mass'
    416                PRINT*,'Grid point ',ig
    417                PRINT*,'Ps = ',pplev(ig,1)
    418                PRINT*,'d Ps = ',pdpsrfc(ig)
    419                STOP
    420             ENDIF
    421          END IF
    422          
    423       ENDDO ! end of ngrid loop.
    424 
    425 
    426 !---------------------------------------------------------------------------------------------
    427 !     Surface albedo and emissivity of the ground below the snow (emisref)
    428 !---------------------------------------------------------------------------------------------
    429 
    430 
    431       DO ig=1,ngrid
    432      
    433 !         IF(latitude(ig).LT.0.) THEN
    434 !            icap=2 ! Southern Hemisphere
    435 !         ELSE
    436 !            icap=1 ! Nortnern hemisphere
    437 !         ENDIF
    438 
    439          if(.not.picen2(ig).ge.0.) THEN
    440             if(picen2(ig).le.-1.e-8) print*,   &
    441               'WARNING : in condense_n2cloud: picen2(',ig,')=', picen2(ig)
    442             picen2(ig)=0.
    443          endif
    444 !         if (picen2(ig) .gt. 1.) then  ! N2 Albedo condition changed to ~1 mm coverage. Change by MT2015.
    445 !           DO nw=1,L_NSPECTV
    446 !               albedo(ig,nw) = albedo_n2_ice_SPECTV(nw)
    447 !           ENDDO
    448 !            emisref(ig)   = emisice(icap)
    449 !         else
    450 !            DO nw=1,L_NSPECTV
    451 !               albedo(ig,nw) = albedo_bareground(ig) ! Note : If you have some water, it will be taken into account in the "hydrol" routine.
    452 !           ENDDO
    453 !            emisref(ig)   = emissiv
    454 !            pemisurf(ig)  = emissiv
    455 !         end if
    456          
    457       END DO
    458 
    459       return
    460      
    461       end subroutine condense_n2
    462 
    463 
    464 
    465 !-------------------------------------------------------------------------
    466 !-------------------------------------------------------------------------
    467 !-------------------------------------------------------------------------
    468 !-------------------------------------------------------------------------
    469 !-------------------------------------------------------------------------
    470 !-------------------------------------------------------------------------
    471 
    472 
    473 
    474       subroutine get_tcond_n2(p,tcond)  ! Calculates the condensation temperature for N2
    475      
    476 
    477       implicit none
    478 
    479       real p, tcond
    480 
    481       IF (p.ge.0.529995) then
    482      ! tcond Fray and Schmitt for N2 phase beta (T>35.6 K) FIT TB
    483          tcond_n2 = (1.)/(0.01583606505-1.211938776e-3*log(7.963685594e-5*p*vmr))
    484       ELSE
    485      ! tcond Fray and Schmitt for N2 phase alpha(T<35.6 K) FIT BT
    486          tcond_n2 = (1.)/(1./35.6-1.089633028e-3*log(1.968275338*p*vmr))
    487       ENDIF
    488       return
    489 
    490       end subroutine get_tcond_n2
    491 
    492 
    493 
    494 !-------------------------------------------------------------------------
    495 !-------------------------------------------------------------------------
    496 !-------------------------------------------------------------------------
    497 !-------------------------------------------------------------------------
    498 !-------------------------------------------------------------------------
    499 !-------------------------------------------------------------------------
    500 
    501 
    502 
    503       subroutine get_tnuc_n2(p,tnuc)
    504       ! Calculates the nucleation temperature for N2, based on a simple super saturation criterion. JL 2011.
    505 
    506       use callkeys_mod, only: n2supsat
    507 
    508       implicit none
    509 
    510       real p, peff, tnuc
    511       real, parameter :: ptriple=518000.0
    512 
    513       peff=p/n2supsat
    514 
    515       !! AF24: this code below is for CO2, needs to be udpated for N2
    516 
    517       ! if(peff.lt.ptriple) then
    518       !    tnuc = (-3167.8)/(log(.01*peff)-23.23) ! Fanale's formula
    519       ! else
    520       !    tnuc = 684.2-92.3*log(peff)+4.32*log(peff)**2
    521       !    ! liquid-vapour transition (based on CRC handbook 2003 data)
    522       ! endif
    523      
    524       return
    525 
    526       end subroutine get_tnuc_n2
     863         if (m.gt.0) then
     864            sigw=(w(l+1)+Mtot)/masse(m)
     865            qm(l+1)= (-1/w(l+1))*(MQtot + (-w(l+1)-Mtot)* &
     866           (q(m)-0.5*(1.+sigw)*dzq(m))  )
     867         else
     868            qm(l+1)= (-1/w(l+1))*(MQtot + (-w(l+1)-Mtot)*qm(1))
     869         end if   
     870      end if
     871   enddo
     872
     873   return
     874   end  subroutine vl1d
     875

  • trunk/LMDZ.PLUTO/libf/phypluto/datafile_mod.F90

    r3184 r3193  
    2121      character(LEN=18),parameter :: aerdir="aerosol_properties"
    2222
     23      ! Data haze properties
     24      character(len=300),save :: hazeprop
     25     
    2326      end module datafile_mod
    2427!-----------------------------------------------------------------------

  • trunk/LMDZ.PLUTO/libf/phypluto/inifis_mod.F90

    r3184 r3193  
    1313  use radcommon_h, only: ini_radcommon_h
    1414  use radii_mod, only: radfixed, Nmix_n2
    15   use datafile_mod, only: datadir
     15  use datafile_mod, only: datadir, hazeprop
    1616  use comdiurn_h, only: sinlat, coslat, sinlon, coslon
    1717  use comgeomfi_h, only: totarea, totarea_planet
    … …  
    136136     if (is_master) write(*,*) trim(rname)//": datadir = ",trim(datadir)
    137137
     138     if (is_master) write(*,*) "Haze optical properties datafile"
     139     hazeprop="../../../datadir"  ! default path
     140     call getin_p("hazeprop",hazeprop)
     141     if (is_master) write(*,*) trim(rname)//" hazeprop = ",trim(hazeprop)
     142
     143
    138144     if (is_master) write(*,*) trim(rname)//&
    139145       ": Run with or without tracer transport ?"
    … …  
    788794     if (is_master) write(*,*)trim(rname)//": generic_condensation = ",generic_condensation
    789795     
    790      if (is_master) write(*,*)trim(rname)//": Generic Condensation of tracers ?"
    791      generic_rain=.false. !default value
    792      call getin_p("generic_rain",generic_rain)
    793      if (is_master) write(*,*)trim(rname)//": generic_rain = ",generic_rain
    794 
    795      if (is_master) write(*,*)trim(rname)//": Compute water cycle ?"
    796      water=.false. ! default value
    797      call getin_p("water",water)
    798      if (is_master) write(*,*)trim(rname)//": water = ",water
    799          
    800 ! Test of incompatibility:
    801 ! if water is true, there should be at least a tracer
    802      if (water.and.(.not.tracer)) then
    803         call abort_physic(rname,'if water is ON, tracer must be ON too!',1)
    804      endif
    805 
    806      if (is_master) write(*,*)trim(rname)//": Include water condensation ?"
    807      watercond=.false. ! default value
    808      call getin_p("watercond",watercond)
    809      if (is_master) write(*,*)trim(rname)//": watercond = ",watercond
    810 
    811796! Test of incompatibility:
    812797     if (is_master) write(*,*)trim(rname)//": Spectral Dependant albedo ?"

  • trunk/LMDZ.PLUTO/libf/phypluto/initracer.F90

    r3184 r3193  
    520520      rho_n2=1000.   ! N2 ice density (kg.m-3)
    521521
     522      lw_co=274000.
     523      lw_ch4=586700.
     524      lw_n2=2.5e5
     525     
     526      if (haze) then
     527        ! the sedimentation radius remains radius(igcm_haze)
     528        if (fractal) then
     529           nmono=nb_monomer
     530        else
     531           nmono=1
     532        endif 
     533
     534        ia=0
     535        if (aerohaze) then
     536           ia=ia+1
     537           iaero_haze=ia
     538           write(*,*) '--- number of haze aerosol = ', iaero_haze
     539
     540           block=0  ! Only one type of haze is active : the first one set in traceur.def
     541           do iq=1,nqmx
     542             tracername=noms(iq)
     543             write(*,*) "--> tracername ",iq,'/',nqmx,' = ',tracername
     544             if (tracername(1:4).eq."haze".and.block.eq.0) then
     545               i_haze=iq
     546               block=1
     547               write(*,*) "i_haze=",i_haze
     548               write(*,*) "Careful: if you set many haze traceurs in
     549    & traceur.def,only ",tracername," will be radiatively active
     550    & (first one in traceur.def)"
     551             endif
     552           enddo
     553        endif
     554     endif
     555
    522556!     Initialization for water vapor !AF24: removed
    523557

  • trunk/LMDZ.PLUTO/libf/phypluto/suaer_corrk.F90

    r3184 r3193  
    1010      use radinc_h,    only: L_NSPECTI,L_NSPECTV,nsizemax,naerkind
    1111      use radcommon_h, only: blamv,blami,lamrefir,lamrefvis
    12       use datafile_mod, only: datadir, aerdir
     12      use datafile_mod, only: datadir, aerdir, hazeprop
    1313
    1414      ! outputs
    … …  
    136136        print*, 'naerkind= 0'
    137137      endif
     138     
     139     
    138140      do iaer=1,naerkind
    139        if (iaer.eq.iaero_n2) then
    140         forgetit=.true.
    141           if (.not.noaero) then
    142               print*, 'naerkind= n2', iaer
    143           end if
    144 !     visible
    145         if(forgetit)then
    146            file_id(iaer,1) = 'optprop_n2_vis_ff.dat' ! Francois' values
    147         else
    148            file_id(iaer,1) = 'optprop_n2ice_vis_n50.dat'
    149         endif
    150         lamrefvis(iaer)=1.5E-6   ! 1.5um OMEGA/MEx ???
    151 
    152 !     infrared
    153         if(forgetit)then
    154            file_id(iaer,2) = 'optprop_n2_ir_ff.dat' ! Francois' values
    155         else
    156            file_id(iaer,2) = 'optprop_n2ice_ir_n50.dat'
    157         endif
    158         lamrefir(iaer)=12.1E-6   ! Dummy in generic phys. (JVO 20)
    159        endif ! N2 aerosols 
    160 ! NOTE: these lamref's are currently for dust, not N2 ice.
    161 ! JB thinks this shouldn't matter too much, but it needs to be
    162 ! fixed / decided for the final version.
    163 
    164        if (iaer.eq.iaero_dust) then
    165         print*, 'naerkind= dust', iaer
    166 
    167 !     visible
    168          file_id(iaer,1) = 'optprop_dustvis_n50.dat'
    169          !lamrefvis(iaer)=1.5E-6   ! 1.5um OMEGA/MEx
    170          lamrefvis(iaer)=0.67e-6
    171 !     infrared
    172          file_id(iaer,2) = 'optprop_dustir_n50.dat'
    173          lamrefir(iaer)=9.3E-6     ! Dummy in generic phys. (JVO 20)
    174        endif
    175 
    176        if (iaer.eq.iaero_h2so4) then
    177          print*, 'naerkind= h2so4', iaer
    178 
    179 !     visible
    180          file_id(iaer,1) = 'optprop_h2so4vis_n20.dat'
    181          lamrefvis(iaer)=1.5E-6   ! no idea, must find
    182 !     infrared
    183          file_id(iaer,2) = 'optprop_h2so4ir_n20.dat'
    184          lamrefir(iaer)=9.3E-6 ! ! Dummy in generic phys. (JVO 20)
    185 ! added by LK
    186        endif
    187 
    188        if (iaer.eq.iaero_back2lay) then
    189          print*, 'naerkind= back2lay', iaer
     141         !     naerkind=1: haze
     142         if (iaer.eq.1) then
     143
     144         ! Only one table of optical properties :
     145         write(*,*)'Suaer haze optical properties, using: ', &
     146                           TRIM(hazeprop)
     147
     148         ! visible
     149         file_id(iaer,1)=TRIM(hazeprop)
     150         ! infrared
     151         file_id(iaer,2)=file_id(iaer,1)
    190152         
    191 !     visible
    192          file_id(iaer,1) = TRIM(optprop_back2lay_vis)
    193          lamrefvis(iaer)=0.8E-6  ! This is the important one.
    194 !     infrared
    195          file_id(iaer,2) = TRIM(optprop_back2lay_ir)
    196          lamrefir(iaer)=6.E-6    ! This is dummy.
    197 ! added by SG
    198        endif
    199      
    200        if (iaer.eq.iaero_nh3) then
    201          print*, 'naerkind= nh3', iaer
    202 
    203 !     visible
    204          file_id(iaer,1) = 'optprop_nh3ice_vis.dat'
    205          lamrefvis(iaer)=0.756E-6  !
    206 !     infrared
    207          file_id(iaer,2) = 'optprop_nh3ice_ir.dat'
    208          lamrefir(iaer)=6.E-6  ! dummy
    209 ! added by SG
    210        endif
    211 
    212        do ia=1,nlayaero
    213          if (iaer.eq.iaero_nlay(ia)) then
    214            print*, 'naerkind= nlay', iaer
    215            
    216 !       visible
    217            file_id(iaer,1) = TRIM(optprop_aeronlay_vis(ia))
    218            lamrefvis(iaer)=aeronlay_lamref(ia)
    219 !       infrared
    220            file_id(iaer,2) = TRIM(optprop_aeronlay_ir(ia))
    221            lamrefir(iaer)=6.E-6 ! Dummy value
     153         lamrefvis(iaer)=0.607E-6   ! reference wavelength for opacity vis pivot wavelength Cheng et al 2017
     154         lamrefir(iaer)=2.E-6   !  reference wavelength for opacity IR (in the LEISA range)
     155
    222156         endif
    223        enddo
    224 ! added by JVO
    225      
    226        if (iaer.eq.iaero_aurora) then
    227          print*, 'naerkind= aurora', iaer
    228 
    229 !     visible
    230          file_id(iaer,1) = 'optprop_aurora_vis.dat'
    231          lamrefvis(iaer)=0.3E-6  !
    232 !     infrared
    233          file_id(iaer,2) = 'optprop_aurora_ir.dat'
    234          lamrefir(iaer)=6.E-6  ! dummy
    235 ! added by SG
    236        endif
    237 
    238        
    239 ! the following was added by LT
    240        do ia=1,aerogeneric ! Read Radiative Generic Condensable Species data
    241          if (iaer .eq. iaero_generic(ia)) then
    242             if (index(noms(i_rgcs_ice(ia)),'Fe') .ne. 0) then
    243                print*,"Reading Fe file"
    244                file_id(iaer,1)='Fe.ocst.txt'
    245                file_id(iaer,2)='Fe.ocst.txt'
    246                lamrefvis(iaer) = 1.0E-6 !random pick
    247                lamrefir(iaer) = 1.0E-6 !dummy but random pick
    248             else if (index(noms(i_rgcs_ice(ia)),'Mn') .ne. 0) then
    249                print*,"Reading MnS file"
    250                file_id(iaer,1)='MnS.ocst.txt'
    251                file_id(iaer,2)='MnS.ocst.txt'
    252                lamrefvis(iaer) = 1.0E-6 !random pick
    253                lamrefir(iaer) = 1.0E-6 !dummy but random pick   
    254             else if (index(noms(i_rgcs_ice(ia)),'Mg') .ne. 0) then 
    255                print*,"Reading Mg2SiO4 file"
    256                file_id(iaer,1)='Mg2SiO4.ocst.txt'
    257                file_id(iaer,2)='Mg2SiO4.ocst.txt'
    258                lamrefvis(iaer) = 1.0E-6 !random pick
    259                lamrefir(iaer) = 1.0E-6 !dummy but random pick 
    260             else if (index(noms(i_rgcs_ice(ia)),'Cr') .ne. 0) then
    261                print*,"Reading Cr file"
    262                file_id(iaer,1)='Cr.ocst.txt'
    263                file_id(iaer,2)='Cr.ocst.txt'
    264                lamrefvis(iaer) = 1.0E-6 !random pick
    265                lamrefir(iaer) = 1.0E-6 !dummy but random pick
    266             else
    267 ! If you want to add another specie, copy,paste & adapt the elseif block up here to your new specie (LT 2022)
    268                call abort_physic("suaer_corrk", "Unknown specie in radiative generic condensable species",1)
    269             endif
    270          endif
    271        enddo ! ia=1,aerogeneric
    272       enddo ! of do iaer=1,naerkind
     157      enddo   
     158
     159      IF (naerkind .gt. 1) THEN
     160         print*,'naerkind = ',naerkind
     161         print*,'but we only have data for 1 type, exiting.'
     162         call abort
     163      ENDIF
    273164     
    274165!------------------------------------------------------------------
Note: See TracChangeset for help on using the changeset viewer.