Changeset 3195 for trunk/LMDZ.PLUTO

Timestamp:

Jan 31, 2024, 4:36:51 PM (10 months ago)

Author:

afalco

Message:

Pluto PCM:
Imported condense n2 from pluto.old.
Aerosol data from Pluto.old not yet working.
AF

Location:

trunk/LMDZ.PLUTO

Files:

• changelog.txt (modified) (1 diff)
• deftank/callphys.def (modified) (4 diffs)
• deftank/traceur.def (modified) (1 diff)
• libf/dynphy_lonlat/phypluto/lect_start_archive.F (modified) (44 diffs)
• libf/phypluto/aeropacity.F90 (modified) (6 diffs)
• libf/phypluto/aerosol_mod.F90 (modified) (2 diffs)
• libf/phypluto/callcorrk.F90 (modified) (40 diffs)
• libf/phypluto/callkeys_mod.F90 (modified) (4 diffs)
• libf/phypluto/callsedim.F (modified) (7 diffs)
• libf/phypluto/condense_n2.F90 (modified) (49 diffs)
• libf/phypluto/convadj.F (modified) (15 diffs)
• libf/phypluto/inifis_mod.F90 (modified) (4 diffs)
• libf/phypluto/initracer.F90 (modified) (23 diffs)
• libf/phypluto/phys_state_var_mod.F90 (modified) (4 diffs)
• libf/phypluto/physiq_mod.F90 (modified) (93 diffs)
• libf/phypluto/radii_mod.F90 (modified) (5 diffs)
• libf/phypluto/suaer_corrk.F90 (modified) (17 diffs)
• libf/phypluto/surfdat_h.F90 (modified) (1 diff)
• libf/phypluto/surfini.F (modified) (3 diffs)
• libf/phypluto/tracer_h.F90 (modified) (3 diffs)
• libf/phypluto/vdifc_mod.F (modified) (22 diffs)

trunk/LMDZ.PLUTO/changelog.txt

@@ -1838 +1838 @@
 == 25/01/2024 == AF
 Added LMDZ.PLUTO, a copy of the generic model cleaned of some 
-unnecessary  
+unnecessary modules

trunk/LMDZ.PLUTO/deftank/callphys.def

@@ -6 +6 @@
 diurnal   = .true.
 # Seasonal cycle ? if season=false, Ls stays constant, to value set in "start"
-season    = .true. 
+season    = .true.
 # Tidally resonant orbit ? must have diurnal=false, correct rotation rate in newstart
 tlocked   = .false.
@@ -82 +82 @@
 Fat1AU = 1366.0
 
-## Tracer and aerosol options 
+## Tracer and aerosol options
 ## ~~~~~~~~~~~~~~~~~~~~~~~~~~
 # Gravitational sedimentation of tracers (just H2O ice for now) ?
 sedimentation = .true.
 
-naerkind = 1
+naerkind = 0
 hazeprop = optprop_tholins_fractal100nm.dat
 
@@ -105 +105 @@
 ## extra non-standard definitions for Earth
 #########################################################################
- 
-## Tracer and aerosol options 
+
+## Tracer and aerosol options
 ## ~~~~~~~~~~~~~~~~~~~~~~~~~~
 # atm mass update due to tracer evaporation/condensation?
@@ -123 +123 @@
 
 
-## Condensation 
+## Condensation
 ## ~~~~~~~~~~~
 # call N2 condensation ?

trunk/LMDZ.PLUTO/deftank/traceur.def

@@ -1 @@
-1
-n2
+0

trunk/LMDZ.PLUTO/libf/dynphy_lonlat/phypluto/lect_start_archive.F

@@ -7 +7 @@
 
       USE comsoil_h, ONLY: nsoilmx, layer, mlayer, volcapa, inertiedat
-      USE tracer_h, ONLY: igcm_n2_ice
+      USE tracer_h, ONLY: igcm_haze
       USE infotrac, ONLY: tname, nqtot
 !      USE slab_ice_h, ONLY: noceanmx
@@ -39 +39 @@
 
 
-c Variables pour les lectures des fichiers "ini" 
+c Variables pour les lectures des fichiers "ini"
 c--------------------------------------------------
 !      INTEGER sizei,
@@ -53 +53 @@
 !      character (len=50) :: tmpname
 
-c Variable histoire 
+c Variable histoire
 c------------------
       REAL,INTENT(OUT) :: vcov(iip1,jjm,llm),ucov(iip1,jjp1,llm) ! vents covariants
@@ -59 +59 @@
       REAL,INTENT(OUT) :: q(iip1,jjp1,llm,nqtot)
 
-c Physique sur grille scalaire 
+c Physique sur grille scalaire
 c----------------------------
 
@@ -112 +112 @@
 
 
-c Variable de l'ancienne grille 
+c Variable de l'ancienne grille
 c---------------------------------------------------------
 
@@ -150 +150 @@
       logical :: therminertia_3D=.true. ! flag
 ! therminertia_3D=.true. if thermal inertia is 3D and read from datafile
-c Variable intermediaires iutilise pour l'extrapolation verticale 
+c Variable intermediaires iutilise pour l'extrapolation verticale
 c----------------------------------------------------------------
-      real, dimension(:,:,:), allocatable :: var,varp1 
+      real, dimension(:,:,:), allocatable :: var,varp1
       real, dimension(:), allocatable :: oldgrid, oldval
       real, dimension(:), allocatable :: newval
@@ -158 +158 @@
       real,intent(out) :: surfith(iip1,jjp1) ! surface thermal inertia
       ! surface thermal inertia at old horizontal grid resolution
-      real, dimension(:,:), allocatable :: surfithold 
+      real, dimension(:,:), allocatable :: surfithold
 
       character(len=30) :: txt ! to store some text
@@ -172 +172 @@
 !-----------------------------------------------------------------------
 
-! 1.2 Read the various dimension lengths of data in file 
+! 1.2 Read the various dimension lengths of data in file
 
       ierr= NF_INQ_DIMID(nid,"Time",dimid)
@@ -245 +245 @@
 
 ! 1.3 Report dimensions
-      
+
       write(*,*) "Start_archive dimensions:"
       write(*,*) "longitude: ",imold
@@ -259 +259 @@
       endif
       write(*,*) "time lenght: ",timelen
-      write(*,*) 
+      write(*,*)
 
 !-----------------------------------------------------------------------
@@ -311 +311 @@
 
 C-----------------------------------------------------------------------
-c 3.1. Lecture du tableau des parametres du run 
+c 3.1. Lecture du tableau des parametres du run
 c     (pour  la lecture ulterieure de "ptotalold" et "n2icetotalold")
 c-----------------------------------------------------------------------
@@ -355 +355 @@
          PRINT*, "lect_start_archive: Field <rlatu> not found"
          CALL abort
-      ENDIF 
+      ENDIF
 #ifdef NC_DOUBLE
       ierr = NF_GET_VAR_DOUBLE(nid, nvarid, rlatuold)
@@ -440 +440 @@
 c 3.4 Read Soil layers depths
 c-----------------------------------------------------------------------
-     
+
       ierr=NF_INQ_VARID(nid,"soildepth",nvarid)
       if (ierr.ne.NF_NOERR) then
@@ -540 +540 @@
       ptotalold = tab_cntrl(tab0+49)
       n2icetotalold = tab_cntrl(tab0+50)
- 
+
 c-----------------------------------------------------------------------
 c 4. Lecture du temps et choix
 c-----------------------------------------------------------------------
- 
+
 c  lecture du temps
 c
@@ -559 +559 @@
       IF (ierr .NE. NF_NOERR) THEN
          ierr = NF_INQ_VARID (nid, "temps", nvarid)
-      endif 
+      endif
 #ifdef NC_DOUBLE
       ierr = NF_GET_VAR_DOUBLE(nid, nvarid, timelist)
@@ -584 +584 @@
 
    6  FORMAT(i7,i7,f9.3)
- 
+
       write(*,*)
       write(*,*) 'Choice for the date'
@@ -595 +595 @@
         endif
       end do
- 
+
       if (memo.eq.0) then
         write(*,*)
@@ -618 +618 @@
 c 5.1 Lecture des champs 2D (n2ice, emis,ps,tsurf,Tg[10], qsurf)
 c-----------------------------------------------------------------------
- 
+
       start=(/1,1,memo,0/)
       count=(/imold+1,jmold+1,1,0/)
-       
+
       ierr = NF_INQ_VARID (nid, "emis", nvarid)
       IF (ierr .NE. NF_NOERR) THEN
@@ -756 +756 @@
 !          PRINT*, "lect_start_archive: Failed loading <sea_ice>"
 !       ENDIF
-  
+
 !       ENDIF! ok_slab_ocean
 c
@@ -763 +763 @@
       write(*,*)
 
-! Surface tracers:      
+! Surface tracers:
       ! initialize all surface tracers to zero
       qsurfold(1:imold+1,1:jmold+1,1:nqtot)=0
@@ -777 +777 @@
           endif
 
-        
+
         write(*,*) "lect_start_archive: loading tracer ",trim(txt)
         ierr = NF_INQ_VARID (nid,txt,nvarid)
@@ -838 +838 @@
 c
        enddo ! of do isoil=1,nsoilold
-      
+
       ! reset 'start' and 'count' to "3D" behaviour
       start=(/1,1,1,memo/)
       count=(/imold+1,jmold+1,nsoilold,1/)
-      
+
       else
        write(*,*) "lect_start_archive: loading tsoil "
@@ -856 +856 @@
 #endif
        endif ! of if (ierr.ne.NF_NOERR)
-       
+
       endif ! of if (olddepthdef)
 
@@ -928 +928 @@
 c
 
-! Tracers:      
+! Tracers:
       qold(1:imold+1,1:jmold+1,1:lmold,1:nqtot)=0.
 
@@ -1040 +1040 @@
 c   INTERPOLATION DANS LA NOUVELLE GRILLE et initialisation des variables
 c=======================================================================
-c  Interpolation horizontale puis passage dans la grille physique pour 
-c  les variables physique 
+c  Interpolation horizontale puis passage dans la grille physique pour
+c  les variables physique
 c  Interpolation verticale puis horizontale pour chaque variable 3D
 c=======================================================================
@@ -1048 +1048 @@
 c 6.1   Variable 2d :
 c-----------------------------------------------------------------------
-c Relief 
+c Relief
       call interp_horiz (phisold,phisold_newgrid,imold,jmold,iim,jjm,1,
      &                   rlonuold,rlatvold,rlonu,rlatv)
 
-! N2 ice is now in qsurf(igcm_n2_ice)
+! N2 ice is now in qsurf(igcm_haze)
 !      call interp_horiz (n2iceold,n2ices,imold,jmold,iim,jjm,1,
 !     &                   rlonuold,rlatvold,rlonu,rlatv)
@@ -1088 +1088 @@
 c       On assure la conservation de la masse de l'atmosphere + calottes
 c-----------------------------------------------------------------------
+      !AF: TODO: mass conservation: check this. haze?
 
       ptotal =  0.
@@ -1096 +1097 @@
       END DO
       n2icetotal = 0.
-      if (igcm_n2_ice.ne.0) then
+      if (igcm_haze.ne.0) then
         ! recast surface N2 ice on new grid
-        call interp_horiz(qsurfold(1,1,igcm_n2_ice),
-     &                  qsurfs(1,1,igcm_n2_ice),
+        call interp_horiz(qsurfold(1,1,igcm_haze),
+     &                  qsurfs(1,1,igcm_haze),
      &                  imold,jmold,iim,jjm,1,
      &                  rlonuold,rlatvold,rlonu,rlatv)
@@ -1105 +1106 @@
          DO i=1,iim
            !n2icetotal = n2icetotal + n2iceS(i,j)*aire(i,j)
-           n2icetotal=n2icetotal+qsurfS(i,j,igcm_n2_ice)*aire(i,j)
+           n2icetotal=n2icetotal+qsurfS(i,j,igcm_haze)*aire(i,j)
          END DO
        END DO
@@ -1115 +1116 @@
       write(*,*)'Old grid: atmospheric mass :',ptotalold
       write(*,*)'New grid: atmospheric mass :',ptotal
-      write (*,*) 'Ratio new atm./ old atm =', ptotal/ptotalold 
+      write (*,*) 'Ratio new atm./ old atm =', ptotal/ptotalold
       write(*,*)
       write(*,*)'Old grid: mass of N2 ice:',n2icetotalold
@@ -1131 +1132 @@
       END DO
 
-      if ( n2icetotalold.gt.0.) then 
+      if ( n2icetotalold.gt.0.) then
 !         DO j=1,jjp1
 !            DO i=1,iip1
@@ -1157 +1158 @@
      &f soil thermal inertia; might be wiser to reset it.'
         write(*,*)
-       
+
         do i=1,imold+1
          do j=1,jmold+1
@@ -1177 +1178 @@
         ! We have inertiedatold
        if((imold.ne.iim).or.(jmold.ne.jjm)) then
-       write(*,*)'lect_start_archive: WARNING: horizontal interpolation 
+       write(*,*)'lect_start_archive: WARNING: horizontal interpolation
      &of thermal inertia; might be better to reset it.'
        write(*,*)
        endif
-       
+
         ! Do horizontal interpolation
         if (depthinterpol) then
@@ -1210 +1211 @@
       call gr_dyn_fi (nsoilmx,iim+1,jjm+1,ngrid,
      &                  inertiedatS,inertiedat)
-      
+
 c-----------------------------------------------------------------------
 c 6.2.2 Soil temperature
@@ -1282 +1283 @@
         deallocate(oldval)
         deallocate(newval)
-       
+
        else
         tsoiloldnew(:,:,:)=tsoilold(:,:,:)
@@ -1303 +1304 @@
 c 6.4 Variable 3d :
 c-----------------------------------------------------------------------
-      
+
 c temperatures atmospheriques
       write (*,*) 'lect_start_archive: told ', told (1,jmold+1,1)  ! INFO
@@ -1321 +1322 @@
      &                   rlonuold,rlatvold,rlonu,rlatv)
       call gr_dyn_fi(llm,iim+1,jjm+1,ngrid,du_nonoro_gwdS,du_nonoro_gwd)
-      
+
       call interp_vert
      &    (dv_nonoro_gwdold,var,lmold,llm,apsold,bpsold,aps,bps,
@@ -1328 +1329 @@
      &                   rlonuold,rlatvold,rlonu,rlatv)
       call gr_dyn_fi(llm,iim+1,jjm+1,ngrid,dv_nonoro_gwdS,dv_nonoro_gwd)
-      
+
       call interp_vert
      &    (east_gwstressold,var,lmold,llm,apsold,bpsold,aps,bps,
@@ -1335 +1336 @@
      &                   rlonuold,rlatvold,rlonu,rlatv)
       call gr_dyn_fi(llm,iim+1,jjm+1,ngrid,east_gwstressS,east_gwstress)
-      
+
       call interp_vert
      &    (west_gwstressold,var,lmold,llm,apsold,bpsold,aps,bps,
@@ -1379 +1380 @@
           end do
         end do
-      end do 
+      end do
       write (*,*) 'lect_start_archive: ucov ', ucov (1,2,1)  ! INFO
 c     write(48,*) 'ucov',ucov
@@ -1410 +1411 @@
      &                  rlonuold,rlatvold,rlonu,rlatv)
       enddo
-cccccccccccccccccccccccccccccc      
-c  make sure that sum of q = 1      
-c dominent species is = 1 - sum(all other species)      
-cccccccccccccccccccccccccccccc      
+cccccccccccccccccccccccccccccc
+c  make sure that sum of q = 1
+c dominent species is = 1 - sum(all other species)
+cccccccccccccccccccccccccccccc
 c      iqmax=1
-c      
+c
 c      if (nqold.gt.10) then
 c       do l=1,llm
@@ -1428 +1429 @@
 c           qtot(i,j,l)=0
 c           do iq=1,nqold
-c            if (iq.ne.iqmax) then        
-c              q(i,j,l,iqmax)=q(i,j,l,iqmax)-q(i,j,l,iq)        
+c            if (iq.ne.iqmax) then
+c              q(i,j,l,iqmax)=q(i,j,l,iqmax)-q(i,j,l,iq)
 c            endif
 c           enddo !iq
@@ -1437 +1438 @@
 c     $    qtot(i,j,l)
 c           enddo !iq
-c          enddo !i   
-c         enddo !j   
-c       enddo !l  
+c          enddo !i
+c         enddo !j
+c       enddo !l
 c      endif
 ccccccccccccccccccccccccccccccc
@@ -1451 +1452 @@
          end do
       enddo
-      
+
 !      call gr_dyn_fi (1,iim+1,jjm+1,ngrid,n2ices,n2ice)
-! no need to transfer "n2ice" any more; it is in qsurf(igcm_n2_ice)
+! no need to transfer "n2ice" any more; it is in qsurf(igcm_haze)
 
       endif !! if nqtot .ne. 0

trunk/LMDZ.PLUTO/libf/phypluto/aeropacity.F90

@@ -10 +10 @@
 
        use radinc_h, only : L_TAUMAX,naerkind
-       use aerosol_mod, only: iaero_nlay, iaero_generic, &
-                              iaero_aurora, iaero_back2lay, iaero_n2, &
-                              iaero_dust, iaero_h2so4, &
-                              iaero_nh3, i_rgcs_ice, noaero
+       use aerosol_mod, only: iaero_haze, i_haze, iaero_generic, &
+                              noaero
        USE tracer_h, only: noms,rho_n2,rho_ice,rho_q,mmol
        use comcstfi_mod, only: g, pi, mugaz, avocado
        use geometry_mod, only: latitude
        use callkeys_mod, only: aerofixn2,kastprof,cloudlvl,     &
-                dusttau,                                &
                 pres_bottom_tropo,pres_top_tropo,obs_tau_col_tropo,     &
                 pres_bottom_strato,pres_top_strato,obs_tau_col_strato,  &
-                tau_nh3_cloud, pres_nh3_cloud,                          &
-                nlayaero, aeronlay_tauref, aeronlay_choice,             & 
+                nlayaero, aeronlay_tauref, aeronlay_choice,             &
                 aeronlay_pbot, aeronlay_ptop, aeronlay_sclhght,         &
                 aerogeneric
@@ -29 +25 @@
 
 !==================================================================
-!     
+!
 !     Purpose
 !     -------
 !     Compute aerosol optical depth in each gridbox.
-!     
+!
 !     Authors
-!     ------- 
+!     -------
 !     F. Forget
-!     F. Montmessin (water ice scheme) 
+!     F. Montmessin (water ice scheme)
 !     update J.-B. Madeleine (2008)
 !     dust removal, simplification by Robin Wordsworth (2009)
@@ -44 +40 @@
 !
 !     Input
-!     ----- 
+!     -----
 !     ngrid             Number of horizontal gridpoints
 !     nlayer            Number of layers
@@ -97 +93 @@
       CHARACTER(LEN=20) :: tracername ! to temporarily store text
 
-      ! for fixed dust profiles
-      real topdust, expfactor, zp
-      REAL taudusttmp(ngrid) ! Temporary dust opacity used before scaling
-      REAL tauh2so4tmp(ngrid) ! Temporary h2so4 opacity used before scaling
-
       real CLFtot
       integer igen_ice,igen_vap ! to store the index of generic tracer
@@ -131 +122 @@
         endif
 
-        if ((iaero_n2.ne.0).and.(.not.noaero)) then
-          print*, 'iaero_n2=  ',iaero_n2
-        endif
-        if (iaero_dust.ne.0) then
-          print*,'iaero_dust= ',iaero_dust
-        endif
-        if (iaero_h2so4.ne.0) then
-          print*,'iaero_h2so4= ',iaero_h2so4
-        endif
-        if (iaero_back2lay.ne.0) then
-          print*,'iaero_back2lay= ',iaero_back2lay
-        endif
-        if (iaero_nh3.ne.0) then
-          print*,'iaero_nh3= ',iaero_nh3
-        endif
-        if (iaero_nlay(1).ne.0) then
-          print*,'iaero_nlay= ',iaero_nlay(:)
-        endif
-        if (iaero_aurora.ne.0) then
-          print*,'iaero_aurora= ',iaero_aurora
-        endif
-        if (aerogeneric .ne. 0) then 
-          print*,"iaero_generic= ",iaero_generic(:)
-        endif
         firstcall=.false.
       ENDIF ! of IF (firstcall)
@@ -161 +128 @@
 !     ---------------------------------------------------------
 !==================================================================
-!    N2 ice aerosols
+!    Haze aerosols
 !==================================================================
 
-      if (iaero_n2.ne.0) then
-           iaer=iaero_n2
+      if (iaero_haze.ne.0) then
+        iaer=iaero_haze
 !       1. Initialization
-            aerosol(1:ngrid,1:nlayer,iaer)=0.0
+         aerosol(1:ngrid,1:nlayer,iaer)=0.0
 !       2. Opacity calculation
-            if (noaero) then ! aerosol set to zero
-             aerosol(1:ngrid,1:nlayer,iaer)=0.0
-            elseif (aerofixn2.or.(i_n2ice.eq.0)) then !  N2 ice cloud prescribed
-               aerosol(1:ngrid,1:nlayer,iaer)=1.e-9
-               !aerosol(1:ngrid,12,iaer)=4.0 ! single cloud layer option
-            else
-               DO ig=1, ngrid
-                  DO l=1,nlayer-1 ! to stop the rad tran bug
-
-                     aerosol0 =                         &
-                          (  0.75 * QREFvis3d(ig,l,iaer) /        &
-                          ( rho_n2 * reffrad(ig,l,iaer) )  ) *   &
-                          ( pq(ig,l,i_n2ice) + 1.E-9 ) *         &
-                          ( pplev(ig,l) - pplev(ig,l+1) ) / g
-                     aerosol0           = max(aerosol0,1.e-9)
-                     aerosol0           = min(aerosol0,L_TAUMAX)
-                     aerosol(ig,l,iaer) = aerosol0
-!                     aerosol(ig,l,iaer) = 0.0
-!                     print*, aerosol(ig,l,iaer)
-!        using cloud fraction
-!                     aerosol(ig,l,iaer) = -log(1 - CLF + CLF*exp(-aerosol0/CLF))
-!                     aerosol(ig,l,iaer) = min(aerosol(ig,l,iaer),L_TAUMAX)
-
-
-                  ENDDO
+         DO ig=1, ngrid
+               DO l=1,nlayer-1 ! to stop the rad tran bug
+                  ! if fractal, radius doit etre equivalent sphere radius
+                  aerosol0 =                         &
+                       (  0.75 * QREFvis3d(ig,l,iaer) /        &
+                       ( rho_q(i_haze) * reffrad(ig,l,iaer) )  ) *   &
+                       ( pq(ig,l,i_haze) + 1.E-10 ) *         &
+                       ( pplev(ig,l) - pplev(ig,l+1) ) / g
+                  aerosol0           = max(aerosol0,1.e-10)
+                  aerosol0           = min(aerosol0,L_TAUMAX)
+                  aerosol(ig,l,iaer) = aerosol0
                ENDDO
-            end if ! if fixed or varying
-      end if ! if N2 aerosols   
-!==================================================================
-!     Water ice / liquid (AF24: removed)
-!==================================================================
-
-!==================================================================
-!             Dust (AF24: removed)
-!==================================================================
-
-!==================================================================
-!           H2SO4 (AF24: removed)
-!==================================================================
-!==================================================================
-!    Two-layer aerosols (unknown composition)
-!    S. Guerlet (2013) - Modif by J. Vatant d'Ollone (2020)
-!    
-!    This scheme is deprecated and left for retrocompatibility
-!    You should use the n-layer scheme below !
-!
-!==================================================================
-
-!==================================================================
-!    Saturn/Jupiter ammonia cloud = thin cloud (scale height 0.2 hard coded...)
-!    S. Guerlet (2013)
-!    JVO 20 : You should now use the generic n-layer scheme below
-!==================================================================
-
-!=========================================================================================================
-!    Generic N-layers aerosols/clouds
-!    Author : J. Vatant d'Ollone (2020)
-!    
-!    Purpose: Replaces and extents the former buggy 2-layer scheme as well as hard-coded NH3 cloud
-!    
-!    + Each layer can have different optical properties, size of particle ...
-!    + Enables up to n=4 layers as we apparently cannot run with more scatterers (could be worth checking...)
-!    + You have different choices for vertical profile of the aerosol layers :
-!           * aeronlay_choice = 1 : Layer tau is spread between ptop and pbot following atm scale height.
-!           * aeronlay_choice = 2 : Layer tau follows its own scale height above cloud deck (pbot).
-!                                   In this case ptop is dummy and sclhght gives the ratio H_cl/H_atm.
-!           * aeronlay_choice = ... feel free to add more cases  !
-!    + Layers can overlap if needed (if you want a 'transition layer' as in the 2-scheme, just add it)
-!
-!=========================================================================================================
-
-      if (iaero_nlay(1) .ne.0) then
-
-        DO ia=1,nlayaero
-           iaer=iaero_nlay(ia)
-
-!          a. Initialization
-           aerosol(1:ngrid,1:nlayer,iaer)=0.D0
-
-!          b. Opacity calculation
-           
-           ! Case 1 : Follows atmospheric scale height between boundaries pressures
-           ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-           IF (aeronlay_choice(ia).eq.1) THEN
-
-             dp_layer(:)=0.D0
-             DO ig=1,ngrid
-               DO l=1,nlayer-1
-                 !! i. Opacity follows scale height
-                 IF ( pplev(ig,l).le.aeronlay_pbot(ia)   .AND.          &
-                      pplev(ig,l).ge.aeronlay_ptop(ia) ) THEN
-                   aerosol(ig,l,iaer) = ( pplev(ig,l) - pplev(ig,l+1) )
-                   dp_layer(ig) = dp_layer(ig) + aerosol(ig,l,iaer)
-                 !! ii. Outside aerosol layer boundaries: no aerosols
-                 ELSE
-                   aerosol(ig,l,iaer) = 0.D0
-                 ENDIF
-               ENDDO
-             ENDDO
-             ! iii. Re-normalize to required total opacity
-             DO ig=1,ngrid
-               DO l=1,nlayer-1
-                 IF ( pplev(ig,l).le.aeronlay_pbot(ia)   .AND.          &
-                      pplev(ig,l).ge.aeronlay_ptop(ia) ) THEN
-                  aerosol(ig,l,iaer) = aerosol(ig,l,iaer) / dp_layer(ig) &
-                                     * aeronlay_tauref(ia)
-                 ENDIF
-               ENDDO
-             ENDDO
-
-           ! Case 2 : Follows input scale height
-           ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-           ELSE IF (aeronlay_choice(ia).eq.2) THEN
-           
-             cloud_slope  = 1.D0/aeronlay_sclhght(ia)
-             pcloud_deck  = aeronlay_pbot(ia)
-             dp_layer(:)  = 0.D0
-
-             DO ig=1,ngrid
-               DO l=1,nlayer-1
-                 !! i. Below cloud layer: no opacity
-                 IF (pplev(ig,l) .gt. pcloud_deck) THEN
-                   aerosol(ig,l,iaer) = 0.D0            
-                 !! ii. Follows scale height above cloud deck
-                 ELSEIF (pplev(ig,l) .le. pcloud_deck) THEN 
-                   aerosol(ig,l,iaer) = ((pplev(ig,l)/pcloud_deck)**(cloud_slope))
-                   dp_layer(ig) = dp_layer(ig) + aerosol(ig,l,iaer)
-                 ENDIF
-               ENDDO
-             ENDDO
-             ! iii. Re-normalize to required total opacity
-             DO ig=1,ngrid
-               DO l=1,nlayer-1
-                 IF (pplev(ig,l) .le. pcloud_deck) THEN 
-                  aerosol(ig,l,iaer) = aerosol(ig,l,iaer) / dp_layer(ig) &
-                                     * aeronlay_tauref(ia)
-                 ENDIF
-               ENDDO
-             ENDDO
-
-           ENDIF ! aeronlay_choice
-
-          ENDDO ! loop on n aerosol layers
-
-      end if ! if N-layer aerosols
-  
-!==================================================================
-!    Jovian auroral aerosols (unknown composition) NON-GENERIC: vertical and meridional profile tuned to observations
-!    S. Guerlet (2015)
-!==================================================================
-
-
-      if (iaero_aurora .ne.0) then
-           iaer=iaero_aurora
-!       1. Initialization
-            aerosol(1:ngrid,1:nlayer,iaer)=0.D0 
-         pm = 2000. !!case study: maxi aerosols at 20 hPa
-!       2. Opacity calculation
-          DO ig=1,ngrid
-
-          !! Test Jupiter (based on Zhang et al 2013 observations, but a bit different), decembre 2015
-              DO l=1,nlayer
-                aerosol(ig,l,iaer) = (pplev(ig,l)/pm)**2 * exp(-(pplev(ig,l)/pm)**2)
-              ENDDO
-          ENDDO
-         
- !       3. Meridional distribution, and re-normalize to observed total column
-         dp_layer(:)=0.D0
-         DO ig=1,ngrid
-          !!Jupiter
-          !!Hem sud:
-          IF (latitude(ig)*180.D0/pi .lt. -45.D0 .and. latitude(ig)*180.D0/pi .gt. -70.) THEN
-          obs_tau_col_aurora= 10.D0**(-0.06D0*latitude(ig)*180.D0/pi-3.4D0) 
-          ELSEIF (latitude(ig)*180.D0/pi .lt. -37.D0 .and. latitude(ig)*180.D0/pi .ge. -45.) THEN
-          obs_tau_col_aurora= 10.D0**(-0.3D0*latitude(ig)*180.D0/pi-14.3D0) 
-           ELSEIF (latitude(ig)*180./pi .le. -70. ) THEN
-          obs_tau_col_aurora= 10**(0.06*70.-3.4D0) 
-          !!Hem Nord:  
-          ELSEIF (latitude(ig)*180.D0/pi .gt. 30.D0 .and. latitude(ig)*180.D0/pi .lt. 70.) THEN
-          obs_tau_col_aurora= 10.D0**(0.03D0*latitude(ig)*180.D0/pi-1.17D0)  
-          ELSEIF (latitude(ig)*180.D0/pi .gt. 22.D0 .and. latitude(ig)*180.D0/pi .le. 30.) THEN
-          obs_tau_col_aurora= 10.D0**(0.3D0*latitude(ig)*180.D0/pi-9.4D0)  
-          ELSEIF (latitude(ig)*180.D0/pi .ge. 70.) THEN
-          obs_tau_col_aurora= 10**(0.03*70.-1.17D0)  
-          ELSEIF (latitude(ig)*180.D0/pi .ge. -37. .and. latitude(ig)*180.D0/pi .le. 22.) THEN
-         obs_tau_col_aurora = 0.001D0    !!Jupiter: mini pas a zero
-          ENDIF
-
-          DO l=1,nlayer  
-               dp_layer(ig) = dp_layer(ig) + aerosol(ig,l,iaer)/obs_tau_col_aurora
-          ENDDO
          ENDDO
-
-         DO ig=1,ngrid
-           DO l=1,nlayer-1
-                aerosol(ig,l,iaer)=aerosol(ig,l,iaer)/dp_layer(ig)
-           ENDDO
-         ENDDO
-
-
-      end if ! if Auroral aerosols  
-!===========================================================================
-!    Radiative Generic Condensable aerosols scheme
-!    Only used when we give aerogeneric != 0 in callphys.def
-!    Computes the generic aerosols' opacity in the same fashion as water of 
-!    dust, using the QREFvis3d of the concerned specie
-!    Lucas Teinturier (2022)
-!===========================================================================
-      if (aerogeneric .ne. 0) then ! we enter the scheme
-        do ia=1,aerogeneric
-          iaer = iaero_generic(ia)
-          ! Initialization
-          aerosol(1:ngrid,1:nlayer,iaer) = 0.D0
-          igen_ice = i_rgcs_ice(ia)
-          ! Let's loop on the horizontal and vertical grid
-          do ig=1,ngrid
-            do l=1,nlayer
-              aerosol(ig,l,iaer) = ( 0.75*QREFvis3d(ig,l,iaer)  / &
-                                  (rho_q(igen_ice) * reffrad(ig,l,iaer)) ) * &
-                                  (pq(ig,l,igen_ice)+1E-9 ) *                &
-                                  (pplev(ig,l) - pplev(ig,l+1)) /g 
-            enddo !l=1,nlayer
-          enddo !ig=1,ngrid
-        enddo !ia=1,aerogeneric
-      endif !aerogeneric .ne. 0
+         !QREF est le meme dans toute la colonne par def si size uniforme
+         !print*, 'TB17: QREFvis3d=',QREFvis3d(1,:,1)
+         !print*, 'TB17: rho_q=',rho_q(i_haze)
+         !print*, 'TB17: reffrad=',reffrad(1,:,1)
+         !print*, 'TB17: pq=',pq(1,:,i_haze)
+         !print*, 'TB17: deltap=',pplev(1,1) - pplev(1,nlayer)
+   end if ! if haze aerosols
 
 ! --------------------------------------------------------------------------
 ! Column integrated visible optical depth in each point (used for diagnostic)
 
-      tau_col(:)=0.0
-      do iaer = 1, naerkind
-         do l=1,nlayer
-            do ig=1,ngrid
-               tau_col(ig) = tau_col(ig) + aerosol(ig,l,iaer)
-            end do
+   tau_col(:)=0.0
+   do iaer = 1, naerkind
+      do l=1,nlayer
+         do ig=1,ngrid
+            tau_col(ig) = tau_col(ig) + aerosol(ig,l,iaer)
          end do
       end do
+   end do
 
-      ! do ig=1,ngrid
-      !    do l=1,nlayer
-      !       do iaer = 1, naerkind
-      !          if(aerosol(ig,l,iaer).gt.1.e3)then
-      !             print*,'WARNING: aerosol=',aerosol(ig,l,iaer)
-      !             print*,'at ig=',ig,',  l=',l,', iaer=',iaer
-      !             print*,'QREFvis3d=',QREFvis3d(ig,l,iaer)
-      !             print*,'reffrad=',reffrad(ig,l,iaer)
-      !          endif
-      !       end do
-      !    end do
-      ! end do
+   do ig=1,ngrid
+      do l=1,nlayer
+         do iaer = 1, naerkind
+            if(aerosol(ig,l,iaer).gt.1.e3)then
+               print*,'WARNING: aerosol=',aerosol(ig,l,iaer)
+               print*,'at ig=',ig,',  l=',l,', iaer=',iaer
+               print*,'QREFvis3d=',QREFvis3d(ig,l,iaer)
+               print*,'reffrad=',reffrad(ig,l,iaer)
+            endif
+         end do
+      end do
+   end do
 
-      ! do ig=1,ngrid
-      !    if(tau_col(ig).gt.1.e3)then
-      !       print*,'WARNING: tau_col=',tau_col(ig)
-      !       print*,'at ig=',ig
-      !       print*,'aerosol=',aerosol(ig,:,:)
-      !       print*,'QREFvis3d=',QREFvis3d(ig,:,:)
-      !       print*,'reffrad=',reffrad(ig,:,:)
-      !    endif
-      ! end do
+   do ig=1,ngrid
+      if(tau_col(ig).gt.1.e3)then
+         print*,'WARNING: tau_col=',tau_col(ig)
+         print*,'at ig=',ig
+         print*,'aerosol=',aerosol(ig,:,:)
+         print*,'QREFvis3d=',QREFvis3d(ig,:,:)
+         print*,'reffrad=',reffrad(ig,:,:)
+      endif
+   end do
+  end subroutine aeropacity
 
-    end subroutine aeropacity
-      
 end module aeropacity_mod

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

@@ -8 +8 @@
 !                 corresponding aerosol was not activated in callphys.def
 !                 -- otherwise a value is set via iniaerosol
-      integer, save, protected :: iaero_n2 = 0 
-      integer, save, protected :: iaero_dust = 0
-      integer, save, protected :: iaero_h2so4 = 0
+      integer, save :: iaero_haze = 0
+      integer, save :: i_haze = 0
       logical, save, protected :: noaero = .false.
-!$OMP THREADPRIVATE(iaero_n2,iaero_h2o,iaero_dust,iaero_h2so4,noaero)
+!$OMP THREADPRIVATE(iaero_haze,i_haze,noaero)
 
 ! two-layer simple aerosol model
       integer, save, protected :: iaero_back2lay = 0
- ! NH3 cloud
-      integer, save, protected :: iaero_nh3 = 0
 ! N-layer aerosol model (replaces the 2-layer and hard-coded clouds)
       integer,dimension(:), allocatable, save, protected :: iaero_nlay
-! Auroral aerosols
-      integer, save, protected :: iaero_aurora = 0
-!$OMP THREADPRIVATE(iaero_back2lay,iaero_nh3,iaero_nlay,iaero_aurora)
+!$OMP THREADPRIVATE(iaero_back2lay,iaero_nlay)
 
 ! Generic aerosols
@@ -33 +28 @@
 contains
 
-  SUBROUTINE iniaerosol
+  !==================================================================
+   subroutine haze_prof(ngrid,nlayer,zzlay,pplay,pt,reffrad,profmmr)
+!==================================================================
+!     Purpose
+!     -------
+!     Get fixed haze properties
+!     profile of haze (from txt file) and fixed radius profile
+!
+!==================================================================
+      use radinc_h, only: naerkind
+      use datafile_mod
+      use tracer_h
+      use comcstfi_mod, only: r, pi
 
-  use mod_phys_lmdz_para, only : is_master
-  use radinc_h, only: naerkind
-  use tracer_h, only: n_rgcs, nqtot, is_rgcs
-  use callkeys_mod, only: aeron2, dusttau, &
-                          aeroback2lay, aeronh3, nlayaero, aeronlay, &
-                          aeroaurora, aerogeneric
+!-----------------------------------------------------------------------
+!     Arguments
+      Implicit none
 
-  IMPLICIT NONE
-!=======================================================================
-!   subject:
-!   --------
-!   Initialization related to aerosols 
-!   (N2 aerosols, dust, water, chemical species, ice...)   
-!
-!   author: Laura Kerber, S. Guerlet
-!   ------
-!        
-!=======================================================================
+      integer,intent(in) :: ngrid
+      integer,intent(in) :: nlayer
+      real,intent(in) :: zzlay(ngrid,nlayer)
+      real,intent(in) :: pplay(ngrid,nlayer)
+      real,intent(in) :: pt(ngrid,nlayer)
+      real, intent(in) :: reffrad(ngrid,nlayer,naerkind)      ! haze particles radii (m)
 
-  integer :: i, ia, iq
+      real, intent(out) :: profmmr(ngrid,nlayer)              ! mmr haze kg/kg
 
-  ! Special case, dyn. allocation for n-layer depending on callphys.def
-  IF(.NOT.ALLOCATED(iaero_nlay)) ALLOCATE(iaero_nlay(nlayaero))
-  iaero_nlay(:) = 0
-  ! Do the same for iaero_generic and i_rgcs_ice
-  IF (.not. allocated(iaero_generic)) allocate(iaero_generic(aerogeneric))
-  if (.not. allocated(i_rgcs_ice)) allocate(i_rgcs_ice(aerogeneric)) 
+!     Local variables
+      integer :: iaer,l,ig,ifine
 
-  ! Init of i_rgcs_ice
-  i_rgcs_ice(:) =0
-  ia = 1
-  do iq=1,nqtot
-    if (is_rgcs(iq) .eq. 1) then 
-        i_rgcs_ice(ia)=iq
-        ia = ia+1
-     endif
-  enddo
+      LOGICAL firstcall
+      SAVE firstcall
+      DATA firstcall/.true./
 
-  iaero_generic(:)=0
-  ia=0
-  if (aeron2) then
-     ia=ia+1
-     iaero_n2=ia
-  endif
-  if (is_master) write(*,*) '--- N2 aerosol = ', iaero_n2
+      !!read altitudes and haze mmrs
+      integer Nfine
+      !parameter(Nfine=21)
+      parameter(Nfine=701)
+      character(len=100) :: file_path
+      real,save :: levdat(Nfine),densdat(Nfine)
 
-  if (dusttau.gt.0) then
-     ia=ia+1
-     iaero_dust=ia
-  endif
-  if (is_master) write(*,*) '--- Dust aerosol = ', iaero_dust
+!---------------- INPUT ------------------------------------------------
 
-     
-  if (aeroback2lay) then
-     ia=ia+1
-     iaero_back2lay=ia
-  endif
-  if (is_master) write(*,*) '--- Two-layer aerosol = ', iaero_back2lay
+      !! Read data
+      IF (firstcall) then
+        firstcall=.false.
+        file_path=trim(datadir)//'/haze_prop/hazemmr.txt'
+        open(224,file=file_path,form='formatted')
+        do ifine=1,Nfine
+           read(224,*) levdat(ifine), densdat(ifine)
+        enddo
+        close(224)
+        print*, 'TB22 read Haze MMR profile'
+      ENDIF
 
-  if (aeronh3) then
-     ia=ia+1
-     iaero_nh3=ia
-  endif
-  if (is_master) write(*,*) '--- NH3 Cloud = ', iaero_nh3
+      !! Interpolate on the model vertical grid
+      do ig=1,ngrid
+        CALL interp_line(levdat,densdat,Nfine,zzlay(ig,:)/1000.,profmmr(ig,:),nlayer)
+      enddo
 
-  if (aeronlay) then
-     do i=1,nlayaero
-       ia=ia+1
-       iaero_nlay(i)=ia
-     enddo
-  endif
-  if (is_master) write(*,*) '--- N-layer aerosol = ', iaero_nlay
+      !! Get profile Mass mixing ratio from number density:    part.cm-3 --> m-3 --> m3 m-3
+      !                                --> kg m-3 --> kg/kg
+      do iaer=1,naerkind
+            if(iaer.eq.iaero_haze.and.1.eq.2) then !TB22 activate/deactivate mmr or part density
+              !print*, 'Haze profile is fixed'
+              do ig=1,ngrid
+                 do l=1,nlayer
+                    !from number density in cm-3
+                    profmmr(ig,l)=profmmr(ig,l)*1.e6*4./3.*pi*reffrad(ig,l,iaer)**3*rho_q(i_haze)/(pplay(ig,l)/(r*pt(ig,l)))
+                 enddo
+              enddo
+            endif
+      enddo
+   end subroutine haze_prof
+!==================================================================
 
-  if (aeroaurora) then
-     ia=ia+1
-     iaero_aurora=ia
-  endif
-  if (is_master) write(*,*) '--- Auroral aerosols = ', iaero_aurora
-
-  if (aerogeneric .ne. 0) then 
-     do i=1,aerogeneric
-        ia = ia+1
-        iaero_generic(i) = ia
-     enddo
-  endif 
-      
-  if (is_master) then
-    write(*,*)'--- Radiative Generic Condensable Species = ',iaero_generic
-
-    write(*,*) '=== Number of aerosols= ', ia
-  endif ! of is_master
-
-! For the zero aerosol case, we currently make a dummy n2 aerosol which is zero everywhere.
-! (See aeropacity.F90 for how this works). A better solution would be to turn off the 
-! aerosol machinery in the no aerosol case, but this would be complicated. LK
-
-  if (ia.eq.0) then  !For the zero aerosol case. 
-     ia = 0
-     noaero = .true.
-  endif
-
-  if (.not.noaero .and. ia.ne.naerkind) then
-    if (is_master) then
-      print*, 'Aerosols counted not equal to naerkind'
-      print*, 'set correct value for nearkind in callphys.def'
-      print*, 'which should be ',ia
-      print*, 'according to current options in callphys.def'
-      print*, 'or change/correct incompatible options there'
-      print*, 'Abort in iniaerosol'
-    endif
-    call abort_physic("iniaerosl",'wrong number of aerosols',1)
-  endif ! of if (ia.ne.naerkind)
-
-  END SUBROUTINE iniaerosol
 
 end module aerosol_mod

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

@@ -6 +6 @@
 
       subroutine callcorrk(ngrid,nlayer,pq,nq,qsurf,           &
-          albedo,albedo_equivalent,emis,mu0,pplev,pplay,pt,    & 
+          albedo,albedo_equivalent,emis,mu0,pplev,pplay,pt,    &
           tsurf,fract,dist_star,aerosol,muvar,                 &
           dtlw,dtsw,fluxsurf_lw,                               &
@@ -26 +26 @@
       use ioipsl_getin_p_mod, only: getin_p
       use gases_h, only: ngasmx
-      use radii_mod, only : su_aer_radii,n2_reffrad,dust_reffrad,h2so4_reffrad,back2lay_reffrad
-      use aerosol_mod, only : iaero_n2,iaero_dust,iaero_h2so4, &
-                              iaero_back2lay, iaero_aurora
+      use radii_mod, only : su_aer_radii
+      use aerosol_mod, only : iaero_haze
       use aeropacity_mod, only: aeropacity
       use aeroptproperties_mod, only: aeroptproperties
-      use tracer_h, only: igcm_n2_ice
       use tracer_h, only: constants_epsi_generic
       use comcstfi_mod, only: pi, mugaz, cpp
@@ -37 +35 @@
                               diagdtau,kastprof,strictboundcorrk,specOLR, &
                               tplanckmin,tplanckmax,global1d, &
-                              generic_condensation
+                              generic_condensation, aerohaze, haze_radproffix
       use optcv_mod, only: optcv
       use optci_mod, only: optci
@@ -57 +55 @@
 !
 !     Authors
-!     ------- 
+!     -------
 !     Emmanuel 01/2001, Forget 09/2001
 !     Robin Wordsworth (2009)
@@ -65 +63 @@
 !-----------------------------------------------------------------------
 !     Declaration of the arguments (INPUT - OUTPUT) on the LMD GCM grid
-!     Layer #1 is the layer near the ground. 
+!     Layer #1 is the layer near the ground.
 !     Layer #nlayer is the layer at the top.
 !-----------------------------------------------------------------------
@@ -90 +88 @@
       logical,intent(in) :: firstcall              ! Signals first call to physics.
       logical,intent(in) :: lastcall               ! Signals last call to physics.
-      
+
       ! OUTPUT
       REAL,INTENT(OUT) :: aerosol(ngrid,nlayer,naerkind) ! Aerosol tau at reference wavelenght.
@@ -109 +107 @@
       REAL,INTENT(OUT) :: int_dtaui(ngrid,nlayer,L_NSPECTI) ! VI optical thickness of layers within narrowbands for diags ().
       REAL,INTENT(OUT) :: int_dtauv(ngrid,nlayer,L_NSPECTV) ! IR optical thickness of layers within narrowbands for diags ().
-      
-      
-      
-      
-
-! Globally varying aerosol optical properties on GCM grid ; not needed everywhere so not in radcommon_h.   
+
+
+
+
+
+! Globally varying aerosol optical properties on GCM grid ; not needed everywhere so not in radcommon_h.
 ! made "save" variables so they are allocated once in for all, not because
 ! the values need be saved from a time step to the next
@@ -208 +206 @@
       real vtmp(nlayer)
       REAL*8 muvarrad(L_LEVELS)
-      
-      ! Included by MT for albedo calculations.      
+
+      ! Included by MT for albedo calculations.
       REAL*8 albedo_temp(L_NSPECTV) ! For equivalent albedo calculation.
       REAL*8 surface_stellar_flux   ! Stellar flux reaching the surface. Useful for equivalent albedo calculation.
-     
+
       ! local variable
       integer ok ! status (returned by NetCDF functions)
@@ -254 +252 @@
           allocate(tauaero(L_LEVELS,naerkind))
         endif
-        
-        if(.not.allocated(QXVAER)) then 
+
+        if(.not.allocated(QXVAER)) then
           allocate(QXVAER(L_LEVELS,L_NSPECTV,naerkind), stat=ok)
           if (ok /= 0) then
@@ -334 +332 @@
 
          call su_aer_radii(ngrid,nlayer,reffrad,nueffrad)
-         
-         
+
+
 !--------------------------------------------------
 !             Set up correlated k
@@ -353 +351 @@
          ! call sugas_corrk       ! Set up gaseous absorption properties.
          ! call suaer_corrk       ! Set up aerosol optical properties.
-        
+
 
          ! now that L_NGAUSS has been initialized (by sugas_corrk)
@@ -435 +433 @@
          endif
 
-         
+
          if(varfixed .and. generic_condensation)then
             write(*,*) "Deep generic tracer vapor mixing ratio ? (no effect if negative) "
@@ -449 +447 @@
 
 !=======================================================================
-!          I.b  Initialization on every call   
+!          I.b  Initialization on every call
 !=======================================================================
- 
+
       qxvaer(:,:,:)=0.0
       qsvaer(:,:,:)=0.0
@@ -467 +465 @@
 !     Effective radius and variance of the aerosols
 !--------------------------------------------------
-
-      do iaer=1,naerkind
-
-         if ((iaer.eq.iaero_n2).and.tracer.and.(igcm_n2_ice.gt.0)) then ! Treat condensed n2 particles.
-            call n2_reffrad(ngrid,nlayer,nq,pq,reffrad(1,1,iaero_n2))
-
-            call planetwide_maxval(reffrad(:,:,iaero_n2),maxvalue)
-            call planetwide_minval(reffrad(:,:,iaero_n2),minvalue)
-            if (is_master) then
-              print*,'Max. N2 ice particle size = ',maxvalue/1.e-6,' um'
-               print*,'Min. N2 ice particle size = ',minvalue/1.e-6,' um'
-            end if
-         end if
-         ! Trear other aerosols here
-
-         if(iaer.eq.iaero_dust)then
-            call dust_reffrad(ngrid,nlayer,reffrad(1,1,iaero_dust))
-            if (is_master) then
-               print*,'Dust particle size = ',reffrad(1,1,iaer)/1.e-6,' um'
-            end if
-         endif
-         
-         if(iaer.eq.iaero_h2so4)then
-            call h2so4_reffrad(ngrid,nlayer,reffrad(1,1,iaero_h2so4))
-            if (is_master) then
-               print*,'H2SO4 particle size =',reffrad(1,1,iaer)/1.e-6,' um'
-            end if
-         endif
-         
-          if(iaer.eq.iaero_back2lay)then
-            call back2lay_reffrad(ngrid,reffrad(1,1,iaero_back2lay),nlayer,pplev)
-         endif
-
-         !  For n-layer aerosol size set once for all at firstcall in su_aer_radii
-
-!         if(iaer.eq.iaero_aurora)then
-!           call aurora_reffrad(ngrid,nlayer,reffrad(1,1,iaero_aurora))
-!         endif
-        
-     end do !iaer=1,naerkind.
+      if (aerohaze) then
+         do iaer=1,naerkind
+            if ((iaer.eq.iaero_haze)) then
+               call su_aer_radii(ngrid,nlayer,reffrad(1,1,iaer), &
+                            nueffrad(1,1,iaer))
+            endif
+         end do !iaer=1,naerkind.
+         if (haze_radproffix) then
+            print*, 'haze_radproffix=T : fixed profile for haze rad'
+         else
+            print*,'reffrad haze:',reffrad(1,1,iaero_haze)
+            print*,'nueff haze',nueffrad(1,1,iaero_haze)
+         endif
+      endif
 
 
@@ -518 +490 @@
            QVISsQREF3d,omegaVIS3d,gVIS3d,                          &
            QIRsQREF3d,omegaIR3d,gIR3d,                             &
-           QREFvis3d,QREFir3d)                                     
+           QREFvis3d,QREFir3d)
 
       ! Get aerosol optical depths.
       call aeropacity(ngrid,nlayer,nq,pplay,pplev, pt,pq,aerosol,      &
-           reffrad,nueffrad,QREFvis3d,QREFir3d,                             & 
-           tau_col,cloudfrac,totcloudfrac,clearsky)                
- 
-!-----------------------------------------------------------------------    
+           reffrad,nueffrad,QREFvis3d,QREFir3d,                             &
+           tau_col,cloudfrac,totcloudfrac,clearsky)
+
+!-----------------------------------------------------------------------
       do ig=1,ngrid ! Starting Big Loop over every GCM column
 !-----------------------------------------------------------------------
@@ -539 +511 @@
 !    The transformation in the vertical is the same as for temperature.
 !-----------------------------------------------------------------------
-           
-           
+
+
             do iaer=1,naerkind
                ! Shortwave.
-               do nw=1,L_NSPECTV 
-               
+               do nw=1,L_NSPECTV
+
                   do l=1,nlayer
 
@@ -580 +552 @@
 
                end do ! L_NSPECTV
-             
+
                do nw=1,L_NSPECTI
                   ! Longwave
@@ -618 +590 @@
 
                end do ! L_NSPECTI
-               
+
             end do ! naerkind
 
@@ -627 +599 @@
                   do nw=1,L_NSPECTV
                      if(qsvaer(k,nw,iaer).gt.1.05*qxvaer(k,nw,iaer))then
-                        message='Serious problems with qsvaer values' 
+                        message='Serious problems with qsvaer values'
                         call abort_physic(subname,message,1)
                      endif
@@ -635 +607 @@
                   end do
 
-                  do nw=1,L_NSPECTI 
+                  do nw=1,L_NSPECTI
                      if(qsiaer(k,nw,iaer).gt.1.05*qxiaer(k,nw,iaer))then
-                        message='Serious problems with qsvaer values' 
+                        message='Serious problems with qsvaer values'
                         call abort_physic(subname,message,1)
                      endif
@@ -651 +623 @@
 !     Aerosol optical depths
 !-----------------------------------------------------------------------
-            
-         do iaer=1,naerkind     ! a bug was here           
+
+         do iaer=1,naerkind     ! a bug was here
             do k=0,nlayer-1
-               
+
                pweight=(pplay(ig,L_NLAYRAD-k)-pplev(ig,L_NLAYRAD-k+1))/   &
                        (pplev(ig,L_NLAYRAD-k)-pplev(ig,L_NLAYRAD-k+1))
@@ -667 +639 @@
             tauaero(1,iaer)          = tauaero(2,iaer)
             !tauaero(1,iaer)          = 0.
-            !JL18 at time of testing, the two above conditions gave the same results bit for bit. 
-            
+            !JL18 at time of testing, the two above conditions gave the same results bit for bit.
+
          end do ! naerkind
 
@@ -696 +668 @@
 
             call generic_tracer_index(nq,iq,igcm_generic_vap,igcm_generic_ice,call_ice_vap_generic)
-            
+
             if (call_ice_vap_generic) then ! to call only one time the ice/vap pair of a tracer
 
@@ -704 +676 @@
                   do l=1,nlayer
                      qvar(2*l)   = pq(ig,nlayer+1-l,i_var)
-                     qvar(2*l+1) = pq(ig,nlayer+1-l,i_var)    
-                     !JL13index            qvar(2*l+1) = (pq(ig,nlayer+1-l,i_var)+pq(ig,max(nlayer-l,1),i_var))/2    
+                     qvar(2*l+1) = pq(ig,nlayer+1-l,i_var)
+                     !JL13index            qvar(2*l+1) = (pq(ig,nlayer+1-l,i_var)+pq(ig,max(nlayer-l,1),i_var))/2
                      !JL13index            ! Average approximation as for temperature...
                   end do
@@ -713 +685 @@
 
                   do l=1,nlayer ! Here we will assign fixed water vapour profiles globally.
-                                          
+
                      call Psat_generic(pt(ig,l),pplay(ig,l),metallicity,psat,qsat)
 
                      if (qsat .lt. qvap_deep) then
                         pq_temp(l) = qsat      ! fully saturated everywhere
-                     else 
+                     else
                         pq_temp(l) = qvap_deep
                      end if
 
                   end do
-                  
+
                   do l=1,nlayer
                      qvar(2*l)   = pq_temp(nlayer+1-l)
                      qvar(2*l+1) = (pq_temp(nlayer+1-l)+pq_temp(max(nlayer-l,1)))/2
                   end do
-                  
+
                   qvar(1)=qvar(2)
-         
+
                else
                   do k=1,L_LEVELS
@@ -739 +711 @@
             endif
 
-         end do ! do iq=1,nq loop on tracers 
+         end do ! do iq=1,nq loop on tracers
 
       end if ! if (generic_condensation)
@@ -780 +752 @@
          muvarrad(2*l+1) = (muvar(ig,nlayer+2-l)+muvar(ig,max(nlayer+1-l,1)))/2
       END DO
-   
+
       muvarrad(1) = muvarrad(2)
-      muvarrad(2*nlayer+1)=muvar(ig,1)         
-      
+      muvarrad(2*nlayer+1)=muvar(ig,1)
+
       ! Keep values inside limits for which we have radiative transfer coefficients !!!
       if(L_REFVAR.gt.1)then ! (there was a bug here)
@@ -805 +777 @@
          tlevrad(2*l+1) = (pt(ig,nlayer+1-l)+pt(ig,max(nlayer-l,1)))/2
       END DO
-      
+
       plevrad(1) = 0.
-!      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. 
-      
+!      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.
+
       tlevrad(1) = tlevrad(2)
       tlevrad(2*nlayer+1)=tsurf(ig)
-      
-      pmid(1) = pplay(ig,nlayer)/scalep   
+
+      pmid(1) = pplay(ig,nlayer)/scalep
       pmid(2) =  pmid(1)
 
       tmid(1) = tlevrad(2)
       tmid(2) = tmid(1)
-    
+
       DO l=1,L_NLAYRAD-1
          tmid(2*l+1) = tlevrad(2*l+1)
@@ -829 +801 @@
 !!Alternative interpolation:
 !         pmid(3) = pmid(1)
-!         pmid(4) = pmid(1) 
+!         pmid(4) = pmid(1)
 !         tmid(3) = tmid(1)
 !         tmid(4) = tmid(1)
@@ -870 +842 @@
             else
               print*,'***********************************************'
-              print*,'we allow model to continue with tlevrad<tgasmin' 
+              print*,'we allow model to continue with tlevrad<tgasmin'
               print*,'  ... we assume we know what you are doing ... '
               print*,'  ... but do not let this happen too often ... '
@@ -886 +858 @@
             else
               print*,'***********************************************'
-              print*,'we allow model to continue with tlevrad>tgasmax'  
+              print*,'we allow model to continue with tlevrad>tgasmax'
               print*,'  ... we assume we know what you are doing ... '
               print*,'  ... but do not let this happen too often ... '
@@ -1000 +972 @@
 
          ! Equivalent Albedo Calculation (for OUTPUT). MT2015
-         if(fract(ig) .ge. 1.0e-4) then ! equivalent albedo makes sense only during daylight.       
-            surface_stellar_flux=sum(nfluxgndv_nu(1:L_NSPECTV))      
+         if(fract(ig) .ge. 1.0e-4) then ! equivalent albedo makes sense only during daylight.
+            surface_stellar_flux=sum(nfluxgndv_nu(1:L_NSPECTV))
             if(surface_stellar_flux .gt. 1.0e-3) then ! equivalent albedo makes sense only if the stellar flux received by the surface is positive.
-               DO nw=1,L_NSPECTV                  
+               DO nw=1,L_NSPECTV
                   albedo_temp(nw)=albedo(ig,nw)*nfluxgndv_nu(nw)
                ENDDO
@@ -1025 +997 @@
 
          call sfluxi(plevrad,tlevrad,dtaui,taucumi,ubari,albi,      &
-              wnoi,dwni,cosbi,wbari,nfluxtopi,nfluxtopi_nu,         & 
+              wnoi,dwni,cosbi,wbari,nfluxtopi,nfluxtopi_nu,         &
               fmneti,fluxupi,fluxdni,fluxupi_nu,fzeroi,taugsurfi)
 
@@ -1033 +1005 @@
 
 !     Flux incident at the top of the atmosphere
-         fluxtop_dn(ig)=fluxtopvdn 
+         fluxtop_dn(ig)=fluxtopvdn
 
          fluxtop_lw(ig)  = real(nfluxtopi)
@@ -1039 +1011 @@
          fluxsurf_lw(ig) = real(fluxdni(L_NLAYRAD))
          fluxsurf_sw(ig) = real(fluxdnv(L_NLAYRAD))
-         
-!        Flux absorbed by the surface. By MT2015.         
+
+!        Flux absorbed by the surface. By MT2015.
          fluxsurfabs_sw(ig) = fluxsurf_sw(ig)*(1.-albedo_equivalent(ig))
 
@@ -1056 +1028 @@
 !     Spectral output, for exoplanet observational comparison
          if(specOLR)then
-            do nw=1,L_NSPECTI 
+            do nw=1,L_NSPECTI
                OLR_nu(ig,nw)=nfluxtopi_nu(nw)/DWNI(nw) !JL Normalize to the bandwidth
             end do
-            do nw=1,L_NSPECTV 
+            do nw=1,L_NSPECTV
                GSR_nu(ig,nw)=nfluxgndv_nu(nw)/DWNV(nw)
                OSR_nu(ig,nw)=nfluxoutv_nu(nw)/DWNV(nw) !JL Normalize to the bandwidth
@@ -1072 +1044 @@
             dtlw(ig,L_NLAYRAD+1-l)=(fmneti(l)-fmneti(l-1))  &
                 *glat(ig)/(cpp*scalep*(plevrad(2*l+1)-plevrad(2*l-1)))
-         END DO      
+         END DO
 
 !     These are values at top of atmosphere
@@ -1098 +1070 @@
           enddo
         enddo
-      endif        
-
-
-!-----------------------------------------------------------------------    
+      endif
+
+
+!-----------------------------------------------------------------------
       end do ! End of big loop over every GCM column.
 !-----------------------------------------------------------------------
@@ -1118 +1090 @@
          open(116,file='surf_vals.out')
          write(116,*) tsurf(1),pplev(1,1),fluxtop_dn(1),         &
-                      real(-nfluxtopv),real(nfluxtopi) 
+                      real(-nfluxtopv),real(nfluxtopi)
          close(116)
 
@@ -1147 +1119 @@
                write(118,*) plevrad(2*l)
                do nw=1,L_NSPECTI
-                  write(119,*) fluxupi_nu(l,nw) 
+                  write(119,*) fluxupi_nu(l,nw)
                enddo
-            enddo 
+            enddo
             close(118)
             close(119)

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

@@ -1 +1 @@
 MODULE callkeys_mod
-IMPLICIT NONE  
+IMPLICIT NONE
 
       logical,save :: callrad,corrk,calldifv,UseTurbDiff
@@ -10 +10 @@
       logical,save :: callgasvis,continuum,graybody
 !$OMP THREADPRIVATE(callgasvis,continuum,graybody)
-      logical,save :: strictboundcorrk                                     
+      logical,save :: strictboundcorrk
 !$OMP THREADPRIVATE(strictboundcorrk)
-      logical,save :: strictboundcia                                     
+      logical,save :: strictboundcia
 !$OMP THREADPRIVATE(strictboundcia)
 
@@ -57 +57 @@
 
 !! Pluto-specific variables
-      logical,save :: haze_radproffix,haze_proffix 
-!$OMP THREADPRIVATE(haze_radproffix,haze_proffix)      
-      logical,save :: haze,fasthaze,changeti,changetid,aerohaze         
-!$OMP THREADPRIVATE(haze,fasthaze,changeti,changetid,aerohaze)      
-      logical,save :: fast,metcloud,monoxcloud,glaflow,triton,paleo,nlte,strobel
-!$OMP THREADPRIVATE(fast,metcloud,monoxcloud,glaflow,triton,paleo,nlte,strobel      
-      logical,save :: kbo,nbsub,cooling,mode_n2,thres_n2ice,thres_coice 
-!$OMP THREADPRIVATE(kbo,nbsub,cooling,mode_n2,thres_n2ice,thres_coice)      
-      logical,save :: deltab,nb_monomer                                 
-!$OMP THREADPRIVATE(deltab,nb_monomer)      
-      logical,save :: metlateq,tholateq,metls1,metls2                   
-!$OMP THREADPRIVATE(metlateq,tholateq,metls1,metls2)      
-      logical,save :: mode_ch4,mode_tholins                             
-!$OMP THREADPRIVATE(mode_ch4,mode_tholins)      
-      logical,save :: source_haze,ch4lag,latlag,vmrlag,kfix,hazeconservch4      
-!$OMP THREADPRIVATE(source_haze,ch4lag,latlag,vmrlag,kfix,hazeconservch4)      
-      logical,save :: fracsource,latsource,lonsource                    
-!$OMP THREADPRIVATE(fracsource,latsource,lonsource)      
-      logical,save :: spelon1,spelon2,spelat1,spelat2,specalb           
-!$OMP THREADPRIVATE(spelon1,spelon2,spelat1,spelat2,specalb)      
-      logical,save :: assymflux,mode_hs,tsurfmax,albmin_ch4
-!$OMP THREADPRIVATE(assymflux,mode_hs,tsurfmax,albmin_ch4)      
-      logical,save :: feedback_met,thres_ch4ice,fdch4_latn,fdch4_lats   
-!$OMP THREADPRIVATE(feedback_met,thres_ch4ice,fdch4_latn,fdch4_lats)      
-      logical,save :: globmean1d,kmix_proffix,rad_haze          
-!$OMP THREADPRIVATE(globmean1d,kmix_proffix,rad_haze)      
-      logical,save :: fdch4_lone,fdch4_lonw,fdch4_maxalb,fdch4_depalb,fdch4_finalb
-!$OMP THREADPRIVATE(fdch4_lone,fdch4_lonw,fdch4_maxalb,fdch4_depalb,fdch4_finalb      
-      logical,save :: tholatn,tholats,tholone,tholonw                   
-!$OMP THREADPRIVATE(tholatn,tholats,tholone,tholonw)      
-      logical,save :: fdch4_ampl,fdch4_maxice,condmetsurf,condcosurf,vertdiff   
-!$OMP THREADPRIVATE(fdch4_ampl,fdch4_maxice,condmetsurf,condcosurf,vertdiff)      
-      logical,save :: convergeps,conservn2,patchflux,condensn2,no_n2frost
-!$OMP THREADPRIVATE(convergeps,conservn2,patchflux,condensn2,no_n2frost      
-
+      logical,save :: methane,carbox
+      !$OMP THREADPRIVATE(methane,carbox)
+      logical,save :: haze,haze_proffix,haze_radproffix
+!$OMP THREADPRIVATE(haze,haze_proffix,haze_radproffix)
+      logical,save :: fasthaze,changeti,changetid,aerohaze,fractal
+!$OMP THREADPRIVATE(fasthaze,changeti,changetid,aerohaze,fractal)
+      logical,save :: fast,metcloud,monoxcloud,glaflow,triton,paleo
+!$OMP THREADPRIVATE(fast,metcloud,monoxcloud,glaflow,triton,paleo)
+      logical,save :: nlte,strobel
+!$OMP THREADPRIVATE(nlte,strobel)
+      logical,save :: kbo
+!$OMP THREADPRIVATE(kbo)
+      logical,save :: cooling
+!$OMP THREADPRIVATE(cooling)
+      logical,save :: source_haze,hazeconservch4
+!$OMP THREADPRIVATE(source_haze,hazeconservch4)
+      logical,save :: ch4lag,tsurfmax
+!$OMP THREADPRIVATE(ch4lag,tsurfmax)
+      logical,save :: specalb
+!$OMP THREADPRIVATE(specalb)
+      logical,save :: assymflux
+!$OMP THREADPRIVATE(assymflux)
+      logical,save :: condmetsurf,condcosurf,vertdiff
+!$OMP THREADPRIVATE(condmetsurf,condcosurf,vertdiff)
+      logical,save :: convergeps,conservn2,condensn2,no_n2frost
+!$OMP THREADPRIVATE(convergeps,conservn2,condensn2,no_n2frost)
+      logical,save :: globmean1d,kmix_proffix
+!$OMP THREADPRIVATE(globmean1d,kmix_proffix)
+      integer,save :: nbsub
+!$OMP THREADPRIVATE(nbsub)
+      integer,save :: mode_n2
+!$OMP THREADPRIVATE(mode_n2)
+      integer,save :: mode_ch4
+!$OMP THREADPRIVATE(mode_ch4)
+      integer,save :: mode_tholins
+!$OMP THREADPRIVATE(mode_tholins)
+      integer,save :: kfix
+!$OMP THREADPRIVATE(kfix)
+      integer,save :: mode_hs
+!$OMP THREADPRIVATE(mode_hs)
+      integer,save :: feedback_met
+!$OMP THREADPRIVATE(feedback_met)
+      integer,save :: patchflux
+!$OMP THREADPRIVATE(patchflux)
+      integer,save :: nb_monomer
+!$OMP THREADPRIVATE(nb_monomer)
+      real,save    :: Nmix_co
+!$OMP THREADPRIVATE(Nmix_co)
+      real,save    :: Nmix_ch4
+!$OMP THREADPRIVATE(Nmix_ch4)
+      real,save    :: tau_n2
+!$OMP THREADPRIVATE(tau_n2)
+      real,save    :: tau_ch4
+!$OMP THREADPRIVATE(tau_ch4)
+      real,save    :: tau_co
+!$OMP THREADPRIVATE(tau_co)
+      real,save    :: tau_prechaze
+!$OMP THREADPRIVATE(tau_prechaze)
+      real,save    :: paleoyears
+!$OMP THREADPRIVATE(paleoyears)
+      real,save    :: dayfrac
+!$OMP THREADPRIVATE(dayfrac)
+      real,save    :: thresh_non2
+!$OMP THREADPRIVATE(thresh_non2)
+      real,save    :: vmrch4fix
+!$OMP THREADPRIVATE(vmrch4fix)
+      real,save    :: fluxgeo
+!$OMP THREADPRIVATE(fluxgeo)
+      real,save    :: fluxgeo2
+!$OMP THREADPRIVATE(fluxgeo2)
+      real,save    :: deltab
+!$OMP THREADPRIVATE(deltab)
+      real,save    :: metlateq
+!$OMP THREADPRIVATE(metlateq)
+      real,save    :: metls1
+!$OMP THREADPRIVATE(metls1)
+      real,save    :: metls2
+!$OMP THREADPRIVATE(metls2)
+      real,save    :: tholateq
+!$OMP THREADPRIVATE(tholateq)
+      real,save    :: tholatn
+!$OMP THREADPRIVATE(tholatn)
+      real,save    :: tholats
+!$OMP THREADPRIVATE(tholats)
+      real,save    :: tholone
+!$OMP THREADPRIVATE(tholone)
+      real,save    :: tholonw
+!$OMP THREADPRIVATE(tholonw)
+      real,save    :: spelon1
+!$OMP THREADPRIVATE(spelon1)
+      real,save    :: spelat2
+!$OMP THREADPRIVATE(spelat2)
+      real,save    :: spelat1
+!$OMP THREADPRIVATE(spelat1)
+      real,save    :: spelon2
+!$OMP THREADPRIVATE(spelon2)
+      real,save    :: latlag
+!$OMP THREADPRIVATE(latlag)
+      real,save    :: vmrlag
+!$OMP THREADPRIVATE(vmrlag)
+      real,save    :: albmin_ch4
+!$OMP THREADPRIVATE(albmin_ch4)
+      real,save    :: fracsource
+!$OMP THREADPRIVATE(fracsource)
+      real,save    :: latsource
+!$OMP THREADPRIVATE(latsource)
+      real,save    :: lonsource
+!$OMP THREADPRIVATE(lonsource)
+      real,save    :: thres_ch4ice
+!$OMP THREADPRIVATE(thres_ch4ice)
+      real,save    :: thres_n2ice
+!$OMP THREADPRIVATE(thres_n2ice)
+      real,save    :: thres_coice
+!$OMP THREADPRIVATE(thres_coice)
+      real,save    :: fdch4_latn
+!$OMP THREADPRIVATE(fdch4_latn)
+      real,save    :: fdch4_lats
+!$OMP THREADPRIVATE(fdch4_lats)
+      real,save    :: fdch4_lone
+!$OMP THREADPRIVATE(fdch4_lone)
+      real,save    :: fdch4_lonw
+!$OMP THREADPRIVATE(fdch4_lonw)
+      real,save    :: fdch4_maxice
+!$OMP THREADPRIVATE(fdch4_maxice)
+      real,save    :: fdch4_maxalb
+!$OMP THREADPRIVATE(fdch4_maxalb)
+      real,save    :: fdch4_ampl
+!$OMP THREADPRIVATE(fdch4_ampl)
+      real,save    :: fdch4_depalb
+!$OMP THREADPRIVATE(fdch4_depalb)
+      real,save    :: fdch4_finalb
+!$OMP THREADPRIVATE(fdch4_finalb)
+      real,save    :: rad_haze
+!$OMP THREADPRIVATE(rad_haze)
 
       logical,save :: global1d
@@ -169 +268 @@
       real,save :: kmixmin
 !$OMP THREADPRIVATE(kmixmin)
-      
+
       logical,save :: iscallphys=.false.!existence of callphys.def
 !$OMP THREADPRIVATE(iscallphys)
 
       ! do we read a startphy.nc file (default=.true.)
-      logical,save :: startphy_file=.true. 
+      logical,save :: startphy_file=.true.
 !$OMP THREADPRIVATE(startphy_file)
 

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

@@ -6 +6 @@
       ! use radii_mod, only: h2o_reffrad
       ! use aerosol_mod, only : iaero_h2o
-      USE tracer_h, only : igcm_n2_ice,radius,rho_q
+      USE tracer_h, only : igcm_haze,radius,rho_q
       use comcstfi_mod, only: g
 
@@ -12 +12 @@
 
 !==================================================================
-!     
+!
 !     Purpose
 !     -------
 !     Calculates sedimentation of aerosols depending on their
 !     density and radius.
-!     
+!
 !     Authors
 !     -------
 !     F. Forget (1999)
 !     Tracer generalisation by E. Millour (2009)
-!     
+!
 !==================================================================
 
@@ -43 +43 @@
       real,intent(in) :: pdqfi(ngrid,nlay,nq)  ! tendency on tracers before
                                                ! sedimentation (kg/kg.s-1)
-      
+
       real,intent(out) :: pdqsed(ngrid,nlay,nq) ! tendency due to sedimentation (kg/kg.s-1)
       real,intent(out) :: pdqs_sed(ngrid,nq)    ! flux at surface (kg.m-2.s-1)
@@ -53 +53 @@
 
       ! for particles with varying radii:
-      real,allocatable,save :: reffrad(:,:,:) ! particle radius (m) 
+      real,allocatable,save :: reffrad(:,:,:) ! particle radius (m)
       real,allocatable,save :: nueffrad(:,:,:) ! aerosol effective radius variance
 !$OMP THREADPRIVATE(reffrad,nueffrad)
@@ -80 +80 @@
         allocate(nueffrad(ngrid,nlay,naerkind))
       ENDIF ! of IF (firstcall)
-      
+
 !=======================================================================
 !     Preliminary calculation of the layer characteristics
@@ -87 +87 @@
       do  l=1,nlay
         do ig=1, ngrid
-          masse(ig,l)=(pplev(ig,l) - pplev(ig,l+1)) /g 
+          masse(ig,l)=(pplev(ig,l) - pplev(ig,l+1)) /g
           epaisseur(ig,l)= zlev(ig,l+1) - zlev(ig,l)
           zt(ig,l)=pt(ig,l)+pdt(ig,l)*ptimestep
@@ -97 +97 @@
 ! [This has been rearranged by L. Kerber to allow the sedimentation
 !  of general tracers.]
- 
+
       do iq=1,nq
-       if((radius(iq).gt.1.e-9).and.(iq.ne.igcm_n2_ice)) then ! JVO 08/2017 : be careful radius was tested uninitialized (fixed) ... 
-       
-!         (no sedim for gases, and n2_ice sedim is done in condense_n2)      
+       if((radius(iq).gt.1.e-9).and.(iq.ne.igcm_haze)) then ! JVO 08/2017 : be careful radius was tested uninitialized (fixed) ...
+
+!         (no sedim for gases, and n2_ice sedim is done in condense_n2)
 
 ! store locally updated tracers
 
-          do l=1,nlay 
+          do l=1,nlay
             do ig=1, ngrid
               zqi(ig,l,iq)=pq(ig,l,iq)+pdqfi(ig,l,iq)*ptimestep
             enddo
           enddo ! of do l=1,nlay
-        
+
 !======================================================================
-! Sedimentation 
+! Sedimentation
 !======================================================================
 ! Water      !AF24: deleted
 
 ! ! General Case
-!        else 
+!        else
              call newsedim(ngrid,nlay,1,ptimestep,
      &            pplev,masse,epaisseur,zt,radius(iq),rho_q(iq),

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

@@ -1 @@
-subroutine condens_n2(klon,nlayer,nq,ptimestep, &
+subroutine condense_n2(klon,klev,nq,ptimestep, &
       pcapcal,pplay,pplev,ptsrf,pt, &
       pphi,pdt,pdu,pdv,pdtsrf,pu,pv,pq,pdq, &
@@ -8 +8 @@
   use radinc_h, only : naerkind
   use comgeomfi_h
-  use comcstfi_mod, only: g, r, cpp
-  USE surfdat_h, only: phisfi
-  USE tracer_h, only: noms, igcm_n2
-  USE callkeys_mod, only: fast,ch4lag,latlag,lw_n2,nbsub,no_n2frost,tsurfmax,kmixmin,source_haze,vmrlag
-  USE dimphy, only : klon,klev
+  use comcstfi_mod, only: g, r, cpp, pi
+  USE surfdat_h, only: phisfi,kp,p00
+  USE tracer_h, only: noms, igcm_n2, lw_n2
+  USE callkeys_mod, only: fast,ch4lag,latlag,nbsub,no_n2frost,tsurfmax,kmixmin,source_haze,vmrlag
+  USE comvert_mod, ONLY: ap,bp
+  use geometry_mod, only: latitude
 
 
@@ -22 +23 @@
 !     Condense and/or sublime N2 ice on the ground and in the
 !     atmosphere, and sediment the ice.
-!  
+!
 !     Inputs
-!     ------ 
+!     ------
 !     klon                 Number of vertical columns
-!     nlayer                Number of layers
-!     pplay(klon,nlayer)   Pressure layers
-!     pplev(klon,nlayer+1) Pressure levels 
-!     pt(klon,nlayer)      Temperature (in K)
+!     klev                Number of layers
+!     pplay(klon,klev)   Pressure layers
+!     pplev(klon,klev+1) Pressure levels
+!     pt(klon,klev)      Temperature (in K)
 !     ptsrf(klon)          Surface temperature
-!     
-!     pdt(klon,nlayermx)   Time derivative before condensation/sublimation of pt
+!
+!     pdt(klon,klev)   Time derivative before condensation/sublimation of pt
 !     pdtsrf(klon)         Time derivative before condensation/sublimation of ptsrf
 !     picen2(klon)         n2 ice at the surface (kg/m2)
-!     
+!
 !     Outputs
 !     -------
 !     pdpsrf(klon)         \  Contribution of condensation/sublimation
-!     pdtc(klon,nlayermx)  /  to the time derivatives of Ps, pt, and ptsrf
+!     pdtc(klon,klev)  /  to the time derivatives of Ps, pt, and ptsrf
 !     pdtsrfc(klon)       /
 !     pdicen2(klon)         Tendancy n2 ice at the surface (kg/m2)
-!     
+!
 !     Both
 !     ----
@@ -49 +50 @@
 !
 !     Authors
-!     ------- 
+!     -------
 !     Francois Forget (1996,2013)
 !     Converted to Fortran 90 and slightly modified by R. Wordsworth (2009)
-!     
-!     
+!
+!
 !==================================================================
 
@@ -59 +60 @@
 !     Arguments
 
-  INTEGER klon, nlayer, nq
-
-  REAL ptimestep 
-  REAL pplay(klon,nlayer),pplev(klon,nlayer+1)
+  INTEGER klon, klev, nq
+
+  REAL ptimestep
+  REAL pplay(klon,klev),pplev(klon,klev+1)
   REAL pcapcal(klon)
-  REAL pt(klon,nlayer)
+  REAL pt(klon,klev)
   REAL ptsrf(klon),flu1(klon),flu2(klon),flu3(klon)
-  REAL pphi(klon,nlayer)
-  REAL pdt(klon,nlayer),pdtsrf(klon),pdtc(klon,nlayer)
+  REAL pphi(klon,klev)
+  REAL pdt(klon,klev),pdtsrf(klon),pdtc(klon,klev)
   REAL pdtsrfc(klon),pdpsrf(klon)
   REAL picen2(klon),psolaralb(klon),pemisurf(klon)
- 
-
-  REAL pu(klon,nlayer) ,  pv(klon,nlayer)
-  REAL pdu(klon,nlayer) , pdv(klon,nlayer)
-  REAL pduc(klon,nlayer) , pdvc(klon,nlayer)
-  REAL pq(ngridmx,nlayer,nq),pdq(klon,nlayer,nq)
-  REAL pdqc(klon,nlayer,nq)
+
+
+  REAL pu(klon,klev) ,  pv(klon,klev)
+  REAL pdu(klon,klev) , pdv(klon,klev)
+  REAL pduc(klon,klev) , pdvc(klon,klev)
+  REAL pq(klon,klev,nq),pdq(klon,klev,nq)
+  REAL pdqc(klon,klev,nq)
 
 !-----------------------------------------------------------------------
@@ -83 +84 @@
   INTEGER l,ig,ilay,it,iq,ich4_gas
 
-  REAL*8 zt(ngridmx,nlayermx)
+  REAL*8 zt(klon,klev)
   REAL tcond_n2
   REAL pcond_n2
   REAL glob_average2d         ! function
-  REAL zqn2(ngridmx,nlayermx) ! N2 MMR used to compute Tcond/zqn2
-  REAL ztcond (ngridmx,nlayermx)
-  REAL ztcondsol(ngridmx),zfallheat
-  REAL pdicen2(ngridmx)
-  REAL zcondicea(ngridmx,nlayermx), zcondices(ngridmx)
-  REAL zfallice(ngridmx,nlayermx+1)
-  REAL zmflux(nlayermx+1)
-  REAL zu(nlayermx),zv(nlayermx)
-  REAL zq(nlayermx,nqmx),zq1(nlayermx)
-  REAL ztsrf(ngridmx)
-  REAL ztc(nlayermx), ztm(nlayermx+1) 
-  REAL zum(nlayermx+1) , zvm(nlayermx+1)
-  REAL zqm(nlayermx+1,nqmx),zqm1(nlayermx+1)
-  LOGICAL condsub(ngridmx)
+  REAL zqn2(klon,klev) ! N2 MMR used to compute Tcond/zqn2
+  REAL ztcond (klon,klev)
+  REAL ztcondsol(klon),zfallheat
+  REAL pdicen2(klon)
+  REAL zcondicea(klon,klev), zcondices(klon)
+  REAL zfallice(klon,klev+1)
+  REAL zmflux(klev+1)
+  REAL zu(klev),zv(klev)
+  REAL zq(klev,nq),zq1(klev)
+  REAL ztsrf(klon)
+  REAL ztc(klev), ztm(klev+1)
+  REAL zum(klev+1) , zvm(klev+1)
+  REAL zqm(klev+1,nq),zqm1(klev+1)
+  LOGICAL condsub(klon)
   REAL subptimestep
   Integer Ntime
-  real masse (nlayermx),w(nlayermx+1)
-  real wq(ngridmx,nlayermx+1)
+  real masse (klev),w(klev+1)
+  real wq(klon,klev+1)
   real vstokes,reff
   real dWtotsn2
-  real condnconsn2(ngridmx)
+  real condnconsn2(klon)
   real nconsMAXn2
 !     Special diagnostic variables
-  real tconda1(ngridmx,nlayermx)
-  real tconda2(ngridmx,nlayermx)
-  real zdtsig (ngridmx,nlayermx)
-  real zdtlatent (ngridmx,nlayermx)
-  real zdt (ngridmx,nlayermx)
-  REAL albediceF(ngridmx)
-  SAVE albediceF
+  real tconda1(klon,klev)
+  real tconda2(klon,klev)
+  real zdtsig (klon,klev)
+  real zdtlatent (klon,klev)
+  real zdt (klon,klev)
+  REAL albediceF(klon)
+!   SAVE albediceF
   INTEGER nsubtimestep,itsub    !number of subtimestep when calling vl1d
   REAL subtimestep        !ptimestep/nsubtimestep
-  REAL dtmax              
-
-  REAL zplevhist(ngridmx)
-  REAL zplevnew(ngridmx)
-  REAL zplev(ngridmx)
-  REAL zpicen2(ngridmx)
-  REAL ztsrfhist(ngridmx)
-  REAL zdtsrf(ngridmx)
+  REAL dtmax
+
+  REAL zplevhist(klon)
+  REAL zplevnew(klon)
+  REAL zplev(klon)
+  REAL zpicen2(klon)
+  REAL ztsrfhist(klon)
+  REAL zdtsrf(klon)
   REAL globzplevnew
 
-  REAL vmrn2(ngridmx)
-  SAVE vmrn2
-  REAL stephan   
+  REAL vmrn2(klon)
+!   SAVE vmrn2
+  REAL stephan
   DATA stephan/5.67e-08/  ! Stephan Boltzman constant
   SAVE stephan
@@ -167 +168 @@
   IF (firstcall) THEN
      ccond=cpp/(g*lw_n2)
-     print*,'In condens_n2cloud: ccond=',ccond,' latcond=',lw_n2
+     print*,'In condense_n2cloud: ccond=',ccond,' latcond=',lw_n2
 
      ! calculate global mean surface pressure for the fast mode
      IF (fast) THEN
-        DO ig=1,ngridmx
+        DO ig=1,klon
            kp(ig) = exp(-phisfi(ig)/(r*38.))
         ENDDO
@@ -179 +180 @@
      vmrn2(:) = 1.
      IF (ch4lag) then
-        DO ig=1,ngridmx
-           if (lati(ig)*180./pi.ge.latlag) then
+        DO ig=1,klon
+           if (latitude(ig)*180./pi.ge.latlag) then
               vmrn2(ig) = vmrlag
            endif
         ENDDO
-     ENDIF   
+     ENDIF
      IF (no_n2frost) then
-        DO ig=1,ngridmx
+        DO ig=1,klon
            if (picen2(ig).eq.0.) then
               vmrn2(ig) = 1.e-15
            endif
         ENDDO
-     ENDIF   
+     ENDIF
      firstcall=.false.
   ENDIF
 
 !-----------------------------------------------------------------------
-!     Calculation of n2 condensation / sublimation 
+!     Calculation of n2 condensation / sublimation
 
 !     Variables used:
 !     picen2(klon)            : amount of n2 ice on the ground     (kg/m2)
-!     zcondicea(klon,nlayermx): condensation rate in layer l       (kg/m2/s)
+!     zcondicea(klon,klev): condensation rate in layer l       (kg/m2/s)
 !     zcondices(klon)         : condensation rate on the ground    (kg/m2/s)
-!     zfallice(klon,nlayermx) : amount of ice falling from layer l (kg/m2/s)
-!     zdtlatent(klon,nlayermx): dT/dt due to phase changes         (K/s)
+!     zfallice(klon,klev) : amount of ice falling from layer l (kg/m2/s)
+!     zdtlatent(klon,klev): dT/dt due to phase changes         (K/s)
 
 !     Tendencies initially set to 0
@@ -213 +214 @@
   pdicen2(1:klon) = 0.
   zfallheat=0
-  pdqc(1:klon,1:nlayer,1:nq)=0
-  pdtc(1:klon,1:nlayer)=0
-  pduc(1:klon,1:nlayer)=0
-  pdvc(1:klon,1:nlayer)=0
-  zfallice(1:klon,1:nlayer+1)=0
-  zcondicea(1:klon,1:nlayer)=0
-  zdtlatent(1:klon,1:nlayer)=0
-  zt(1:klon,1:nlayer)=0.
+  pdqc(1:klon,1:klev,1:nq)=0
+  pdtc(1:klon,1:klev)=0
+  pduc(1:klon,1:klev)=0
+  pdvc(1:klon,1:klev)=0
+  zfallice(1:klon,1:klev+1)=0
+  zcondicea(1:klon,1:klev)=0
+  zdtlatent(1:klon,1:klev)=0
+  zt(1:klon,1:klev)=0.
 
 !-----------------------------------------------------------------------
@@ -229 +230 @@
 !      (calcul of zcondicea , zfallice and pdtc)
 
-  zt(1:klon,1:nlayer)=pt(1:klon,1:nlayer)+ pdt(1:klon,1:nlayer)*ptimestep
+  zt(1:klon,1:klev)=pt(1:klon,1:klev)+ pdt(1:klon,1:klev)*ptimestep
   if (igcm_n2.ne.0) then
-     zqn2(1:klon,1:nlayer) = 1. ! & temporaire
-!        zqn2(1:klon,1:nlayer)= pq(1:klon,1:nlayer,igcm_n2) + pdq(1:klon,1:nlayer,igcm_n2)*ptimestep
+     zqn2(1:klon,1:klev) = 1. ! & temporaire
+!        zqn2(1:klon,1:klev)= pq(1:klon,1:klev,igcm_n2) + pdq(1:klon,1:klev,igcm_n2)*ptimestep
   else
-     zqn2(1:klon,1:nlayer) = 1.
+     zqn2(1:klon,1:klev) = 1.
   end if
-  
+
   if (.not.fast) then
 !     Forecast the atmospheric frost temperature 'ztcond' with function tcond_n2
-    DO l=1,nlayer
+    DO l=1,klev
       DO ig=1,klon
-          ztcond (ig,l) = tcond_n2(pplay(ig,l),zqn2(ig,l)) 
+          ztcond (ig,l) = tcond_n2(pplay(ig,l),zqn2(ig,l))
       ENDDO
     ENDDO
 
-    DO l=nlayer,1,-1
+    DO l=klev,1,-1
       DO ig=1,klon
        pdtc(ig,l)=0.  ! final tendancy temperature set to 0
@@ -251 +252 @@
        IF((zt(ig,l).LT.ztcond(ig,l)).or.(zfallice(ig,l+1).gt.0))THEN
            condsub(ig)=.true.  !condensation at level l
-           IF (zfallice(ig,l+1).gt.0) then 
-             zfallheat=zfallice(ig,l+1)*&  
+           IF (zfallice(ig,l+1).gt.0) then
+             zfallheat=zfallice(ig,l+1)*&
             (pphi(ig,l+1)-pphi(ig,l) +&
            cpice*(ztcond(ig,l+1)-ztcond(ig,l)))/lw_n2
@@ -266 +267 @@
                 .AND. (zfallice(ig,l+1).gt.0)) THEN
 
-              zdtlatent(ig,l)=(-zfallice(ig,l+1)+zfallheat)/&     
+              zdtlatent(ig,l)=(-zfallice(ig,l+1)+zfallheat)/&
                (ccond*(pplev(ig,l)-pplev(ig,l+1)))
               zcondicea(ig,l)= -zfallice(ig,l+1)
@@ -272 +273 @@
 
            zfallice(ig,l) = zcondicea(ig,l)+zfallice(ig,l+1)
-      
+
         END IF
-     
+
       ENDDO
     ENDDO
@@ -285 +286 @@
 !     Adding subtimesteps : in fast version, pressures too low lead to
 !     instabilities.
-  IF (fast) THEN 
+  IF (fast) THEN
      IF (pplev(1,1).gt.0.3) THEN
-         nsubtimestep= 1  
+         nsubtimestep= 1
      ELSE
-         nsubtimestep= nbsub !max(nint(ptimestep/dtmax),1) 
+         nsubtimestep= nbsub !max(nint(ptimestep/dtmax),1)
      ENDIF
   ELSE
@@ -307 +308 @@
        ENDDO
      ELSE
-     ! additional loop : 
+     ! additional loop :
      ! 1)  pressure updated
      ! 2)  direct redistribution of pressure over the globe
@@ -339 +340 @@
        ENDDO
      ENDIF  ! (itsub=1)
-   
+
    DO ig=1,klon
      ! forecast of frost temperature ztcondsol
@@ -353 +354 @@
 !     Condensation or partial sublimation of N2 ice
         if (ztsrf(ig) .LT. ztcondsol(ig)) then  ! condensation
-!             Include a correction to account for the cooling of air near 
+!             Include a correction to account for the cooling of air near
 !             the surface before condensing:
          zcondices(ig)=pcapcal(ig)*(ztcondsol(ig)-ztsrf(ig)) &
@@ -359 +360 @@
         else                                    ! sublimation
           zcondices(ig)=pcapcal(ig)*(ztcondsol(ig)-ztsrf(ig)) &
-          /(lw_n2*subtimestep) 
+          /(lw_n2*subtimestep)
         end if
 
@@ -369 +370 @@
         IF((zpicen2(ig)/subtimestep).LE. &
             -zcondices(ig))THEN
-           zcondices(ig) = -zpicen2(ig)/subtimestep 
+           zcondices(ig) = -zpicen2(ig)/subtimestep
            pdtsrfc(ig)=(lw_n2/pcapcal(ig))*       &
                (zcondices(ig))
@@ -376 +377 @@
 !     Changing N2 ice amount and pressure
 
-        pdicen2(ig) = zcondices(ig) 
+        pdicen2(ig) = zcondices(ig)
         pdpsrf(ig)   = -pdicen2(ig)*g
     !    pdpsrf(ig)   = 0. ! OPTION to check impact N2 sub/cond
@@ -411 +412 @@
      if(.not.picen2(ig).ge.0.) THEN
 !           if(picen2(ig) + pdicen2(ig)*ptimestep.le.-1.e-8) then
-          print*, 'WARNING NEG RESERVOIR in condense_n2: picen2(',ig,')=', picen2(ig) + pdicen2(ig)*ptimestep  
+          print*, 'WARNING NEG RESERVOIR in condense_n2: picen2(',ig,')=', picen2(ig) + pdicen2(ig)*ptimestep
 !              pdicen2(ig)= -picen2(ig)/ptimestep
 !           else
@@ -420 +421 @@
 
 ! ***************************************************************
-!  Correction to account for redistribution between sigma or hybrid 
+!  Correction to account for redistribution between sigma or hybrid
 !  layers when changing surface pressure (and warming/cooling
 !  of the n2 currently changing phase).
@@ -431 +432 @@
 !  """"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
         zmflux(1) = -zcondices(ig)
-        DO l=1,nlayer
+        DO l=1,klev
          zmflux(l+1) = zmflux(l) -zcondicea(ig,l)  &
-         + (bp(l)-bp(l+1))*(zfallice(ig,1)-zmflux(1))     
-! zmflux set to 0 if very low to avoid: top layer is disappearing in v1ld  
+         + (bp(l)-bp(l+1))*(zfallice(ig,1)-zmflux(1))
+! zmflux set to 0 if very low to avoid: top layer is disappearing in v1ld
       if (abs(zmflux(l+1)).lt.1E-13.OR.bp(l+1).eq.0.) zmflux(l+1)=0.
         END DO
 
 ! Mass of each layer
-! ------------------ 
-       DO l=1,nlayer
+! ------------------
+       DO l=1,klev
          masse(l)=(pplev(ig,l) - pplev(ig,l+1))/g
        END DO
@@ -447 +448 @@
 !  Corresponding fluxes for T,U,V,Q
 !  """"""""""""""""""""""""""""""""
-!           averaging operator for TRANSPORT  
+!           averaging operator for TRANSPORT
 !           """"""""""""""""""""""""""""""""
 
 !           Subtimestep loop to perform the redistribution gently and simultaneously with
-!           the other tendencies  
+!           the other tendencies
 !           Estimation of subtimestep (using only  the first layer, the most critical)
         dtmax=ptimestep
-        if (zmflux(1).gt.1.e-20) then 
+        if (zmflux(1).gt.1.e-20) then
             dtmax=min(dtmax,masse(1)*zqn2(ig,1)/abs(zmflux(1)))
         endif
-        nsubtimestep= max(nint(ptimestep/dtmax),nint(2.)) 
+        nsubtimestep= max(nint(ptimestep/dtmax),nint(2.))
         subtimestep=ptimestep/float(nsubtimestep)
 
 !          New flux for each subtimestep
-       do l=1,nlayer+1
+       do l=1,klev+1
             w(l)=-zmflux(l)*subtimestep
        enddo
 !          initializing variables that will vary during subtimestep:
-       do l=1,nlayer
-           ztc(l)  =pt(ig,l)   
-           zu(l)   =pu(ig,l) 
-           zv(l)   =pv(ig,l)  
-           do iq=1,nqmx
-              zq(l,iq) = pq(ig,l,iq) 
+       do l=1,klev
+           ztc(l)  =pt(ig,l)
+           zu(l)   =pu(ig,l)
+           zv(l)   =pv(ig,l)
+           do iq=1,nq
+              zq(l,iq) = pq(ig,l,iq)
            enddo
        end do
@@ -477 +478 @@
 !          """"""""""""""""""""""
        do itsub=1,nsubtimestep
-!             Progressively adding tendancies from other processes. 
-          do l=1,nlayer
+!             Progressively adding tendancies from other processes.
+          do l=1,klev
              ztc(l)  =ztc(l)   +(pdt(ig,l) + zdtlatent(ig,l))*subtimestep
              zu(l)   =zu(l)   +pdu( ig,l) * subtimestep
              zv(l)   =zv(l)   +pdv( ig,l) * subtimestep
-             do iq=1,nqmx
+             do iq=1,nq
                 zq(l,iq) = zq(l,iq) + pdq(ig,l,iq)* subtimestep
              enddo
@@ -488 +489 @@
 
 !             Change of mass in each layer
-          do l=1,nlayer
+          do l=1,klev
              masse(l)=masse(l)+pdpsrf(ig)*subtimestep*(pplev(ig,l) - pplev(ig,l+1))&
                                              /(g*pplev(ig,1))
           end do
 
-           ztm(2:nlayermx+1)=0.
-           zum(2:nlayermx+1)=0.
-           zvm(2:nlayermx+1)=0.
-           zqm1(1:nlayermx+1)=0.
+           ztm(2:klev+1)=0.
+           zum(2:klev+1)=0.
+           zvm(2:klev+1)=0.
+           zqm1(1:klev+1)=0.
 
 !             Van Leer scheme:
-          call vl1d(nlayer,ztc,2.,masse,w,ztm) 
-          call vl1d(nlayer,zu ,2.,masse,w,zum) 
-          call vl1d(nlayer,zv ,2.,masse,w,zvm) 
-         do iq=1,nqmx
-           do l=1,nlayer
+          call vl1d(klev,ztc,2.,masse,w,ztm)
+          call vl1d(klev,zu ,2.,masse,w,zum)
+          call vl1d(klev,zv ,2.,masse,w,zvm)
+         do iq=1,nq
+           do l=1,klev
             zq1(l)=zq(l,iq)
            enddo
            zqm1(1)=zqm(1,iq)
-           call vl1d(nlayer,zq1,2.,masse,w,zqm1)
-           do l=2,nlayer
+           call vl1d(klev,zq1,2.,masse,w,zqm1)
+           do l=2,klev
             zqm(l,iq)=zqm1(l)
            enddo
@@ -516 +517 @@
           if (igcm_n2.ne.0) then
             zqm(1,igcm_n2)=1.
-            do l=1,nlayer-1
-              if (w(l)*zqm(l,igcm_n2).gt.zq(l,igcm_n2)*masse(l)) then 
+            do l=1,klev-1
+              if (w(l)*zqm(l,igcm_n2).gt.zq(l,igcm_n2)*masse(l)) then
                  zqm(l+1,igcm_n2)=max(zqm(l+1,igcm_n2),    &
                  (zqm(l,igcm_n2)*w(l) -zq(l,igcm_n2)*masse(l))/w(l+1) )
-              else 
+              else
                 exit
               endif
@@ -528 +529 @@
 !             Value transfert at the surface interface when condensation sublimation:
 
-          if (zmflux(1).lt.0) then 
+          if (zmflux(1).lt.0) then
 !               Surface condensation
             zum(1)= zu(1)
-            zvm(1)= zv(1) 
+            zvm(1)= zv(1)
             ztm(1) = ztc(1)
-          else 
+          else
 !               Surface sublimation:
             ztm(1) = ztsrf(ig) + pdtsrfc(ig)*ptimestep
-            zum(1) = 0  
-            zvm(1) = 0  
+            zum(1) = 0
+            zvm(1) = 0
           end if
-          do iq=1,nqmx
+          do iq=1,nq
              zqm(1,iq)=0. ! most tracer do not condense !
           enddo
@@ -550 +551 @@
             DO iq=1,nq
              tname=noms(iq)
-             if (tname(1:4).eq."haze") then 
+             if (tname(1:4).eq."haze") then
                    !zqm(1,iq)=0.02
                    !zqm(1,iq)=-pdicen2(ig)*0.02
@@ -559 +560 @@
              endif
             ENDDO
-           ENDIF      
-          ENDIF      
-          ztm(nlayer+1)= ztc(nlayer) ! should not be used, but... 
-          zum(nlayer+1)= zu(nlayer)  ! should not be used, but...
-          zvm(nlayer+1)= zv(nlayer)  ! should not be used, but...
-          do iq=1,nqmx
-           zqm(nlayer+1,iq)= zq(nlayer,iq)
+           ENDIF
+          ENDIF
+          ztm(klev+1)= ztc(klev) ! should not be used, but...
+          zum(klev+1)= zu(klev)  ! should not be used, but...
+          zvm(klev+1)= zv(klev)  ! should not be used, but...
+          do iq=1,nq
+           zqm(klev+1,iq)= zq(klev,iq)
           enddo
 
-!             Tendencies on T, U, V, Q 
+!             Tendencies on T, U, V, Q
 !             """""""""""""""""""""""
-          DO l=1,nlayer
+          DO l=1,klev
 
 !               Tendencies on T
@@ -591 +592 @@
 
 !             Tendencies on Q
-          do iq=1,nqmx
+          do iq=1,nq
             if (iq.eq.igcm_n2) then
 !                 SPECIAL Case when the tracer IS N2 :
-              DO l=1,nlayer
+              DO l=1,klev
               pdqc(ig,l,iq)= (1/masse(l)) *        &
                ( zmflux(l)*(zqm(l,iq) - zq(l,iq))    &
@@ -601 +602 @@
               END DO
             else
-              DO l=1,nlayer
+              DO l=1,klev
                 pdqc(ig,l,iq)= (1/masse(l)) *         &
                ( zmflux(l)*(zqm(l,iq) - zq(l,iq))     &
@@ -610 +611 @@
           enddo
 !             Update variables at the end of each subtimestep.
-          do l=1,nlayer
+          do l=1,klev
             ztc(l)  =ztc(l)  + zdtsig(ig,l)  *subtimestep
             zu(l)   =zu(l)   + pduc(ig,l)  *subtimestep
             zv(l)   =zv(l)   + pdvc(ig,l)  *subtimestep
-            do iq=1,nqmx
+            do iq=1,nq
               zq(l,iq) = zq(l,iq) + pdqc(ig,l,iq)*subtimestep
             enddo
@@ -620 +621 @@
        enddo   ! loop on nsubtimestep
 !          Recomputing Total tendencies
-       do l=1,nlayer
+       do l=1,klev
          pdtc(ig,l)   = (ztc(l) - zt(ig,l) )/ptimestep
          pduc(ig,l)   = (zu(l) - (pu(ig,l) + pdu(ig,l)*ptimestep))/ptimestep
          pdvc(ig,l)   = (zv(l) - (pv(ig,l) + pdv(ig,l)*ptimestep))/ptimestep
-         do iq=1,nqmx
+         do iq=1,nq
            pdqc(ig,l,iq) = (zq(l,iq) - (pq(ig,l,iq) + pdq(ig,l,iq)*ptimestep))/ptimestep
 
@@ -632 +633 @@
              if((pq(ig,l,iq) +(pdqc(ig,l,iq)+ pdq(ig,l,iq))*ptimestep) &
                      .lt.0.01) then ! if n2 < 1 %  !
-               pdqc(ig,l,iq)=(0.01-pq(ig,l,iq))/ptimestep-pdq(ig,l,iq)  
+               pdqc(ig,l,iq)=(0.01-pq(ig,l,iq))/ptimestep-pdq(ig,l,iq)
              end if
             end if
@@ -639 +640 @@
        end do
 ! *******************************TEMPORAIRE ******************
-       if (klon.eq.1) then 
+       if (klon.eq.1) then
               write(*,*) 'nsubtimestep=' ,nsubtimestep
            write(*,*) 'masse avant' , (pplev(ig,1) - pplev(ig,2))/g
-              write(*,*) 'masse apres' , masse(1) 
+              write(*,*) 'masse apres' , masse(1)
               write(*,*) 'zmflux*DT, l=1' ,  zmflux(1)*ptimestep
               write(*,*) 'zmflux*DT, l=2' ,  zmflux(2)*ptimestep
@@ -653 +654 @@
 ! ***********************************************************
     end if ! if (condsub)
-   END DO  ! loop on ig 
+   END DO  ! loop on ig
   endif ! not fast
 
@@ -660 +661 @@
   if (olkin) then
   DO ig=1,klon
-     if (lati(ig).lt.0) then
+     if (latitude(ig).lt.0) then
        pdtsrfc(ig) = max(0.,pdtsrfc(ig))
        pdpsrf(ig) = 0.
        pdicen2(ig)  = 0.
-       do l=1,nlayer
+       do l=1,klev
          pdtc(ig,l)  =  max(0.,zdtlatent(ig,l))
          pduc(ig,l)  = 0.
          pdvc(ig,l)  = 0.
-         do iq=1,nqmx
+         do iq=1,nq
            pdqc(ig,l,iq) = 0
          enddo
@@ -681 +682 @@
 ! ***************************************************************
 
-!     DO l=1,nlayer
+!     DO l=1,klev
 !        DO ig=1,klon
 !         Taux de cond en kg.m-2.pa-1.s-1
@@ -689 +690 @@
 !        END DO
 !     END DO
-!     call WRITEDIAGFI(ngridmx,'tconda1',              &
+!     call WRITEDIAGFI(klon,'tconda1',              &
 !     'Taux de condensation N2 atmospherique /Pa',    &
 !      'kg.m-2.Pa-1.s-1',3,tconda1)
-!     call WRITEDIAGFI(ngridmx,'tconda2',              &
+!     call WRITEDIAGFI(klon,'tconda2',              &
 !     'Taux de condensation N2 atmospherique /m',     &
 !      'kg.m-3.s-1',3,tconda2)
@@ -698 +699 @@
 
   return
-  end subroutine condens_n2
+  end subroutine condense_n2
 
 !-------------------------------------------------------------------------
@@ -729 +730 @@
 !       tcond_n2 = (1.)/(0.026315-0.0011877*log(.7143*p*vmr))
   IF (t.ge.35.6) then
-!     tcond Fray and Schmitt for N2 phase beta (T>35.6 K) FIT TB 
+!     tcond Fray and Schmitt for N2 phase beta (T>35.6 K) FIT TB
   ! pcond_n2 = 0.125570*1.e5/vmr*exp((2.5e5*0.98)/296.925*(1./63.1470-1./t))
     pcond_n2 = 0.125570e5/vmr*exp(825.1241896*(1./63.147-1./t))
@@ -743 +744 @@
 
   real function glob_average2d(var)
-!   Calculates the global average of variable var     
+!   Calculates the global average of variable var
   use comgeomfi_h
+  use dimphy, only: klon
+  USE mod_grid_phy_lmdz, ONLY: nbp_lon, nbp_lat
+  use geometry_mod, only: cell_area, latitude
+
   implicit none
 
 !      INTEGER klon
-  REAL var(ngridmx)
+  REAL var(klon)
   INTEGER ig
 
   glob_average2d = 0.
-  DO ig=2,ngridmx-1
-       glob_average2d = glob_average2d + var(ig)*area(ig)
+  DO ig=2,klon-1
+       glob_average2d = glob_average2d + var(ig)*cell_area(ig)
   END DO
-  glob_average2d = glob_average2d + var(1)*area(1)*iim
-  glob_average2d = glob_average2d + var(ngridmx)*area(ngridmx)*iim
-  glob_average2d = glob_average2d/(totarea+(area(1)+area(ngridmx))*(iim-1))
+  glob_average2d = glob_average2d + var(1)*cell_area(1)*nbp_lon
+  glob_average2d = glob_average2d + var(klon)*cell_area(klon)*nbp_lon
+  glob_average2d = glob_average2d/(totarea+(cell_area(1)+cell_area(klon))*(nbp_lon-1))
 
   end function glob_average2d
@@ -763 +768 @@
 ! *****************************************************************
 
-  subroutine vl1d(nlayer,q,pente_max,masse,w,qm)
-!    
+  subroutine vl1d(klev,q,pente_max,masse,w,qm)
+!
 !     Operateur de moyenne inter-couche pour calcul de transport type
 !     Van-Leer " pseudo amont " dans la verticale
@@ -776 +781 @@
 !   Arguments:
 !   ----------
-  integer nlayer
-  real masse(nlayer),pente_max
-  REAL q(nlayer),qm(nlayer+1)
-  REAL w(nlayer+1)
-!
-!      Local 
+  integer klev
+  real masse(klev),pente_max
+  REAL q(klev),qm(klev+1)
+  REAL w(klev+1)
+!
+!      Local
 !   ---------
 !
   INTEGER l
 !
-  real dzq(nlayer),dzqw(nlayer),adzqw(nlayer),dzqmax
+  real dzq(klev),dzqw(klev),adzqw(klev),dzqmax
   real sigw, Mtot, MQtot
-  integer m 
-
-
-!    On oriente tout dans le sens de la pression 
+  integer m
+
+
+!    On oriente tout dans le sens de la pression
 !     W > 0 WHEN DOWN !!!!!!!!!!!!!
 
-  do l=2,nlayer
+  do l=2,klev
         dzqw(l)=q(l-1)-q(l)
         adzqw(l)=abs(dzqw(l))
   enddo
 
-  do l=2,nlayer-1
+  do l=2,klev-1
         if(dzqw(l)*dzqw(l+1).gt.0.) then
             dzq(l)=0.5*(dzqw(l)+dzqw(l+1))
@@ -810 +815 @@
 
      dzq(1)=0.
-     dzq(nlayer)=0.
-
-   do l = 1,nlayer-1
+     dzq(klev)=0.
+
+   do l = 1,klev-1
 
 !         Regular scheme (transfered mass < layer mass)
@@ -829 +834 @@
          Mtot = masse(m)
          MQtot = masse(m)*q(m)
-         !if (m.lt.nlayer) then ! because some compilers will have problems
-         !                      ! evaluating masse(nlayer+1)
-         do while ((m.lt.nlayer).and.(w(l+1).gt.(Mtot+masse(m+1))))
+         !if (m.lt.klev) then ! because some compilers will have problems
+         !                      ! evaluating masse(klev+1)
+         do while ((m.lt.klev).and.(w(l+1).gt.(Mtot+masse(m+1))))
             m=m+1
             Mtot = Mtot + masse(m)
             MQtot = MQtot + masse(m)*q(m)
-         !   if (m.eq.nlayer) exit
+         !   if (m.eq.klev) exit
          end do
          !endif
-         if (m.lt.nlayer) then
+         if (m.lt.klev) then
             sigw=(w(l+1)-Mtot)/masse(m+1)
             qm(l+1)= (1/w(l+1))*(MQtot + (w(l+1)-Mtot)* &
@@ -848 +853 @@
             stop
          end if
-      else      ! if(w(l+1).lt.0) 
-         m = l-1 
+      else      ! if(w(l+1).lt.0)
+         m = l-1
          Mtot = masse(m+1)
          MQtot = masse(m+1)*q(m+1)
@@ -863 +868 @@
          if (m.gt.0) then
             sigw=(w(l+1)+Mtot)/masse(m)
-            qm(l+1)= (-1/w(l+1))*(MQtot + (-w(l+1)-Mtot)* & 
+            qm(l+1)= (-1/w(l+1))*(MQtot + (-w(l+1)-Mtot)* &
            (q(m)-0.5*(1.+sigw)*dzq(m))  )
          else
             qm(l+1)= (-1/w(l+1))*(MQtot + (-w(l+1)-Mtot)*qm(1))
-         end if   
+         end if
       end if
    enddo

trunk/LMDZ.PLUTO/libf/phypluto/convadj.F

@@ -7 +7 @@
       USE tracer_h
       use comcstfi_mod, only: g
-      use callkeys_mod, only: tracer,water
+      use callkeys_mod, only: tracer
 
       implicit none
 
 !==================================================================
-!     
+!
 !     Purpose
 !     -------
@@ -23 +23 @@
 !     Original author unknown.
 !     Modif. 2005 by F. Forget.
-!     
+!
 !==================================================================
 
@@ -84 +84 @@
       ENDDO
 
-      if(tracer) then      
+      if(tracer) then
         DO iq =1, nq
          DO l=1,nlay
@@ -118 +118 @@
         ENDDO
       ENDDO
-      
+
 !     Make a list of them
       jcnt=0
@@ -136 +136 @@
 
           i = jadrs(jj)
- 
+
 !     Calculate sigma in this column
           DO l=1,nlay+1
             sig(l)=pplev(i,l)/pplev(i,1)
-        
+
           ENDDO
          DO l=1,nlay
@@ -154 +154 @@
             IF (l2 .GT. nlay) EXIT
             IF (zhc(i, l2).LT.zhc(i, l2-1)) THEN
- 
+
 !     l2 is the highest level of the unstable column
- 
+
               l1 = l2 - 1
               l  = l1
@@ -170 +170 @@
                 zhm = zhm + sdsig(l) * (zh2(i, l) - zhm) / zsm
                 zhmc = zhm
- 
+
 !     do we have to extend the column downwards?
- 
+
                 down = .false.
                 IF (l1 .ne. 1) then    !--  and then
@@ -179 +179 @@
                   END IF
                 END IF
- 
+
                 ! this could be a problem...
 
                 if (down) then
- 
+
                   l1 = l1 - 1
                   l  = l1
- 
+
                 else
- 
+
 !     can we extend the column upwards?
- 
+
                   if (l2 .eq. nlay) exit
- 
+
                   if (zhc(i, l2+1) .ge. zhmc) exit
 
@@ -224 +224 @@
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
-!     to conserve tracers/ KE, we must calculate zum, zvm and zqm using 
-!     the up-to-date column values. If we do not do this, there are cases 
+!     to conserve tracers/ KE, we must calculate zum, zvm and zqm using
+!     the up-to-date column values. If we do not do this, there are cases
 !     where convection stops at one level and starts at the next where we
 !     can break conservation of stuff (particularly tracers) significantly.
@@ -249 +249 @@
                  zv2(i,l)=zv2(i,l)+zalpha*(zvm-zv2(i,l))
                  do iq=1,nq
-!                  zq2(i,l,iq)=zq2(i,l,iq)+zalpha*(zqm(iq)-zq2(i,l,iq)) 
+!                  zq2(i,l,iq)=zq2(i,l,iq)+zalpha*(zqm(iq)-zq2(i,l,iq))
                    zq2(i,l,iq)=zqm(iq)
                  end do
@@ -266 +266 @@
 !     check conservation
          cadjncons=0.0
-         if(water)then
-         do l = 1, nlay
-            masse = (pplev(i,l) - pplev(i,l+1))/g
-            iq    = igcm_h2o_vap
-            cadjncons = cadjncons + 
-     &           masse*(zq2(i,l,iq)-zq(i,l,iq))/ptimestep 
-         end do
-         endif
 
          if(cadjncons.lt.-1.e-6)then
@@ -283 +275 @@
          do l = 1, nlay
             print*,'dsig         = ',dsig(l)
-         end do         
+         end do
          do l = 1, nlay
             print*,'dsig         = ',dsig(l)
@@ -311 +303 @@
             print*,'zh2(ig,:)        = ',zh2(i,l)
          end do
-         do l = 1, nlay
-            print*,'zq(ig,:,vap)     = ',zq(i,l,igcm_h2o_vap)
-         end do
-         do l = 1, nlay
-            print*,'zq2(ig,:,vap)    = ',zq2(i,l,igcm_h2o_vap)
-         end do
-            print*,'zqm(vap)         = ',zqm(igcm_h2o_vap)
+         ! do l = 1, nlay
+         !    print*,'zq(ig,:,vap)     = ',zq(i,l,igcm_h2o_vap)
+         ! end do
+         ! do l = 1, nlay
+         !    print*,'zq2(ig,:,vap)    = ',zq2(i,l,igcm_h2o_vap)
+         ! end do
+         !    print*,'zqm(vap)         = ',zqm(igcm_h2o_vap)
             print*,'jadrs=',jadrs
 
@@ -335 +327 @@
       ENDDO
 
-      if(tracer) then 
+      if(tracer) then
         do iq=1, nq
           do  l=1,nlay
             DO ig=1,ngrid
-              pdqadj(ig,l,iq)=(zq2(ig,l,iq)-zq(ig,l,iq))/ptimestep 
-            end do 
-          end do 
-        end do 
+              pdqadj(ig,l,iq)=(zq2(ig,l,iq)-zq(ig,l,iq))/ptimestep
+            end do
+          end do
+        end do
       end if
 
@@ -349 +341 @@
 !      if (ngrid.eq.1) then
 !         ig=1
-!         iq =1 
-!         write(*,*)'**** l, pq(ig,l,iq),zq(ig,l,iq),zq2(ig,l,iq)'  
+!         iq =1
+!         write(*,*)'**** l, pq(ig,l,iq),zq(ig,l,iq),zq2(ig,l,iq)'
 !         do l=nlay,1,-1
 !           write(*,*) l, pq(ig,l,iq),zq(ig,l,iq),zq2(ig,l,iq)

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

@@ -24 +24 @@
   use planetwide_mod, only: planetwide_sumval
   use callkeys_mod
+  use surfdat_h
   use wstats_mod, only: callstats
   use ioipsl_getin_p_mod, only : getin_p
@@ -277 +278 @@
      if (is_master) write(*,*) trim(rname)//": calladj = ",calladj
 
-     if (is_master) write(*,*) trim(rname)//": call N2 condensation ?"
-     n2cond=.false. ! default value
-     call getin_p("n2cond",n2cond)
-     if (is_master) write(*,*) trim(rname)//": n2cond = ",n2cond
-! Test of incompatibility
-     if (n2cond.and.(.not.tracer)) then
-        call abort_physic(rname,"Error: We need a N2 ice tracer to condense N2!",1)
-     endif 
- 
-     if (is_master) write(*,*) trim(rname)//": N2 supersaturation level ?"
-     n2supsat=1.0 ! default value
-     call getin_p("n2supsat",n2supsat)
-     if (is_master) write(*,*) trim(rname)//": n2supsat = ",n2supsat
-
      if (is_master) write(*,*) trim(rname)//&
      ": Radiative timescale for Newtonian cooling (in Earth days)?"
@@ -479 +466 @@
      endif  
 
-     if (is_master) write(*,*)trim(rname)//&
-       ": Number mixing ratio of N2 ice particles:"
-     Nmix_n2=1.e6 ! default value
-     call getin_p("Nmix_n2",Nmix_n2)
-     if (is_master) write(*,*)trim(rname)//": Nmix_n2 = ",Nmix_n2
-
-     if (is_master) write(*,*)trim(rname)//&
+      if (is_master) write(*,*)trim(rname)//&
        "Number of radiatively active aerosols:"
      naerkind=0 ! default value
@@ -491 +472 @@
      if (is_master) write(*,*)trim(rname)//": naerkind = ",naerkind
 
-     if (is_master) write(*,*)trim(rname)//": Opacity of dust (if used):"
-     dusttau=0. ! default value
-     call getin_p("dusttau",dusttau)
-     if (is_master) write(*,*)trim(rname)//": dusttau = ",dusttau
-
-     if (is_master) write(*,*)trim(rname)//": Radiatively active N2 aerosols?"
-     aeron2=.false.     ! default value
-     call getin_p("aeron2",aeron2)
-     if (is_master) write(*,*)trim(rname)//": aeron2 = ",aeron2
-
-     if (is_master) write(*,*)trim(rname)//": Fixed N2 aerosol distribution?"
-     aerofixn2=.false.     ! default value
-     call getin_p("aerofixn2",aerofixn2)
-     if (is_master) write(*,*)trim(rname)//": aerofixn2 = ",aerofixn2
-
-     if (is_master) write(*,*)trim(rname)//&
-       ": Radiatively active sulfuric acid aerosols?"
-     aeroh2so4=.false.     ! default value
-     call getin_p("aeroh2so4",aeroh2so4)
-     if (is_master) write(*,*)trim(rname)//": aeroh2so4 = ",aeroh2so4
- 
-     if (is_master) write(*,*)trim(rname)//&
-       ": Radiatively active auroral aerosols?"
-     aeroaurora=.false.     ! default value
-     call getin_p("aeroaurora",aeroaurora)
-     if (is_master) write(*,*)trim(rname)//": aeroaurora = ",aeroaurora
-
-     if (is_master) write(*,*)trim(rname)//&
-     ": Radiatively active Generic Condensable aerosols?"
-     aerogeneric=0     ! default value
-     call getin_p("aerogeneric",aerogeneric)
-     if (is_master) write(*,*)trim(rname)//":aerogeneric",aerogeneric
-
+     
+!***************************************************************
+         !! TRACERS options
+
+     if (is_master)write(*,*)trim(rname)//&
+      "call N2 condensation ?"
+     N2cond=.true. ! default value
+     call getin_p("N2cond",N2cond)
+     if (is_master)write(*,*)trim(rname)//&
+      " N2cond = ",N2cond
+
+     if (is_master)write(*,*)trim(rname)//&
+      "call N2 cloud condensation ?"
+     condensn2=.false. ! default value
+     call getin_p("condensn2",condensn2)
+     if (is_master)write(*,*)trim(rname)//&
+      "condensn2 = ",condensn2
+
+     if (is_master)write(*,*)trim(rname)//&
+      "call no N2 frost formation?"
+     no_n2frost=.false. ! default value
+     call getin_p("no_n2frost",no_n2frost)
+     if (is_master)write(*,*)trim(rname)//&
+      "no_n2frost = ",no_n2frost
+
+     if (is_master)write(*,*)trim(rname)//&
+      "N2 condensation subtimestep?"
+     nbsub=20 ! default value
+     call getin_p("nbsub",nbsub)
+     if (is_master)write(*,*)trim(rname)//&
+      " nbsub = ",nbsub
+
+     if (is_master)write(*,*)trim(rname)//&
+      "Gravitationnal sedimentation ?"
+     sedimentation=.true. ! default value
+     call getin_p("sedimentation",sedimentation)
+     if (is_master)write(*,*)trim(rname)//&
+      " sedimentation = ",sedimentation
+
+     if (is_master)write(*,*)trim(rname)//&
+      "Compute glacial flow ?"
+     glaflow=.false. ! default value
+     call getin_p("glaflow",glaflow)
+     if (is_master)write(*,*)trim(rname)//&
+      " glaflow = ",glaflow
+
+     if (is_master)write(*,*)trim(rname)//&
+      "Compute methane cycle ?"
+     methane=.false. ! default value
+     call getin_p("methane",methane)
+     if (is_master)write(*,*)trim(rname)//&
+      " methane = ",methane
+     condmetsurf=.true. ! default value
+     call getin_p("condmetsurf",condmetsurf)
+     if (is_master)write(*,*)trim(rname)//&
+      " condmetsurf = ",condmetsurf
+
+     if (is_master)write(*,*)trim(rname)//&
+      "Compute methane clouds ?"
+     metcloud=.false. ! default value
+     call getin_p("metcloud",metcloud)
+     if (is_master)write(*,*)trim(rname)//&
+      " metcloud = ",metcloud
+
+     if (is_master)write(*,*)trim(rname)//&
+      "Compute CO cycle ?"
+     carbox=.false. ! default value
+     call getin_p("carbox",carbox)
+     if (is_master)write(*,*)trim(rname)//&
+      " carbox = ",carbox
+     condcosurf=.true. ! default value
+     call getin_p("condcosurf",condcosurf)
+     if (is_master)write(*,*)trim(rname)//&
+      " condcosurf = ",condcosurf
+
+     if (is_master)write(*,*)trim(rname)//&
+      "Compute CO clouds ?"
+     monoxcloud=.false. ! default value
+     call getin_p("monoxcloud",monoxcloud)
+     if (is_master)write(*,*)trim(rname)//&
+      " monoxcloud = ",monoxcloud
+
+     if (is_master)write(*,*)trim(rname)//&
+     "atmospheric redistribution (s):"
+     tau_n2=1. ! default value
+     call getin_p("tau_n2",tau_n2)
+     if (is_master)write(*,*)trim(rname)//&
+     " tau_n2 = ",tau_n2
+     tau_ch4=1.E7 ! default value
+     call getin_p("tau_ch4",tau_ch4)
+     if (is_master)write(*,*)trim(rname)//&
+     " tau_ch4 = ",tau_ch4
+     tau_co=1. ! default value
+     call getin_p("tau_co",tau_co)
+     if (is_master)write(*,*)trim(rname)//&
+     " tau_co = ",tau_co
+     tau_prechaze=1. ! default value
+     call getin_p("tau_prechaze",tau_prechaze)
+     if (is_master)write(*,*)trim(rname)//&
+     " tau_prechaze = ",tau_prechaze
+
+     if (is_master)write(*,*)trim(rname)//&
+      "Day fraction for limited cold trap in SP?"
+     dayfrac=0. ! default value
+     call getin_p("dayfrac",dayfrac)
+     if (is_master)write(*,*)trim(rname)//&
+      " dayfrac = ",dayfrac
+     thresh_non2=0. ! default value
+     call getin_p("thresh_non2",thresh_non2)
+     if (is_master)write(*,*)trim(rname)//&
+      " thresh_non2 = ",thresh_non2
+
+!***************************************************************
+     !! Haze options
+
+     if (is_master)write(*,*)trim(rname)//&
+     "Production of haze ?"
+     haze=.false. ! default value
+     call getin_p("haze",haze)
+     if (is_master)write(*,*)trim(rname)//&
+     " haze = ",haze
+     hazeconservch4=.false. ! conservch4, by default value ch4 is photolyzed
+     call getin_p("hazeconservch4",hazeconservch4)
+     if (is_master)write(*,*)trim(rname)//&
+     "hazconservch4 = ",hazeconservch4
+     if (is_master)write(*,*)trim(rname)//&
+     "Production of haze (fast version) ?"
+     fasthaze=.false. ! default value
+     call getin_p("fasthaze",fasthaze)
+     if (is_master)write(*,*)trim(rname)//&
+     "fasthaze = ",fasthaze
+
+     if (is_master)write(*,*)trim(rname)//&
+     "Add source of haze ?"
+     source_haze=.false. ! default value
+     call getin_p("source_haze",source_haze)
+     if (is_master)write(*,*)trim(rname)//&
+     " source_haze = ",source_haze
+     mode_hs=0 ! mode haze source default value
+     call getin_p("mode_hs",mode_hs)
+     if (is_master)write(*,*)trim(rname)//&
+     " mode_hs",mode_hs
+     kfix=1 ! default value
+     call getin_p("kfix",kfix)
+     if (is_master)write(*,*)trim(rname)//&
+     "levels kfix",kfix
+     fracsource=0.1 ! default value
+     call getin_p("fracsource",fracsource)
+     if (is_master)write(*,*)trim(rname)//&
+     " fracsource",fracsource
+     latsource=30. ! default value
+     call getin_p("latsource",latsource)
+     if (is_master)write(*,*)trim(rname)//&
+     " latsource",latsource
+     lonsource=180. ! default value
+     call getin_p("lonsource",lonsource)
+     if (is_master)write(*,*)trim(rname)//&
+     " lonsource",lonsource
+
+     if (is_master)write(*,*)trim(rname)//&
+     "Radiatively active haze ?"
+     aerohaze=.false. ! default value
+     call getin_p("aerohaze",aerohaze)
+     if (is_master)write(*,*)trim(rname)//&
+     "aerohaze = ",aerohaze
+
+     if (is_master)write(*,*)trim(rname)//&
+     "Haze monomer radius ?"
+     rad_haze=10.e-9 ! default value
+     call getin_p("rad_haze",rad_haze)
+     if (is_master)write(*,*)trim(rname)//&
+     "rad_haze = ",rad_haze
+
+     if (is_master)write(*,*)trim(rname)//&
+     "fractal particle ?"
+     fractal=.false. ! default value
+     call getin_p("fractal",fractal)
+     if (is_master)write(*,*)trim(rname)//&
+     "fractal = ",fractal
+     nb_monomer=10 ! default value
+     call getin_p("nb_monomer",nb_monomer)
+     if (is_master)write(*,*)trim(rname)//&
+     "nb_monomer = ",nb_monomer
+     
+     if (is_master)write(*,*)trim(rname)//&
+     "fixed radius profile from txt file ?"
+     haze_radproffix=.false. ! default value
+     call getin_p("haze_radproffix",haze_radproffix)
+     if (is_master)write(*,*)trim(rname)//&
+     "haze_radproffix = ",haze_radproffix
+     if (is_master)write(*,*)trim(rname)//&
+     "fixed MMR profile from txt file ?"
+     haze_proffix=.false. ! default value
+     call getin_p("haze_proffix",haze_proffix)
+     if (is_master)write(*,*)trim(rname)//&
+     "haze_proffix = ",haze_proffix
+
+     if (is_master)write(*,*)trim(rname)//&
+     "Number mix ratio of haze particles for co clouds:"
+     Nmix_co=100000. ! default value
+     call getin_p("Nmix_co",Nmix_co)
+     if (is_master)write(*,*)trim(rname)//&
+     " Nmix_co = ",Nmix_co
+
+     if (is_master)write(*,*)trim(rname)//&
+     "Number mix ratio of haze particles for ch4 clouds:"
+     Nmix_ch4=100000. ! default value
+     call getin_p("Nmix_ch4",Nmix_ch4)
+     if (is_master)write(*,*)trim(rname)//&
+     " Nmix_ch4 = ",Nmix_ch4
+
+!***************************************************************
+     !! Surface properties
+
+!*********** N2 *********************************
+
+     if (is_master)write(*,*)trim(rname)//&
+      "Mode for changing N2 albedo"
+     mode_n2=0 ! default value
+     call getin_p("mode_n2",mode_n2)
+     if (is_master)write(*,*)trim(rname)//&
+      " mode_n2 = ",mode_n2
+     thres_n2ice=1. ! default value
+     call getin_p("thres_n2ice",thres_n2ice)
+     if (is_master)write(*,*)trim(rname)//&
+      " thres_n2ice = ",thres_n2ice
+
+     if (is_master)write(*,*)trim(rname)//&
+      "Diff of N2 albedo with thickness"
+     deltab=0. ! default value
+     call getin_p("deltab",deltab)
+     if (is_master)write(*,*)trim(rname)//&
+      " deltab = ",deltab
+
+     if (is_master)write(*,*)trim(rname)//&
+      "albedo N2 beta "
+     alb_n2b=0.7 ! default value
+     call getin_p("alb_n2b",alb_n2b)
+     if (is_master)write(*,*)trim(rname)//&
+      " alb_n2b = ",alb_n2b
+
+     if (is_master)write(*,*)trim(rname)//&
+      "albedo N2 alpha "
+     alb_n2a=0.7 ! default value
+     call getin_p("alb_n2a",alb_n2a)
+     if (is_master)write(*,*)trim(rname)//&
+      " alb_n2a = ",alb_n2a
+
+     if (is_master)write(*,*)trim(rname)//&
+      "emis N2 beta "
+     emis_n2b=0.7 ! default value
+     call getin_p("emis_n2b",emis_n2b)
+     if (is_master)write(*,*)trim(rname)//&
+      " emis_n2b = ",emis_n2b
+
+     if (is_master)write(*,*)trim(rname)//&
+      "emis N2 alpha "
+     emis_n2a=0.7 ! default value
+     call getin_p("emis_n2a",emis_n2a)
+     if (is_master)write(*,*)trim(rname)//&
+      " emis_n2a = ",emis_n2a
+
+!*********** CH4 *********************************
+
+     if (is_master)write(*,*)trim(rname)//&
+      "Mode for changing CH4 albedo"
+     mode_ch4=0 ! default value
+     call getin_p("mode_ch4",mode_ch4)
+     if (is_master)write(*,*)trim(rname)//&
+      " mode_ch4 = ",mode_ch4
+     feedback_met=0 ! default value
+     call getin_p("feedback_met",feedback_met)
+     if (is_master)write(*,*)trim(rname)//&
+      " feedback_met = ",feedback_met
+     thres_ch4ice=1. ! default value, mm
+     call getin_p("thres_ch4ice",thres_ch4ice)
+     if (is_master)write(*,*)trim(rname)//&
+      " thres_ch4ice = ",thres_ch4ice
+     fdch4_latn=-1. 
+     fdch4_lats=0. 
+     fdch4_lone=0. 
+     fdch4_lonw=-1. 
+     fdch4_depalb=0.5 
+     fdch4_finalb=0.95 
+     fdch4_maxalb=0.99 
+     fdch4_ampl=200. 
+     fdch4_maxice=100. 
+     call getin_p("  fdch4_latn",fdch4_latn)
+     call getin_p("  fdch4_lats",fdch4_lats)
+     call getin_p("  fdch4_lone",fdch4_lone)
+     call getin_p("  fdch4_lonw",fdch4_lonw)
+     call getin_p("  fdch4_depalb",fdch4_depalb)
+     call getin_p("  fdch4_finalb",fdch4_finalb)
+     call getin_p("  fdch4_maxalb",fdch4_maxalb)
+     call getin_p("  fdch4_maxice",fdch4_maxice)
+     call getin_p("  fdch4_ampl",fdch4_ampl)
+     if (is_master)write(*,*)trim(rname)//&
+      " Values for albedo feedback = ",fdch4_latn,&
+     fdch4_lats,fdch4_lone,fdch4_lonw,fdch4_depalb,&
+     fdch4_finalb,fdch4_maxalb,fdch4_maxice,fdch4_ampl
+
+     if (is_master)write(*,*)trim(rname)//&
+      "Latitude for diff albedo methane"
+     metlateq=25. ! default value
+     call getin_p("metlateq",metlateq)
+     if (is_master)write(*,*)trim(rname)//&
+      " metlateq = ",metlateq
+
+     if (is_master)write(*,*)trim(rname)//&
+      "Ls1 and Ls2 for change of ch4 albedo"
+     metls1=-1. ! default value
+     metls2=-2. ! default value
+     call getin_p("metls1",metls1)
+     call getin_p("metls2",metls2)
+
+     if (is_master)write(*,*)trim(rname)//&
+      "albedo CH4 "
+     alb_ch4=0.5 ! default value
+     call getin_p("alb_ch4",alb_ch4)
+     if (is_master)write(*,*)trim(rname)//&
+      " alb_ch4 = ",alb_ch4
+
+     if (is_master)write(*,*)trim(rname)//&
+      "albedo equatorial CH4 "
+     alb_ch4_eq=alb_ch4 ! default value
+     call getin_p("alb_ch4_eq",alb_ch4_eq)
+     if (is_master)write(*,*)trim(rname)//&
+      " alb_ch4_eq = ",alb_ch4_eq
+
+     if (is_master)write(*,*)trim(rname)//&
+      "albedo south hemis CH4 "
+     alb_ch4_s=alb_ch4 ! default value
+     call getin_p("alb_ch4_s",alb_ch4_s)
+     if (is_master)write(*,*)trim(rname)//&
+      " alb_ch4_s = ",alb_ch4_s
+
+     if (is_master)write(*,*)trim(rname)//&
+      "emis CH4 "
+     emis_ch4=0.5 ! default value
+     call getin_p("emis_ch4",emis_ch4)
+     if (is_master)write(*,*)trim(rname)//&
+      " emis_ch4 = ",emis_ch4
+
+     if (is_master)write(*,*)trim(rname)//&
+      "CH4 lag for n2 sublimation limitation"
+     ch4lag=.false. ! default value
+     latlag=-90. ! default value
+     vmrlag=1. ! default value
+     call getin_p("ch4lag",ch4lag)
+     call getin_p("latlag",latlag)
+     call getin_p("vmrlag",vmrlag)
+     if (is_master)write(*,*)trim(rname)//&
+      " ch4lag = ",ch4lag
+     if (is_master)write(*,*)trim(rname)//&
+      " latlag = ",latlag
+     if (is_master)write(*,*)trim(rname)//&
+      " vmrlag = ",vmrlag
+
+     if (is_master)write(*,*)trim(rname)//&
+      "Max temperature for surface ?"
+     tsurfmax=.false. ! default value
+     albmin_ch4=0.3 ! default value
+     call getin_p("tsurfmax",tsurfmax)
+     call getin_p("albmin_ch4",albmin_ch4)
+     if (is_master)write(*,*)trim(rname)//&
+      " tsurfmax = ",tsurfmax
+     if (is_master)write(*,*)trim(rname)//&
+      " albmin_ch4 = ",albmin_ch4
+
+!*********** CO *********************************
+
+     if (is_master)write(*,*)trim(rname)//&
+      "albedo CO "
+     alb_co=0.4 ! default value
+     call getin_p("alb_co",alb_co)
+     if (is_master)write(*,*)trim(rname)//&
+      " alb_co = ",alb_co
+     thres_coice=1. ! default value, mm
+     call getin_p("thres_coice",thres_coice)
+     if (is_master)write(*,*)trim(rname)//&
+      " thres_coice = ",thres_coice
+
+     if (is_master)write(*,*)trim(rname)//&
+      "emis CO "
+     emis_co=0.4 ! default value
+     call getin_p("emis_co",emis_co)
+     if (is_master)write(*,*)trim(rname)//&
+      " emis_co = ",emis_co
+
+!*********** THOLINS *********************************
+     if (is_master)write(*,*)trim(rname)//&
+      "Mode for changing tholins albedo/emis"
+     mode_tholins=0 ! default value
+     call getin_p("mode_tholins",mode_tholins)
+     if (is_master)write(*,*)trim(rname)//&
+      " mode_tholins = ",mode_tholins
+
+     if (is_master)write(*,*)trim(rname)//&
+      "albedo tho "
+     alb_tho=0.1 ! default value
+     call getin_p("alb_tho",alb_tho)
+     if (is_master)write(*,*)trim(rname)//&
+      " alb_tho = ",alb_tho
+
+     if (is_master)write(*,*)trim(rname)//&
+      "albedo tho eq"
+     alb_tho_eq=0.1 ! default value
+     call getin_p("alb_tho_eq",alb_tho_eq)
+     if (is_master)write(*,*)trim(rname)//&
+      " alb_tho_eq = ",alb_tho_eq
+
+     if (is_master)write(*,*)trim(rname)//&
+      "emis tholins "
+     emis_tho=1. ! default value
+     call getin_p("emis_tho",emis_tho)
+     if (is_master)write(*,*)trim(rname)//&
+      " emis_tho = ",emis_tho
+
+     if (is_master)write(*,*)trim(rname)//&
+      "emis tholins eq"
+     emis_tho_eq=1. ! default value
+     call getin_p("emis_tho_eq",emis_tho_eq)
+     if (is_master)write(*,*)trim(rname)//&
+      " emis_tho_eq = ",emis_tho_eq
+
+     if (is_master)write(*,*)trim(rname)//&
+      "Latitude for diff albedo tholins"
+     tholateq=25. ! default value
+     call getin_p("tholateq",tholateq)
+     if (is_master)write(*,*)trim(rname)//&
+      " tholateq = ",tholateq
+     tholatn=-1. 
+     tholats=0. 
+     tholone=0. 
+     tholonw=-1. 
+     alb_tho_spe=0.1 ! default value
+     emis_tho_spe=1. ! default value
+     call getin_p("  tholatn",tholatn)
+     call getin_p("  tholats",tholats)
+     call getin_p("  tholone",tholone)
+     call getin_p("  tholonw",tholonw)
+     if (is_master)write(*,*)trim(rname)//&
+      " Values for special tholins albedo = ",tholatn,&
+        tholats,tholone,tholonw,alb_tho_spe,emis_tho_spe
+
+     if (is_master)write(*,*)trim(rname)//&
+      "Specific albedo"
+     spelon1=-180. ! default value
+     spelon2=180. ! default value
+     spelat1=-90. ! default value
+     spelat2=90. ! default value
+     specalb=.false. ! default value
+     if (is_master)write(*,*)trim(rname)//&
+      "albedo/emis spe "
+     albspe=0.1 ! default value
+     emispe=1. ! default value
+     call getin_p("spelon1",spelon1)
+     call getin_p("spelon2",spelon2)
+     call getin_p("spelat1",spelat1)
+     call getin_p("spelat2",spelat2)
+     call getin_p("specalb",specalb)
+     call getin_p("albspe",albspe)
+     call getin_p("emispe",emispe)
+
+     if (is_master)write(*,*)trim(rname)//&
+      " specific = ",specalb
+
+!********************** TI *****************************
+
+     if (is_master)write(*,*)trim(rname)//&
+      "Change TI with time"
+     changeti=.false. ! default value
+     call getin_p("changeti",changeti)
+     if (is_master)write(*,*)trim(rname)//&
+      " changeti = ",changeti
+     changetid=.false. ! default value for diurnal TI
+     call getin_p("changetid",changetid)
+     if (is_master)write(*,*)trim(rname)//&
+      " changetid = ",changetid
+
+     if (is_master)write(*,*)trim(rname)//&
+      "IT N2 "
+     ITN2=800. ! default value
+     call getin_p("ITN2",ITN2)
+     if (is_master)write(*,*)trim(rname)//&
+      " ITN2 = ",ITN2
+     ITN2d=20. ! default value
+     call getin_p("ITN2d",ITN2d)
+     if (is_master)write(*,*)trim(rname)//&
+      " ITN2d = ",ITN2d
+
+     if (is_master)write(*,*)trim(rname)//&
+      "IT CH4"
+     ITCH4=800. ! default value
+     call getin_p("ITCH4",ITCH4)
+     if (is_master)write(*,*)trim(rname)//&
+      " ITCH4 = ",ITCH4
+     ITCH4d=20. ! default value
+     call getin_p("ITCH4d",ITCH4d)
+     if (is_master)write(*,*)trim(rname)//&
+      " ITCH4d = ",ITCH4d
+
+     if (is_master)write(*,*)trim(rname)//&
+      "IT H2O"
+     ITH2O=800. ! default value
+     call getin_p("ITH2O",ITH2O)
+     if (is_master)write(*,*)trim(rname)//&
+      " ITH2O = ",ITH2O
+     ITH2Od=20. ! default value
+     call getin_p("ITH2Od",ITH2Od)
+     if (is_master)write(*,*)trim(rname)//&
+      " ITH2Od = ",ITH2Od
 
 !=================================

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

@@ -3 +3 @@
       use surfdat_h, ONLY: dryness
       USE tracer_h
+      USE callkeys_mod, only: aerohaze,nb_monomer,haze,fractal,fasthaze,rad_haze
       USE recombin_corrk_mod, ONLY: ini_recombin
       USE mod_phys_lmdz_para, only: is_master, bcast
       use generic_cloud_common_h
+      use aerosol_mod, only: iaero_haze,i_haze
       IMPLICIT NONE
 !=======================================================================
 !   subject:
 !   --------
-!   Initialization related to tracer 
+!   Initialization related to tracer
 !   (transported dust, water, chemical species, ice...)
 !
@@ -23 +25 @@
 !            Ehouarn Millour (oct. 2008) identify tracers by their names
 !            Y Jaziri & J. Vatant d'Ollone (2020) : Modern traceur.def
-!            L Teinturier (2022): Tracer names are now read here instead of 
-!                                  inside interfaces 
+!            L Teinturier (2022): Tracer names are now read here instead of
+!                                  inside interfaces
 !=======================================================================
 
@@ -30 +32 @@
 
       character(len=500) :: tracline   ! to read traceur.def lines
-      ! character(len=30)  :: txt        ! to store some text
       integer :: blank      !to store the index of 1st blank when reading tracers names
       integer iq,ig,count,ierr
+      real r0_lift , reff_lift, rho_haze
+      integer nqhaze(nq)               ! to store haze tracers
+      integer i, ia, block
+      character(len=20) :: txt ! to store some text
+      CHARACTER(LEN=20) :: tracername ! to temporarily store text
 
 !-----------------------------------------------------------------------
@@ -86 +92 @@
          ENDDO
         ENDIF ! if modern or standard traceur.def
-       
+
        endif ! of if (is_master)
-       
+
        ! share the information with other procs/threads (if any)
        CALL bcast(nqtot)
        CALL bcast(moderntracdef)
-       
+
 !! -----------------------------------------------------------------------
        !! For the moment number of tracers in dynamics and physics are equal
@@ -108 +114 @@
        IF (.NOT.ALLOCATED(alpha_devil))    ALLOCATE(alpha_devil(nq))
        IF (.NOT.ALLOCATED(qextrhor))       ALLOCATE(qextrhor(nq))
-       IF (.NOT.ALLOCATED(igcm_dustbin))   ALLOCATE(igcm_dustbin(nq))
+      !  IF (.NOT.ALLOCATED(igcm_dustbin))   ALLOCATE(igcm_dustbin(nq))
        IF (.NOT.ALLOCATED(is_chim))        ALLOCATE(is_chim(nqtot))
        IF (.NOT.ALLOCATED(is_rad))         ALLOCATE(is_rad(nqtot))
@@ -140 +146 @@
        is_recomb_qset(:) = 0
        is_recomb_qotf(:) = 0
-       
+
        ! Added by JVO 2017 : these arrays are handled later
        ! -> initialization is the least we can do, please !!!
@@ -163 +169 @@
 
 ! Initialization: Read tracers names from traceur.def
-        do iq=1,nq 
+        do iq=1,nq
           if (is_master) read(407,'(A)') tracline
           call bcast(tracline)
@@ -173 +179 @@
       ! 0. initialize tracer indexes to zero:
       ! NB: igcm_* indexes are commons in 'tracer.h'
-      do iq=1,nq
-        igcm_dustbin(iq)=0
-      enddo
-      igcm_dust_mass=0
-      igcm_dust_number=0
-      igcm_h2o_vap=0
-      igcm_h2o_ice=0
-      igcm_co2=0
-      igcm_co=0
-      igcm_o=0
-      igcm_o1d=0
-      igcm_o2=0
-      igcm_o3=0
-      igcm_h=0
-      igcm_h2=0
-      igcm_oh=0
-      igcm_ho2=0
-      igcm_h2o2=0
       igcm_n2=0
-      igcm_n=0
-      igcm_n2d=0
-      igcm_no=0
-      igcm_no2=0
-      igcm_ar=0
-      igcm_ar_n2=0
-      igcm_n2_ice=0
-
-      igcm_ch4=0
-      igcm_ch3=0
-      igcm_ch=0
-      igcm_3ch2=0
-      igcm_1ch2=0
-      igcm_cho=0
-      igcm_ch2o=0
-      igcm_ch3o=0
-      igcm_c=0
-      igcm_c2=0
-      igcm_c2h=0
-      igcm_c2h2=0
-      igcm_c2h3=0
-      igcm_c2h4=0
-      igcm_c2h6=0
-      igcm_ch2co=0
-      igcm_ch3co=0
-      igcm_hcaer=0
+      igcm_ch4_gas=0
+      igcm_ch4_ice=0
+      igcm_prec_haze=0
+      igcm_co_gas=0
+      igcm_co_ice=0
+
+      nqhaze(:)=0
+      i=0
+      DO iq=1,nq
+         txt=noms(iq)
+         IF (txt(1:4).eq."haze") THEN
+            i=i+1
+            nqhaze(i)=iq
+         ENDIF
+      ENDDO
+      if ((haze.or.fasthaze).and.i==0) then
+         print*, 'Haze active but no haze tracer in traceur.def'
+         stop
+      endif
+
+      igcm_haze=0
+      igcm_haze10=0
+      igcm_haze30=0
+      igcm_haze50=0
+      igcm_haze100=0
+
+!     Eddy diffusion tracers
+      igcm_eddy1e6=0
+      igcm_eddy1e7=0
+      igcm_eddy5e7=0
+      igcm_eddy1e8=0
+      igcm_eddy5e8=0
+      write(*,*) 'initracer: noms() ', noms
+      rho_n2=1000        ! n2 ice
+      rho_ch4_ice=520.       ! rho ch4 ice
+      rho_co_ice=520.       ! rho ch4 ice
+      rho_haze=800.     ! haze
 
       ! 1. find dust tracers
@@ -224 +223 @@
       ! 2. find chemistry and water tracers
       do iq=1,nq
-        if (noms(iq).eq."co2") then
-          igcm_co2=iq
-          mmol(igcm_co2)=44.
-          count=count+1
-!          write(*,*) 'co2: count=',count
-        endif
-!         if (noms(iq).eq."co2_ice") then
-!           igcm_co2_ice=iq
-!           mmol(igcm_co2_ice)=44.
-!           count=count+1
-! !          write(*,*) 'co2_ice: count=',count
-!         endif
-        if (noms(iq).eq."h2o_vap") then
-          igcm_h2o_vap=iq
-          mmol(igcm_h2o_vap)=18.
-          count=count+1
-!          write(*,*) 'h2o_vap: count=',count
-        endif
-        if (noms(iq).eq."h2o_ice") then
-          igcm_h2o_ice=iq
-          mmol(igcm_h2o_ice)=18.
-          count=count+1
-!          write(*,*) 'h2o_ice: count=',count
-        endif
-        if (noms(iq).eq."co") then
-          igcm_co=iq
-          mmol(igcm_co)=28.
-          count=count+1
-        endif
-        if (noms(iq).eq."o") then
-          igcm_o=iq
-          mmol(igcm_o)=16.
-          count=count+1
-        endif
-        if (noms(iq).eq."o1d") then
-          igcm_o1d=iq
-          mmol(igcm_o1d)=16.
-          count=count+1
-        endif
-        if (noms(iq).eq."o2") then
-          igcm_o2=iq
-          mmol(igcm_o2)=32.
-          count=count+1
-        endif
-        if (noms(iq).eq."o3") then
-          igcm_o3=iq
-          mmol(igcm_o3)=48.
-          count=count+1
-        endif
-        if (noms(iq).eq."h") then
-          igcm_h=iq
-          mmol(igcm_h)=1.
-          count=count+1
-        endif
-        if (noms(iq).eq."h2") then
-          igcm_h2=iq
-          mmol(igcm_h2)=2.
-          count=count+1
-        endif
-        if (noms(iq).eq."oh") then
-          igcm_oh=iq
-          mmol(igcm_oh)=17.
-          count=count+1
-        endif
-        if (noms(iq).eq."ho2") then
-          igcm_ho2=iq
-          mmol(igcm_ho2)=33.
-          count=count+1
-        endif
-        if (noms(iq).eq."h2o2") then
-          igcm_h2o2=iq
-          mmol(igcm_h2o2)=34.
-          count=count+1
-        endif
         if (noms(iq).eq."n2") then
           igcm_n2=iq
           mmol(igcm_n2)=28.
           count=count+1
-        endif
-        if (noms(iq).eq."n2_ice") then
-          igcm_n2_ice=iq
-          mmol(igcm_n2_ice)=28.
-          count=count+1
-!          write(*,*) 'n2_ice: count=',count
-        endif
-        if (noms(iq).eq."ch4") then
-          igcm_ch4=iq
-          mmol(igcm_ch4)=16.
-          count=count+1
-        endif
-        if (noms(iq).eq."ar") then
-          igcm_ar=iq
-          mmol(igcm_ar)=40.
-          count=count+1
-        endif
-        if (noms(iq).eq."n") then
-          igcm_n=iq
-          mmol(igcm_n)=14.
-          count=count+1
-        endif
-        if (noms(iq).eq."no") then
-          igcm_no=iq
-          mmol(igcm_no)=30.
-          count=count+1
-        endif
-        if (noms(iq).eq."no2") then
-          igcm_no2=iq
-          mmol(igcm_no2)=46.
-          count=count+1
-        endif
-        if (noms(iq).eq."n2d") then
-          igcm_n2d=iq
-          mmol(igcm_n2d)=28.
-          count=count+1
-        endif
-        if (noms(iq).eq."ch3") then
-          igcm_ch3=iq
-          mmol(igcm_ch3)=15.
-          count=count+1
-        endif
-        if (noms(iq).eq."ch") then
-          igcm_ch=iq
-          mmol(igcm_ch)=13.
-          count=count+1
-        endif
-        if (noms(iq).eq."3ch2") then
-          igcm_3ch2=iq
-          mmol(igcm_3ch2)=14.
-          count=count+1
-        endif
-        if (noms(iq).eq."1ch2") then
-          igcm_1ch2=iq
-          mmol(igcm_1ch2)=14.
-          count=count+1
-        endif
-        if (noms(iq).eq."cho") then
-          igcm_cho=iq
-          mmol(igcm_cho)=29.
-          count=count+1
-        endif
-        if (noms(iq).eq."ch2o") then
-          igcm_ch2o=iq
-          mmol(igcm_ch2o)=30.
-          count=count+1
-        endif
-        if (noms(iq).eq."ch3o") then
-          igcm_ch3o=iq
-          mmol(igcm_ch3o)=31.
-          count=count+1
-        endif
-        if (noms(iq).eq."c") then
-          igcm_c=iq
-          mmol(igcm_c)=12.
-          count=count+1
-        endif
-        if (noms(iq).eq."c2") then
-          igcm_c2=iq
-          mmol(igcm_c2)=24.
-          count=count+1
-        endif
-        if (noms(iq).eq."c2h") then
-          igcm_c2h=iq
-          mmol(igcm_c2h)=25.
-          count=count+1
-        endif
-        if (noms(iq).eq."c2h2") then
-          igcm_c2h2=iq
-          mmol(igcm_c2h2)=26.
-          count=count+1
-        endif
-        if (noms(iq).eq."c2h3") then
-          igcm_c2h3=iq
-          mmol(igcm_c2h3)=27.
-          count=count+1
-        endif
-        if (noms(iq).eq."c2h4") then
-          igcm_c2h4=iq
-          mmol(igcm_c2h4)=28.
-          count=count+1
-        endif
-        if (noms(iq).eq."c2h6") then
-          igcm_c2h6=iq
-          mmol(igcm_c2h6)=30.
-          count=count+1
-        endif
-        if (noms(iq).eq."ch2co") then
-          igcm_ch2co=iq
-          mmol(igcm_ch2co)=42.
-          count=count+1
-        endif
-        if (noms(iq).eq."ch3co") then
-          igcm_ch3co=iq
-          mmol(igcm_ch3co)=43.
-          count=count+1
-        endif
-        if (noms(iq).eq."hcaer") then
-          igcm_hcaer=iq
-          mmol(igcm_hcaer)=50.
-          count=count+1
-        endif
-
+          write(*,*) 'Tracer ',count,' = n2'
+        endif
+        if (noms(iq).eq."ch4_gas") then
+          igcm_ch4_gas=iq
+          mmol(igcm_ch4_gas)=16.
+          count=count+1
+          write(*,*) 'Tracer ',count,' = ch4 gas'
+        endif
+        if (noms(iq).eq."ch4_ice") then
+          igcm_ch4_ice=iq
+          mmol(igcm_ch4_ice)=16.
+          radius(iq)=3.e-6
+          rho_q(iq)=rho_ch4_ice
+          count=count+1
+          write(*,*) 'Tracer ',count,' = ch4 ice'
+        endif
+        if (noms(iq).eq."co_gas") then
+          igcm_co_gas=iq
+          mmol(igcm_co_gas)=28.
+          count=count+1
+          write(*,*) 'Tracer ',count,' = co gas'
+        endif
+        if (noms(iq).eq."co_ice") then
+          igcm_co_ice=iq
+          mmol(igcm_co_ice)=28.
+          radius(iq)=3.e-6
+          rho_q(iq)=rho_co_ice
+          count=count+1
+          write(*,*) 'Tracer ',count,' = co ice'
+        endif
+        if (noms(iq).eq."prec_haze") then
+          igcm_prec_haze=iq
+          count=count+1
+          write(*,*) 'Tracer ',count,' = prec haze'
+        endif
+        if (noms(iq).eq."haze") then
+          igcm_haze=iq
+          count=count+1
+          radius(iq)=rad_haze
+          rho_q(iq)=rho_haze
+          write(*,*) 'Tracer ',count,' = haze'
+        endif
+        if (noms(iq).eq."haze10") then
+          igcm_haze10=iq
+          count=count+1
+          radius(iq)=10.e-9
+          rho_q(iq)=rho_haze
+          write(*,*) 'Tracer ',count,' = haze10'
+        endif
+        if (noms(iq).eq."haze30") then
+          igcm_haze30=iq
+          count=count+1
+          radius(iq)=30.e-9
+          rho_q(iq)=rho_haze
+          write(*,*) 'Tracer ',count,' = haze30'
+        endif
+        if (noms(iq).eq."haze50") then
+          igcm_haze50=iq
+          count=count+1
+          radius(iq)=50.e-9
+          rho_q(iq)=rho_haze
+          write(*,*) 'Tracer ',count,' = haze50'
+        endif
+        if (noms(iq).eq."haze100") then
+          igcm_haze100=iq
+          count=count+1
+          radius(iq)=100.e-9
+          rho_q(iq)=rho_haze
+          write(*,*) 'Tracer ',count,' = haze100'
+        endif
+!       Eddy diffusion tracers
+        if (noms(iq).eq."eddy1e6") then
+          igcm_eddy1e6=iq
+          count=count+1
+          write(*,*) 'Tracer ',count,' = eddy1e6'
+        endif
+        if (noms(iq).eq."eddy1e7") then
+          igcm_eddy1e7=iq
+          count=count+1
+          write(*,*) 'Tracer ',count,' = eddy1e7'
+        endif
+        if (noms(iq).eq."eddy5e7") then
+          igcm_eddy5e7=iq
+          count=count+1
+          write(*,*) 'Tracer ',count,' = eddy5e7'
+        endif
+        if (noms(iq).eq."eddy1e8") then
+          igcm_eddy1e8=iq
+          count=count+1
+          write(*,*) 'Tracer ',count,' = eddy1e8'
+        endif
+        if (noms(iq).eq."eddy5e8") then
+          igcm_eddy5e8=iq
+          count=count+1
+          write(*,*) 'Tracer ',count,' = eddy5e8'
+        endif
       enddo ! of do iq=1,nq
 
       ! 3. find condensable traceurs different from h2o and n2
       do iq=1,nq
-        if ((index(noms(iq),"vap") .ne. 0) .and. (index(noms(iq),"h2o") .eq. 0) .and. (index(noms(iq),"n2") .eq. 0)) then
-          count=count+1
-        endif
-        if ((index(noms(iq),"ice") .ne. 0) .and. (index(noms(iq),"h2o") .eq. 0) .and. (index(noms(iq),"n2") .eq. 0)) then
+        if ((index(noms(iq),"vap") .ne. 0) .and. (index(noms(iq),"n2") .eq. 0)) then
+          count=count+1
+        endif
+        if ((index(noms(iq),"ice") .ne. 0) .and. (index(noms(iq),"n2") .eq. 0)) then
           count=count+1
         endif
 
       enddo ! of do iq=1,nq
-      
+
       ! check that we identified all tracers:
       if (count.ne.nq) then
@@ -456 +354 @@
       if (is_master) then
        rewind(407)
-       do 
+       do
         read(407,'(A)') tracline
-        if (index(tracline,"#") .ne. 1) then 
+        if (index(tracline,"#") .ne. 1) then
           exit
         endif
@@ -502 +400 @@
       n_rgcs = sum(is_rgcs)
       write(*,*)'Number of Radiative Generic Condensable Species is n_rgcs = ',n_rgcs
-      if (n_rgcs> ngt/2) then 
+      if (n_rgcs> ngt/2) then
         write(*,*) 'You have more Radiative Generic Condensable Species than Generic Condensable Species'
         write(*,*)'This is not possible: check your Modern traceur.def'
@@ -523 +421 @@
       lw_ch4=586700.
       lw_n2=2.5e5
-      
+
       if (haze) then
         ! the sedimentation radius remains radius(igcm_haze)
@@ -530 +428 @@
         else
            nmono=1
-        endif  
+        endif
 
         ia=0
@@ -539 +437 @@
 
            block=0  ! Only one type of haze is active : the first one set in traceur.def
-           do iq=1,nqmx
+           do iq=1,nq
              tracername=noms(iq)
-             write(*,*) "--> tracername ",iq,'/',nqmx,' = ',tracername
+             write(*,*) "--> tracername ",iq,'/',nq,' = ',tracername
              if (tracername(1:4).eq."haze".and.block.eq.0) then
                i_haze=iq
                block=1
                write(*,*) "i_haze=",i_haze
-               write(*,*) "Careful: if you set many haze traceurs in 
-    & traceur.def,only ",tracername," will be radiatively active 
-    & (first one in traceur.def)"
+               write(*,*) "Careful: if you set many haze traceurs in&
+     traceur.def,only ",tracername," will be radiatively active&
+     (first one in traceur.def)"
              endif
            enddo
@@ -563 +461 @@
       write(*,*) 'alpha_devil = ', alpha_devil
       write(*,*) 'radius  = ', radius
-      write(*,*) 'Qext  = ', qext 
+      write(*,*) 'Qext  = ', qext
       write(*,*)
 
@@ -579 +477 @@
           integer,           intent(in) :: iq       ! tracer index
           character(len=500),intent(in) :: tracline ! traceur.def lines with parameters
-  
+
           read(tracline,*) noms(iq)
           ! JVO 20 : We should add a sanity check aborting when duplicates in names !
@@ -593 +491 @@
                  mmol(iq)
           end if
-          ! option aki 
+          ! option aki
           if (index(tracline,'aki='   ) /= 0) then
               read(tracline(index(tracline,'aki=')+len('aki='):),*) &
@@ -603 +501 @@
                   aki(iq)
           end if
-          ! option cpi 
+          ! option cpi
           if (index(tracline,'cpi='   ) /= 0) then
               read(tracline(index(tracline,'cpi=')+len('cpi='):),*) &
@@ -613 +511 @@
                   cpi(iq)
           end if
-          ! option is_chim 
+          ! option is_chim
           if (index(tracline,'is_chim='   ) /= 0) then
           read(tracline(index(tracline,'is_chim=')+len('is_chim='):),*) &
@@ -623 +521 @@
                   is_chim(iq)
           end if
-          ! option is_rad 
+          ! option is_rad
           if (index(tracline,'is_rad=') /= 0) then
               read(tracline(index(tracline,'is_rad=') &
@@ -668 +566 @@
           end if
           !option is_condensable (LT)
-          if (index(tracline,'is_condensable=') /=0) then 
+          if (index(tracline,'is_condensable=') /=0) then
             read(tracline(index(tracline,'is_condensable=') &
               +len('is_condensable='):),*) is_condensable(iq)
             write(*,*) ' Parameter value (traceur.def) :'// &
-              ' is_condensable=', is_condensable(iq)       
+              ' is_condensable=', is_condensable(iq)
           else
               write(*,*) ' Parameter value (default)     :'// &
-              ' is_condensable=', is_condensable(iq)       
+              ' is_condensable=', is_condensable(iq)
           endif
           !option radius
-          if (index(tracline,'radius=') .ne. 0) then 
+          if (index(tracline,'radius=') .ne. 0) then
             read(tracline(index(tracline,'radius=') &
               +len('radius='):),*) radius(iq)
@@ -685 +583 @@
           else
             write(*,*) ' Parameter value (default)     :'// &
-            ' radius=', radius(iq)," m"     
+            ' radius=', radius(iq)," m"
           endif
           !option rho
-          if (index(tracline,'rho=') .ne. 0) then 
+          if (index(tracline,'rho=') .ne. 0) then
             read(tracline(index(tracline,'rho=') &
               +len('rho='):),*) rho_q(iq)
@@ -695 +593 @@
           else
             write(*,*) ' Parameter value (default)     :'// &
-              ' rho=', rho_q(iq)       
+              ' rho=', rho_q(iq)
           endif
-          !option is_rgcs 
-          if (index(tracline,'is_rgcs') .ne. 0) then 
+          !option is_rgcs
+          if (index(tracline,'is_rgcs') .ne. 0) then
             read(tracline(index(tracline,'is_rgcs=') &
               +len('is_rgcs='):),*) is_rgcs(iq)

trunk/LMDZ.PLUTO/libf/phypluto/phys_state_var_mod.F90

@@ -13 +13 @@
       use comsaison_h, only: mu0, fract
       use radcommon_h, only: glat
-!      use slab_ice_h, only : noceanmx
-      ! USE ocean_slab_mod, ONLY: nslay
       USE radinc_h, only : L_NSPECTI, L_NSPECTV,naerkind
       use surfdat_h, only: phisfi, albedodat,  &
@@ -73 +71 @@
       real,dimension(:,:,:),allocatable,save :: int_dtauv   ! VI optical thickness of layers within narrowbands for diags ().
       real,dimension(:,:,:),allocatable,save :: int_dtaui   ! IR optical thickness of layers within narrowbands for diags ().
-!$OMP THREADPRIVATE(int_dtaui,int_dtauv) 
+!$OMP THREADPRIVATE(int_dtaui,int_dtauv)
 
       real,allocatable,dimension(:),save :: tau_col ! Total Aerosol Optical Depth.
@@ -100 +98 @@
       real,dimension(:,:),allocatable,save :: lscaledEz ! Heat from largescale (W.m-2)
 !$OMP THREADPRIVATE(madjdE,madjdEz,lscaledE,lscaledEz)
-      
+
       real,dimension(:),allocatable,save :: genericconddE ! Heat from generic condensation (W.m-2)
 !$OMP THREADPRIVATE(genericconddE)
@@ -122 +120 @@
 !zpic(:)   ! Maximum de l'OESM
 !zval(:)   ! Minimum de l'OESM
-!rugoro(:) ! longueur de rugosite de l'OESM 
+!rugoro(:) ! longueur de rugosite de l'OESM
 
         ALLOCATE(phisfi(klon))

trunk/LMDZ.PLUTO/libf/phypluto/physiq_mod.F90

@@ -1 +1 @@
       module physiq_mod
-      
+
       implicit none
-      
+
       contains
-      
+
       subroutine physiq(ngrid,nlayer,nq,   &
                   firstcall,lastcall,      &
@@ -14 +14 @@
 
 !!
-      use write_field_phy, only: Writefield_phy  
+      use write_field_phy, only: Writefield_phy
 !!
-      use ioipsl_getin_p_mod, only: getin_p 
+      use ioipsl_getin_p_mod, only: getin_p
       use radinc_h, only : L_NSPECTI,L_NSPECTV,naerkind, corrkdir, banddir
       use generic_cloud_common_h, only : epsi_generic, Psat_generic
@@ -23 +23 @@
       use radcommon_h, only: sigma, glat, grav, BWNV, WNOI, DWNI, DWNV, WNOV
       use suaer_corrk_mod, only: suaer_corrk
-      use radii_mod, only: n2_reffrad
-      use aerosol_mod, only: iniaerosol, iaero_n2 
+      use radii_mod, only: su_aer_radii,haze_reffrad_fix
+      use aerosol_mod, only: iaero_haze, i_haze, haze_prof
       use surfdat_h, only: phisfi, zmea, zstd, zsig, zgam, zthe, &
                            dryness
@@ -33 +33 @@
       USE comgeomfi_h, only: totarea, totarea_planet
       USE tracer_h, only: noms, mmol, radius, rho_q, qext, &
+                          igcm_n2,&
                           alpha_lift, alpha_devil, qextrhor, &
-                          igcm_dustbin, &
-                          igcm_n2_ice, nesp, is_chim, is_condensable,constants_epsi_generic
+                          nesp, is_chim, is_condensable,constants_epsi_generic
       use time_phylmdz_mod, only: ecritphy, iphysiq, nday
       use phyetat0_mod, only: phyetat0
@@ -48 +48 @@
       use callkeys_mod, only: albedo_spectral_mode, calladj, calldifv, &
                               callrad, callsoil, nosurf, &
-                              n2cond, corrk, diagdtau, &
+                              aerohaze, corrk, diagdtau,&
                               diurnal, enertest, fat1au, &
                               icetstep, intheat, iradia, kastprof, &
                               lwrite, mass_redistrib, meanOLR, &
-                              nearn2cond, noseason_day, &
+                              n2cond,nearn2cond, noseason_day, &
                               season, sedimentation,generic_condensation, &
+                              aerohaze, haze_proffix, &
                               specOLR, &
                               startphy_file, testradtimes, &
@@ -67 +68 @@
       use turb_mod, only : q2,sensibFlux,turb_resolved
       use mass_redistribution_mod, only: mass_redistribution
-      use condensation_generic_mod, only: condensation_generic 
+      use condensation_generic_mod, only: condensation_generic
       use datafile_mod, only: datadir
 #ifndef MESOSCALE
@@ -85 +86 @@
 #endif
 
-#ifdef CPP_XIOS      
+#ifdef CPP_XIOS
       use xios_output_mod, only: initialize_xios_output, &
                                  update_xios_timestep, &
@@ -96 +97 @@
 
 !==================================================================
-!     
+!
 !     Purpose
 !     -------
@@ -186 +187 @@
 !           Frederic Hourdin        15/10/93
 !           Francois Forget        1994
-!           Christophe Hourdin        02/1997 
+!           Christophe Hourdin        02/1997
 !           Subroutine completely rewritten by F. Forget (01/2000)
 !           Water ice clouds: Franck Montmessin (update 06/2003)
 !           Radiatively active tracers: J.-B. Madeleine (10/2008-06/2009)
 !           New correlated-k radiative scheme: R. Wordsworth (2009)
-!           Many specifically Martian subroutines removed: R. Wordsworth (2009)       
+!           Many specifically Martian subroutines removed: R. Wordsworth (2009)
 !           Improved water cycle: R. Wordsworth / B. Charnay (2010)
 !           To F90: R. Wordsworth (2010)
 !           New turbulent diffusion scheme: J. Leconte (2012)
 !           Loops converted to F90 matrix format: J. Leconte (2012)
-!           No more ngridmx/nqmx, F90 commons and adaptation to parallel: A. Spiga (2012)
+!           No more ngrid/nq, F90 commons and adaptation to parallel: A. Spiga (2012)
 !           Purge of the code : M. Turbet (2015)
 !           Photochemical core developped by F. Lefevre: B. Charnay (2017)
@@ -217 +218 @@
       integer,intent(in) :: nlayer            ! Number of atmospheric layers.
       integer,intent(in) :: nq                ! Number of tracers.
-      
+
       logical,intent(in) :: firstcall ! Signals first call to physics.
       logical,intent(in) :: lastcall  ! Signals last call to physics.
-      
+
       real,intent(in) :: pday                  ! Number of elapsed sols since reference Ls=0.
       real,intent(in) :: ptime                 ! "Universal time", given as fraction of sol (e.g.: 0.5 for noon).
@@ -253 +254 @@
 ! -----------------
 
-!     Aerosol (dust or ice) extinction optical depth  at reference wavelength 
+!     Aerosol (dust or ice) extinction optical depth  at reference wavelength
 !     for the "naerkind" optically active aerosols:
 
@@ -264 +265 @@
       integer l,ig,ierr,iq,nw,isoil,iesp, igcm_generic_vap, igcm_generic_ice
       logical call_ice_vap_generic ! to call only one time the ice/vap pair of a tracer
-      
+
       real zls                       ! Solar longitude (radians).
       real zlss                      ! Sub solar point longitude (radians).
@@ -272 +273 @@
 
 ! VARIABLES for the thermal plume model (AF24: deleted)
-      
+
 ! TENDENCIES due to various processes :
 
@@ -281 +282 @@
       real zdtsdif(ngrid)     ! Turbdiff/vdifc routines.
       ! real zdtsurf_hyd(ngrid) ! Hydrol routine.
-            
-      ! For Atmospheric Temperatures : (K/s)    
+
+      ! For Atmospheric Temperatures : (K/s)
       real dtlscale(ngrid,nlayer)                             ! Largescale routine.
-      real dt_generic_condensation(ngrid,nlayer)              ! condensation_generic routine. 
+      real dt_generic_condensation(ngrid,nlayer)              ! condensation_generic routine.
       real zdtc(ngrid,nlayer)                                 ! Condense_n2 routine.
       real zdtdif(ngrid,nlayer)                               ! Turbdiff/vdifc routines.
@@ -290 +291 @@
       real zdtsw1(ngrid,nlayer), zdtlw1(ngrid,nlayer)         ! Callcorrk routine.
       real zdtchim(ngrid,nlayer)                              ! Calchim routine.
-                              
+
       ! For Surface Tracers : (kg/m2/s)
       real dqsurf(ngrid,nq)                 ! Cumulated tendencies.
       real zdqsurfc(ngrid)                  ! Condense_n2 routine.
+      REAL zdqsc(ngrid,nq)                  ! Condense_n2 routine.
       real zdqsdif(ngrid,nq)                ! Turbdiff/vdifc routines.
       real zdqssed(ngrid,nq)                ! Callsedim routine.
       real zdqsurfmr(ngrid,nq)              ! Mass_redistribution routine.
-                  
+
       ! For Tracers : (kg/kg_of_air/s)
       real zdqc(ngrid,nlayer,nq)      ! Condense_n2 routine.
@@ -311 +313 @@
 !$OMP THREADPRIVATE(zdqchim,zdqschim)
 
+
+      !! PLUTO variables
+      REAL zdqch4cloud(ngrid,nlayer,nq)
+      REAL zdqsch4cloud(ngrid,nq)
+      REAL zdtch4cloud(ngrid,nlayer)
+      REAL zdqcocloud(ngrid,nlayer,nq)
+      REAL zdqscocloud(ngrid,nq)
+      REAL zdtcocloud(ngrid,nlayer)
+      REAL rice_ch4(ngrid,nlayer)    ! Methane ice geometric mean radius (m)
+      REAL rice_co(ngrid,nlayer)     ! CO ice geometric mean radius (m)
+
+      REAL zdqsch4fast(ngrid)    ! used only for fast model nogcm
+      REAL zdqch4fast(ngrid)    ! used only for fast model nogcm
+      REAL zdqscofast(ngrid)    ! used only for fast model nogcm
+      REAL zdqcofast(ngrid)    ! used only for fast model nogcm
+      REAL zdqflow(ngrid,nq)
+
+      REAL zdteuv(ngrid,nlayer)    ! (K/s)
+      REAL zdtconduc(ngrid,nlayer) ! (K/s)
+      REAL zdumolvis(ngrid,nlayer)
+      REAL zdvmolvis(ngrid,nlayer)
+      real zdqmoldiff(ngrid,nlayer,nq)
+
+      ! Haze relatated tendancies
+      REAL zdqhaze(ngrid,nlayer,nq)
+      REAL zdqprodhaze(ngrid,nq)
+      REAL zdqsprodhaze(ngrid)
+      REAL zdqhaze_col(ngrid)
+      REAL zdqphot_prec(ngrid,nlayer)
+      REAL zdqphot_ch4(ngrid,nlayer)
+      REAL zdqconv_prec(ngrid,nlayer)
+      REAL zdq_source(ngrid,nlayer,nq)
+      ! Fast Haze relatated tendancies
+      REAL fluxbot(ngrid)
+      REAL gradflux(ngrid)
+      REAL fluxlym_sol_bot(ngrid)      ! Solar flux Lyman alpha ph.m-2.s-1 reaching the surface
+      REAL fluxlym_ipm_bot(ngrid)      ! IPM (Interplanetary) flux Lyman alpha ph.m-2.s-1 reaching the surface
+      REAL flym_sol(ngrid)      ! Incident Solar flux Lyman alpha ph.m-2.s-1
+      REAL flym_ipm(ngrid)      ! Incident IPM (Interplanetary) flux Lyman alpha ph.m-2.s-1
+
+
+
       REAL array_zero1(ngrid)
       REAL array_zero2(ngrid,nlayer)
-                        
+
       ! For Winds : (m/s/s)
       real zdvadj(ngrid,nlayer), zduadj(ngrid,nlayer)       ! Convadj routine.
@@ -320 +364 @@
       real zdhdif(ngrid,nlayer)                             ! Turbdiff/vdifc routines.
       real zdhadj(ngrid,nlayer)                             ! Convadj routine.
-     
+      REAL zdvc(ngrid,nlayer),zduc(ngrid,nlayer)            ! condense_n2 routine.
+
       ! For Pressure and Mass :
       real zdmassmr(ngrid,nlayer) ! Atmospheric Mass tendency for mass_redistribution (kg_of_air/m2/s).
@@ -326 +371 @@
       real zdpsrfmr(ngrid)        ! Pressure tendency for mass_redistribution routine (Pa/s).
 
-      
-    
+
+
 ! Local variables for LOCAL CALCULATIONS:
 ! ---------------------------------------
@@ -352 +397 @@
       real vmr(ngrid,nlayer)                        ! volume mixing ratio
       real time_phys
-      
+
       real ISR,ASR,OLR,GND,DYN,GSR,Ts1,Ts2,Ts3,TsS ! for Diagnostic.
-      
+
       real qcol(ngrid,nq) ! Tracer Column Mass (kg/m2).
+
+      !     Pluto adding variables
+      real vmr_ch4(ngrid)  ! vmr ch4
+      real vmr_co(ngrid)  ! vmr co
+      real rho(ngrid,nlayer) ! density
+      real zrho_ch4(ngrid,nlayer) ! density methane kg.m-3
+      real zrho_co(ngrid,nlayer) ! density CO kg.m-3
+      real zrho_haze(ngrid,nlayer) ! density haze kg.m-3
+      real zdqrho_photprec(ngrid,nlayer) !photolysis rate kg.m-3.s-1
+      real zq1temp_ch4(ngrid) !
+      real qsat_ch4(ngrid) !
+      real qsat_ch4_l1(ngrid) !
+!      CHARACTER(LEN=20) :: txt ! to temporarly store text for eddy tracers
+      real profmmr(ngrid,nlayer) ! fixed profile of haze if haze_proffix
+      real sensiblehf1(ngrid) ! sensible heat flux
+      real sensiblehf2(ngrid) ! sensible heat flux
 
 !     included by RW for H2O Manabe scheme
@@ -367 +428 @@
       real,save :: dEtotSW, dEtotsSW, dEtotLW, dEtotsLW
 !$OMP THREADPRIVATE(dEtotSW, dEtotsSW, dEtotLW, dEtotsLW)
-      
+
 !JL12 conservation test for mean flow kinetic energy has been disabled temporarily
 
-      real dtmoist_max,dtmoist_min     
+      real dtmoist_max,dtmoist_min
       real dItot, dItot_tmp, dVtot, dVtot_tmp
 
@@ -408 +469 @@
       real :: flux_sens_lat(ngrid)
       real :: qsurfint(ngrid,nq)
-#ifdef MESOSCALE 
+#ifdef MESOSCALE
       REAL :: lsf_dt(nlayer)
       REAL :: lsf_dq(nlayer)
@@ -416 +477 @@
       logical, save ::  check_physics_inputs=.false.
       logical, save ::  check_physics_outputs=.false.
-!$OPM THREADPRIVATE(check_physics_inputs,check_physics_outputs)     
+!$OPM THREADPRIVATE(check_physics_inputs,check_physics_outputs)
 
       ! Misc
       character*2 :: str2
-      character(len=10) :: tmp1 
+      character(len=10) :: tmp1
       character(len=10) :: tmp2
 !==================================================================================================
@@ -458 +519 @@
 !        Initialize aerosol indexes.
 !        ~~~~~~~~~~~~~~~~~~~~~~~~~~~
-         call iniaerosol
+         ! call iniaerosol
          ! allocate related local arrays
          ! (need be allocated instead of automatic because of "naerkind")
@@ -472 +533 @@
 !        Read 'startfi.nc' file.
 !        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-#ifndef MESOSCALE 
+#ifndef MESOSCALE
          call phyetat0(startphy_file,                                 &
                        ngrid,nlayer,"startfi.nc",0,0,nsoilmx,nq,      &
                        day_ini,time_phys,tsurf,tsoil,emis,q2,qsurf)
-                       
-#else   
+
+#else
 
          day_ini = pday
@@ -505 +566 @@
          write (*,*) 'In physiq day_ini =', day_ini
 
-!        Initialize albedo calculation. 
+!        Initialize albedo calculation.
 !        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
          albedo(:,:)=0.0
@@ -511 +572 @@
          albedo_snow_SPECTV(:)=0.0
          albedo_n2_ice_SPECTV(:)=0.0
-         call surfini(ngrid,nq,qsurf,albedo,albedo_bareground,albedo_snow_SPECTV,albedo_n2_ice_SPECTV) 
-
-!        Initialize orbital calculation. 
+         call surfini(ngrid,nq,qsurf,albedo,albedo_bareground,albedo_snow_SPECTV,albedo_n2_ice_SPECTV)
+
+!        Initialize orbital calculation.
 !        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
          call iniorbit(apoastr,periastr,year_day,peri_day,obliquit)
@@ -555 +616 @@
          endif
 
-         if(meanOLR)then          
-            call system('rm -f rad_bal.out') ! to record global radiative balance.           
-            call system('rm -f tem_bal.out') ! to record global mean/max/min temperatures.            
+         if(meanOLR)then
+            call system('rm -f rad_bal.out') ! to record global radiative balance.
+            call system('rm -f tem_bal.out') ! to record global mean/max/min temperatures.
             call system('rm -f h2o_bal.out') ! to record global hydrological balance.
          endif
@@ -563 +624 @@
 
          !! Initialize variables for water cycle ! AF24: removed
-            
+
 !        Set metallicity for GCS
 !        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
          metallicity=0.0 ! default value --- is not used here but necessary to call function Psat_generic
          call getin_p("metallicity",metallicity) ! --- is not used here but necessary to call function Psat_generic
-         
+
 !        Set some parameters for the thermal plume model !AF24: removed
 !        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-         
+
 #ifndef MESOSCALE
          if (ngrid.ne.1) then ! Note : no need to create a restart file in 1d.
@@ -580 +641 @@
 
 !!         call WriteField_phy("post_physdem_qsurf",qsurf(1:ngrid,igcm_h2o_vap),1)
-#endif         
-         if (corrk) then 
-            ! We initialise the spectral grid here instead of 
-            ! at firstcall of callcorrk so we can output XspecIR, XspecVI 
+#endif
+         if (corrk) then
+            ! We initialise the spectral grid here instead of
+            ! at firstcall of callcorrk so we can output XspecIR, XspecVI
             ! when using Dynamico
             print*, "physiq_mod: Correlated-k data base folder:",trim(datadir)
@@ -595 +656 @@
             call setspv ! Basic visible properties.
             call sugas_corrk       ! Set up gaseous absorption properties.
-            call suaer_corrk       ! Set up aerosol optical properties.
+            if (aerohaze) then
+               call suaer_corrk       ! Set up aerosol optical properties.
+            endif
          endif
 
@@ -606 +669 @@
                                      year_day,presnivs,pseudoalt,WNOI,WNOV)
 #endif
-         
+
 !!         call WriteField_phy("post_xios_qsurf",qsurf(1:ngrid,igcm_h2o_vap),1)
 
@@ -613 +676 @@
 
 !!      call WriteField_phy("post_firstcall_qsurf",qsurf(1:ngrid,igcm_h2o_vap),1)
-!!      call writediagfi(ngrid,"firstcall_post_qsurf"," "," ",2,qsurf(1:ngrid,igcm_h2o_vap))  
+!!      call writediagfi(ngrid,"firstcall_post_qsurf"," "," ",2,qsurf(1:ngrid,igcm_h2o_vap))
 
       if (check_physics_inputs) then
@@ -627 +690 @@
 ! ------------------------------------------------------
 
-#ifdef CPP_XIOS      
+#ifdef CPP_XIOS
       ! update XIOS time/calendar
       call update_xios_timestep
-#endif      
+#endif
 
       ! Initialize various variables
       ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-      
+
       if ( .not.nearn2cond ) then
          pdt(1:ngrid,1:nlayer) = 0.0
-      endif      
+      endif
       zdtsurf(1:ngrid)      = 0.0
       pdq(1:ngrid,1:nlayer,1:nq) = 0.0
@@ -644 +707 @@
       pdv(1:ngrid,1:nlayer) = 0.0
       pdpsrf(1:ngrid)       = 0.0
-      zflubid(1:ngrid)      = 0.0      
+      zflubid(1:ngrid)      = 0.0
       flux_sens_lat(1:ngrid) = 0.0
       taux(1:ngrid) = 0.0
@@ -660 +723 @@
 
       call orbite(zls,dist_star,declin,right_ascen)
-      
+
       if (diurnal) then
               zlss=-2.*pi*(zday-.5)
       else if(diurnal .eqv. .false.) then
               zlss=9999.
-      endif 
+      endif
 
 
@@ -673 +736 @@
       ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
       zzlay(1:ngrid,1:nlayer)=pphi(1:ngrid,1:nlayer)
-      do l=1,nlayer         
+      do l=1,nlayer
          zzlay(1:ngrid,l)= zzlay(1:ngrid,l)/glat(1:ngrid)
       enddo
@@ -685 +748 @@
             zzlev(ig,l)=(z1*zzlay(ig,l-1)+z2*zzlay(ig,l))/(z1+z2)
          enddo
-      enddo     
+      enddo
 
       !Altitude of top interface (nlayer+1), using the thicknesss of the level below the top one. LT22
-      
+
       zzlev(1:ngrid,nlayer+1) = 2*zzlev(1:ngrid,nlayer)-zzlev(1:ngrid,nlayer-1)
 
@@ -694 +757 @@
       ! Note : Potential temperature calculation may not be the same in physiq and dynamic...
       ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-      do l=1,nlayer         
+      do l=1,nlayer
          do ig=1,ngrid
             zpopsk(ig,l)=(pplay(ig,l)/pplev(ig,1))**rcp
@@ -745 +808 @@
          endif
 
-      endif 
+      endif
 
       if (callrad) then
          if( mod(icount-1,iradia).eq.0.or.lastcall) then
-           
+
             ! Eclipse incoming sunlight !AF24: removed
 
 !!            call writediagfi(ngrid,"corrk_pre_dqsurf"," "," ",2,dqsurf(1:ngrid,igcm_h2o_vap))
-!!            call writediagfi(ngrid,"corrk_pre_qsurf"," "," ",2,qsurf(1:ngrid,igcm_h2o_vap)) 
+!!            call writediagfi(ngrid,"corrk_pre_qsurf"," "," ",2,qsurf(1:ngrid,igcm_h2o_vap))
 
 
             if (corrk) then
-            
+
 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 ! II.a Call correlated-k radiative transfer scheme
@@ -766 +829 @@
                endif
                ! if(water) then !AF24: removed
-               
+
                if(generic_condensation) then
                   do iq=1,nq
 
                      call generic_tracer_index(nq,iq,igcm_generic_vap,igcm_generic_ice,call_ice_vap_generic)
-                     
+
                      if (call_ice_vap_generic) then ! to call only one time the ice/vap pair of a tracer
 
                         epsi_generic=constants_epsi_generic(iq)
 
-                        muvar(1:ngrid,1:nlayer)=mugaz/(1.e0+(1.e0/epsi_generic-1.e0)*pq(1:ngrid,1:nlayer,igcm_generic_vap)) 
-                        muvar(1:ngrid,nlayer+1)=mugaz/(1.e0+(1.e0/epsi_generic-1.e0)*pq(1:ngrid,nlayer,igcm_generic_vap))  
+                        muvar(1:ngrid,1:nlayer)=mugaz/(1.e0+(1.e0/epsi_generic-1.e0)*pq(1:ngrid,1:nlayer,igcm_generic_vap))
+                        muvar(1:ngrid,nlayer+1)=mugaz/(1.e0+(1.e0/epsi_generic-1.e0)*pq(1:ngrid,nlayer,igcm_generic_vap))
 
                      endif
-                  end do ! do iq=1,nq loop on tracers 
+                  end do ! do iq=1,nq loop on tracers
                ! take into account generic condensable specie (GCS) effect on mean molecular weight
-               
+
                else
                   muvar(1:ngrid,1:nlayer+1)=mugaz
-               endif 
-      
+               endif
+
                ! if(ok_slab_ocean) then !AF24: removed
-               
+
                ! standard callcorrk
                ! clearsky=.false.
@@ -799 +862 @@
                               int_dtaui,int_dtauv,                                &
                               tau_col,cloudfrac,totcloudfrac,                     &
-                              .false.,firstcall,lastcall)     
+                              .false.,firstcall,lastcall)
 
                 !GG (feb2021): Option to "artificially" decrease the raditive time scale in
                 !the deep atmosphere  press > 0.1 bar. Suggested by MT.
-                !! COEFF_RAD to be "tuned" to facilitate convergence of tendency 
-        
+                !! COEFF_RAD to be "tuned" to facilitate convergence of tendency
+
                 !coeff_rad=0.   ! 0 values, it doesn't accelerate the convergence
                 !coeff_rad=0.5
-                !coeff_rad=1.                 
+                !coeff_rad=1.
                 !do l=1, nlayer
                 !  do ig=1,ngrid
@@ -836 +899 @@
                ! Net atmospheric radiative heating rate (K.s-1)
                dtrad(1:ngrid,1:nlayer)=zdtsw(1:ngrid,1:nlayer)+zdtlw(1:ngrid,1:nlayer)
-               
-               ! Late initialization of the Ice Spectral Albedo. We needed the visible bands to do that ! 
+
+               ! Late initialization of the Ice Spectral Albedo. We needed the visible bands to do that !
                if (firstcall .and. albedo_spectral_mode) then
                   call spectral_albedo_calc(albedo_snow_SPECTV)
@@ -847 +910 @@
             else
 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-! II.c Atmosphere has no radiative effect 
+! II.c Atmosphere has no radiative effect
 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
                fluxtop_dn(1:ngrid)  = fract(1:ngrid)*mu0(1:ngrid)*Fat1AU/dist_star**2
@@ -856 +919 @@
                print*,'------------WARNING---WARNING------------' ! by MT2015.
                print*,'You are in corrk=false mode, '
-               print*,'and the surface albedo is taken equal to the first visible spectral value'               
-               
+               print*,'and the surface albedo is taken equal to the first visible spectral value'
+
                fluxsurfabs_sw(1:ngrid) = fluxtop_dn(1:ngrid)*(1.-albedo(1:ngrid,1))
                fluxrad_sky(1:ngrid)    = fluxsurfabs_sw(1:ngrid)
@@ -867 +930 @@
 
          endif ! of if(mod(icount-1,iradia).eq.0)
-       
+
 
          ! Transformation of the radiative tendencies
@@ -875 +938 @@
          fluxrad(1:ngrid)=fluxrad_sky(1:ngrid)-zplanck(1:ngrid)
          pdt(1:ngrid,1:nlayer)=pdt(1:ngrid,1:nlayer)+dtrad(1:ngrid,1:nlayer)
-         
+
          ! Test of energy conservation
          !----------------------------
@@ -908 +971 @@
 
       if (calldifv) then
-      
+
          zflubid(1:ngrid)=fluxrad(1:ngrid)+fluxgrd(1:ngrid)
 
          ! JL12 the following if test is temporarily there to allow us to compare the old vdifc with turbdiff.
          if (UseTurbDiff) then
-         
+
             call turbdiff(ngrid,nlayer,nq,                  &
                           ptimestep,capcal,                      &
@@ -924 +987 @@
 
          else
-         
+
             zdh(1:ngrid,1:nlayer)=pdt(1:ngrid,1:nlayer)/zpopsk(1:ngrid,1:nlayer)
- 
+
             call vdifc(ngrid,nlayer,nq,zpopsk,           &
                        ptimestep,capcal,lwrite,               &
@@ -955 +1018 @@
 !!         call writediagfi(ngrid,"vdifc_post_zdqsdif"," "," ",2,zdqsdif(1:ngrid,igcm_h2o_vap))
 
-         if (tracer) then 
+         if (tracer) then
            pdq(1:ngrid,1:nlayer,1:nq)=pdq(1:ngrid,1:nlayer,1:nq)+ zdqdif(1:ngrid,1:nlayer,1:nq)
            dqsurf(1:ngrid,1:nq)=dqsurf(1:ngrid,1:nq) + zdqsdif(1:ngrid,1:nq)
@@ -966 +1029 @@
          !-------------------------
          if(enertest)then
-         
+
             dEzdiff(:,:)=cpp*mass(:,:)*zdtdif(:,:)
             do ig = 1, ngrid
@@ -972 +1035 @@
                dEzdiff(ig,1)= dEzdiff(ig,1)+ sensibFlux(ig)! subtract flux to the ground
             enddo
-            
+
             call planetwide_sumval(dEdiff(:)*cell_area(:)/totarea_planet,dEtot)
             dEdiffs(:)=capcal(:)*zdtsdif(:)-zflubid(:)-sensibFlux(:)
             call planetwide_sumval(dEdiffs(:)*cell_area(:)/totarea_planet,dEtots)
             call planetwide_sumval(sensibFlux(:)*cell_area(:)/totarea_planet,AtmToSurf_TurbFlux)
-            
+
             if (is_master) then
-            
+
                if (UseTurbDiff) then
                          print*,'In TurbDiff sensible flux (atm=>surf) =',AtmToSurf_TurbFlux,' W m-2'
@@ -990 +1053 @@
                end if
             endif ! end of 'is_master'
-            
+
          ! JL12 : note that the black body radiative flux emitted by the surface has been updated by the implicit scheme but not given back elsewhere.
          endif ! end of 'enertest'
@@ -1008 +1071 @@
 !   IV. Convection :
 !-------------------
-      
+
 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~
 ! IV.a Thermal plume model : AF24: removed
 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~
-      
+
 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 ! IV.b Dry convective adjustment :
@@ -1037 +1100 @@
          zdtadj(1:ngrid,1:nlayer) = zdhadj(1:ngrid,1:nlayer)*zpopsk(1:ngrid,1:nlayer) ! for diagnostic only
 
-         if(tracer) then 
-            pdq(1:ngrid,1:nlayer,1:nq) = pdq(1:ngrid,1:nlayer,1:nq) + zdqadj(1:ngrid,1:nlayer,1:nq) 
+         if(tracer) then
+            pdq(1:ngrid,1:nlayer,1:nq) = pdq(1:ngrid,1:nlayer,1:nq) + zdqadj(1:ngrid,1:nlayer,1:nq)
          end if
 
@@ -1048 +1111 @@
 
          ! ! Test water conservation !AF24: removed
-         
+
       endif ! end of 'calladj'
 !----------------------------------------------
 !   Non orographic Gravity Waves: AF24: removed
 !---------------------------------------------
-      
+
 !-----------------------------------------------
-!   V. Carbon dioxide condensation-sublimation :
+!   V. Nitrogen condensation-sublimation :
 !-----------------------------------------------
 
       if (n2cond) then
-      
+
          if (.not.tracer) then
             print*,'We need a N2 ice tracer to condense N2'
             call abort
          endif
+
          call condense_n2(ngrid,nlayer,nq,ptimestep,                    &
-                           capcal,pplay,pplev,tsurf,pt,                  &
-                           pdt,zdtsurf,pq,pdq,                           &
-                           qsurf,zdqsurfc,albedo,emis,                   &
-                           albedo_bareground,albedo_n2_ice_SPECTV,      &
-                           zdtc,zdtsurfc,pdpsrf,zdqc)
+                           capcal,pplay,pplev,tsurf,pt,                 &
+                           pphi,pdt,pdu,pdv,zdtsurf,pu,pv,pq,pdq,       &
+                           qsurf(1,igcm_n2),albedo,emis,                &
+                           zdtc,zdtsurfc,pdpsrf,zduc,zdvc,              &
+                           zdqc,zdqsc(1,igcm_n2))
 
          pdt(1:ngrid,1:nlayer) = pdt(1:ngrid,1:nlayer)+zdtc(1:ngrid,1:nlayer)
@@ -1075 +1139 @@
 
          pdq(1:ngrid,1:nlayer,1:nq)   = pdq(1:ngrid,1:nlayer,1:nq)+ zdqc(1:ngrid,1:nlayer,1:nq)
-         dqsurf(1:ngrid,igcm_n2_ice) = dqsurf(1:ngrid,igcm_n2_ice) + zdqsurfc(1:ngrid)
+         ! dqsurf(1:ngrid,igcm_n2_ice) = dqsurf(1:ngrid,igcm_n2_ice) + zdqsurfc(1:ngrid)
 
 !!         call writediagfi(ngrid,"condense_n2_post_dqsurf"," "," ",2,dqsurf(1:ngrid,igcm_h2o_vap))
@@ -1094 +1158 @@
 
 !---------------------------------------------
-!   VI. Specific parameterizations for tracers 
+!   VI. Specific parameterizations for tracers
 !---------------------------------------------
 
       if (tracer) then
-      
+
   ! ---------------------
   !   VI.1. Water and ice !AF24: removed
@@ -1110 +1174 @@
 
          !Generic Condensation
-         if (generic_condensation) then 
+         if (generic_condensation) then
             call condensation_generic(ngrid,nlayer,nq,ptimestep,pplev,pplay,   &
                                           pt,pq,pdt,pdq,dt_generic_condensation, &
@@ -1128 +1192 @@
             end if
 
-            ! if (.not. water) then 
+            ! if (.not. water) then
                ! Compute GCS (Generic Condensable Specie) cloud fraction. For now we can not have both water cloud fraction and GCS cloud fraction
                ! Water is the priority
@@ -1150 +1214 @@
                      enddo
                   endif
-            end do ! do iq=1,nq loop on tracers 
+            end do ! do iq=1,nq loop on tracers
             ! endif ! .not. water
 
@@ -1164 +1228 @@
 
             ! if(watertest)then !AF24: removed
-           
+
             call callsedim(ngrid,nlayer,ptimestep,           &
                           pplev,zzlev,pt,pdt,pq,pdq,        &
@@ -1170 +1234 @@
 
             ! if(watertest)then !AF24: removed
-            
+
             ! Whether it falls as rain or snow depends only on the surface temperature
             pdq(1:ngrid,1:nlayer,1:nq) = pdq(1:ngrid,1:nlayer,1:nq) + zdqsed(1:ngrid,1:nlayer,1:nq)
@@ -1178 +1242 @@
 
             ! ! Test water conservation !AF24: removed
-            
+
          end if ! end of 'sedimentation'
 
@@ -1200 +1264 @@
                zdmassmr_col(ig)=SUM(zdmassmr(ig,1:nlayer))
             enddo
-            
+
             ! call writediagfi(ngrid,"mass_evap","mass gain"," ",3,zdmassmr)
             ! call writediagfi(ngrid,"mass_evap_col","mass gain col"," ",2,zdmassmr_col)
@@ -1209 +1273 @@
                                      pu,pv,pdt,zdtsurf,pdq,pdu,pdv,zdmassmr,        &
                                      zdtmr,zdtsurfmr,zdpsrfmr,zdumr,zdvmr,zdqmr,zdqsurfmr)
-         
+
             pdq(1:ngrid,1:nlayer,1:nq) = pdq(1:ngrid,1:nlayer,1:nq) + zdqmr(1:ngrid,1:nlayer,1:nq)
             dqsurf(1:ngrid,1:nq)       = dqsurf(1:ngrid,1:nq) + zdqsurfmr(1:ngrid,1:nq)
@@ -1217 +1281 @@
             pdpsrf(1:ngrid)            = pdpsrf(1:ngrid) + zdpsrfmr(1:ngrid)
             zdtsurf(1:ngrid)           = zdtsurf(1:ngrid) + zdtsurfmr(1:ngrid)
-            
+
          endif
 
@@ -1238 +1302 @@
          qsurf(1:ngrid,1:nq) = qsurf(1:ngrid,1:nq) + ptimestep*dqsurf(1:ngrid,1:nq)
 
-         ! Add qsurf to qsurf_hist, which is what we save in diagfi.nc. At the same time, we set the water 
+         ! Add qsurf to qsurf_hist, which is what we save in diagfi.nc. At the same time, we set the water
          ! content of ocean gridpoints back to zero, in order to avoid rounding errors in vdifc, rain.
          qsurf_hist(:,:) = qsurf(:,:)
@@ -1250 +1314 @@
 
 !------------------------------------------------
-!   VII. Surface and sub-surface soil temperature 
+!   VII. Surface and sub-surface soil temperature
 !------------------------------------------------
 
@@ -1256 +1320 @@
       ! ! Increment surface temperature
       ! if(ok_slab_ocean)then  !AF24: removed
-      
-      tsurf(1:ngrid)=tsurf(1:ngrid)+ptimestep*zdtsurf(1:ngrid) 
+
+      tsurf(1:ngrid)=tsurf(1:ngrid)+ptimestep*zdtsurf(1:ngrid)
 
       ! Compute soil temperatures and subsurface heat flux.
       if (callsoil) then
          call soil(ngrid,nsoilmx,.false.,lastcall,inertiedat,   &
-                   ptimestep,tsurf,tsoil,capcal,fluxgrd)     
+                   ptimestep,tsurf,tsoil,capcal,fluxgrd)
       endif
 
 
 !       if (ok_slab_ocean) then !AF24: removed
-      
+
       ! Test energy conservation
       if(enertest)then
-         call planetwide_sumval(cell_area(:)*capcal(:)*zdtsurf(:)/totarea_planet,dEtots)         
+         call planetwide_sumval(cell_area(:)*capcal(:)*zdtsurf(:)/totarea_planet,dEtots)
          if (is_master) print*,'Surface energy change                 =',dEtots,' W m-2'
       endif
@@ -1286 +1350 @@
       zu(1:ngrid,1:nlayer) = pu(1:ngrid,1:nlayer) + pdu(1:ngrid,1:nlayer)*ptimestep
       zv(1:ngrid,1:nlayer) = pv(1:ngrid,1:nlayer) + pdv(1:ngrid,1:nlayer)*ptimestep
-      
+
       !! Recast thermal plume vertical velocity array for outputs
       !! AF24: removed
-      
+
       ! Diagnostic.
       zdtdyn(1:ngrid,1:nlayer)     = (pt(1:ngrid,1:nlayer)-ztprevious(1:ngrid,1:nlayer)) / ptimestep
@@ -1352 +1416 @@
             endif
          end do
-                     
+
          if(ngrid.eq.1)then
             DYN=0.0
          endif
-         
+
          if (is_master) then
             print*,'  ISR            ASR            OLR            GND            DYN [W m^-2]'
@@ -1412 +1476 @@
          ! Generalised for arbitrary aerosols now. By LK
          reffcol(1:ngrid,1:naerkind)=0.0
-         if(n2cond.and.(iaero_n2.ne.0))then
-            call n2_reffrad(ngrid,nlayer,nq,zq,reffrad(1,1,iaero_n2))
-            do ig=1,ngrid
-               reffcol(ig,iaero_n2) = SUM(zq(ig,1:nlayer,igcm_n2_ice)*reffrad(ig,1:nlayer,iaero_n2)*mass(ig,1:nlayer))
-            enddo
-         endif
-         ! if(water.and.(iaero_h2o.ne.0))then !AF24: removed
-         
+            ! call n2_reffrad(ngrid,nlayer,nq,zq,reffrad(1,1,iaero_haze))
+         if (haze_proffix.and.i_haze.gt.0.) then
+               call haze_prof(ngrid,nlayer,zzlay,pplay,pt, &
+                                        reffrad,profmmr)
+               zdqhaze(:,:,i_haze)=(profmmr(:,:)-pq(:,:,igcm_n2)) & ! AF: TODO: replace by igcm_haze?
+                                                        /ptimestep
+              else
+               !! AF: TODO import from pluto.old?
+               ! call hazecloud(ngrid,nlayer,nq,ptimestep,&
+               !    pplay,pplev,pq,pdq,dist_sol,mu0,zfluxuv,zdqhaze,&
+               !    zdqphot_prec,zdqphot_ch4,zdqconv_prec,declin)
+              endif
+
+         DO iq=1, nq ! should be updated
+            DO l=1,nlayer
+              DO ig=1,ngrid
+                   pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqhaze(ig,l,iq)
+              ENDDO
+            ENDDO
+         ENDDO
       endif ! end of 'tracer'
 
@@ -1431 +1507 @@
 
             call generic_tracer_index(nq,iq,igcm_generic_vap,igcm_generic_ice,call_ice_vap_generic)
-            
+
             if (call_ice_vap_generic) then ! to call only one time the ice/vap pair of a tracer
-      
+
                do l = 1, nlayer
                   do ig=1,ngrid
@@ -1442 +1518 @@
 
             end if
-         
+
          end do ! iq=1,nq
 
@@ -1449 +1525 @@
 
       if (is_master) print*,'--> Ls =',zls*180./pi
-      
-      
+
+
 !----------------------------------------------------------------------
 !        Writing NetCDF file  "RESTARTFI" at the end of the run
@@ -1468 +1544 @@
          !! Update surface ice distribution to iterate to steady state if requested
          !! AF24: removed
-         
+
          ! endif
 #ifndef MESOSCALE
-         
+
          if (ngrid.ne.1) then
             write(*,*)'PHYSIQ: for physdem ztime_fin =',ztime_fin
-            
+
             call physdem1("restartfi.nc",nsoilmx,ngrid,nlayer,nq, &
                           ptimestep,ztime_fin,                    &
@@ -1493 +1569 @@
 !  which can later be used to make the statistic files of the run:
 !  if flag "callstats" from callphys.def is .true.)
-         
+
 
          call wstats(ngrid,"ps","Surface pressure","Pa",2,ps)
@@ -1503 +1579 @@
                      "Thermal IR radiative flux to space","W.m-2",2,    &
                      fluxtop_lw)
-                     
+
 !            call wstats(ngrid,"fluxsurf_sw",                               &
 !                        "Solar radiative flux to surface","W.m-2",2,       &
-!                         fluxsurf_sw_tot)                     
+!                         fluxsurf_sw_tot)
 !            call wstats(ngrid,"fluxtop_sw",                                &
 !                        "Solar radiative flux to space","W.m-2",2,         &
@@ -1527 +1603 @@
             do iq=1,nq
                call wstats(ngrid,noms(iq),noms(iq),'kg/kg',3,zq(1,1,iq))
-               call wstats(ngrid,trim(noms(iq))//'_surf',trim(noms(iq))//'_surf',  & 
+               call wstats(ngrid,trim(noms(iq))//'_surf',trim(noms(iq))//'_surf',  &
                            'kg m^-2',2,qsurf(1,iq) )
-               call wstats(ngrid,trim(noms(iq))//'_col',trim(noms(iq))//'_col',    & 
+               call wstats(ngrid,trim(noms(iq))//'_col',trim(noms(iq))//'_col',    &
                           'kg m^-2',2,qcol(1,iq) )
-                          
-!              call wstats(ngrid,trim(noms(iq))//'_reff',                          & 
-!                          trim(noms(iq))//'_reff',                                   & 
+
+!              call wstats(ngrid,trim(noms(iq))//'_reff',                          &
+!                          trim(noms(iq))//'_reff',                                   &
 !                          'm',3,reffrad(1,1,iq))
 
             end do
-            
-         endif ! end of 'tracer' 
+
+         endif ! end of 'tracer'
 
          !AF24: deleted slab ocean and water
@@ -1546 +1622 @@
             call mkstats(ierr)
          endif
-         
+
 
 #ifndef MESOSCALE
-       
+
 !-----------------------------------------------------------------------------------------------------
 !           OUTPUT in netcdf file "DIAGFI.NC", containing any variable for diagnostic
@@ -1577 +1653 @@
 
 !       ! Oceanic layers !AF24: removed
-      
+
       ! ! Thermal plume model  !AF24: removed
-      
+
 !        GW non-oro outputs  !AF24: removed
-      
+
       ! Total energy balance diagnostics
       if(callrad)then
-      
+
          !call writediagfi(ngrid,"ALB","Surface albedo"," ",2,albedo_equivalent)
          !call writediagfi(ngrid,"ALB_1st","First Band Surface albedo"," ",2,albedo(:,1))
@@ -1591 +1667 @@
          call writediagfi(ngrid,"OLR","outgoing longwave rad.","W m-2",2,fluxtop_lw)
          call writediagfi(ngrid,"shad","rings"," ", 2, fract)
-         
+
 !           call writediagfi(ngrid,"ASRcs","absorbed stellar rad (cs).","W m-2",2,fluxabs_sw1)
 !           call writediagfi(ngrid,"OLRcs","outgoing longwave rad (cs).","W m-2",2,fluxtop_lw1)
@@ -1604 +1680 @@
          call writediagfi(ngrid,"GND","heat flux from ground","W m-2",2,fluxgrd)
          ! endif
-         
+
          call writediagfi(ngrid,"DYN","dynamical heat input","W m-2",2,fluxdyn)
-         
+
       endif ! end of 'callrad'
-        
+
       if(enertest) then
-      
+
          if (calldifv) then
-         
+
             call writediagfi(ngrid,"q2","turbulent kinetic energy","J.kg^-1",3,q2)
             call writediagfi(ngrid,"sensibFlux","sensible heat flux","w.m^-2",2,sensibFlux)
-            
+
 !             call writediagfi(ngrid,"dEzdiff","turbulent diffusion heating (-sensible flux)","w.m^-2",3,dEzdiff)
 !             call writediagfi(ngrid,"dEdiff","integrated turbulent diffusion heating (-sensible flux)","w.m^-2",2,dEdiff)
 !             call writediagfi(ngrid,"dEdiffs","In TurbDiff (correc rad+latent heat) surf nrj change","w.m^-2",2,dEdiffs)
-             
-         endif
-          
+
+         endif
+
          if (corrk) then
             call writediagfi(ngrid,"dEzradsw","radiative heating","w.m^-2",3,dEzradsw)
             call writediagfi(ngrid,"dEzradlw","radiative heating","w.m^-2",3,dEzradlw)
          endif
-          
+
 !          if(watercond) then  !AF24: removed
 
@@ -1633 +1709 @@
             call writediagfi(ngrid,"genericconddE","heat from generic condensation","W m-2",2,genericconddE)
             call writediagfi(ngrid,"dt_generic_condensation","heating from generic condensation","K s-1",3,dt_generic_condensation)
-         
-         endif
-          
+
+         endif
+
       endif ! end of 'enertest'
 
         ! Diagnostics of optical thickness
         ! Warning this is exp(-tau), I let you postproc with -log to have tau itself - JVO 19
-        if (diagdtau) then                
+        if (diagdtau) then
           do nw=1,L_NSPECTV
             write(str2,'(i2.2)') nw
@@ -1657 +1733 @@
         call writediagfi(ngrid,"dtrad","radiative heating","K s-1",3,dtrad)
         call writediagfi(ngrid,"zdtdyn","Dyn. heating","T s-1",3,zdtdyn)
-        
+
         ! For Debugging.
         !call writediagfi(ngrid,'rnat','Terrain type',' ',2,real(rnat))
         !call writediagfi(ngrid,'pphi','Geopotential',' ',3,pphi)
-        
-
-      ! Output aerosols.
-      if (igcm_n2_ice.ne.0.and.iaero_n2.ne.0) &
-         call writediagfi(ngrid,'N2ice_reff','N2ice_reff','m',3,reffrad(1,1,iaero_n2))
-      if (igcm_n2_ice.ne.0.and.iaero_n2.ne.0) &
-         call writediagfi(ngrid,'N2ice_reffcol','N2ice_reffcol','um kg m^-2',2,reffcol(1,iaero_n2))
+
+
+      ! Output aerosols.!AF: TODO: write haze aerosols
+      ! if (igcm_n2_ice.ne.0.and.iaero_haze.ne.0) &
+      !    call writediagfi(ngrid,'N2ice_reff','N2ice_reff','m',3,reffrad(1,1,iaero_haze))
+      ! if (igcm_n2_ice.ne.0.and.iaero_haze.ne.0) &
+      !    call writediagfi(ngrid,'N2ice_reffcol','N2ice_reffcol','um kg m^-2',2,reffcol(1,iaero_haze))
       ! if (igcm_h2o_ice.ne.0.and.iaero_h2o.ne.0) &  !AF24: removed
-      
+
       ! Output tracers.
       if (tracer) then
-      
+
          do iq=1,nq
             call writediagfi(ngrid,noms(iq),noms(iq),'kg/kg',3,zq(1,1,iq))
-            call writediagfi(ngrid,trim(noms(iq))//'_surf',trim(noms(iq))//'_surf',  & 
+            call writediagfi(ngrid,trim(noms(iq))//'_surf',trim(noms(iq))//'_surf',  &
                              'kg m^-2',2,qsurf_hist(1,iq) )
-            call writediagfi(ngrid,trim(noms(iq))//'_col',trim(noms(iq))//'_col',    & 
+            call writediagfi(ngrid,trim(noms(iq))//'_col',trim(noms(iq))//'_col',    &
                             'kg m^-2',2,qcol(1,iq) )
-!          call writediagfi(ngrid,trim(noms(iq))//'_surf',trim(noms(iq))//'_surf',  & 
-!                          'kg m^-2',2,qsurf(1,iq) )                          
+!          call writediagfi(ngrid,trim(noms(iq))//'_surf',trim(noms(iq))//'_surf',  &
+!                          'kg m^-2',2,qsurf(1,iq) )
 
             ! if(watercond.or.CLFvarying)then  !AF24: removed
-            
+
             if(generic_condensation)then
                call writediagfi(ngrid,"rneb_generic","GCS cloud fraction (generic condensation)"," ",3,rneb_generic)
@@ -1692 +1768 @@
             ! if(generic_rain)then  !AF24: removed
             ! if((hydrology).and.(.not.ok_slab_ocean))then  !AF24: removed
-            
+
             call writediagfi(ngrid,"tau_col","Total aerosol optical depth","[]",2,tau_col)
 
          enddo ! end of 'nq' loop
-         
+
        endif ! end of 'tracer'
 
 
       ! Output spectrum.
-      if(specOLR.and.corrk)then      
+      if(specOLR.and.corrk)then
          call writediagspecIR(ngrid,"OLR3D","OLR(lon,lat,band)","W/m^2/cm^-1",3,OLR_nu)
          call writediagspecVI(ngrid,"OSR3D","OSR(lon,lat,band)","W/m^2/cm^-1",3,OSR_nu)
@@ -1713 +1789 @@
    if (.not.calldifv) comm_LATENT_HF(:)=0.0
    ! if ((tracer).and.(water)) then  !AF24: removed
-   
+
    if ((tracer).and.(generic_condensation)) then
    ! .and.(.not. water)
@@ -1724 +1800 @@
       do iq=1,nq
          call generic_tracer_index(nq,iq,igcm_generic_vap,igcm_generic_ice,call_ice_vap_generic)
-                  
+
          if (call_ice_vap_generic) then ! to call only one time the ice/vap pair of a tracer
 
@@ -1741 +1817 @@
             comm_RH(1:ngrid,1:nlayer) = RH_generic(1:ngrid,1:nlayer,iq)
 
-         endif 
-      end do ! do iq=1,nq loop on tracers 
+         endif
+      end do ! do iq=1,nq loop on tracers
 
    else
@@ -1785 +1861 @@
 
 ! XIOS outputs
-#ifdef CPP_XIOS      
+#ifdef CPP_XIOS
       ! Send fields to XIOS: (NB these fields must also be defined as
       ! <field id="..." /> in context_lmdz_physics.xml to be correctly used)
       CALL send_xios_field("ls",zls)
-      
+
       CALL send_xios_field("ps",ps)
       CALL send_xios_field("area",cell_area)
@@ -1797 +1873 @@
       CALL send_xios_field("v",zv)
       CALL send_xios_field("omega",omega)
-      
+
 !       IF (calltherm) THEN  !AF24: removed
       ! IF (water) THEN  !AF24: removed
-      
+
       CALL send_xios_field("ISR",fluxtop_dn)
       CALL send_xios_field("OLR",fluxtop_lw)
       CALL send_xios_field("ASR",fluxabs_sw)
 
-      if (specOLR .and. corrk) then 
+      if (specOLR .and. corrk) then
          call send_xios_field("OLR3D",OLR_nu)
          call send_xios_field("IR_Bandwidth",DWNI)
@@ -1818 +1894 @@
       endif
 #endif
-      
+
       if (check_physics_outputs) then
          ! Check the validity of updated fields at the end of the physics step
@@ -1825 +1901 @@
 
       icount=icount+1
-      
+
       end subroutine physiq
-      
+
    end module physiq_mod

trunk/LMDZ.PLUTO/libf/phypluto/radii_mod.F90

@@ -4 +4 @@
 !  module to centralize the radii calculations for aerosols
 !==================================================================
-      
+
 !     N2 cloud properties (initialized in inifis)
       real,save :: Nmix_n2 ! Number mixing ratio of N2 ice particles
@@ -21 +21 @@
 !     Purpose
 !     -------
-!     Compute the effective radii of liquid and icy water particles
-!     Jeremy Leconte (2012)
-!     Extended to dust, N2, NH3, 2-lay,Nlay,auroral aerosols by ??
-!     Added Radiative Generic Condensable Species effective radii
-!     calculations  (Lucas Teinturier 2022)
+!     Compute the effective radii of haze particles  (TB)
 !
-!     Authors
-!     -------
-!     Jeremy Leconte (2012)
 !
 !==================================================================
@@ -35 +28 @@
       use ioipsl_getin_p_mod, only: getin_p
       use radinc_h, only: naerkind
-      use aerosol_mod, only: iaero_back2lay, iaero_n2, iaero_dust, &
-                             iaero_h2so4, iaero_nh3, iaero_nlay, &
-                             iaero_aurora, iaero_generic, i_rgcs_ice
-      use callkeys_mod, only: size_nh3_cloud, nlayaero, aeronlay_size, &
+      use aerosol_mod, only: iaero_haze, i_haze, &
+                              iaero_generic, i_rgcs_ice
+      use callkeys_mod, only: nlayaero, aeronlay_size, &
                               aeronlay_nueff,aerogeneric
-      use tracer_h, only: radius, nqtot, is_rgcs
+      use tracer_h, only: radius, nqtot, is_rgcs, nmono
       Implicit none
 
@@ -47 +39 @@
 
       real, intent(out) :: reffrad(ngrid,nlayer,naerkind)      !aerosols radii (K)
-      real, intent(out) :: nueffrad(ngrid,nlayer,naerkind)     !variance     
+      real, intent(out) :: nueffrad(ngrid,nlayer,naerkind)     !variance
 
-      logical, save :: firstcall=.true.
-!$OMP THREADPRIVATE(firstcall)
-      integer :: iaer, ia , iq, i_rad  
+      integer :: iaer, ia , iq, i_rad
 
       do iaer=1,naerkind
-!     these values will change once the microphysics gets to work
-!     UNLESS tracer=.false., in which case we should be working with
-!     a fixed aerosol layer, and be able to define reffrad in a 
-!     .def file. To be improved!
-!                |-> Done in th n-layer aerosol case (JVO 20)
-
-         if(iaer.eq.iaero_n2)then ! N2 ice
-            reffrad(1:ngrid,1:nlayer,iaer) = 1.e-4
-            nueffrad(1:ngrid,1:nlayer,iaer) = 0.1 
+         print*, 'TB22 iaer',iaer,iaero_haze
+         if(iaer.eq.iaero_haze)then
+           ! Equivalent sphere radius
+           reffrad(1:ngrid,1:nlayer,iaer)=radius(i_haze)*nmono**(1./3.)
+           !reffrad(1:ngrid,1:nlayer,iaer) = 2.e-6 ! haze
+           nueffrad(1:ngrid,1:nlayer,iaer) = 0.02 ! haze
+           print*, 'TB22 Hello2',radius(i_haze)*nmono**(1./3.)
          endif
-
-         if(iaer.eq.iaero_dust)then ! dust
-            reffrad(1:ngrid,1:nlayer,iaer) = 1.e-5
-            nueffrad(1:ngrid,1:nlayer,iaer) = 0.1 
-         endif
- 
-         if(iaer.eq.iaero_h2so4)then ! H2SO4 ice
-            reffrad(1:ngrid,1:nlayer,iaer) = 1.e-6
-            nueffrad(1:ngrid,1:nlayer,iaer) = 0.1 
-         endif
-            
-         if(iaer.eq.iaero_back2lay)then ! Two-layer aerosols
-            reffrad(1:ngrid,1:nlayer,iaer) = 2.e-6
-            nueffrad(1:ngrid,1:nlayer,iaer) = 0.1 
-         endif
-
-
-         if(iaer.eq.iaero_nh3)then ! Nh3 cloud
-            reffrad(1:ngrid,1:nlayer,iaer) = size_nh3_cloud
-            nueffrad(1:ngrid,1:nlayer,iaer) = 0.1 
-         endif
-
-         do ia=1,nlayaero
-            if(iaer.eq.iaero_nlay(ia))then ! N-layer aerosols
-               reffrad(1:ngrid,1:nlayer,iaer) = aeronlay_size(ia)
-               nueffrad(1:ngrid,1:nlayer,iaer) = aeronlay_nueff(ia) 
-            endif
-         enddo
-
-         if(iaer.eq.iaero_aurora)then ! Auroral aerosols
-            reffrad(1:ngrid,1:nlayer,iaer) = 3.e-7
-            nueffrad(1:ngrid,1:nlayer,iaer) = 0.1 
-         endif
-
-         do ia=1,aerogeneric     ! Radiative Generic Condensable Species
-            if (iaer .eq. iaero_generic(ia)) then 
-               i_rad = i_rgcs_ice(ia)
-               reffrad(1:ngrid,1:nlayer,iaer)=radius(i_rad)
-               nueffrad(1:ngrid,1:nlayer,iaer) = 0.1
-            endif 
-         enddo  ! generic radiative condensable aerosols
-         
-      enddo ! iaer=1,naerkind
+      enddo
 
    end subroutine su_aer_radii
@@ -113 +59 @@
 
 !==================================================================
-   subroutine n2_reffrad(ngrid,nlayer,nq,pq,reffrad)
-!==================================================================
+   subroutine haze_reffrad_fix(ngrid,nlayer,zzlay,reffrad,nueffrad)
+      !==================================================================
 !     Purpose
 !     -------
-!     Compute the effective radii of n2 ice particles !AF24: copied from CO2
-!
-!     Authors
-!     -------
-!     Jeremy Leconte (2012)
+!     Compute the effective radii of haze particles
+!     fixed profile of radius (TB)
 !
 !==================================================================
-      USE tracer_h, only:igcm_n2_ice,rho_n2
+      use radinc_h, only: naerkind
+      USE tracer_h, only:rho_n2,nmono
       use comcstfi_mod, only: pi
+      use aerosol_mod, only: iaero_haze, i_haze
+      use datafile_mod
       Implicit none
 
-      integer,intent(in) :: ngrid,nlayer,nq
+      integer,intent(in) :: ngrid,nlayer
+      real,intent(in) :: zzlay(ngrid,nlayer)
+      real, intent(out) :: reffrad(ngrid,nlayer,naerkind)      ! haze particles radii (m)
+      real, intent(out) :: nueffrad(ngrid,nlayer,naerkind)     !
 
-      real, intent(in) :: pq(ngrid,nlayer,nq) !tracer mixing ratios (kg/kg)
-      real, intent(out) :: reffrad(ngrid,nlayer)      !n2 ice particles radii (m)
+      real :: zrad
+      real,external :: CBRT
 
-      integer :: ig,l
-      real :: zrad   
-      real,external :: CBRT            
-            
-      
+      logical, save :: firstcall=.true.
+      !$OMP THREADPRIVATE(firstcall)
 
-      if (radfixed) then
-         reffrad(1:ngrid,1:nlayer) = 5.e-5 ! N2 ice
-      else
-         do l=1,nlayer
-            do ig=1,ngrid
-               zrad = CBRT( 3*pq(ig,l,igcm_n2_ice)/(4*Nmix_n2*pi*rho_n2) )
-               reffrad(ig,l) = min(max(zrad,1.e-6),100.e-6)
-            enddo
-         enddo      
-      end if
+!     Local variables
+      integer :: iaer,l,ifine,ig
+      real :: radvec(ngrid,nlayer)
 
-   end subroutine n2_reffrad
+      !!read altitudes and radius
+      integer Nfine
+      !parameter(Nfine=21)
+      parameter(Nfine=701)
+      character(len=100) :: file_path
+      real,save :: levdat(Nfine),raddat(Nfine)
+
+!---------------- INPUT ------------------------------------------------
+
+      IF (firstcall) then
+         firstcall=.false.
+         file_path=trim(datadir)//'/haze_prop/hazerad.txt'
+         open(223,file=file_path,form='formatted')
+         do ifine=1,Nfine
+            read(223,*) levdat(ifine), raddat(ifine)
+         enddo
+         close(223)
+         print*, 'TB22 READ HAZERAD'
+       ENDIF
+
+       ! in radii mod levs has been put in km
+       DO ig=1,ngrid
+         CALL interp_line(levdat,raddat,Nfine,zzlay(ig,:)/1000,radvec(ig,:),nlayer)
+       enddo
+
+       do iaer=1,naerkind
+             if(iaer.eq.iaero_haze)then
+                  ! spherical radius or eq spherical radius
+                  ! TB22: fractal has no impact on reffrad if haze_radproffix
+                  do ig=1,ngrid
+                     do l=1,nlayer
+                        reffrad(ig,l,iaer)=radvec(ig,l)*1.e-9    !  nm => m
+                     enddo
+                  enddo
+                  nueffrad(1:ngrid,1:nlayer,iaer) = 0.02 ! haze
+             endif
+       enddo
+
+   end subroutine haze_reffrad_fix
 !==================================================================
 
 
-
-!==================================================================
-   subroutine dust_reffrad(ngrid,nlayer,reffrad)
-!==================================================================
-!     Purpose
-!     -------
-!     Compute the effective radii of dust particles
-!
-!     Authors
-!     -------
-!     Jeremy Leconte (2012)
-!
-!==================================================================
-      Implicit none
-
-      integer,intent(in) :: ngrid
-      integer,intent(in) :: nlayer
-
-      real, intent(out) :: reffrad(ngrid,nlayer)      !dust particles radii (m)
-            
-      reffrad(1:ngrid,1:nlayer) = 2.e-6 ! dust
-
-   end subroutine dust_reffrad
-!==================================================================
-
-
-!==================================================================
-   subroutine h2so4_reffrad(ngrid,nlayer,reffrad)
-!==================================================================
-!     Purpose
-!     -------
-!     Compute the effective radii of h2so4 particles
-!
-!     Authors
-!     -------
-!     Jeremy Leconte (2012)
-!
-!==================================================================
-      Implicit none
-
-      integer,intent(in) :: ngrid
-      integer,intent(in) :: nlayer
-
-      real, intent(out) :: reffrad(ngrid,nlayer)      !h2so4 particle radii (m)
-                
-      reffrad(1:ngrid,1:nlayer) = 1.e-6 ! h2so4
-
-   end subroutine h2so4_reffrad
-!==================================================================
-
-!==================================================================
-   subroutine back2lay_reffrad(ngrid,reffrad,nlayer,pplev)
-!==================================================================
-!     Purpose
-!     -------
-!     Compute the effective radii of particles in a 2-layer model
-!
-!     Authors
-!     -------
-!     Sandrine Guerlet (2013)
-!
-!==================================================================
-      use callkeys_mod, only: pres_bottom_tropo,pres_top_tropo,size_tropo,  &
-                              pres_bottom_strato,size_strato
- 
-      Implicit none
-
-      integer,intent(in) :: ngrid
-
-      real, intent(out) :: reffrad(ngrid,nlayer)      ! particle radii (m)
-      REAL,INTENT(IN) :: pplev(ngrid,nlayer+1) ! inter-layer pressure (Pa)
-      INTEGER,INTENT(IN) :: nlayer ! number of atmospheric layers
-      REAL :: expfactor
-      INTEGER l,ig
-            
-      reffrad(:,:)=1e-6  !!initialization, not important
-          DO ig=1,ngrid
-            DO l=1,nlayer-1
-              IF (pplev(ig,l) .le. pres_bottom_tropo .and. pplev(ig,l) .ge. pres_top_tropo) THEN
-                reffrad(ig,l) = size_tropo
-              ELSEIF (pplev(ig,l) .lt. pres_top_tropo .and. pplev(ig,l) .gt. pres_bottom_strato) THEN
-                expfactor=log(size_strato/size_tropo) / log(pres_bottom_strato/pres_top_tropo)
-                reffrad(ig,l)= size_tropo*((pplev(ig,l)/pres_top_tropo)**expfactor)
-              ELSEIF (pplev(ig,l) .le. pres_bottom_strato) then
-                reffrad(ig,l) = size_strato
-              ENDIF
-            ENDDO
-          ENDDO
-
-   end subroutine back2lay_reffrad
-!==================================================================
-
 end module radii_mod
 !==================================================================

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

@@ -16 +16 @@
       use radcommon_h, only: radiustab,nsize,tstellar
       use radcommon_h, only: qrefvis,qrefir,omegarefir !,omegarefvis
-      use aerosol_mod, only: noaero,iaero_n2,iaero_dust,iaero_h2so4
-      use aerosol_mod, only: iaero_back2lay,iaero_nh3,iaero_nlay,iaero_aurora
-      use aerosol_mod, only: iaero_generic,i_rgcs_ice
+      use aerosol_mod, only: noaero,iaero_haze
       use callkeys_mod, only: tplanet, optprop_back2lay_vis, optprop_back2lay_ir, &
                               optprop_aeronlay_vis, optprop_aeronlay_ir,          &
                               aeronlay_lamref, nlayaero,aerogeneric
       use tracer_h, only: noms
-      
+
       use mod_phys_lmdz_para, only : is_master, bcast
 
@@ -41 +39 @@
 !     omeg(nwvl)                     omegav   omegavis(L_NSPECTV)
 !     gfactor(nwvl)                  gav      gvis(L_NSPECTV)
-!     
+!
 !     Authors
-!     ------- 
+!     -------
 !     Richard Fournier (1996) Francois Forget (1996)
 !     Frederic Hourdin
@@ -54 +52 @@
 !     - Add water ice clouds
 !     MODIF R. Wordsworth (2009)
-!     - generalisation to arbitrary spectral bands 
+!     - generalisation to arbitrary spectral bands
 !
 !==================================================================
@@ -100 +98 @@
       REAL omegavIR(L_NSPECTI)          ! Average single scattering albedo
       REAL gavIR(L_NSPECTI)             ! Average assymetry parameter
-      
+
       logical forgetit                  ! use francois' data?
       integer iwvl, ia
@@ -109 +107 @@
 
 !---- Please indicate the names of the optical property files below
-!     Please also choose the reference wavelengths of each aerosol      
+!     Please also choose the reference wavelengths of each aerosol
 
 !--------- README TO UNDERSTAND WHAT FOLLOWS (JVO 20) -------
@@ -136 +134 @@
         print*, 'naerkind= 0'
       endif
-      
-      
+
+
       do iaer=1,naerkind
          !     naerkind=1: haze
          if (iaer.eq.1) then
 
-         ! Only one table of optical properties : 
+         ! Only one table of optical properties :
          write(*,*)'Suaer haze optical properties, using: ', &
                            TRIM(hazeprop)
@@ -150 +148 @@
          ! infrared
          file_id(iaer,2)=file_id(iaer,1)
-         
+
          lamrefvis(iaer)=0.607E-6   ! reference wavelength for opacity vis pivot wavelength Cheng et al 2017
          lamrefir(iaer)=2.E-6   !  reference wavelength for opacity IR (in the LEISA range)
 
          endif
-      enddo   
+      enddo
 
       IF (naerkind .gt. 1) THEN
@@ -162 +160 @@
          call abort
       ENDIF
-      
+
 !------------------------------------------------------------------
 
@@ -185 +183 @@
 !!!!$OMP MASTER
       if (is_master) then
-          
+
             INQUIRE(FILE=TRIM(datadir)//'/'//TRIM(aerdir)//&
                     '/'//TRIM(file_id(iaer,idomain)),&
@@ -205 +203 @@
                           '/'//TRIM(file_id(iaer,idomain)),&
                      FORM='formatted')
-              ENDIF              
+              ENDIF
             ENDIF
             IF(.NOT.file_ok) THEN
@@ -220 +218 @@
                write(*,*)' http://www.lmd.jussieu.fr/',&
                '~lmdz/planets/LMDZ.GENERIC/datagcm/'
-               CALL ABORT 
+               CALL ABORT
             ENDIF
 
@@ -239 +237 @@
                read(file_unit,*) nsize(iaer,idomain)
                endwhile = .true.
-            CASE DEFAULT reading1_seq ! ---------------------------- 
-               WRITE(*,*) 'readoptprop: ',& 
-               'Error while loading optical properties.' 
-               CALL ABORT 
+            CASE DEFAULT reading1_seq ! ----------------------------
+               WRITE(*,*) 'readoptprop: ',&
+               'Error while loading optical properties.'
+               CALL ABORT
             END SELECT reading1_seq ! ==============================
          ENDIF
@@ -256 +254 @@
       ALLOCATE(radiusdyn(nsize(iaer,idomain))) ! radiusdyn
       ALLOCATE(ep(nwvl,nsize(iaer,idomain))) ! ep
-      ALLOCATE(omeg(nwvl,nsize(iaer,idomain))) ! omeg 
+      ALLOCATE(omeg(nwvl,nsize(iaer,idomain))) ! omeg
       ALLOCATE(gfactor(nwvl,nsize(iaer,idomain))) ! g
 
@@ -316 +314 @@
 
       jfile = jfile+1
-      IF ((idomain.NE.iaero_n2).OR.(iaer.NE.iaero_n2)) THEN
+      IF ((idomain.NE.iaero_haze).OR.(iaer.NE.iaero_haze)) THEN
          endwhile = .true.
       ENDIF
@@ -334 +332 @@
 
 !     1.5 If Francois' data, convert wvl to metres
-       if(iaer.eq.iaero_n2)then
+       if(iaer.eq.iaero_haze)then
          if(forgetit)then
          !   print*,'please sort out forgetit for naerkind>1'
@@ -440 +438 @@
       DEALLOCATE(wvl)             ! wvl
       DEALLOCATE(radiusdyn)       ! radiusdyn
-      DEALLOCATE(ep)              ! ep 
-      DEALLOCATE(omeg)            ! omeg 
+      DEALLOCATE(ep)              ! ep
+      DEALLOCATE(omeg)            ! omeg
       DEALLOCATE(gfactor)         ! g
 
@@ -447 +445 @@
   END DO                    ! Loop on idomain
 !========================================================================
- 
+
   ! cleanup
   deallocate(file_id)
 
 end subroutine suaer_corrk
-      
+
 end module suaer_corrk_mod

trunk/LMDZ.PLUTO/libf/phypluto/surfdat_h.F90

@@ -16 +16 @@
        real,allocatable,dimension(:) :: zmea,zstd,zsig,zgam,zthe
 !$OMP THREADPRIVATE(zmea,zstd,zsig,zgam,zthe)
+       real ttop
+       real,allocatable,dimension(:) :: kp ! TB ref pressure 
+       real p00
+!$OMP THREADPRIVATE(ttop,kp,p00)
+! surface properties ! TB16
+       real alb_n2b,alb_n2a,alb_ch4,alb_co,alb_tho,emis_n2b,emis_n2a
+!$OMP THREADPRIVATE(alb_n2b,alb_n2a,alb_ch4,alb_co,alb_tho,emis_n2b,emis_n2a)       
+       real emis_ch4,emis_co,emis_tho,emis_tho_eq,alb_tho_eq,alb_ch4_eq
+!$OMP THREADPRIVATE(emis_ch4,emis_co,emis_tho,emis_tho_eq,alb_tho_eq,alb_ch4_eq)       
+       real ITN2,ITCH4,ITH2O,ITN2d,ITCH4d,ITH2Od,alb_ch4_s,albspe,emispe
+!$OMP THREADPRIVATE(ITN2,ITCH4,ITH2O,ITN2d,ITCH4d,ITH2Od,alb_ch4_s,albspe,emispe)       
+       real alb_tho_spe,emis_tho_spe
+!$OMP THREADPRIVATE(alb_tho_spe,emis_tho_spe)       
 
        real,allocatable,dimension(:) :: dryness  !"Dryness coefficient" for grnd water ice sublimation

trunk/LMDZ.PLUTO/libf/phypluto/surfini.F

@@ -3 +3 @@
 
       USE surfdat_h, only: albedodat
-      USE tracer_h, only: igcm_n2_ice
+      USE tracer_h, only: igcm_haze
       use planetwide_mod, only: planetwide_maxval, planetwide_minval
       use radinc_h, only : L_NSPECTV
@@ -9 +9 @@
 
       IMPLICIT NONE
-      
-      
+
+
 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 cccccccccccccc                                                                 cccccccccccccc
 cccccccccccccc   Spectral Albedo Initialisation - Routine modified by MT2015.  cccccccccccccc
-cccccccccccccc                                                                 cccccccccccccc 
+cccccccccccccc                                                                 cccccccccccccc
 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 
@@ -56 +56 @@
 
       ! Step 3 : We modify the albedo considering some N2 at the surface. We dont take into account water ice (this is made in hydrol after the first timestep) ...
-      if (igcm_n2_ice.ne.0) then
+      if (igcm_haze.ne.0) then
          DO ig=1,ngrid
-            IF (qsurf(ig,igcm_n2_ice) .GT. 1.) THEN ! This was changed by MT2015. Condition for ~1mm of N2 ice deposit.
+            IF (qsurf(ig,igcm_haze) .GT. 1.) THEN ! This was changed by MT2015. Condition for ~1mm of N2 ice deposit.
                DO nw=1,L_NSPECTV
                   albedo(ig,nw)=albedo_n2_ice_SPECTV(nw)
                ENDDO
-            END IF   
-         ENDDO   
+            END IF
+         ENDDO
       else
          write(*,*) "surfini: No N2 ice tracer on surface  ..."
          write(*,*) "         and therefore no albedo change."
-      endif     
+      endif
       call planetwide_minval(albedo,min_albedo)
       call planetwide_maxval(albedo,max_albedo)

trunk/LMDZ.PLUTO/libf/phypluto/tracer_h.F90

@@ -14 +14 @@
 
        character*30, save, allocatable :: noms(:)   ! name of the tracer
-       real, save, allocatable :: mmol(:)     ! mole mass of tracer (g/mol) 
+       real, save, allocatable :: mmol(:)     ! mole mass of tracer (g/mol)
        real, save, allocatable :: aki(:)      ! to compute coefficient of thermal concduction if photochem
        real, save, allocatable :: cpi(:)      ! to compute cpnew in concentration.F if photochem
@@ -28 +28 @@
        real,save :: rho_dust     ! Mars dust density (kg.m-3)
        real,save :: rho_ice     ! Water ice density (kg.m-3)
+       real,save :: rho_ch4_ice     ! ch4 ice density (kg.m-3)
+       real,save :: rho_co_ice     ! co ice density (kg.m-3)
        real,save :: rho_n2     ! N2 ice density (kg.m-3)
+       real,save :: lw_ch4     ! Latent heat CH4 gas -> solid
+       real,save :: lw_co      ! Latent heat CO gas -> solid
+       real,save :: lw_n2      ! Latent heat N2 gas -> solid
+       integer,save :: nmono
        real,save :: ref_r0        ! for computing reff=ref_r0*r0 (in log.n. distribution)
 !$OMP THREADPRIVATE(noms,mmol,aki,cpi,radius,rho_q,qext,alpha_lift,alpha_devil,qextrhor, &
-        !$OMP varian,r3n_q,rho_dust,rho_ice,rho_n2,ref_r0)
+        !$OMP varian,r3n_q,rho_dust,rho_ice,rho_n2,lw_n2,ref_r0)
 
        integer, save, allocatable :: is_chim(:) ! 1 if tracer used in chemistry, else 0
@@ -59 +65 @@
 ! tracer indexes: these are initialized in initracer and should be 0 if the
 !                 corresponding tracer does not exist
-       ! dust
-       integer,save,allocatable :: igcm_dustbin(:) ! for dustbin 'dust' tracers
-       ! dust, special doubleq case
-       integer,save :: igcm_dust_mass   ! dust mass mixing ratio (for transported dust)
-       integer,save :: igcm_dust_number ! dust number mixing ratio (transported dust)
-       ! water
-       integer,save :: igcm_h2o_vap ! water vapour
-       integer,save :: igcm_h2o_ice ! water ice
-       ! chemistry:
-       integer,save :: igcm_co2
-       integer,save :: igcm_co
-       integer,save :: igcm_o
-       integer,save :: igcm_o1d
-       integer,save :: igcm_o2
-       integer,save :: igcm_o3
-       integer,save :: igcm_h
-       integer,save :: igcm_h2
-       integer,save :: igcm_oh
-       integer,save :: igcm_ho2
-       integer,save :: igcm_h2o2
+
+       !Pluto chemistry
+       integer,save :: igcm_co_gas
        integer,save :: igcm_n2
        integer,save :: igcm_ar
-       integer,save :: igcm_n
-       integer,save :: igcm_no
-       integer,save :: igcm_no2
-       integer,save :: igcm_n2d
-       integer,save :: igcm_ch4
+       integer,save :: igcm_ch4_gas ! methane gas
+       ! other tracers
+       integer,save :: igcm_ar_n2 ! for simulations using co2 +neutral gaz
+       integer,save :: igcm_ch4_ice ! methane ice
+       integer,save :: igcm_co_ice ! methane ice
+       integer,save :: igcm_prec_haze
+       integer,save :: igcm_haze
+       integer,save :: igcm_haze10
+       integer,save :: igcm_haze30
+       integer,save :: igcm_haze50
+       integer,save :: igcm_haze100
+       integer,save :: igcm_eddy1e6
+       integer,save :: igcm_eddy1e7
+       integer,save :: igcm_eddy5e7
+       integer,save :: igcm_eddy1e8
+       integer,save :: igcm_eddy5e8
 
-       integer,save :: igcm_ch3
-       integer,save :: igcm_ch
-       integer,save :: igcm_3ch2
-       integer,save :: igcm_1ch2
-       integer,save :: igcm_cho
-       integer,save :: igcm_ch2o
-       integer,save :: igcm_ch3o
-       integer,save :: igcm_c
-       integer,save :: igcm_c2
-       integer,save :: igcm_c2h
-       integer,save :: igcm_c2h2
-       integer,save :: igcm_c2h3
-       integer,save :: igcm_c2h4
-       integer,save :: igcm_c2h6
-       integer,save :: igcm_ch2co
-       integer,save :: igcm_ch3co
-       integer,save :: igcm_hcaer
-
-       ! other tracers
-       integer,save :: igcm_ar_n2 ! for simulations using n2 +neutral gaz
-       integer,save :: igcm_n2_ice ! N2 ice 
-!$OMP THREADPRIVATE(igcm_dustbin,igcm_dust_mass,igcm_dust_number,igcm_h2o_vap,igcm_h2o_ice,      &
-        !$OMP igcm_co2,igcm_co,igcm_o,igcm_o1d,igcm_o2,igcm_o3,igcm_h,igcm_h2,igcm_oh,        &
-        !$OMP igcm_ho2,igcm_h2o2,igcm_n2,igcm_ar,igcm_ar_n2,igcm_n2_ice,                         &
-       !$OMP igcm_n,igcm_no,igcm_no2,igcm_n2d,igcm_ch4,igcm_ch3,igcm_ch,igcm_3ch2,               &
-       !$OMP igcm_1ch2,igcm_cho,igcm_ch2o,igcm_ch3o,igcm_c,igcm_c2,igcm_c2h,igcm_c2h2,           &
-       !$OMP igcm_c2h3,igcm_c2h4,igcm_c2h6,igcm_ch2co,igcm_ch3co,igcm_hcaer)
+!$OMP THREADPRIVATE(igcm_co_gas,igcm_n2,igcm_ar,igcm_ch4_gas,igcm_ar_n2,&igcm_ch4_ice,igcm_co_ice,igcm_prec_haze,igcm_haze,igcm_haze10,igcm_haze30,igcm_haze50,igcm_haze100,igcm_eddy1e6,igcm_eddy1e7,igcm_eddy5e7,igcm_eddy1e8,igcm_eddy5e8)
 
        end module tracer_h

trunk/LMDZ.PLUTO/libf/phypluto/vdifc_mod.F

@@ -1 +1 @@
       module vdifc_mod
-      
+
       implicit none
-      
+
       contains
-      
-      subroutine vdifc(ngrid,nlay,nq,ppopsk,         
-     &     ptimestep,pcapcal,lecrit,                        
+
+      subroutine vdifc(ngrid,nlay,nq,ppopsk,
+     &     ptimestep,pcapcal,lecrit,
      &     pplay,pplev,pzlay,pzlev,pz0,
      &     pu,pv,ph,pq,ptsrf,pemis,pqsurf,
@@ -13 +13 @@
      &     pdqdif,pdqsdif)
 
-!      use watercommon_h, only : RLVTT, T_h2O_ice_liq, RCPD, mx_eau_sol 
+!      use watercommon_h, only : RLVTT, T_h2O_ice_liq, RCPD, mx_eau_sol
 !     &    ,Psat_water, Lcpdqsat_water
       use radcommon_h, only : sigma
       use surfdat_h, only: dryness
-      use tracer_h, only: igcm_h2o_vap, igcm_h2o_ice
       use comcstfi_mod, only: g, r, cpp, rcp
       use callkeys_mod, only: tracer,nosurf
@@ -24 +23 @@
 
 !==================================================================
-!     
+!
 !     Purpose
 !     -------
 !     Turbulent diffusion (mixing) for pot. T, U, V and tracers
-!     
+!
 !     Implicit scheme
 !     We start by adding to variables x the physical tendencies
@@ -34 +33 @@
 !
 !     x(t+1) =  x(t) + dt * (dx/dt)phys(t)  +  dt * (dx/dt)difv(t+1)
-!     
+!
 !     Authors
-!     ------- 
+!     -------
 !     F. Hourdin, F. Forget, R. Fournier (199X)
 !     R. Wordsworth, B. Charnay (2010)
-!     
+!
 !==================================================================
 
@@ -62 +61 @@
       REAL,INTENT(IN) :: pcapcal(ngrid)
       REAL,INTENT(INOUT) :: pq2(ngrid,nlay+1)
-      
+
 !     Arguments added for condensation
       REAL,INTENT(IN) :: ppopsk(ngrid,nlay)
@@ -70 +69 @@
 !     Tracers
 !     --------
-      integer,intent(in) :: nq 
+      integer,intent(in) :: nq
       real,intent(in) :: pqsurf(ngrid,nq)
-      real,intent(in) :: pq(ngrid,nlay,nq), pdqfi(ngrid,nlay,nq) 
-      real,intent(out) :: pdqdif(ngrid,nlay,nq) 
-      real,intent(out) :: pdqsdif(ngrid,nq) 
-      
+      real,intent(in) :: pq(ngrid,nlay,nq), pdqfi(ngrid,nlay,nq)
+      real,intent(out) :: pdqdif(ngrid,nlay,nq)
+      real,intent(out) :: pdqsdif(ngrid,nq)
+
 !     local
 !     -----
@@ -101 +100 @@
       SAVE firstcall
 !$OMP THREADPRIVATE(firstcall)
-      
+
 !     variables added for N2 condensation
 !     ------------------------------------
@@ -134 +133 @@
       real masse, Wtot, Wdiff
 
-      real dqsdif_total(ngrid) 
-      real zq0(ngrid) 
+      real dqsdif_total(ngrid)
+      real zq0(ngrid)
 
       forceWC=.true.
@@ -151 +150 @@
 !     PRINT*,'          acond,bcond',acond,bcond
 
-!         if(water)then
-!                iliq=igcm_h2o_vap
-!                ivap=igcm_h2o_vap
-!                iice=igcm_h2o_ice ! simply to make the code legible               
-!                                  ! to be generalised later
-!         endif
 
          firstcall=.false.
@@ -208 +201 @@
             DO ilev=1,nlay
                DO ig=1,ngrid
-                  zq(ig,ilev,iq)=pq(ig,ilev,iq) + 
+                  zq(ig,ilev,iq)=pq(ig,ilev,iq) +
      &                 pdqfi(ig,ilev,iq)*ptimestep
                ENDDO
@@ -220 +213 @@
 !
 !     Source of turbulent kinetic energy at the surface
-!     ------------------------------------------------- 
+!     -------------------------------------------------
 !     Formula is Cd_0 = (karman / log[1+z1/z0])^2
 
@@ -242 +235 @@
 
 !     Turbulent diffusion coefficients in the boundary layer
-!     ------------------------------------------------------ 
+!     ------------------------------------------------------
 
       call vdif_kc(ngrid,nlay,ptimestep,g,pzlev,pzlay
@@ -272 +265 @@
 !     donc les entrees sont /zcdv/ pour la condition a la limite sol
 !     et /zkv/ = Ku
-      
+
       CALL multipl((nlay-1)*ngrid,zkv(1,2),zb0(1,2),zb(1,2))
       CALL multipl(ngrid,zcdv,zb0,zb)
@@ -405 +398 @@
 !             z1(ig) = pcapcal(ig)*ptsrf(ig) + cpp*zb(ig,1)*zc(ig,1)
 !      &           + zdplanck(ig)*ptsrf(ig) + pfluxsrf(ig)*ptimestep
-!             z2(ig) = pcapcal(ig) + cpp*zb(ig,1)*(1.-zd(ig,1)) 
+!             z2(ig) = pcapcal(ig) + cpp*zb(ig,1)*(1.-zd(ig,1))
 !      &           +zdplanck(ig)
 !             ztsrf2(ig) = z1(ig) / z2(ig)
@@ -432 +425 @@
 !     Calculate vertical flux from the bottom to the first layer (dust)
 !     -----------------------------------------------------------------
-         do ig=1,ngrid  
+         do ig=1,ngrid
             rho(ig) = zb0(ig,1) /ptimestep
          end do
@@ -440 +433 @@
 !     Implicit inversion of q
 !     -----------------------
-         do iq=1,nq 
-
-            if (iq.ne.igcm_h2o_vap) then
+         do iq=1,nq
+
+            ! if (iq.ne.igcm_h2o_vap) then
 
                DO ig=1,ngrid
@@ -448 +441 @@
                   zcq(ig,nlay)=za(ig,nlay)*zq(ig,nlay,iq)*z1(ig)
                   zdq(ig,nlay)=zb(ig,nlay)*z1(ig)
-               ENDDO 
-            
+               ENDDO
+
                DO ilay=nlay-1,2,-1
                   DO ig=1,ngrid
@@ -460 +453 @@
                ENDDO
 
-   !             if ((water).and.(iq.eq.iice)) then
-   !                ! special case for water ice tracer: do not include
-   !                ! h2o ice tracer from surface (which is set when handling
-   !                ! h2o vapour case (see further down).
-   !                ! zb(ig,1)=0 if iq ne ivap
-   !                DO ig=1,ngrid
-   !                   z1(ig)=1./(za(ig,1)+
-   !   &                    zb(ig,2)*(1.-zdq(ig,2)))
-   !                   zcq(ig,1)=(za(ig,1)*zq(ig,1,iq)+
-   !   &                    zb(ig,2)*zcq(ig,2))*z1(ig)
-   !                ENDDO
-   !             else             ! general case
                   DO ig=1,ngrid
                      z1(ig)=1./(za(ig,1)+
@@ -497 +478 @@
                end do
 
-            endif               ! if (iq.ne.igcm_h2o_vap)
-
 !     Removed (AF24): Calculate temperature tendency including latent heat term
 !     and assuming an infinite source of water on the ground
@@ -517 +496 @@
             pdvdif(ig,ilev)=(zv(ig,ilev)-
      &           (pv(ig,ilev)))/ptimestep
-            hh = ph(ig,ilev)+pdhfi(ig,ilev)*ptimestep 
+            hh = ph(ig,ilev)+pdhfi(ig,ilev)*ptimestep
 
             pdhdif(ig,ilev)=( zh(ig,ilev)- hh )/ptimestep
@@ -525 +504 @@
       DO ig=1,ngrid  ! computing sensible heat flux (atm => surface)
          sensibFlux(ig)=cpp*zb(ig,1)/ptimestep*(zh(ig,1)-ztsrf2(ig))
-      ENDDO      
+      ENDDO
 
       if (tracer) then
@@ -561 +540 @@
 !          endif
 
-      endif 
+      endif
 
       ! if(water)then