Changeset 520 for trunk


Ignore:
Timestamp:
Feb 10, 2012, 12:33:14 PM (13 years ago)
Author:
tnavarro
Message:

10/02/12 == TN

Major update on watercycle: a smaller integration timestep is now used

in watercloud.F, sedimentation of clouds is done in watercloud instead of
callsedim.F

Temperature-dependant contact parameter in nuclea.F
No dust lifting if CO2 ice in vdif.c
Ice integrated column opacity is written in diagfi from physiq.F, instead

of aeropacity.F. Mandatory if iradia is not 1.

New definition of permanent ice in surfini.F and possibility to have an ice

cap in it in 1d.

Update in deftank: callphys.def.outliers,run.def.1d; added traceur.def.scavenging

Location:
trunk/LMDZ.MARS
Files:
14 edited

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.MARS/README

    r519 r520  
    13861386   is either something.F or something.F90. Also added that object files are
    13871387   removed from archive before compiling a new version.
     1388
     1389== 10/02/12 == TN
     1390>> Major update on watercycle: a smaller integration timestep is now used
     1391   in watercloud.F, sedimentation of clouds is done in watercloud instead of
     1392   callsedim.F
     1393>> Temperature-dependant contact parameter in nuclea.F
     1394>> No dust lifting if CO2 ice in vdif.c
     1395>> Ice integrated column opacity is written in diagfi from physiq.F, instead
     1396   of aeropacity.F. Mandatory if iradia is not 1.
     1397>> New definition of permanent ice in surfini.F and possibility to have an ice
     1398   cap in it in 1d.
     1399>> Update in deftank: callphys.def.outliers,run.def.1d; added traceur.def.scavenging
     1400
     1401

  • trunk/LMDZ.MARS/deftank/callphys.def.outliers

    r283 r520  
    8888# WATER: Compute water cycle
    8989water = .true.
     90# WATER: Microphysical scheme for water-ice clouds?
     91microphys = .true.
    9092# WATER: current permanent caps at both poles. True IS RECOMMENDED
    9193#        (with .true., North cap is a source of water and South pole
    … …  
    9597albedo_h2o_ice = 0.4
    9698# WATER: water ice thermal inertia (SI) ? STILL TBD !!
    97 # (je trouve dans les 2400 SI à 173 K)
    9899inert_h2o_ice = 1000.
    99 temptag = .false.
    100100# WATER: water ice frost minimal thickness (ie kg.m^-2, ie 0.92 mm) for albedo
    101101frost_albedo_threshold = 0.005
    102102# PHOTOCHEMISTRY: include chemical species
    103103photochem  = .false.
     104# WATER: parameter contact
     105mteta = 0.95
     106# WATER: Effective variance for sedimentation for the log-normal
     107#        distribution of ice particles ?
     108nuice_sed = 0.5
     109# WATER: water ice cloud formation coupled with sedimentation is computed
     110#        every "imicro" times per physical timestep
     111imicro = 1
    104112
    105113## Thermospheric options (relevant if tracer=T) :

  • trunk/LMDZ.MARS/deftank/run.def.1d

    r38 r520  
    3535# hybrid vertical coordinate ? (.true. for hybrid and .false. for sigma levels)
    3636hybrid=.true.
    37 
     37# Water ice cap on ground ?
     38watercaptag=.false.
    3839
    3940###### Initial atmospheric temperature profile

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

    r420 r520  
    11      SUBROUTINE aeropacity(ngrid,nlayer,nq,zday,pplay,pplev,ls,
    2      &    pq,tauscaling,tauref,tau,aerosol,reffrad,nueffrad,
     2     &    pq,tauscaling,tauref,tau,taucloudtes,aerosol,reffrad,nueffrad,
    33     &    QREFvis3d,QREFir3d,omegaREFvis3d,omegaREFir3d)
    44                                                   
    … …  
    386386          ENDDO
    387387        ENDDO
    388 c       3. Outputs
    389         IF (ngrid.NE.1) THEN
    390           CALL WRITEDIAGFI(ngridmx,'tauVIS','tauext VIS refwvl',
    391      &      ' ',2,taucloudvis)
    392           CALL WRITEDIAGFI(ngridmx,'tauTES','tauabs IR refwvl',
    393      &      ' ',2,taucloudtes)
    394           IF (callstats) THEN
    395             CALL wstats(ngridmx,'tauVIS','tauext VIS refwvl',
    396      &        ' ',2,taucloudvis)
    397             CALL wstats(ngridmx,'tauTES','tauabs IR refwvl',
    398      &        ' ',2,taucloudtes)
    399           ENDIF
    400         ELSE
    401 c         CALL writeg1d(ngrid,1,taucloudtes,'tautes','NU')
    402         ENDIF
     388c       3. Outputs -- Now done in physiq.F
     389!        IF (ngrid.NE.1) THEN
     390!          CALL WRITEDIAGFI(ngridmx,'tauVIS','tauext VIS refwvl',
     391!     &      ' ',2,taucloudvis)
     392!          CALL WRITEDIAGFI(ngridmx,'tauTES','tauabs IR refwvl',
     393!     &      ' ',2,taucloudtes)
     394!          IF (callstats) THEN
     395!            CALL wstats(ngridmx,'tauVIS','tauext VIS refwvl',
     396!     &        ' ',2,taucloudvis)
     397!            CALL wstats(ngridmx,'tauTES','tauabs IR refwvl',
     398!     &        ' ',2,taucloudtes)
     399!          ENDIF
     400!        ELSE
     401!         CALL writeg1d(ngrid,1,taucloudtes,'tautes','NU')
     402!        ENDIF
    403403c==================================================================
    404404        END SELECT aerkind

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

    r370 r520  
    22     $     emis,mu0,pplev,pplay,pt,tsurf,fract,dist_sol,igout,
    33     $     dtlw,dtsw,fluxsurf_lw,fluxsurf_sw,fluxtop_lw,fluxtop_sw,
    4      &     tauref,tau,aerosol,tauscaling,rdust,rice,nuice,co2ice)
     4     &     tauref,tau,aerosol,tauscaling,taucloudtes,rdust,rice,
     5     &     nuice,co2ice)
    56
    67       IMPLICIT NONE
    … …  
    174175
    175176      REAL tauref(ngrid), tau(ngrid,naerkind)
     177      REAL taucloudtes(ngridmx)! Cloud opacity at infrared
     178                               !   reference wavelength using
     179                               !   Qabs instead of Qext
     180                               !   (direct comparison with TES)
    176181      REAL aerosol(ngrid,nlayer,naerkind)
    177182      REAL rdust(ngridmx,nlayermx)  ! Dust geometric mean radius (m)
    … …  
    395400c     Computing aerosol optical depth in each layer:
    396401      CALL aeropacity(ngrid,nlayer,nq,zday,pplay,pplev,ls,
    397      &      pq,tauscaling,tauref,tau,aerosol,reffrad,nueffrad,
    398      &      QREFvis3d,QREFir3d,omegaREFvis3d,omegaREFir3d)
     402     &      pq,tauscaling,tauref,tau,taucloudtes,aerosol,reffrad,
     403     &      nueffrad,QREFvis3d,QREFir3d,omegaREFvis3d,omegaREFir3d)
    399404
    400405c     Starting loop on sub-domain

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

    r472 r520  
    44     &                pq, pdqfi, pdqsed,pdqs_sed,nq,
    55     &                tau, tauscaling)
     6! to use  'getin'
     7      USE ioipsl_getincom
    68      IMPLICIT NONE
    79
    … …  
    2123c   declarations:
    2224c   -------------
    23 
     25     
    2426#include "dimensions.h"
    2527#include "dimphys.h"
    … …  
    6971c     Sedimentation radius of water ice
    7072      real rsedcloud(ngridmx,nlayermx)
     73      real beta      ! correction for the shape of the ice particles (cf. newsedim)
     74      save beta
    7175c     Cloud density (kg.m-3)
    7276      real rhocloud(ngridmx,nlayermx)
    … …  
    118122      SAVE idust_mass,idust_number
    119123      SAVE iccn_mass,iccn_number
     124     
    120125
    121126c     Firstcall:
    … …  
    201206        ENDIF !of if (scavenging)
    202207
    203         IF (water) THEN
    204           write(*,*) "water_param nueff Sedimentation:", nuice_sed
    205           IF (activice) THEN
    206             write(*,*) "water_param nueff Radiative:", nuice_ref
    207           ENDIF
    208         ENDIF
     208!        IF (water) THEN
     209!         write(*,*) "correction for the shape of the ice particles ?"
     210!         beta=0.75 ! default value
     211!         call getin("ice_shape",beta)
     212!         write(*,*) " ice_shape = ",beta
     213!
     214!          write(*,*) "water_param nueff Sedimentation:", nuice_sed
     215!          IF (activice) THEN
     216!            write(*,*) "water_param nueff Radiative:", nuice_ref
     217!          ENDIF
     218!        ENDIF
    209219     
    210220        firstcall=.false.
    … …  
    269279          if ((doubleq.and.
    270280     &        ((iq.eq.idust_mass).or.
    271      &         (iq.eq.idust_number))).or.
    272      &        (scavenging.and.
    273      &        ((iq.eq.iccn_mass).or.
    274      &        (iq.eq.iccn_number)))) then
     281     &         (iq.eq.idust_number)))) then !.or.
     282!     &        (scavenging.and.
     283!     &        ((iq.eq.iccn_mass).or.
     284!     &        (iq.eq.iccn_number)))) then
    275285
    276286c           Computing size distribution:
    277287c           ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    278288
    279             if ((iq.eq.idust_mass).or.(iq.eq.idust_number)) then
     289c            if ((iq.eq.idust_mass).or.(iq.eq.idust_number)) then
    280290              do  l=1,nlay
    281291                do ig=1, ngrid
    … …  
    284294              end do
    285295              sigma0 = varian
    286             else
    287               do  l=1,nlay
    288                 do ig=1, ngrid
    289                   r0(ig,l)=r0ccn(ig,l)
    290                 end do
    291               end do
    292               sigma0 = sqrt(log(1.+nuice_sed))
    293             endif
     296c            else             ! ccn
     297c              do  l=1,nlay
     298c                do ig=1, ngrid
     299c                  r0(ig,l)=r0ccn(ig,l)
     300c                end do
     301c              end do
     302c              sigma0 = sqrt(log(1.+nuice_sed))
     303c            endif
    294304
    295305c        Computing mass mixing ratio for each particle size
    … …  
    297307          IF ((iq.EQ.idust_mass).or.(iq.EQ.iccn_mass)) then
    298308            radpower = 2
    299           ELSE
     309          ELSE  ! number
    300310            radpower = -1
    301311          ENDIF
    … …  
    348358
    349359          do ir=1,nr
    350              IF ((iq.EQ.idust_mass).or.(iq.EQ.idust_number)) then
     360!             IF ((iq.EQ.idust_mass).or.(iq.EQ.idust_number)) then
     361              ! Dust sedimentation
    351362               call newsedim(ngrid,nlay,1,1,ptimestep,
    352363     &         pplev,masse,epaisseur,pt,rd(ir),rho_dust,qr(1,1,ir),
    353364     &         wq,0.5)
    354              ELSE
    355                call newsedim(ngrid,nlay,1,ngrid*nlay,ptimestep,
    356      &         pplev,masse,epaisseur,pt,rd(ir),rhocloud,qr(1,1,ir),
    357      &         wq,1.0)
    358              ENDIF
     365!             ELSE
     366!               ! CCN sedimentation
     367!               call newsedim(ngrid,nlay,1,ngrid*nlay,ptimestep,
     368!     &         pplev,masse,epaisseur,pt,rd(ir),rhocloud,qr(1,1,ir),
     369!     &         wq,beta)
     370!             ENDIF
    359371
    360372c            Tendencies
    … …  
    371383          enddo ! of do ir=1,nr
    372384c -----------------------------------------------------------------
    373 c         WATER CYCLE CASE
    374 c -----------------------------------------------------------------
    375           else if (water.and.(iq.eq.igcm_h2o_ice)) then
    376             if (microphys) then
    377               call newsedim(ngrid,nlay,ngrid*nlay,ngrid*nlay,
    378      &        ptimestep,pplev,masse,epaisseur,pt,rsedcloud,rhocloud,
    379      &        zqi(1,1,iq),wq,1.0)
    380             else
    381               call newsedim(ngrid,nlay,ngrid*nlay,1,
    382      &        ptimestep,pplev,masse,epaisseur,pt,rsedcloud,rho_q(iq),
    383      &        zqi(1,1,iq),wq,1.0)
    384             endif ! of if (microphys)
     385c         WATER CYCLE CASE --- NOW DONE IN WATERCLOUD !!
     386c -----------------------------------------------------------------
     387!          else if (water.and.(iq.eq.igcm_h2o_ice)) then
     388!            if (microphys) then
     389!              ! water ice sedimentation
     390!              call newsedim(ngrid,nlay,ngrid*nlay,ngrid*nlay,
     391!     &        ptimestep,pplev,masse,epaisseur,pt,rsedcloud,rhocloud,
     392!     &        zqi(1,1,iq),wq,beta)
     393!            else
     394!              ! water ice sedimentation
     395!              call newsedim(ngrid,nlay,ngrid*nlay,1,
     396!     &        ptimestep,pplev,masse,epaisseur,pt,rsedcloud,rho_q(iq),
     397!     &        zqi(1,1,iq),wq,beta)
     398!            endif ! of if (microphys)
    385399c           Tendencies
    386400c           ~~~~~~~~~~
    387             do ig=1,ngrid
    388               pdqs_sed(ig,iq)=wq(ig,1)/ptimestep
    389             end do
     401!            do ig=1,ngrid
     402!              pdqs_sed(ig,iq)=wq(ig,1)/ptimestep
     403!            end do
    390404c -----------------------------------------------------------------
    391405c         GENERAL CASE
    392406c -----------------------------------------------------------------
    393           else
     407c          else
     408           else if ((iq .ne. iccn_mass) .and. (iq .ne. iccn_number)
     409     &       .and. (iq .ne. igcm_h2o_ice)) then ! because it is done in watercloud
    394410            call newsedim(ngrid,nlay,1,1,ptimestep,
    395411     &      pplev,masse,epaisseur,pt,radius(iq),rho_q(iq),
    … …  
    418434c     Update the dust particle size "rdust"
    419435c     -------------------------------------
    420       if (doubleq) then
    421        DO l = 1, nlay
     436      DO l = 1, nlay
    422437        DO ig=1,ngrid
    423438          rdust(ig,l)=
    … …  
    426441          rdust(ig,l)=min(max(rdust(ig,l),1.e-10),500.e-6)
    427442        ENDDO
    428        ENDDO
    429       endif !of if (doubleq)
     443      ENDDO
    430444     
    431445c     Update the ice particle size "rice"
    432446c     -------------------------------------
    433       IF(scavenging) THEN
    434         DO l = 1, nlay
    435           DO ig=1,ngrid
    436             Mo   = zqi(ig,l,igcm_h2o_ice) +
    437      &             zqi(ig,l,iccn_mass)* tauscaling(ig) + 1.e-30
    438             No   = zqi(ig,l,iccn_number)* tauscaling(ig)+ 1.e-30
    439             rhocloud(ig,l) = zqi(ig,l,igcm_h2o_ice) / Mo * rho_ice
    440      &                       +zqi(ig,l,iccn_mass)* tauscaling(ig)
    441      &                        / Mo * rho_dust
    442             rhocloud(ig,l) =
    443      &         min(max(rhocloud(ig,l),rho_ice),rho_dust)
    444             rice(ig,l) =
    445      &        ( Mo / No * 0.75 / pi / rhocloud(ig,l) ) **(1./3.)
    446            if ((Mo.lt.1.e-20) .or. (No.le.1)) rice(ig,l) = 1.e-8
     447!      IF(scavenging) THEN
     448!        DO l = 1, nlay
     449!          DO ig=1,ngrid
     450!            Mo   = zqi(ig,l,igcm_h2o_ice) +
     451!     &             zqi(ig,l,iccn_mass)* tauscaling(ig) + 1.e-30
     452!            No   = zqi(ig,l,iccn_number)* tauscaling(ig)+ 1.e-30
     453!            rhocloud(ig,l) = zqi(ig,l,igcm_h2o_ice) / Mo * rho_ice
     454!     &                       +zqi(ig,l,iccn_mass)* tauscaling(ig)
     455!     &                        / Mo * rho_dust
     456!            rhocloud(ig,l) =
     457!     &         min(max(rhocloud(ig,l),rho_ice),rho_dust)
     458!            rice(ig,l) =
     459!     &        ( Mo / No * 0.75 / pi / rhocloud(ig,l) ) **(1./3.)
     460!           if ((Mo.lt.1.e-15) .or. (No.le.50)) rice(ig,l) = 1.e-8
    447461c           print*, "Mice,Mo, No",zqi(ig,l,igcm_h2o_ice),Mo, No
    448462c           print*, "rice, rho apres", rice(ig,l), rhocloud(ig,l)
    449            
    450           ENDDO
    451         ENDDO
    452       ELSE
    453        if (water) then
    454         DO l = 1, nlay
    455           DO ig=1,ngrid
    456             ccntyp =
    457      &     1.3e+8*max(tau(ig,1),0.001)/0.1*exp(-zlay(ig,l)/10000.)
    458             ccntyp = ccntyp /ccn_factor
    459             rice(ig,l)=max( CBRT ( (zqi(ig,l,igcm_h2o_ice)/rho_ice
    460      &      +ccntyp*(4./3.)*pi*rdust(ig,l)**3.)
    461      &      /(ccntyp*4./3.*pi) ), rdust(ig,l))
    462           ENDDO
    463         ENDDO
    464        endif ! of if (water)
    465       ENDIF ! of IF (scavenging)
     463!           
     464!          ENDDO
     465!        ENDDO
     466!      ELSE
     467!        DO l = 1, nlay
     468!          DO ig=1,ngrid
     469!            ccntyp =
     470!     &     1.3e+8*max(tau(ig,1),0.001)/0.1*exp(-zlay(ig,l)/10000.)
     471!            ccntyp = ccntyp /ccn_factor
     472!            rice(ig,l)=max( CBRT ( (zqi(ig,l,igcm_h2o_ice)/rho_ice
     473!     &      +ccntyp*(4./3.)*pi*rdust(ig,l)**3.)
     474!     &      /(ccntyp*4./3.*pi) ), rdust(ig,l))
     475!          ENDDO
     476!        ENDDO
     477!      ENDIF
    466478
    467479      RETURN

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

    r455 r520  
    11      subroutine improvedclouds(ngrid,nlay,ptimestep,
    2      &             pplev,pplay,pt,pdt,
    3      &             pq,pdq,pdqcloud,pdqscloud,pdtcloud,
     2     &             pplev,pplay,zlev,pt,pdt,
     3     &             pq,pdq,pdqcloud,pdtcloud,
    44     &             nq,tauscaling,rdust,rice,nuice,
    55     &             rsedcloud,rhocloud)
     6! to use  'getin'
     7      USE ioipsl_getincom
    68      implicit none
    79c------------------------------------------------------------------
    … …  
    2224c  Authors: J.-B. Madeleine, based on the work by Franck Montmessin
    2325c           (October 2011)
     26c           T. Navarro, debug & correction (October-December 2011)
    2427c------------------------------------------------------------------
    2528#include "dimensions.h"
    … …  
    3740      integer nq                 ! nombre de traceurs
    3841      REAL ptimestep             ! pas de temps physique (s)
    39       REAL pplev(ngrid,nlay+1),pplay(ngrid,nlay)     ! pression au milieu des couches (Pa)
     42      REAL pplev(ngrid,nlay+1)   ! pression aux inter-couches (Pa)
     43      REAL pplay(ngrid,nlay)     ! pression au milieu des couches (Pa)
     44      REAL zlev(ngrid,nlay+1)    ! altitude at layer boundaries
     45           
    4046      REAL pt(ngrid,nlay)        ! temperature at the middle of the
    4147                                 !   layers (K)
    … …  
    5763      REAL pdqcloud(ngrid,nlay,nq) ! tendance de la condensation
    5864                                   !   H2O(kg/kg.s-1)
    59       REAL pdqscloud(ngrid,nq)     ! flux en surface (kg.m-2.s-1)
    6065      REAL pdtcloud(ngrid,nlay)    ! tendance temperature due
    6166                                   !   a la chaleur latente
    … …  
    7580
    7681      INTEGER ig,l,i
     82     
    7783
    7884      REAL zq(ngridmx,nlayermx,nqmx)  ! local value of tracers
    … …  
    106112      DOUBLE PRECISION rb_cld(nbin_cld+1)! boundary values of each rad_cld bin (m)
    107113      SAVE rb_cld
    108       DOUBLE PRECISION dr_cld(nbin_cld)! width of each rad_cld bin (m)
    109       DOUBLE PRECISION vol_cld(nbin_cld)   ! particle volume for each bin (m3)
     114      DOUBLE PRECISION dr_cld(nbin_cld)   ! width of each rad_cld bin (m)
     115      DOUBLE PRECISION vol_cld(nbin_cld)  ! particle volume for each bin (m3)
    110116
    111117      REAL sigma_ice ! Variance of the ice and CCN distributions
    112118      SAVE sigma_ice
    113 
     119     
    114120c----------------------------------     
    115121c some outputs for 1D -- TEMPORARY
    … …  
    126132      REAL Mccn_col(ngridmx)         ! total column ccn mass
    127133      REAL Nccn_col(ngridmx)         ! total column ccn mass
     134      REAL dMice_col(ngridmx)        ! total column ice mass change
     135      REAL drice_col(ngridmx)        ! total column ice radius change
     136      REAL icetot(ngridmx)        ! total column ice mass
     137      REAL rice_out(ngridmx,nlayermx) ! ice radius change
    128138      REAL rate_out(ngridmx,nlayermx) ! nucleation rate
    129139      INTEGER count
    130140     
    131141      LOGICAL output_sca     ! scavenging outputs flag for tests
     142     
    132143      output_sca = .false.
    133144c----------------------------------     
    … …  
    202213c-----------------------------------------------------------------------
    203214c     Initialize the tendencies
    204       pdqscloud(1:ngrid,1:nq)=0
    205215      pdqcloud(1:ngrid,1:nlay,1:nq)=0
    206216      pdtcloud(1:ngrid,1:nlay)=0
    … …  
    249259      enddo
    250260
    251 
    252261c------------------------------------------------------------------
    253262c     Cloud microphysical scheme
    … …  
    269278       
    270279        IF ((satu .ge. 1)!        ) THEN             ! if we have condensation
    271      &      .or. ( zq(ig,l,igcm_ccn_number)
    272      &             .ge. 10)    ) THEN    ! or sublimation
     280     &      .or. ( zq(ig,l,igcm_ccn_number)*tauscaling(ig)
     281     &             .ge. 2)    ) THEN    ! or sublimation
    273282
    274283
    … …  
    287296     &                           -derf( dlog(rb_cld(i) * Rm) * dev2 ) )
    288297        enddo
    289        
     298               
    290299!        sumcheck = 0
    291300!        do i = 1, nbin_cld
    … …  
    320329        rate_out(ig,l) = 0
    321330        do i = 1, nbin_cld
     331        ! schema simple
    322332          n_aer(i) = n_aer(i) / ( 1. + rate(i)*ptimestep )
    323333          m_aer(i) = m_aer(i) / ( 1. + rate(i)*ptimestep )
    324334          dN       = dN + n_aer(i) * rate(i) * ptimestep
    325335          dM       = dM + m_aer(i) * rate(i) * ptimestep
    326           rate_out(ig,l)=rate_out(ig,l)+rate(i)
     336          !rate_out(ig,l)=rate_out(ig,l)+rate(i)
    327337        enddo
    328        
    329 !        dN = min( max(dN,-zq(ig,l,igcm_ccn_number) ),
    330 !     &                 zq(ig,l,igcm_dust_number) )
    331 !     
    332 !        dM = min( max(dM,-zq(ig,l,igcm_ccn_mass) ),
    333 !     &                 zq(ig,l,igcm_dust_mass) )
    334338
    335339          Mo   = zq0(ig,l,igcm_h2o_ice) +
    … …  
    382386          zq(ig,l,igcm_h2o_vap) = zq0(ig,l,igcm_h2o_vap) -
    383387     &        dMice
    384 
     388     
    385389c         If all the ice particles sublimate, all the condensation
    386390c           nuclei are released:
    … …  
    411415        zq(ig,l,igcm_ccn_mass) = zq(ig,l,igcm_ccn_mass) + dM
    412416        zq(ig,l,igcm_ccn_number) = zq(ig,l,igcm_ccn_number) + dN
    413        
     417
     418               
    414419        pdqcloud(ig,l,igcm_dust_mass)=(zq(ig,l,igcm_dust_mass)
    415420     &                         -zq0(ig,l,igcm_dust_mass))/ptimestep
    … …  
    425430     &                         -zq0(ig,l,igcm_h2o_ice))/ptimestep
    426431     
     432     
    427433        count = count +1
    428434        ELSE ! for coherence (rdust, rice computations etc ..)
    … …  
    447453          lw=(2834.3-0.28*(zt(ig,l)-To)-0.004*(zt(ig,l)-To)**2)*1.e+3
    448454          pdtcloud(ig,l)=-pdqcloud(ig,l,igcm_h2o_vap)*lw/cpp
     455c          pdtcloud(ig,l)=pdt(ig,l)-pdqcloud(ig,l,igcm_h2o_vap)*lw/cpp
    449456         
    450 c----- update rice & rhocloud AFTER scavenging
     457c----- update rice & rhocloud AFTER microphysic
    451458          Mo   = zq(ig,l,igcm_h2o_ice) +
    452459     &           zq(ig,l,igcm_ccn_mass)* tauscaling(ig) + 1.e-30
    … …  
    456463     &                      / Mo * rho_dust
    457464          rhocloud(ig,l) = min(max(rhocloud(ig,l),rho_ice),rho_dust)
     465         
     466          rice_out(ig,l)=rice(ig,l)
    458467          rice(ig,l) =
    459468     &      ( Mo / No * 0.75 / pi / rhocloud(ig,l) ) **(1./3.)
    460469          if ((Mo.lt.1.e-20) .or. (No.le.1)) rice(ig,l) = 1.e-8
     470          rice_out(ig,l)=rice(ig,l)-rice_out(ig,l)
    461471         
    462472          nuice(ig,l)=nuice_ref ! used for rad. transfer calculations
    … …  
    474484        ENDDO
    475485      ENDDO
    476      
     486     
     487!      print*, 'SATU'
     488!      print*, satu_out(1,:)
    477489     
    478490c------------------------------------------------------------------
    … …  
    494506      Mccn_col(:)  = 0
    495507      Nccn_col(:)  = 0
     508      dMice_col(:) = 0
     509      drice_col(:) = 0
     510      icetot(:)    = 0
    496511      do l=1, nlay
    497512        do ig=1,ngrid
    … …  
    512527     &       zq(ig,l,igcm_ccn_number)*tauscaling(ig)
    513528     &       *(pplev(ig,l) - pplev(ig,l+1)) / g
     529          dMice_col(ig) = dMice_col(ig) +
     530     &       Mcon_out(ig,l)
     531     &       *(pplev(ig,l) - pplev(ig,l+1)) / g
     532          drice_col(ig) = drice_col(ig) +
     533     &       rice_out(ig,l)*zq(ig,l,igcm_h2o_ice)
     534     &       *(pplev(ig,l) - pplev(ig,l+1)) / g
     535          icetot(ig) = icetot(ig) +
     536     &       zq(ig,l,igcm_h2o_ice)
     537     &       *(pplev(ig,l) - pplev(ig,l+1)) / g
    514538        enddo ! of do ig=1,ngrid
    515539      enddo ! of do l=1,nlay
     540     
     541      drice_col(ig) = drice_col(ig)/icetot(ig)
    516542
    517543
    … …  
    549575         call WRITEDIAGFI(ngrid,"dN","ccn variation","#",1,
    550576     &                    dN_out)
     577         call WRITEDIAGFI(ngrid,"dicetot","ice col var","kg/m2",0,
     578     &                    dMice_col)
     579         call WRITEDIAGFI(ngrid,"dricetot","ice col var","m",0,
     580     &                    drice_col)
    551581         call WRITEDIAGFI(ngrid,"mcond","h2o condensed mass","kg",1,
    552582     &                    Mcon_out)

  • trunk/LMDZ.MARS/libf/phymars/microphys.h

    r358 r520  
    55
    66!     Number of bins
    7       INTEGER, PARAMETER :: nbin_cld = 50
     7      INTEGER, PARAMETER :: nbin_cld = 5
    88
    99!     Reference temperature, T=273.15 K

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

    r358 r520  
    4444
    4545      integer i
     46     
     47      LOGICAL firstcall
     48      DATA firstcall/.true./
     49      SAVE firstcall
    4650
    4751c     *************************************************
    4852
    4953      mtetalocal = mteta
     54
     55cccccccccccccccccccccccccccccccccccccccccccccccccc
     56ccccccccccc ESSAIS TN MTETA = F (T) cccccccccccccc
     57c      if (temp .gt. 200) then
     58c         mtetalocal = mtetalocal
     59c      else if (temp .lt. 190) then
     60c         mtetalocal = mtetalocal-0.05
     61c      else
     62c         mtetalocal = mtetalocal - (190-temp)*0.005
     63c      endif
     64c----------------exp law, see Trainer 2008, J. Phys. Chem. C 2009, 113, 2036\u20132040
     65       !mtetalocal = max(mtetalocal - 6005*exp(-0.065*temp),0.1)
     66       !mtetalocal = max(mtetalocal - 6005*exp(-0.068*temp),0.1)
     67               !print*, mtetalocal, temp
     68cccccccccccccccccccccccccccccccccccccccccccccccccc
     69cccccccccccccccccccccccccccccccccccccccccccccccccc
     70      IF (firstcall) THEN
     71          print*, ' ' 
     72          print*, 'dear user, please keep in mind that'
     73          print*, 'contact parameter IS constant'
     74          !print*, 'contact parameter IS NOT constant:'
     75          !print*, 'max(mteta - 6005*exp(-0.065*temp),0.1)'
     76          !print*, 'max(mteta - 6005*exp(-0.068*temp),0.1)'
     77          print*, ' ' 
     78         firstcall=.false.
     79      END IF
     80cccccccccccccccccccccccccccccccccccccccccccccccccc
     81cccccccccccccccccccccccccccccccccccccccccccccccccc
     82   
    5083
    5184      if (sat .gt. 1.) then    ! minimum condition to activate nucleation

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

    r506 r520  
    201201      real rsedcloud(ngridmx,nlayermx) ! Cloud sedimentation radius
    202202      real rhocloud(ngridmx,nlayermx)  ! Cloud density (kg.m-3)
    203       REAL surfdust(ngridmx,nlayermx) ! dust surface area (m2/m3, if photochemistry)
    204       REAL surfice(ngridmx,nlayermx)  ! ice surface area  (m2/m3, if photochemistry)
     203      REAL surfdust(ngridmx,nlayermx) ! dust surface area (micron2/cm3, if photochemistry)
     204      REAL surfice(ngridmx,nlayermx)  !  ice surface area (micron2/cm3, if photochemistry)
    205205
    206206c     Variables used by the slope model
    … …  
    308308      real rho(ngridmx,nlayermx)  ! density
    309309      real vmr(ngridmx,nlayermx)  ! volume mixing ratio
    310       REAL colden(ngridmx,nqmx)   ! vertical column of tracers
     310      !real colden(ngridmx,nqmx)   ! vertical column                      !FL
    311311      REAL mtot(ngridmx)          ! Total mass of water vapor (kg/m2)
    312312      REAL icetot(ngridmx)        ! Total mass of water ice (kg/m2)
     313      REAL ccntot(ngridmx)        ! Total number of ccn (nbr/m2)
    313314      REAL rave(ngridmx)          ! Mean water ice effective radius (m)
    314315      REAL opTES(ngridmx,nlayermx)! abs optical depth at 825 cm-1
    315316      REAL tauTES(ngridmx)        ! column optical depth at 825 cm-1
    316317      REAL Qabsice                ! Water ice absorption coefficient
     318      REAL taucloudtes(ngridmx)! Cloud opacity at infrared
     319                               !   reference wavelength using
     320                               !   Qabs instead of Qext
     321                               !   (direct comparison with TES)
     322
    317323
    318324c Test 1d/3d scavenging
    319       real h2o_tot
     325      real h2otot(ngridmx)
     326      REAL satu(ngridmx,nlayermx) ! satu ratio for output
     327      REAL zqsat(ngridmx,nlayermx)    ! saturation
    320328
    321329      REAL time_phys
    … …  
    332340
    333341      REAL lmax_th_out(ngridmx),zmax_th(ngridmx)
    334       REAL, SAVE :: wstar(ngridmx)
    335       REAL, SAVE :: hfmax_th(ngridmx)
     342      REAL, SAVE :: wmax_th(ngridmx)
     343      REAL hfmax_th(ngridmx)
    336344      REAL pdu_th(ngridmx,nlayermx),pdv_th(ngridmx,nlayermx)
    337345      REAL pdt_th(ngridmx,nlayermx),pdq_th(ngridmx,nlayermx,nqmx)
    … …  
    342350      REAL z_out                          ! height of interpolation between z0 and z1 [meters]
    343351      REAL ustar(ngridmx),tstar(ngridmx)  ! friction velocity and friction potential temp
    344       REAL L_mo(ngridmx),wstarpbl(ngridmx),vhf(ngridmx),vvv(ngridmx)
     352      REAL L_mo(ngridmx)
    345353      REAL zu2(ngridmx)
    346354c=======================================================================
    … …  
    359367         fluxrad(:)=0
    360368
    361          wstar(:)=0.
     369         wmax_th(:)=0.
    362370
    363371c        read startfi
    … …  
    569577                 enddo
    570578              endif
    571            else
    572               zdtnlte(:,:)=0.
    573579           endif
    574580
    … …  
    589595     $     emis,mu0,pplev,pplay,pt,tsurf,fract,dist_sol,igout,
    590596     $     zdtlw,zdtsw,fluxsurf_lw,fluxsurf_sw,fluxtop_lw,fluxtop_sw,
    591      &     tauref,tau,aerosol,tauscaling,rdust,rice,nuice,co2ice)
     597     $     tauref,tau,aerosol,tauscaling,taucloudtes,rdust,rice,
     598     $     nuice,co2ice)
    592599
    593600c          Outputs for basic check (middle of domain)
    … …  
    747754          DO ig=1, ngridmx
    748755             IF (zh(ig,1) .lt. tsurf(ig)) THEN
    749                wstar(ig)=1.
    750                hfmax_th(ig)=0.2
    751              ELSE
    752                wstar(ig)=0.
    753                hfmax_th(ig)=0.
    754              ENDIF     
     756               wmax_th(ig)=1.
     757             ENDIF       
    755758          ENDDO
    756759        ENDIF
    … …  
    768771     $        zdum1,zdum2,zdh,pdq,zflubid,
    769772     $        zdudif,zdvdif,zdhdif,zdtsdif,q2,
    770      &        zdqdif,zdqsdif,wstar,zcdv,zcdh,hfmax_th)
     773     &        zdqdif,zdqsdif,wmax_th,zcdv,zcdh)
    771774
    772775#ifdef MESOSCALE
    … …  
    829832     $ pplay,pplev,pphi,zpopsk,
    830833     $ pdu_th,pdv_th,pdt_th,pdq_th,lmax_th,zmax_th,
    831      $ dtke_th,hfmax_th,wstar)
     834     $ dtke_th,hfmax_th,wmax_th)
    832835 
    833836         DO l=1,nlayer
    … …  
    858861        else   !of if calltherm
    859862        lmax_th(:)=0
    860         wstar(:)=0.
    861         hfmax_th(:)=0.
     863        wmax_th(:)=0.
    862864        lmax_th_out(:)=0.
    863865        end if
    … …  
    911913c   -------------------------------------------
    912914
    913       IF (tituscap) THEN
    914          !!! get the actual co2 seasonal cap from Titus observations
    915          CALL geticecover( ngrid, 180.*zls/pi,
     915#ifdef MESOSCALE
     916      !!! get the actual co2 seasonal cap from Titus observations
     917      CALL geticecover( ngrid, 180.*zls/pi,
    916918     .                  180.*long/pi, 180.*lati/pi, co2ice )
    917          co2ice = co2ice * 10000.
    918       ENDIF
     919      co2ice = co2ice * 10000.
     920#endif
    919921
    920922      IF (callcond) THEN
    … …  
    962964c        ----------------------------------------
    963965         IF (water) THEN
     966         
    964967
    965968           call watercloud(ngrid,nlayer,ptimestep,
    … …  
    976979           ENDDO
    977980           endif
    978 
     981           
    979982! increment water vapour and ice atmospheric tracers tendencies
    980983           IF (water) THEN
    … …  
    987990     &                    pdq(ig,l,igcm_h2o_ice)+
    988991     &                    zdqcloud(ig,l,igcm_h2o_ice)
     992                 if (((pq(ig,l,igcm_h2o_ice) +
     993     &                 pdq(ig,l,igcm_h2o_ice)*ptimestep)) .le. 0)
     994     &           then
     995                   pdq(ig,l,igcm_h2o_ice) =
     996     &                 - pq(ig,l,igcm_h2o_ice)/ptimestep + 1e-30
     997                 endif
    989998                 IF (scavenging) THEN
    990999                   pdq(ig,l,igcm_ccn_mass)=
    … …  
    9941003     &                      pdq(ig,l,igcm_ccn_number)+
    9951004     &                      zdqcloud(ig,l,igcm_ccn_number)
     1005                   pdq(ig,l,igcm_dust_mass)=
     1006     &                      pdq(ig,l,igcm_dust_mass)+
     1007     &                      zdqcloud(ig,l,igcm_dust_mass)
     1008                   pdq(ig,l,igcm_dust_number)=
     1009     &                      pdq(ig,l,igcm_dust_number)+
     1010     &                      zdqcloud(ig,l,igcm_dust_number)
    9961011!!!!!!!!!!!!!!!!!!!!! We need to check that we have Nccn & Ndust > 0
    9971012!!!!!!!!!!!!!!!!!!!!! This is due to single precision rounding problems 
    … …  
    10011016                       pdq(ig,l,igcm_ccn_number) =
    10021017     &                 - pq(ig,l,igcm_ccn_number)/ptimestep + 1e-30
     1018                       pdq(ig,l,igcm_ccn_mass) =
     1019     &                 - pq(ig,l,igcm_ccn_mass)/ptimestep + 1e-30
    10031020                   endif
    10041021                   if (((pq(ig,l,igcm_dust_number) +
    … …  
    10071024                       pdq(ig,l,igcm_dust_number) =
    10081025     &                 - pq(ig,l,igcm_dust_number)/ptimestep + 1e-30
     1026                       pdq(ig,l,igcm_dust_mass) =
     1027     &                 - pq(ig,l,igcm_dust_mass)/ptimestep + 1e-30
    10091028                   endif
    10101029!!!!!!!!!!!!!!!!!!!!!
    10111030!!!!!!!!!!!!!!!!!!!!!
    1012                    pdq(ig,l,igcm_dust_mass)=
    1013      &                      pdq(ig,l,igcm_dust_mass)+
    1014      &                      zdqcloud(ig,l,igcm_dust_mass)
    1015                    pdq(ig,l,igcm_dust_number)=
    1016      &                      pdq(ig,l,igcm_dust_number)+
    1017      &                      zdqcloud(ig,l,igcm_dust_number)
    10181031                 ENDIF
    10191032               ENDDO
    … …  
    10291042         END IF  ! of IF (water)
    10301043
    1031 c   7b. Aerosol particles
     1044
     1045c   7b. Chemical species
     1046c     ------------------
     1047
     1048#ifndef MESOSCALE
     1049c        --------------
     1050c        photochemistry :
     1051c        --------------
     1052         IF (photochem .or. thermochem) then
     1053!NB: Photochemistry includes condensation of H2O2
     1054            PRINT*, 'SURFDUST,SURFICE TO BE IMPLEMENTED. YAAAAAARG.'
     1055            call calchim(ptimestep,pplay,pplev,pt,pdt,dist_sol,mu0,
     1056     $                   zzlev,zzlay,zday,pq,pdq,zdqchim,zdqschim,
     1057     $                   zdqcloud,zdqscloud,tauref,co2ice,
     1058     $                   pu,pdu,pv,pdv,surfdust,surfice)
     1059
     1060           ! increment values of tracers:
     1061           DO iq=1,nq ! loop on all tracers; tendencies for non-chemistry
     1062                      ! tracers is zero anyways
     1063             DO l=1,nlayer
     1064               DO ig=1,ngrid
     1065                 pdq(ig,l,iq)=pdq(ig,l,iq)+zdqchim(ig,l,iq)
     1066               ENDDO
     1067             ENDDO
     1068           ENDDO ! of DO iq=1,nq
     1069           ! add condensation tendency for H2O2
     1070           if (igcm_h2o2.ne.0) then
     1071             DO l=1,nlayer
     1072               DO ig=1,ngrid
     1073                 pdq(ig,l,igcm_h2o2)=pdq(ig,l,igcm_h2o2)
     1074     &                                +zdqcloud(ig,l,igcm_h2o2)
     1075               ENDDO
     1076             ENDDO
     1077           endif
     1078
     1079           ! increment surface values of tracers:
     1080           DO iq=1,nq ! loop on all tracers; tendencies for non-chemistry
     1081                      ! tracers is zero anyways
     1082             DO ig=1,ngrid
     1083               dqsurf(ig,iq)=dqsurf(ig,iq)+zdqschim(ig,iq)
     1084             ENDDO
     1085           ENDDO ! of DO iq=1,nq
     1086           ! add condensation tendency for H2O2
     1087           if (igcm_h2o2.ne.0) then
     1088             DO ig=1,ngrid
     1089               dqsurf(ig,igcm_h2o2)=dqsurf(ig,igcm_h2o2)
     1090     &                                +zdqscloud(ig,igcm_h2o2)
     1091             ENDDO
     1092           endif
     1093
     1094         END IF  ! of IF (photochem.or.thermochem)
     1095#endif
     1096
     1097c   7c. Aerosol particles
    10321098c     -------------------
    10331099
    … …  
    10701136     &                pq, pdq, zdqsed, zdqssed,nq,
    10711137     &                tau,tauscaling)
     1138     
     1139     
     1140      !print*, 'h2o_ice zdqsed ds physiq', zdqsed(1,:,igcm_h2o_ice)
    10721141
    10731142           DO iq=1, nq
    … …  
    10841153           ENDDO
    10851154         END IF   ! of IF (sedimentation)
    1086 
    1087 c   7c. Chemical species
    1088 c     ------------------
    1089 
    1090 #ifndef MESOSCALE
    1091 c        --------------
    1092 c        photochemistry :
    1093 c        --------------
    1094          IF (photochem .or. thermochem) then
    1095 
    1096 !           dust and ice surface area
    1097             call surfacearea(ngrid, nlayer, ptimestep, pplay, zzlay,
    1098      $                       pt, pq, pdq, nq,
    1099      $                       rdust, rice, tau, tauscaling,
    1100      $                       surfdust, surfice)
    1101 !           call photochemistry
    1102             call calchim(ptimestep,pplay,pplev,pt,pdt,dist_sol,mu0,
    1103      $                   zzlev,zzlay,zday,pq,pdq,zdqchim,zdqschim,
    1104      $                   zdqcloud,zdqscloud,tauref,co2ice,
    1105      $                   pu,pdu,pv,pdv,surfdust,surfice)
    1106 
    1107            ! increment values of tracers:
    1108            DO iq=1,nq ! loop on all tracers; tendencies for non-chemistry
    1109                       ! tracers is zero anyways
    1110              DO l=1,nlayer
    1111                DO ig=1,ngrid
    1112                  pdq(ig,l,iq)=pdq(ig,l,iq)+zdqchim(ig,l,iq)
    1113                ENDDO
    1114              ENDDO
    1115            ENDDO ! of DO iq=1,nq
    1116            
    1117            ! add condensation tendency for H2O2
    1118            if (igcm_h2o2.ne.0) then
    1119              DO l=1,nlayer
    1120                DO ig=1,ngrid
    1121                  pdq(ig,l,igcm_h2o2)=pdq(ig,l,igcm_h2o2)
    1122      &                                +zdqcloud(ig,l,igcm_h2o2)
    1123                ENDDO
    1124              ENDDO
    1125            endif
    1126 
    1127            ! increment surface values of tracers:
    1128            DO iq=1,nq ! loop on all tracers; tendencies for non-chemistry
    1129                       ! tracers is zero anyways
    1130              DO ig=1,ngrid
    1131                dqsurf(ig,iq)=dqsurf(ig,iq)+zdqschim(ig,iq)
    1132              ENDDO
    1133            ENDDO ! of DO iq=1,nq
    1134 
    1135            ! add condensation tendency for H2O2
    1136            if (igcm_h2o2.ne.0) then
    1137              DO ig=1,ngrid
    1138                dqsurf(ig,igcm_h2o2)=dqsurf(ig,igcm_h2o2)
    1139      &                                +zdqscloud(ig,igcm_h2o2)
    1140              ENDDO
    1141            endif
    1142            
    1143          END IF  ! of IF (photochem.or.thermochem)
    1144 #endif
     1155         
    11451156
    11461157
    … …  
    11581169          ENDDO  ! (ig)
    11591170        ENDDO    ! (iq)
     1171
    11601172
    11611173      endif !  of if (tracer)
    … …  
    12321244      ENDIF  ! of IF (tracer.AND.water.AND.(ngridmx.NE.1))
    12331245
    1234 c
     1246
    12351247c   9.2 Compute soil temperatures and subsurface heat flux:
    12361248c   ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    … …  
    12911303       DO ig=1,ngridmx
    12921304          DO l=1,nlayermx
    1293               zh(ig,l) = zt(ig,l)*(zplev(ig,1)/zplay(ig,l))**rcp
     1305              zh(ig,l) = zt(ig,l)*(zplay(ig,l)/zplev(ig,1))**rcp
    12941306          ENDDO
    12951307       ENDDO
    … …  
    13741386         if (tracer) then
    13751387           if (water) then
     1388           
     1389             if (scavenging) then
     1390               ccntot(:)= 0
     1391               do ig=1,ngrid
     1392                 do l=1,nlayermx
     1393                   ccntot(ig) = ccntot(ig) +
     1394     &              zq(ig,l,igcm_ccn_number)*tauscaling(ig)
     1395     &              *(pplev(ig,l) - pplev(ig,l+1)) / g
     1396                 enddo
     1397               enddo
     1398             endif
    13761399
    13771400             mtot(:)=0
    … …  
    14201443           call wstats(ngrid,"ps","Surface pressure","Pa",2,ps)
    14211444           call wstats(ngrid,"tsurf","Surface temperature","K",2,tsurf)
    1422            call wstats(ngrid,"co2ice","CO2 ice cover",
    1423      &                "kg.m-2",2,co2ice)
    1424            call wstats(ngrid,"fluxsurf_lw",
    1425      &                "Thermal IR radiative flux to surface","W.m-2",2,
    1426      &                fluxsurf_lw)
    1427            call wstats(ngrid,"fluxsurf_sw",
    1428      &                "Solar radiative flux to surface","W.m-2",2,
    1429      &                fluxsurf_sw_tot)
    1430            call wstats(ngrid,"fluxtop_lw",
    1431      &                "Thermal IR radiative flux to space","W.m-2",2,
    1432      &                fluxtop_lw)
    1433            call wstats(ngrid,"fluxtop_sw",
    1434      &                "Solar radiative flux to space","W.m-2",2,
    1435      &                fluxtop_sw_tot)
    1436            call wstats(ngrid,"temp","Atmospheric temperature","K",3,zt)
    1437            call wstats(ngrid,"u","Zonal (East-West) wind","m.s-1",3,zu)
    1438            call wstats(ngrid,"v","Meridional (North-South) wind",
    1439      &                "m.s-1",3,zv)
    1440            call wstats(ngrid,"w","Vertical (down-up) wind",
    1441      &                "m.s-1",3,pw)
    1442            call wstats(ngrid,"rho","Atmospheric density","none",3,rho)
    1443            call wstats(ngrid,"pressure","Pressure","Pa",3,pplay)
     1445c           call wstats(ngrid,"co2ice","CO2 ice cover",
     1446c     &                "kg.m-2",2,co2ice)
     1447c           call wstats(ngrid,"fluxsurf_lw",
     1448c     &                "Thermal IR radiative flux to surface","W.m-2",2,
     1449c     &                fluxsurf_lw)
     1450c           call wstats(ngrid,"fluxsurf_sw",
     1451c     &                "Solar radiative flux to surface","W.m-2",2,
     1452c     &                fluxsurf_sw_tot)
     1453c           call wstats(ngrid,"fluxtop_lw",
     1454c     &                "Thermal IR radiative flux to space","W.m-2",2,
     1455c     &                fluxtop_lw)
     1456c           call wstats(ngrid,"fluxtop_sw",
     1457c     &                "Solar radiative flux to space","W.m-2",2,
     1458c     &                fluxtop_sw_tot)
     1459c           call wstats(ngrid,"temp","Atmospheric temperature","K",3,zt)
     1460c           call wstats(ngrid,"u","Zonal (East-West) wind","m.s-1",3,zu)
     1461c           call wstats(ngrid,"v","Meridional (North-South) wind",
     1462c     &                "m.s-1",3,zv)
     1463c           call wstats(ngrid,"w","Vertical (down-up) wind",
     1464c     &                "m.s-1",3,pw)
     1465c           call wstats(ngrid,"rho","Atmospheric density","none",3,rho)
     1466c           call wstats(ngrid,"pressure","Pressure","Pa",3,pplay)
    14441467c          call wstats(ngrid,"q2",
    14451468c    &                "Boundary layer eddy kinetic energy",
    … …  
    14561479           if (tracer) then
    14571480             if (water) then
    1458 c               vmr=zq(1:ngridmx,1:nlayermx,igcm_h2o_vap)
    1459 c     &                  *mugaz/mmol(igcm_h2o_vap)
    1460 c               call wstats(ngrid,"vmr_h2ovapor",
    1461 c     &                    "H2O vapor volume mixing ratio","mol/mol",
    1462 c     &                    3,vmr)
    1463 c               vmr=zq(1:ngridmx,1:nlayermx,igcm_h2o_ice)
    1464 c     &                  *mugaz/mmol(igcm_h2o_ice)
    1465 c               call wstats(ngrid,"vmr_h2oice",
    1466 c     &                    "H2O ice volume mixing ratio","mol/mol",
    1467 c     &                    3,vmr)
     1481               vmr=zq(1:ngridmx,1:nlayermx,igcm_h2o_vap)
     1482     &                  *mugaz/mmol(igcm_h2o_vap)
     1483               call wstats(ngrid,"vmr_h2ovapor",
     1484     &                    "H2O vapor volume mixing ratio","mol/mol",
     1485     &                    3,vmr)
     1486               vmr=zq(1:ngridmx,1:nlayermx,igcm_h2o_ice)
     1487     &                  *mugaz/mmol(igcm_h2o_ice)
     1488               call wstats(ngrid,"vmr_h2oice",
     1489     &                    "H2O ice volume mixing ratio","mol/mol",
     1490     &                    3,vmr)
    14681491               call wstats(ngrid,"h2o_ice_s",
    14691492     &                    "surface h2o_ice","kg/m2",
    14701493     &                    2,qsurf(1,igcm_h2o_ice))
    1471 c               call wstats(ngrid,'albedo',
    1472 c     &                         'albedo',
    1473 c     &                         '',2,albedo(1:ngridmx,1))
     1494               call wstats(ngrid,'albedo',
     1495     &                         'albedo',
     1496     &                         '',2,albedo(1:ngridmx,1))
    14741497               call wstats(ngrid,"mtot",
    14751498     &                    "total mass of water vapor","kg/m2",
    … …  
    14781501     &                    "total mass of water ice","kg/m2",
    14791502     &                    2,icetot)
    1480 c               call wstats(ngrid,"reffice",
    1481 c     &                    "Mean reff","m",
    1482 c     &                    2,rave)
    1483 c              call wstats(ngrid,"rice",
    1484 c    &                    "Ice particle size","m",
    1485 c    &                    3,rice)
    1486 c              If activice is true, tauTES is computed in aeropacity.F;
     1503               call wstats(ngrid,"reffice",
     1504     &                    "Mean reff","m",
     1505     &                    2,rave)
     1506              call wstats(ngrid,"ccntot",
     1507     &                    "condensation nuclei","Nbr/m2",
     1508     &                    2,ccntot)
     1509              call wstats(ngrid,"rice",
     1510     &                    "Ice particle size","m",
     1511     &                    3,rice)
    14871512               if (.not.activice) then
    14881513                 call wstats(ngrid,"tauTESap",
    14891514     &                      "tau abs 825 cm-1","",
    14901515     &                      2,tauTES)
     1516               else
     1517                 call wstats(ngridmx,'tauTES',
     1518     &                   'tau abs 825 cm-1',
     1519     &                  '',2,taucloudtes)
    14911520               endif
    14921521
    … …  
    14951524             if (thermochem.or.photochem) then
    14961525                do iq=1,nq
    1497                    if (noms(iq) .ne. "dust_mass" .and.
    1498      $                 noms(iq) .ne. "dust_number" .and.
    1499      $                 noms(iq) .ne. "ccn_mass" .and.
    1500      $                 noms(iq) .ne. "ccn_number") then
    1501                    do l=1,nlayer
    1502                       do ig=1,ngrid
    1503                          vmr(ig,l)=zq(ig,l,iq)*mmean(ig,l)/mmol(iq)
    1504                       end do
    1505                    end do
    1506                    call wstats(ngrid,"vmr_"//trim(noms(iq)),
    1507      $                         "Volume mixing ratio","mol/mol",3,vmr)
    1508                    if ((noms(iq).eq."o") .or. (noms(iq).eq."co2").or.
    1509      $                 (noms(iq).eq."o3")) then                     
    1510                       call writediagfi(ngrid,"vmr_"//trim(noms(iq)),
    1511      $                         "Volume mixing ratio","mol/mol",3,vmr)
    1512                    end if
    1513                    do ig = 1,ngrid
    1514                       colden(ig,iq) = 0.                           
    1515                    end do
    1516                    do l=1,nlayer                                   
    1517                       do ig=1,ngrid                                 
    1518                          colden(ig,iq) = colden(ig,iq) + zq(ig,l,iq)
    1519      $                                  *(pplev(ig,l)-pplev(ig,l+1))
    1520      $                                  *6.022e22/(mmol(iq)*g)       
    1521                       end do                                       
    1522                    end do                                         
    1523                    call wstats(ngrid,"c_"//trim(noms(iq)),           
    1524      $                         "column","mol cm-2",2,colden(1,iq)) 
    1525                    call writediagfi(ngrid,"c_"//trim(noms(iq)), 
    1526      $                             "column","mol cm-2",2,colden(1,iq))
    1527                    end if
     1526                   if ((noms(iq).eq."o").or.(noms(iq).eq."co2").or.
     1527     .                (noms(iq).eq."co").or.(noms(iq).eq."n2").or.
     1528     .                (noms(iq).eq."h2").or.
     1529     .                (noms(iq).eq."o3")) then
     1530                        do l=1,nlayer
     1531                          do ig=1,ngrid
     1532                            vmr(ig,l)=zq(ig,l,iq)*mmean(ig,l)/mmol(iq)
     1533                          end do
     1534                        end do
     1535                        call wstats(ngrid,"vmr_"//trim(noms(iq)),
     1536     .                     "Volume mixing ratio","mol/mol",3,vmr)
     1537                   endif
     1538!                   do ig = 1,ngrid
     1539!                      colden(ig,iq) = 0.                             !FL
     1540!                   end do
     1541!                   do l=1,nlayer                                     !FL
     1542!                      do ig=1,ngrid                                  !FL
     1543!                         colden(ig,iq) = colden(ig,iq) + zq(ig,l,iq) !FL
     1544!     $                                  *(pplev(ig,l)-pplev(ig,l+1)) !FL
     1545!     $                                  *6.022e22/(mmol(iq)*g)       !FL
     1546!                      end do                                         !FL
     1547!                   end do                                            !FL
     1548!                   call wstats(ngrid,"c_"//trim(noms(iq)),           !FL
     1549!     $                         "column","mol cm-2",2,colden(1,iq))   !FL
    15281550                end do
    15291551             end if ! of if (thermochem.or.photochem)
    … …  
    15831605c        call WRITEDIAGFI(ngrid,"emis","Surface emissivity","w.m-1",2,
    15841606c    &                  emis)
    1585          call WRITEDIAGFI(ngrid,"pplay","Pressure","Pa",3,zplay)
    1586 !         call WRITEDIAGFI(ngrid,"pplev","Pressure","Pa",3,zplev)
     1607c         call WRITEDIAGFI(ngrid,"pplay","Pressure","Pa",3,zplay)
     1608c         call WRITEDIAGFI(ngrid,"pplev","Pressure","Pa",3,zplev)
    15871609         call WRITEDIAGFI(ngrid,"tsurf","Surface temperature","K",2,
    15881610     &                  tsurf)
    … …  
    16011623c         call WRITEDIAGFI(ngrid,"fluxtop_sw","fluxtop_sw","W.m-2",2,
    16021624c     &                  fluxtop_sw_tot)
    1603          call WRITEDIAGFI(ngrid,"temp","temperature","K",3,zt)
    1604         call WRITEDIAGFI(ngrid,"u","Zonal wind","m.s-1",3,zu)
    1605         call WRITEDIAGFI(ngrid,"v","Meridional wind","m.s-1",3,zv)
    1606         call WRITEDIAGFI(ngrid,"w","Vertical wind","m.s-1",3,pw)
    1607          call WRITEDIAGFI(ngrid,"rho","density","none",3,rho)
     1625c         call WRITEDIAGFI(ngrid,"temp","temperature","K",3,zt)
     1626c        call WRITEDIAGFI(ngrid,"u","Zonal wind","m.s-1",3,zu)
     1627c        call WRITEDIAGFI(ngrid,"v","Meridional wind","m.s-1",3,zv)
     1628c        call WRITEDIAGFI(ngrid,"w","Vertical wind","m.s-1",3,pw)
     1629c         call WRITEDIAGFI(ngrid,"rho","density","none",3,rho)
    16081630c        call WRITEDIAGFI(ngrid,"q2","q2","kg.m-3",3,q2)
    1609 !        call WRITEDIAGFI(ngrid,'Teta','T potentielle','K',3,zh)
     1631c        call WRITEDIAGFI(ngrid,'Teta','T potentielle','K',3,zh)
    16101632c        call WRITEDIAGFI(ngrid,"pressure","Pressure","Pa",3,pplay)
    16111633c        call WRITEDIAGFI(ngrid,"ssurf","Surface stress","N.m-2",2,
    … …  
    16151637c        call WRITEDIAGFI(ngridmx,'lw_htrt','lw heat. rate',
    16161638c    &                   'w.m-2',3,zdtlw)
    1617 c        CALL WRITEDIAGFI(ngridmx,'tauTESap',
    1618 c     &                         'tau abs 825 cm-1',
    1619 c     &                         '',2,tauTES)
     1639            if (.not.activice) then
     1640               CALL WRITEDIAGFI(ngridmx,'tauTESap',
     1641     &                         'tau abs 825 cm-1',
     1642     &                         '',2,tauTES)
     1643             else
     1644               CALL WRITEDIAGFI(ngridmx,'tauTES',
     1645     &                         'tau abs 825 cm-1',
     1646     &                         '',2,taucloudtes)
     1647             endif
    16201648
    16211649#ifdef MESOINI
    … …  
    16501678         call WRITEDIAGFI(ngrid,"co2col","CO2 column","kg.m-2",2,
    16511679     &                  co2col)
     1680!!!!! FL
     1681!            do iq = 1,nq
     1682!               if (noms(iq) .ne. "dust_mass" .and.
     1683!     $             noms(iq) .ne. "dust_number") then
     1684!               call writediagfi(ngrid,"c_"//trim(noms(iq)),         
     1685!     $                         "column","mol cm-2",2,colden(1,iq))
     1686!               end if
     1687!            end do
     1688!!!!! FL
    16521689         endif ! of if (tracer.and.(igcm_co2.ne.0))
    16531690
    … …  
    16711708     &                      'kg.m-2',2,
    16721709     &                       qsurf(1:ngridmx,igcm_h2o_ice))
    1673             if (photochem) then
    1674               do iq = 1,nq
    1675                call writediagfi( ngrid,trim(noms(iq)),
    1676      $                           'mix rat','units',
    1677      $                           3,zq(1:ngridmx,1:nlayermx,iq) )
    1678               enddo
    1679             endif
    16801710#endif
    16811711
    … …  
    16971727c     &                       'Mean reff',
    16981728c     &                       'm',2,rave)
     1729c             CALL WRITEDIAGFI(ngrid,"ccntot",
     1730c     &                    "condensation nuclei","Nbr/m2",
     1731c     &                    2,ccntot)
    16991732c            call WRITEDIAGFI(ngridmx,'rice','Ice particle size',
    17001733c     &                       'm',3,rice)
    1701 c            If activice is true, tauTES is computed in aeropacity.F;
    1702              if (.not.activice) then
    1703                CALL WRITEDIAGFI(ngridmx,'tauTESap',
    1704      &                         'tau abs 825 cm-1',
    1705      &                         '',2,tauTES)
    1706              endif
    17071734             call WRITEDIAGFI(ngridmx,'h2o_ice_s',
    17081735     &                       'surface h2o_ice',
    17091736     &                       'kg.m-2',2,qsurf(1,igcm_h2o_ice))
    17101737
    1711             if (caps) then
    1712              do ig=1,ngridmx
    1713                 if (watercaptag(ig)) watercapflag(ig) = 1
    1714              enddo
     1738c            if (caps) then
     1739c             do ig=1,ngridmx
     1740c                if (watercaptag(ig)) watercapflag(ig) = 1
     1741c             enddo
    17151742c            CALL WRITEDIAGFI(ngridmx,'watercaptag',
    1716 c    &                         'Ice water caps',
    1717 c    &                         '',2,watercapflag)
    1718             endif
     1743c     &                         'Ice water caps',
     1744c     &                         '',2,watercapflag)
     1745cs            endif
    17191746c           CALL WRITEDIAGFI(ngridmx,'albedo',
    17201747c    &                         'albedo',
    … …  
    17611788c             call WRITEDIAGFI(ngridmx,'reffdust','reffdust',
    17621789c     &                        'm',3,rdust*ref_r0)
    1763 c             call WRITEDIAGFI(ngridmx,'rdust','rdust',
    1764 c    &                        'm',3,rdust)
    17651790c             call WRITEDIAGFI(ngridmx,'dustq','Dust mass mr',
    17661791c     &                        'kg/kg',3,pq(1,1,igcm_dust_mass))
    … …  
    17681793c     &                        'part/kg',3,pq(1,1,igcm_dust_number))
    17691794#ifdef MESOINI
    1770              call WRITEDIAGFI(ngridmx,'dustq','Dust mass mr',
    1771      &           'kg/kg',3,pq(1:ngridmx,1:nlayermx,igcm_dust_mass))
    17721795             call WRITEDIAGFI(ngridmx,'dustN','Dust number',
    1773      &           'part/kg',3,pq(1:ngridmx,1:nlayermx,igcm_dust_number))
     1796     &                        'part/kg',3,pq(1,1,igcm_dust_number))
    17741797#endif
    17751798           else
    … …  
    17841807
    17851808           if (scavenging) then
    1786              call WRITEDIAGFI(ngridmx,'ccnq','CCN mass mr',
    1787      &                        'kg/kg',3,pq(1,1,igcm_ccn_mass))
    1788              call WRITEDIAGFI(ngridmx,'ccnN','CCN number',
    1789      &                        'part/kg',3,pq(1,1,igcm_ccn_number))
     1809c             call WRITEDIAGFI(ngridmx,'ccnq','CCN mass mr',
     1810c     &                        'kg/kg',3,pq(1,1,igcm_ccn_mass))
     1811c             call WRITEDIAGFI(ngridmx,'ccnN','CCN number',
     1812c     &                        'part/kg',3,pq(1,1,igcm_ccn_number))
    17901813           endif ! (scavenging)
    17911814
    … …  
    18101833         z_out=0.
    18111834         if (calltherm .and. (z_out .gt. 0.)) then
    1812 
    1813          call pbl_parameters(ngrid,nlayer,z0,
    1814      & g,zzlay,zu,zv,wstar,hfmax_th,zmax_th,tsurf,zh,z_out,
    1815      & Teta_out,u_out,ustar,tstar,wstarpbl,L_mo,vhf,vvv)
    1816 
     1835         call surflayer_interpol(ngrid,nlayer,z0,g,zzlay,zu,zv,wmax_th
     1836     &              ,tsurf,zh,z_out,Teta_out,u_out,ustar,tstar,L_mo)
     1837
     1838         zu2(:)=sqrt(zu(:,1)*zu(:,1)+zv(:,1)*zv(:,1))
     1839         call WRITEDIAGFI(ngridmx,'sqrt(zu2)',
     1840     &              'horizontal velocity norm','m/s',
     1841     &                         2,zu2)
     1842
     1843         call WRITEDIAGFI(ngridmx,'Teta_out',
     1844     &              'potential temperature at z_out','K',
     1845     &                         2,Teta_out)
     1846         call WRITEDIAGFI(ngridmx,'u_out',
     1847     &              'horizontal velocity norm at z_out','m/s',
     1848     &                         2,u_out)
     1849         call WRITEDIAGFI(ngridmx,'u*',
     1850     &              'friction velocity','m/s',
     1851     &                         2,ustar)
     1852         call WRITEDIAGFI(ngridmx,'teta*',
     1853     &              'friction potential temperature','K',
     1854     &                         2,tstar)
     1855         call WRITEDIAGFI(ngrid,'L',
     1856     &              'Monin Obukhov length','m',
     1857     &                         2,L_mo)
    18171858         else
    18181859           if((.not. calltherm).and.(z_out .gt. 0.)) then
    … …  
    18511892     &              'maximum TH heat flux','K.m/s',
    18521893     &                         2,hfmax_th)
    1853         call WRITEDIAGFI(ngridmx,'wstar',
     1894        call WRITEDIAGFI(ngridmx,'wmax_th',
    18541895     &              'maximum TH vertical velocity','m/s',
    1855      &                         2,wstar)
     1896     &                         2,wmax_th)
    18561897
    18571898         endif
    … …  
    19031944         z_out=0.
    19041945         if (calltherm .and. (z_out .gt. 0.)) then
    1905 
    1906          call pbl_parameters(ngrid,nlayer,z0,
    1907      & g,zzlay,zu,zv,wstar,hfmax_th,zmax_th,tsurf,zh,z_out,
    1908      & Teta_out,u_out,ustar,tstar,wstarpbl,L_mo,vhf,vvv)
    1909 
     1946         call surflayer_interpol(ngrid,nlayer,z0,g,zzlay,zu,zv,wmax_th
     1947     &              ,tsurf,zh,z_out,Teta_out,u_out,ustar,tstar,L_mo)
     1948
     1949         zu2(:)=sqrt(zu(:,1)*zu(:,1)+zv(:,1)*zv(:,1))
     1950         call WRITEDIAGFI(ngridmx,'sqrt(zu2)',
     1951     &              'horizontal velocity norm','m/s',
     1952     &                         0,zu2)
     1953
     1954         call WRITEDIAGFI(ngridmx,'Teta_out',
     1955     &              'potential temperature at z_out','K',
     1956     &                         0,Teta_out)
     1957         call WRITEDIAGFI(ngridmx,'u_out',
     1958     &              'horizontal velocity norm at z_out','m/s',
     1959     &                         0,u_out)
     1960         call WRITEDIAGFI(ngridmx,'u*',
     1961     &              'friction velocity','m/s',
     1962     &                         0,ustar)
     1963         call WRITEDIAGFI(ngridmx,'teta*',
     1964     &              'friction potential temperature','K',
     1965     &                         0,tstar)
    19101966         else
    19111967           if((.not. calltherm).and.(z_out .gt. 0.)) then
    … …  
    19211977     &              'hauteur du thermique','point',
    19221978     &                         0,lmax_th_out)
    1923         call WRITEDIAGFI(ngridmx,'zmax_th',
    1924      &              'hauteur du thermique','m',
    1925      &                         0,zmax_th)
    19261979        call WRITEDIAGFI(ngridmx,'hfmax_th',
    19271980     &              'maximum TH heat flux','K.m/s',
    19281981     &                         0,hfmax_th)
    1929         call WRITEDIAGFI(ngridmx,'wstar',
     1982        call WRITEDIAGFI(ngridmx,'wmax_th',
    19301983     &              'maximum TH vertical velocity','m/s',
    1931      &                         0,wstar)
     1984     &                         0,wmax_th)
    19321985
    19331986         co2col(:)=0.
    … …  
    19732026         end if
    19742027         
    1975 cccccccccccccccccccccc scavenging outputs 1D TN ccccccccccccccccccc
     2028cccccccccccccccccc scavenging & water outputs 1D TN ccccccccccccccc
    19762029ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    1977          IF (scavenging) THEN
    1978              CALL WRITEDIAGFI(ngridmx,'tauTESap',
     2030         IF (water) THEN
     2031           CALL WRITEDIAGFI(ngridmx,'tauTESap',
    19792032     &                         'tau abs 825 cm-1',
    1980      &                         '',1,tauTES)
     2033     &                         '',0,tauTES)
    19812034     
    1982            h2o_tot = qsurf(1,igcm_h2o_ice)
     2035           CALL WRITEDIAGFI(ngridmx,'tauTES',
     2036     &                         'tau abs 825 cm-1',
     2037     &                         '',0,taucloudtes)
     2038     
     2039           mtot = 0
     2040           icetot = 0
     2041           h2otot = qsurf(1,igcm_h2o_ice)
     2042           rave = 0
    19832043           do l=1,nlayer
    1984              h2o_tot = h2o_tot +
    1985      &           (zq(1,l,igcm_h2o_vap) + zq(1,l,igcm_h2o_ice))
    1986      &                     * (pplev(1,l) - pplev(1,l+1)) / g
     2044             mtot = mtot +  zq(1,l,igcm_h2o_vap)
     2045     &                 * (pplev(1,l) - pplev(1,l+1)) / g
     2046             icetot = icetot +  zq(1,l,igcm_h2o_ice)
     2047     &                 * (pplev(1,l) - pplev(1,l+1)) / g
     2048             rave = rave + zq(1,l,igcm_h2o_ice)
     2049     &                 * (pplev(1,l) - pplev(1,l+1)) / g
     2050     &                 *  rice(1,l) * (1.+nuice_ref)
    19872051           end do
     2052             rave=rave/max(icetot,1.e-30)
     2053            h2otot = h2otot+mtot+icetot
     2054           
     2055           
     2056             if (scavenging) then
     2057               ccntot= 0
     2058              call watersat(ngridmx*nlayermx,zt,pplay,zqsat)
     2059               do l=1,nlayermx
     2060                   ccntot = ccntot +
     2061     &              zq(1,l,igcm_ccn_number)*tauscaling(1)
     2062     &              *(pplev(1,l) - pplev(1,l+1)) / g
     2063                   satu(1,l) = zq(1,l,igcm_h2o_vap)/zqsat(1,l)
     2064                   satu(1,l) = (max(satu(1,l),float(1))-1)
     2065!     &                      * zq(1,l,igcm_h2o_vap) *
     2066!     &                      (pplev(1,l) - pplev(1,l+1)) / g
     2067               enddo
     2068
     2069               CALL WRITEDIAGFI(ngridmx,'ccntot',
     2070     &                         'ccntot',
     2071     &                         'nbr/m2',0,ccntot)
     2072             endif
     2073             
     2074             
     2075             CALL WRITEDIAGFI(ngridmx,'h2otot',
     2076     &                         'h2otot',
     2077     &                         'kg/m2',0,h2otot)
     2078             CALL WRITEDIAGFI(ngridmx,'mtot',
     2079     &                         'mtot',
     2080     &                         'kg/m2',0,mtot)
     2081             CALL WRITEDIAGFI(ngridmx,'icetot',
     2082     &                         'icetot',
     2083     &                         'kg/m2',0,icetot)
     2084             CALL WRITEDIAGFI(ngridmx,'reffice',
     2085     &                         'reffice',
     2086     &                         'm',0,rave)
    19882087 
    19892088
    … …  
    19942093
    19952094         
    1996             call WRITEDIAGFI(ngridmx,'dqssed q','sedimentation q',
    1997      &                       'kg.m-2.s-1',0,zdqssed(:,igcm_dust_number))
    1998             call WRITEDIAGFI(ngridmx,'dqsdif q','diffusion q',
    1999      &                       'kg.m-2.s-1',0,zdqsdif(:,igcm_dust_number))
     2095        call WRITEDIAGFI(ngridmx,'zdqsed_dustq','sedimentation q',
     2096     &                   'kg.m-2.s-1',1,zdqsed(1,:,igcm_dust_mass))
     2097        call WRITEDIAGFI(ngridmx,'zdqsed_dustN','sedimentation N',
     2098     &                'Nbr.m-2.s-1',1,zdqsed(1,:,igcm_dust_number))
    20002099     
    2001             call WRITEDIAGFI(ngridmx,'dqssed N','sedimentation N',
    2002      &                       'kg.m-2.s-1',0,zdqssed(:,igcm_dust_number))
    2003             call WRITEDIAGFI(ngridmx,'dqsdif N','diffusion N',
    2004      &                       'kg.m-2.s-1',0,zdqsdif(:,igcm_dust_number))
     2100        call WRITEDIAGFI(ngridmx,'zdqcloud*dt ice','cloud ice',
     2101     &            'kg.m-2.s-1',1,zdqcloud(1,:,igcm_h2o_ice)*ptimestep)
     2102        call WRITEDIAGFI(ngridmx,'zdqcloud*dt vap','cloud vap',
     2103     &            'kg.m-2.s-1',1,zdqcloud(1,:,igcm_h2o_vap)*ptimestep)
     2104        call WRITEDIAGFI(ngridmx,'zdqdif*dt ice','dif ice',
     2105     &            'kg.m-2.s-1',1,zdqdif(1,:,igcm_h2o_ice)*ptimestep)
     2106        call WRITEDIAGFI(ngridmx,'zdqdif*dt vap','dif vap',
     2107     &            'kg.m-2.s-1',1,zdqdif(1,:,igcm_h2o_vap)*ptimestep)
     2108        call WRITEDIAGFI(ngridmx,'zdqdif*dt vap 0','dif vap',
     2109     &            'kg.m-2.s-1',0,zdqdif(1,1,igcm_h2o_vap)*ptimestep)
     2110
     2111!       if (scavenging) then
     2112!         call WRITEDIAGFI(ngridmx,'zdqsed_ccnq','sedimentation q',
     2113!     &                   'kg.m-2.s-1',1,zdqsed(1,:,igcm_ccn_mass))
     2114!         call WRITEDIAGFI(ngridmx,'zdqsed_ccnN','sedimentation N',
     2115!    &                 'Nbr.m-2.s-1',1,zdqsed(1,:,igcm_ccn_number))
     2116!      endif
     2117!       call WRITEDIAGFI(ngridmx,'zdqsed_ice','sedimentation q',
     2118!     &                   'kg.m-2.s-1',1,zdqsed(1,:,igcm_h2o_ice))
     2119!     
    20052120     
    2006               call WRITEDIAGFI(ngrid,"rice","ice radius","m",1,
     2121          call WRITEDIAGFI(ngrid,"rice","ice radius","m",1,
    20072122     &                    rice)
    2008          ENDIF
     2123          call WRITEDIAGFI(ngrid,"satu","vap in satu","kg/kg",1,
     2124     &                    satu)
     2125         ENDIF ! of IF (water)
    20092126         
    20102127ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

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

    r420 r520  
    11      subroutine simpleclouds(ngrid,nlay,ptimestep,
    22     &             pplev,pplay,pzlev,pzlay,pt,pdt,
    3      &             pq,pdq,pdqcloud,pdqscloud,pdtcloud,
     3     &             pq,pdq,pdqcloud,pdtcloud,
    44     &             nq,tau,rice,nuice,rsedcloud)
    55      implicit none
    … …  
    6262      real pdqcloud(ngrid,nlay,nq) ! tendance de la condensation
    6363                                   !   H2O(kg/kg.s-1)
    64       real pdqscloud(ngrid,nq)     ! flux en surface (kg.m-2.s-1)
    6564      REAL pdtcloud(ngrid,nlay)    ! tendance temperature due
    6665                                   !   a la chaleur latente
    … …  
    141140
    142141        enddo ! of do ig=1,ngrid
    143       enddo ! of dol=1,nlay
    144 
    145       pdqscloud(1:ngrid,1:nq)=0
     142      enddo ! of do l=1,nlay
     143
    146144      pdqcloud(1:ngrid,1:nlay,1:nq)=0
    147145      pdtcloud(1:ngrid,1:nlay)=0
    … …  
    211209c------------------------------------------------------------------
    212210c     TEST_JBM
    213       IF (ngrid.eq.1) THEN
    214          call WRITEDIAGFI(ngrid,"mcond","h2o condensed mass","kg",1,
    215      &                    Mcon_out)
    216          call WRITEDIAGFI(ngrid,"rdusttyp","rdusttyp","m",1,
    217      &                    rdusttyp)
    218          call WRITEDIAGFI(ngrid,"ccntyp","ccntyp","kg-1",1,
    219      &                    ccntyp)
    220       ENDIF
     211!      IF (ngrid.eq.1) THEN
     212!         call WRITEDIAGFI(ngrid,"mcond","h2o condensed mass","kg",1,
     213!     &                    Mcon_out)
     214!         call WRITEDIAGFI(ngrid,"rdusttyp","rdusttyp","m",1,
     215!     &                    rdusttyp)
     216!         call WRITEDIAGFI(ngrid,"ccntyp","ccntyp","kg-1",1,
     217!     &                    ccntyp)
     218!      ENDIF
    221219c------------------------------------------------------------------
    222220      return

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

    r317 r520  
    11      SUBROUTINE surfini(ngrid,piceco2,qsurf,psolaralb)
     2   ! to use  'getin'
     3      USE ioipsl_getincom
    24      IMPLICIT NONE
    35c=======================================================================
    … …  
    2022      REAL  piceco2(ngrid),psolaralb(ngrid,2)
    2123      REAL qsurf(ngrid,nqmx) !tracer on surface (kg/m2)
     24      REAL icedryness ! ice dryness
    2225
    2326      EXTERNAL ISMIN,ISMAX
    2427      INTEGER ISMIN,ISMAX
     28     
     29! Choose false to have a somewhat non resolution dependant water ice caps distribution,
     30! i.e based only on lat & lon values of each physical point.
     31! Choose true to get a carefully choosen distribution for GCM resolutions 32x24 or 64x48
     32! For vizualisation : soon ...
     33       LOGICAL,SAVE :: improvedicelocation = .true.
    2534c
    2635c=======================================================================
    … …  
    4655
    4756         alternate = 0
    48 
    49          do ig=1,ngridmx
     57         temptag = .true.
     58         
     59         write(*,*) "Ice dryness ?"
     60         icedryness=1. ! default value
     61         call getin("icedryness",icedryness)
     62         write(*,*) " icedryness = ",icedryness
     63
    5064       
    5165#ifdef MESOSCALE
     66
     67      do ig=1,ngridmx
     68
    5269         !write(*,*) "all qsurf to zero. dirty."
    5370         do iq=1,nqmx
    … …  
    6683                 watercaptag(ig)  = .false.
    6784                 dryness(ig)      = 1.
    68          endif 
     85         endif
     86         
     87      endif
    6988#else
    70 
    71          dryness (ig) = 1
    72 
     89       ! print*,'ngridmx',ngridmx
     90
     91      if (improvedicelocation) then
     92     
     93        if (ngridmx .eq. 738) then ! hopefully 32x24
     94          print*,'water ice caps distribution for 32x24 resolution:'
     95          watercaptag(1:9)    = .true. ! central cap - core
     96          watercaptag(26:33)  = .true. ! central cap
     97           
     98        else if (ngridmx .eq. 3010) then ! hopefully 64x48
     99          print*,'water ice caps distribution for 64x48 resolution:'
     100          watercaptag(1:65)   = .true. ! central cap - core
     101!          watercaptag(85)     = .true. ! central cap
     102          watercaptag(83:85)  = .true. ! central cap - BIGGER
     103          watercaptag(93:114) = .true. ! central cap
     104!          watercaptag(254)    = .true. ! Korolev crater
     105!          watercaptag(255)    = .true. ! Korolev crater --DECALE
     106          watercaptag(242:257)= .true. ! Korolev crater -- VERY BIG
     107!--------------------------------------------------------------------         
     108!          watercaptag(136)    = .true. ! outlier, lat = 78.75
     109!          watercaptag(138)    = .true. ! outlier, lat = 78.75
     110!          watercaptag(140)    = .true. ! outlier, lat = 78.75
     111!          watercaptag(142)    = .true. ! outlier, lat = 78.75
     112!          watercaptag(183)    = .true. ! outlier, lat = 78.75
     113!          watercaptag(185)    = .true. ! outlier, lat = 78.75
     114!          watercaptag(187)    = .true. ! outlier, lat = 78.75
     115!          watercaptag(189)    = .true. ! outlier, lat = 78.75
     116!          watercaptag(191)    = .true. ! outlier, lat = 78.75
     117!          watercaptag(193)    = .true. ! outlier, lat = 78.75
     118!--------------------------------------------------------------------         
     119          watercaptag(135:142)= .true. ! BIG outlier, lat = 78.75
     120          watercaptag(181:193)= .true. ! BIG outlier, lat = 78.75
     121           
     122        else if (ngridmx .ne. 1) then
     123      print*,'No improved ice location for this resolution!'
     124      print*,'Please set improvedicelocation to false in surfini.'
     125      print*,'And check lat and lon conditions for ice caps in code.'
     126      call exit(1)
     127         
     128        endif
     129       
     130        ! print caps locations
     131        print*,'latitude | longitude | ig'
     132        do ig=1,ngridmx
     133          dryness (ig) = icedryness
     134          !!! OU alors tu mets dryness = icedrag sur outliers et 1 ailleurs
     135          !!! OU remettre comme c'était avant Aymeric et pis c'est tout
     136          if (watercaptag(ig)) then
     137             print*,'ice water cap', lati(ig)*180./pi,
     138     .              long(ig)*180./pi, ig
     139          endif
     140        enddo
     141
     142!!------------------------ test -----------------------------
     143!!-----------------------------------------------------------   
     144!!      else if (1) then ! DUMMY !!, CAREFUL !!, TOUT CA ....
     145!!       
     146!!         do ig=1,ngridmx
     147!!          dryness (ig) = 1
     148!!          if (albedodat(ig) .ge. 0.35) then
     149!!             watercaptag(ig)    = .true.
     150!!             print*,' albedo criteria ice water cap', lati(ig)*180./pi,
     151!!     .              long(ig)*180./pi, ig
     152!!          endif
     153!!        enddo
     154!!-----------------------------------------------------------   
     155!!-----------------------------------------------------------   
     156     
     157      else
     158
     159         do ig=1,ngridmx
     160
     161         dryness (ig) = icedryness
     162         
    73163!!c Towards olympia planitia water caps ('relatively' low latitude ones)
    74164!!c---------------- proposition par AS --------------------
    … …  
    98188!!     .           ( long(ig)*180./pi .ge. 155. ) .and.
    99189!!     .           ( long(ig)*180./pi .le. 180. ) )
    100 !!c     .               .or.
    101 !!c     .         ( ( lati(ig)*180./pi .ge. 71.  ) .and. ! cap #4 (Korolev crater)
    102 !!c     .           ( lati(ig)*180./pi .le. 72.  ) .and.
    103 !!c     .           ( long(ig)*180./pi .ge. 163. ) .and.
    104 !!c     .           ( long(ig)*180./pi .le. 168. ) )
     190!!     .                .or.
     191!!     .         ( ( lati(ig)*180./pi .ge. 71.  ) .and. ! cap #4 (Korolev crater)
     192!!     .           ( lati(ig)*180./pi .le. 72.  ) .and.
     193!!     .           ( long(ig)*180./pi .ge. 163. ) .and.
     194!!     .           ( long(ig)*180./pi .le. 168. ) )
    105195!!     .         .or.
    106196!!     .         ( ( lati(ig)*180./pi .ge. 74.9 ) .and. ! cap #5
    … …  
    109199!!     .           ( long(ig)*180./pi .le. -120.) ) )
    110200!!     .         then
    111 !!     
    112 !!             if (temptag) then
    113 !!             
    114 !!               if ((alternate .eq. 0)) then  !!! 1/2 en 64x48 sinon trop large en lat
    115 !!                if (ngridmx.ne.1) watercaptag(ig)=.true.
    116 !!                  write(*,*) "outliers ", lati(ig)*180./pi,
    117 !!     .              long(ig)*180./pi
    118 !!                  !dryness(ig) = 1.
    119 !!                  alternate = 1
    120 !!                else
    121 !!                  alternate = 0
    122 !!                endif !end if alternate = 0
    123 !!               
    124 !!             else
    125 !!             
    126 !!               if (ngridmx.ne.1) watercaptag(ig)=.true.
    127 !!                  write(*,*) "outliers ", lati(ig)*180./pi,
    128 !!     .              long(ig)*180./pi
    129 !!     
    130 !!             endif ! end if temptag
    131 !!             
    132 !!           endif
    133 !!
    134 !!
    135 !!c Opposite olympia planitia water cap
    136 !!c---------------- proposition par AS --------------------
    137 !!c--------------------------------------------------------
    138 !!c           if ( ( lati(ig)*180./pi      .ge.  82 ) .and.
    139 !!c     .          ( lati(ig)*180./pi      .le.  84 ) .and.
    140 !!c     .          ( ( long(ig)*180./pi .gt. -030. ) .and.
    141 !!c     .          ( long(ig)*180./pi .lt.  090. ) ) ) then
    142 !!c---------------- proposition par TN --------------------
    143 !!c--------------------------------------------------------
    144 !!           if ( ( lati(ig)*180./pi     .ge.  80 ) .and.
    145 !!     .          ( lati(ig)*180./pi     .le.  84 ) .and.
    146 !!     .          ( ( long(ig)*180./pi .gt. -030. ) .and.
    147 !!     .            ( long(ig)*180./pi .lt.  090. ) ) ) then
    148 !!              if (ngridmx.ne.1) then
    149 !!                watercaptag(ig)=.true.
    150 !!                write(*,*) "central cap add ", lati(ig)*180./pi,
    151 !!     .            long(ig)*180./pi
    152 !!              endif
    153 !!              !dryness(ig) = 1.
    154 !!           endif
     201       if ( ( ( lati(ig)*180./pi .ge. 77.  ) .and. ! cap #2
     202     .           ( lati(ig)*180./pi .le. 80.  ) .and.
     203     .           ( long(ig)*180./pi .ge. 110. ) .and.
     204     .           ( long(ig)*180./pi .le. 181. ) )
     205     .         .or.
     206
     207     .         ( ( lati(ig)*180./pi .ge. 75.  ) .and. ! cap #4 (Korolev crater)
     208     .           ( lati(ig)*180./pi .le. 76.  ) .and.
     209     .           ( long(ig)*180./pi .ge. 150. ) .and.
     210     .           ( long(ig)*180./pi .le. 168. ) )
     211     .         .or.
     212     .         ( ( lati(ig)*180./pi .ge. 77 ) .and. ! cap #5
     213     .           ( lati(ig)*180./pi .le. 80.  ) .and.
     214     .           ( long(ig)*180./pi .ge. -150.) .and.
     215     .           ( long(ig)*180./pi .le. -110.) ) )
     216     .         then
     217     
     218             if (temptag) then
     219             
     220               if ((alternate .eq. 0)) then  ! 1/2 en 64x48 sinon trop large en lat
     221                if (ngridmx.ne.1) watercaptag(ig)=.true.
     222!                  print*,'ice water cap', lati(ig)*180./pi,
     223!     .              long(ig)*180./pi, ig
     224                  !dryness(ig) = 1.
     225                  alternate = 1
     226                else
     227                  alternate = 0
     228                endif !end if alternate = 0
     229               
     230             else
     231             
     232!               if (ngridmx.ne.1) watercaptag(ig)=.true.
     233!                  write(*,*) "outliers ", lati(ig)*180./pi,
     234!     .              long(ig)*180./pi
     235     
     236             endif ! end if temptag
     237             
     238           endif
     239
     240
     241c Opposite olympia planitia water cap
     242c---------------- proposition par AS --------------------
     243c--------------------------------------------------------
     244c           if ( ( lati(ig)*180./pi      .ge.  82 ) .and.
     245c     .          ( lati(ig)*180./pi      .le.  84 ) .and.
     246c     .          ( ( long(ig)*180./pi .gt. -030. ) .and.
     247c     .          ( long(ig)*180./pi .lt.  090. ) ) ) then
     248c---------------- proposition par TN --------------------
     249c--------------------------------------------------------
     250           if ( ( ( lati(ig)*180./pi     .ge.  80 ) .and.
     251     .            ( lati(ig)*180./pi     .le.  84 ) )
     252     .          .and.
     253     .          ( ( long(ig)*180./pi .lt. -95. ) .or.       !!! 32x24
     254     .            ( long(ig)*180./pi .gt.  85. ) ) ) then   !!! 32x24
     255!     .        ( ( ( long(ig)*180./pi .ge. -29. ) .and.       !!! 64x48
     256!     .            ( long(ig)*180./pi .le.  90. ) ) .or.      !!! 64x48
     257!     .          ( ( long(ig)*180./pi .ge. -77. ) .and.       !!! 64x48
     258!     .            ( long(ig)*180./pi .le. -70. ) ) ) ) then     !!! 64x48
     259              if (ngridmx.ne.1) then
     260                watercaptag(ig)=.true.
     261                 print*,'ice water cap', lati(ig)*180./pi,
     262     .            long(ig)*180./pi, ig
     263              endif
     264              !dryness(ig) = 1.
     265           endif
    155266
    156267c Central cap
    … …  
    159270
    160271           if (lati(ig)*180./pi.gt.84) then
    161              PRINT*,'central cap', lati(ig)*180./pi,
    162      .         long(ig)*180./pi
     272             print*,'ice water cap', lati(ig)*180./pi,
     273     .         long(ig)*180./pi, ig
    163274             if (ngridmx.ne.1) watercaptag(ig)=.true.
    164275             !dryness(ig) = 1.
    … …  
    173284c             endif
    174285           endif  ! (lati>80 deg)
     286           
     287         end do ! (ngridmx)
     288         
     289       endif ! of if (improvedicelocation)
    175290#endif     
    176          end do ! (ngridmx)
    177         ENDIF ! (caps & water)
     291
     292       ENDIF ! (caps & water)
    178293
    179294c ===============================================================
    … …  
    210325           PRINT*,'maximum albedo avec water caps',
    211326     s     psolaralb(ISMAX(ngrid,psolaralb,1),1)
     327
    212328         ENDIF
    213329

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

    r358 r520  
    581581     $   thermo,qsurf)
    582582      endif
    583       watercaptag(ngridmx)=.false.
     583
     584c Regarder si le sol est un reservoir de glace d'eau
     585c --------------------------------------------------
     586      watercaptag(ngridmx)=.false. ! Par defaut il n'y pas de glace au sol
     587      print *,'Water ice cap on ground ?'
     588      call getin("watercaptag",watercaptag)
     589      write(*,*) " watercaptag = ",watercaptag
    584590     
    585591

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

    r501 r520  
    661661           else if (doubleq) then
    662662             do ig=1,ngrid
    663            !!! soulevement constant
     663               if(co2ice(ig).lt.1) then ! soulevement pas constant
    664664                 pdqsdif(ig,igcm_dust_mass) =
    665665     &             -alpha_lift(igcm_dust_mass) 
    666666                 pdqsdif(ig,igcm_dust_number) =
    667      &             -alpha_lift(igcm_dust_number) 
     667     &             -alpha_lift(igcm_dust_number)
     668               end if
    668669             end do
    669670           else if (submicron) then
Note: See TracChangeset for help on using the changeset viewer.