Changeset 3203


Ignore:
Timestamp:
Feb 7, 2024, 12:14:38 PM (11 months ago)
Author:
jbclement
Message:

Mars PCM:

  • Addition of the possibility to use subslopes parametization in 1D. To do so, put 'nslope=x' in "run.def" where x (1, 3, 5 or 7) is the number of slopes you want to, or modify it in an appropriate "startfi.nc". Then, default subslopes definition and distribution are used by the model. The already existing case of using one slope whose inclination and orientation are defined in "run.def" is still possible.
  • Some cleanings throughout the code, in particular related to unused variables.

JBC

Location:
trunk
Files:
10 edited

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.COMMON/libf/evolution/orbit_param_criterion_mod.F90

    r3202 r3203  
    7777do ilask = 1,size(yearlask,1)
    7878    read(73,*) yearlask(ilask), obllask(ilask), ecclask(ilask), lsplask(ilask)
    79     yearlask(ilask) = yearlask(ilask)*10./convert_years ! Conversion from Laskar's kilo Earth years to PEM Martian years
     79    yearlask(ilask) = yearlask(ilask)*1000./convert_years ! Conversion from Laskar's kilo Earth years to PEM Martian years
    8080    if (yearlask(ilask) > Year) last_ilask = ilask + 1
    8181enddo

  • trunk/LMDZ.COMMON/libf/evolution/pem.F90

    r3199 r3203  
    320320    endif
    321321
    322     call init_testphys1d('start1D_evol.txt','startfi_evol.nc',.true.,therestart1D,therestartfi,ngrid,nlayer,610., &
    323                          nq,q,time_0,ps(1),ucov,vcov,teta,ndt,ptif,pks,dtphys,zqsat,dq,dqdyn,day0,day,gru,grv,w,  &
     322    call init_testphys1d('start1D_evol.txt','startfi_evol.nc',therestart1D,therestartfi,ngrid,nlayer,610.,nq,q, &
     323                         time_0,ps(1),ucov,vcov,teta,ndt,ptif,pks,dtphys,zqsat,dq,dqdyn,day0,day,gru,grv,w,     &
    324324                         play,plev,latitude,longitude,cell_area,atm_wat_profile,atm_wat_tau,CO2cond_ps)
    325325    ps(2) = ps(1)

  • trunk/LMDZ.COMMON/libf/evolution/pemetat0.F90

    r3190 r3203  
    253253                tsoil_tmp_yr2(:,:,islope) = tsoil_tmp_yr1(:,:,islope)
    254254                call compute_tsoil_pem(ngrid,nsoil_PEM,.false.,TI_PEM(:,:,islope),timestep,tsurf_avg_yr2(:,islope),tsoil_tmp_yr2(:,:,islope))
    255                
     255
    256256                do iloop = nsoil_PCM+1,nsoil_PEM
    257257                    tsoil_PEM(:,iloop,islope) = tsoil_tmp_yr2(:,iloop,islope)

  • trunk/LMDZ.MARS/changelog.txt

    r3201 r3203  
    44664466== 05/02/2024 == JBC
    44674467Small modification of the change introduced in r3200 to make the code simpler and faster.
     4468
     4469== 07/02/2024 == JBC
     4470- Addition of the possibility to use subslopes parametization in 1D. To do so, put 'nslope=x' in "run.def" where x (1, 3, 5 or 7) is the number of slopes you want to, or modify it in an appropriate "startfi.nc". Then, default subslopes definition and distribution are used by the model. The already existing case of using one slope whose inclination and orientation are defined in "run.def" is still possible.
     4471- Some cleanings throughout the code, in particular related to unused variables.

  • trunk/LMDZ.MARS/libf/dynphy_lonlat/phymars/newstart.F

    r3183 r3203  
    17611761      select case(nslope)
    17621762          case(1)
     1763              default_def_slope(1) = -90.
     1764              default_def_slope(2) = 90.
     1765          case(3)
    17631766              default_def_slope(1) = -50.
    1764               default_def_slope(2) = 50.
     1767              default_def_slope(2) = -3.
     1768              default_def_slope(3) = 3.
     1769              default_def_slope(4) = 50.
    17651770          case(5)
    17661771              default_def_slope(1) = -43.
    … …  
    17811786          case default
    17821787              write(*,*) 'Number of slopes not possible: nslope should
    1783      & be 1, 5 or 7!'
     1788     & be 1, 3, 5 or 7!'
    17841789              call abort
    17851790      end select
    … …  
    18011806       ENDDO
    18021807
    1803        call subslope_mola(ngridmx,nslope_new,def_slope,subslope_dist)
     1808       if (ngridmx /= 1) then
     1809         call subslope_mola(ngridmx,nslope_new,def_slope,subslope_dist)
     1810       endif
    18041811
    18051812! Surfdat related stuff

  • trunk/LMDZ.MARS/libf/phymars/dyn1d/init_testphys1d_mod.F90

    r3200 r3203  
    55contains
    66
    7 SUBROUTINE init_testphys1d(start1Dname,startfiname,startfiles_1D,therestart1D,therestartfi,ngrid,nlayer,odpref, &
     7SUBROUTINE init_testphys1d(start1Dname,startfiname,therestart1D,therestartfi,ngrid,nlayer,odpref, &
    88                           nq,q,time,psurf,u,v,temp,ndt,ptif,pks,dttestphys,zqsat,dq,dqdyn,day0,day,gru,grv,w,  &
    99                           play,plev,latitude,longitude,cell_area,atm_wat_profile,atm_wat_tau,CO2cond_ps)
    … …  
    1515use surfdat_h,                only: albedodat, z0_default, z0, emissiv, emisice, albedice, iceradius, dtemisice,      &
    1616                                    zmea, zstd, zsig, zgam, zthe, hmons, summit, base, phisfi, watercaptag, watercap, &
    17                                     tsurf, emis, qsurf, perennial_co2ice
     17                                    tsurf, emis, qsurf, perennial_co2ice, ini_surfdat_h_slope_var, end_surfdat_h_slope_var
    1818use infotrac,                 only: nqtot, tname, nqperes, nqfils
    1919use read_profile_mod,         only: read_profile
    2020use iostart,                  only: open_startphy, get_var, close_startphy
    2121use physics_distribution_mod, only: init_physics_distribution
    22 use comsoil_h,                only: volcapa, nsoilmx, inertiesoil, inertiedat, layer, mlayer, flux_geo, tsoil, qsoil
     22use comsoil_h,                only: volcapa, nsoilmx, inertiesoil, inertiedat, layer, mlayer, flux_geo, tsoil, qsoil, &
     23                                    ini_comsoil_h_slope_var, end_comsoil_h_slope_var
    2324use comvert_mod,              only: ap, bp, aps, bps, pa, preff, presnivs, pseudoalt, scaleheight
    24 use dimradmars_mod,           only: tauvis, totcloudfrac, albedo
     25use dimradmars_mod,           only: tauvis, totcloudfrac, albedo, ini_dimradmars_mod_slope_var, end_dimradmars_mod_slope_var
    2526use regular_lonlat_mod,       only: init_regular_lonlat
    2627use mod_interface_dyn_phys,   only: init_interface_dyn_phys
    … …  
    4041use nonoro_gwd_ran_mod,       only: du_nonoro_gwd, dv_nonoro_gwd
    4142use conf_phys_mod,            only: conf_phys
    42 use paleoclimate_mod,         only: paleoclimate
     43use paleoclimate_mod,         only: paleoclimate, ini_paleoclimate_h, end_paleoclimate_h
    4344use comslope_mod,             only: nslope, subslope_dist, ini_comslope_h, end_comslope_h
    4445! Mostly for XIOS outputs:
    … …  
    5455!=======================================================================
    5556integer,      intent(in) :: ngrid, nlayer
    56 real,         intent(in) :: odpref                                    ! DOD reference pressure (Pa)
    57 character(*), intent(in) :: start1Dname, startfiname                  ! Name of starting files for 1D
    58 logical,      intent(in) :: startfiles_1D, therestart1D, therestartfi ! Use of starting files for 1D
     57real,         intent(in) :: odpref                     ! DOD reference pressure (Pa)
     58character(*), intent(in) :: start1Dname, startfiname   ! Name of starting files for 1D
     59logical,      intent(in) :: therestart1D, therestartfi ! Use of starting files for 1D
    5960
    6061integer, intent(inout) :: nq
    … …  
    9697character(30)          :: header
    9798real, dimension(100)   :: tab_cntrl
    98 real, dimension(1,2,1) :: albedo_read ! surface albedo
    9999
    100100! New flag to compute paleo orbital configurations + few variables JN
    … …  
    103103
    104104! MVals: isotopes as in the dynamics (CRisi)
    105 integer                                  :: ifils, ipere, generation, ierr0
     105integer                                  :: ipere, ierr0
    106106character(30), dimension(:), allocatable :: tnom_transp ! transporting fluid short name
    107107character(80)                            :: line        ! to store a line of text
    … …  
    112112real    :: inertieice = 2100. ! ice thermal inertia
    113113integer :: iref
    114 
    115 ! LL: Subsurface geothermal flux
    116 real :: flux_geo_tmp
    117114
    118115!=======================================================================
    … …  
    122119! Prescribed constants to be set here
    123120!------------------------------------------------------
    124 
    125121! Mars planetary constants
    126122! ------------------------
    … …  
    147143lmixmin = 30       ! mixing length ~100
    148144
    149 ! cap properties and surface emissivities
     145! Cap properties and surface emissivities
    150146! ---------------------------------------
    151147emissiv = 0.95         ! Bare ground emissivity ~.95
    … …  
    159155dtemisice(2) = 2.      ! time scale for snow metamorphism (south)
    160156
    161 ! mesh surface (not a very usefull quantity in 1D)
     157! Mesh surface (not a very usefull quantity in 1D)
    162158! ------------------------------------------------
    163159cell_area(1) = 1.
    164160
    165 ! check if there is a 'traceur.def' file and process it
    166 ! load tracer names from file 'traceur.def'
     161! Check if there is a 'traceur.def' file and process it
     162! Load tracer names from file 'traceur.def'
     163! -----------------------------------------------------
    167164open(90,file = 'traceur.def',status = 'old',form = 'formatted',iostat = ierr)
    168165if (ierr /= 0) then
    … …  
    188185endif
    189186
    190 ! allocate arrays:
     187! Allocate arrays
    191188allocate(tname(nq),q(1,nlayer,nq),zqsat(nlayer))
    192189allocate(dq(1,nlayer,nq),dqdyn(1,nlayer,nq),tnom_transp(nq))
    193190
    194 ! read tracer names from file traceur.def
     191! Read tracer names from file traceur.def
    195192do iq = 1,nq
    196193    read(90,'(80a)',iostat = ierr) line ! store the line from traceur.def
    … …  
    281278nsoil = nsoilmx
    282279
    283 day_step = 48 ! default value for day_step
     280day_step = 48 ! Default value for day_step
    284281write(*,*)'Number of time steps per sol?'
    285282call getin("day_step",day_step)
    286283write(*,*) " day_step = ",day_step
    287284
    288 ecritphy = day_step ! default value for ecritphy, output every time step
    289 
    290 ndt = 10 ! default value for ndt
     285ecritphy = day_step ! Default value for ecritphy, output every time step
     286
     287ndt = 10 ! Default value for ndt
    291288write(*,*)'Number of sols to run?'
    292289call getin("ndt",ndt)
    … …  
    311308! Coefficient to control the surface pressure change
    312309! To be defined in "callphys.def" so that the PEM can read it asa well
     310! --------------------------------------------------------------------
    313311CO2cond_ps = 1. ! Default value
    314312write(*,*) 'Coefficient to control the surface pressure change?'
    … …  
    320318
    321319! Reference pressures
    322 pa = 20.     ! transition pressure (for hybrid coord.)
    323 preff = 610. ! reference surface pressure
     320pa = 20.     ! Transition pressure (for hybrid coord.)
     321preff = 610. ! Reference surface pressure
    324322
    325323! Aerosol properties
    326324! ------------------
    327 tauvis = 0.2 ! default value for tauvis (dust opacity)
     325tauvis = 0.2 ! Default value for tauvis (dust opacity)
    328326write(*,'("Reference dust opacity at ",f4.0," Pa?")')odpref
    329327call getin("tauvis",tauvis)
    … …  
    408406endif
    409407
    410 ! Initialize tracers (surface and atmosphere) here:
     408! Initialize local slope parameters
     409! ---------------------------------
     410nslope = 1 ! default: one slope
     411if (.not. therestartfi) then
     412    ! Number of subgrid scale slopes
     413    write(*,*)'Number of slopes?'
     414    call getin("nslope",nslope)
     415    write(*,*) " nslope = ", nslope
     416    call end_comslope_h
     417    call ini_comslope_h(1,nslope)
     418    select case (nslope)
     419        case(1)
     420            ! Sub-slopes parameters (minimum/maximun angles)
     421            def_slope(1) = -90.
     422            def_slope(2) = 90.
     423            ! Fraction of subslopes in mesh
     424            subslope_dist(1,1) = 1.
     425        case(3)
     426            ! Sub-slopes parameters (minimum/maximun angles)
     427            def_slope(1) = -50.
     428            def_slope(2) = -3.
     429            def_slope(3) = 3.
     430            def_slope(4) = 50.
     431            ! Fraction of subslopes in mesh
     432            subslope_dist(1,1) = 0.1
     433            subslope_dist(1,2) = 0.8
     434            subslope_dist(1,3) = 0.1
     435        case(5)
     436            ! Sub-slopes parameters (minimum/maximun angles)
     437            def_slope(1) = -43.
     438            def_slope(2) = -9.
     439            def_slope(3) = -3.
     440            def_slope(4) = 3.
     441            def_slope(5) = 9.
     442            def_slope(6) = 43.
     443            ! Fraction of subslopes in mesh
     444            subslope_dist(1,1) = 0.025
     445            subslope_dist(1,2) = 0.075
     446            subslope_dist(1,3) = 0.8
     447            subslope_dist(1,4) = 0.075
     448            subslope_dist(1,5) = 0.025
     449        case(7)
     450            ! Sub-slopes parameters (minimum/maximun angles)
     451            def_slope(1) = -43.
     452            def_slope(2) = -19.
     453            def_slope(3) = -9.
     454            def_slope(4) = -3.
     455            def_slope(5) = 3.
     456            def_slope(6) = 9.
     457            def_slope(7) = 19.
     458            def_slope(8) = 43.
     459            ! Fraction of subslopes in mesh
     460            subslope_dist(1,1) = 0.01
     461            subslope_dist(1,2) = 0.02
     462            subslope_dist(1,3) = 0.06
     463            subslope_dist(1,4) = 0.8
     464            subslope_dist(1,5) = 0.06
     465            subslope_dist(1,6) = 0.02
     466            subslope_dist(1,7) = 0.01
     467        case default
     468            error stop 'Number of slopes not possible: nslope should be 1, 3, 5 or 7!'
     469    end select
     470
     471    call end_surfdat_h_slope_var
     472    call ini_surfdat_h_slope_var(1,nq,nslope)
     473    call end_comsoil_h_slope_var
     474    call ini_comsoil_h_slope_var(1,nslope)
     475    call end_dimradmars_mod_slope_var
     476    call ini_dimradmars_mod_slope_var(1,nslope)
     477    call end_paleoclimate_h
     478    call ini_paleoclimate_h(1,nslope)
     479
     480    ! Fraction of subslopes in mesh
     481    !write(*,*)'What subslope distribution do you want to use?'
     482    !call getin("slope_distribution",subslope_dist)
     483    !write(*,*) " subslope_dist = ", subslope_dist(1,:)
     484endif
     485
     486! Slope inclination angle (deg) 0: horizontal, 90: vertical
     487theta_sl = 0. ! default: no inclination
     488call getin("slope_inclination",theta_sl)
     489
     490! Slope orientation (deg)
     491! 0 == Northward, 90 == Eastward, 180 == Southward, 270 == Westward
     492psi_sl = 0. ! default value
     493call getin("slope_orientation",psi_sl)
     494
     495! Initialize tracers (surface and atmosphere)
     496!--------------------------------------------
    411497write(*,*) "init_testphys1d: initializing tracers"
    412498if (.not. therestart1D) then
    413499    call read_profile(nq,nlayer,qsurf(1,:,1),q)
     500    do iq = 1,nq
     501        qsurf(1,iq,:) = qsurf(1,iq,1)
     502    enddo
    414503else
    415504    do iq = 1,nq
    … …  
    430519    call getin("albedo",albedodat(1))
    431520    write(*,*) " albedo = ",albedodat(1)
    432     albedo(1,:,1) = albedodat(1)
     521    albedo(1,:,:) = albedodat(1)
    433522
    434523    inertiedat(1,1) = 400 ! default value for inertiedat
    … …  
    446535    write(*,*) " z0 = ",z0(1)
    447536endif !(.not. therestartfi)
    448 
    449 ! Initialize local slope parameters (only matters if "callslope"
    450 ! is .true. in callphys.def)
    451 ! slope inclination angle (deg) 0: horizontal, 90: vertical
    452 theta_sl(1) = 0. ! default: no inclination
    453 call getin("slope_inclination",theta_sl(1))
    454 ! slope orientation (deg)
    455 ! 0 == Northward, 90 == Eastward, 180 == Southward, 270 == Westward
    456 psi_sl(1) = 0. ! default value
    457 call getin("slope_orientation",psi_sl(1))
    458 
    459 ! Sub-slopes parameters (assuming no sub-slopes distribution for now).
    460 def_slope(1) = -90 ! minimum slope angle
    461 def_slope(2) = 90  ! maximum slope angle
    462 subslope_dist(1,1) = 1 ! fraction of subslopes in mesh
    463537
    464538! For the gravity wave scheme
    … …  
    508582
    509583! Initialize winds for first time step
     584!-------------------------------------
    510585if (.not. therestart1D) then
    511586    u = gru
    … …  
    518593
    519594! Initialize turbulent kinetic energy
     595!------------------------------------
    520596q2 = 0.
    521597
    522598! CO2 ice on the surface
    523599! ----------------------
    524 
    525600if (.not. therestartfi) then
    526     qsurf(1,igcm_co2,1) = 0. ! default value for co2ice
     601    qsurf(1,igcm_co2,:) = 0. ! default value for co2ice
    527602    write(*,*)'Initial CO2 ice on the surface (kg.m-2)'
    528603    call getin("co2ice",qsurf(1,igcm_co2,1))
     604    qsurf(1,igcm_co2,:) = qsurf(1,igcm_co2,1)
    529605    write(*,*) " co2ice = ",qsurf(1,igcm_co2,1)
    530606endif !(.not. therestartfi)
    … …  
    533609! ----------
    534610if (.not. therestartfi) then
    535     emis(:,1) = emissiv
     611    emis = emissiv
    536612    if (qsurf(1,igcm_co2,1) == 1.) then
    537         emis(:,1) = emisice(1) ! northern hemisphere
    538         if (latitude(1) < 0) emis(:,1) = emisice(2) ! southern hemisphere
     613        emis = emisice(1) ! northern hemisphere
     614        if (latitude(1) < 0) emis = emisice(2) ! southern hemisphere
    539615    endif
    540616endif !(.not. therestartfi)
    … …  
    568644
    569645if (.not. therestart1D) then
    570     tsurf(:,1) = tmp2(0)
     646    tsurf = tmp2(0)
    571647    temp = tmp2(1:)
    572648else
    … …  
    610686    endif ! ice_depth > 0
    611687
    612     inertiesoil(1,:,1) = inertiedat(1,:)
    613     tsoil(:,:,1) = tsurf(1,1) ! soil temperature
     688    do isoil = 1,nsoil
     689        inertiesoil(1,isoil,:) = inertiedat(1,isoil)
     690        tsoil(1,isoil,:) = tsurf(1,:) ! soil temperature
     691    enddo
    614692endif !(.not. therestartfi)
    615693
    616 flux_geo_tmp = 0.
    617 call getin("flux_geo",flux_geo_tmp)
    618 flux_geo(:,:) = flux_geo_tmp
     694! Initialize subsurface geothermal flux
     695! -------------------------------------
     696flux_geo = 0.
     697call getin("flux_geo",flux_geo(1,1))
     698write(*,*) " flux_geo = ",flux_geo(1,1)
     699flux_geo = flux_geo(1,1)
    619700
    620701! Initialize soil content
    621 ! -----------------
     702! -----------------------
    622703if (.not. therestartfi) qsoil = 0.
    623704

  • trunk/LMDZ.MARS/libf/phymars/dyn1d/testphys1d.F90

    r3201 r3203  
    77use watersat_mod,        only: watersat
    88use tracer_mod,          only: igcm_h2o_vap, igcm_h2o_ice, igcm_co2, noms
    9 use comcstfi_h,          only: pi, rad, omeg, g, mugaz, rcp, r, cpp
     9use comcstfi_h,          only: pi, g, rcp, cpp
    1010use time_phylmdz_mod,    only: daysec, day_step
    1111use dimradmars_mod,      only: tauvis, totcloudfrac, albedo
    1212use dust_param_mod,      only: tauscaling
    13 use comvert_mod,         only: ap, bp, aps, bps, pa, preff, sig
     13use comvert_mod,         only: ap, bp, aps, bps
    1414use physiq_mod,          only: physiq
    1515use turb_mod,            only: q2
    … …  
    6363integer, parameter                  :: nlayer = llm
    6464real, parameter                     :: odpref = 610. ! DOD reference pressure (Pa)
    65 integer                             :: unitstart     ! unite d'ecriture de "startfi"
    66 integer                             :: ndt, ilayer, ilevel, isoil, idt, iq
     65integer                             :: ndt, ilayer, idt
    6766logical                             :: firstcall, lastcall
    6867integer                             :: day0          ! initial (sol ; =0 at Ls=0)
    … …  
    8180
    8281! Physical and dynamical tendencies (e.g. m.s-2, K/s, Pa/s)
    83 real, dimension(nlayer)             :: du, dv, dtemp, dudyn, dvdyn, dtempdyn
     82real, dimension(nlayer)             :: du, dv, dtemp
    8483real, dimension(1)                  :: dpsurf
    8584real, dimension(:,:,:), allocatable :: dq, dqdyn
    8685
    8786! Various intermediate variables
    88 integer                 :: aslun
    89 real                    :: zls, pks, ptif, qtotinit, qtot
     87real                    :: zls, pks, ptif
    9088real, dimension(nlayer) :: phi, h, s, w
    9189integer                 :: nq = 1 ! number of tracers
    9290real, dimension(1)      :: latitude, longitude, cell_area
    93 character(2)            :: str2
    94 character(7)            :: str7
    95 character(44)           :: txt
    9691
    9792! RV & JBC: Use of starting files for 1D
    … …  
    148143endif
    149144
    150 call init_testphys1d('start1D.txt','startfi.nc',startfiles_1D,therestart1D,therestartfi,ngrid,nlayer,odpref, &
    151                      nq,q,time,psurf,u,v,temp,ndt,ptif,pks,dttestphys,zqsat,dq,dqdyn,day0,day,gru,grv,w,    &
     145call init_testphys1d('start1D.txt','startfi.nc',therestart1D,therestartfi,ngrid,nlayer,odpref,nq,q, &
     146                     time,psurf,u,v,temp,ndt,ptif,pks,dttestphys,zqsat,dq,dqdyn,day0,day,gru,grv,w, &
    152147                     play,plev,latitude,longitude,cell_area,atm_wat_profile,atm_wat_tau,CO2cond_ps)
    153148

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

    r3185 r3203  
    3939     &                      igcm_ccnco2_meteor_mass,
    4040     &                      igcm_ccnco2_meteor_number,
    41      &                      rho_ice_co2,nuiceco2_sed,nuiceco2_ref,
    4241     &                      igcm_dust_mass, igcm_dust_number, igcm_h2o2,
    4342     &                      nuice_ref, rho_ice, rho_dust, ref_r0,
    … …  
    5352     &                        cell_area_for_lonlat_outputs,longitude_deg
    5453      use comgeomfi_h, only: sinlon, coslon, sinlat, coslat
    55       use surfdat_h, only: phisfi, albedodat, zmea, zstd, zsig, zgam,
    56      &                     zthe, z0, albedo_h2o_cap,albedo_h2o_frost,
    57      &                     frost_albedo_threshold,frost_metam_threshold,
    58      &                     tsurf, emis,
    59      &                     capcal, fluxgrd, qsurf,
    60      &                     hmons,summit,base,watercap,watercaptag,
     54      use surfdat_h, only: phisfi, albedodat, z0, albedo_h2o_cap,
     55     &                     albedo_h2o_frost, frost_albedo_threshold,
     56     &                     frost_metam_threshold, tsurf, emis, capcal,
     57     &                     fluxgrd, qsurf, watercap, watercaptag,
    6158     &                     perennial_co2ice
    62       use comsaison_h, only: dist_sol, declin, zls, 
     59      use comsaison_h, only: dist_sol, declin, zls,
    6360     &                       mu0, fract, local_time
    6461      use solarlong_mod, only: solarlong
    … …  
    6663      use nirco2abs_mod, only: nirco2abs
    6764      use slope_mod, only: theta_sl, psi_sl, getslopes, param_slope
    68       use conc_mod, only: init_r_cp_mu, update_r_cp_mu_ak, rnew, 
     65      use conc_mod, only: init_r_cp_mu, update_r_cp_mu_ak, rnew,
    6966     &                    cpnew, mmean
    7067      use time_phylmdz_mod, only: iphysiq, day_step, ecritstart, daysec
    7168      use dimradmars_mod, only: aerosol, totcloudfrac,
    7269     &                          dtrad, fluxrad_sky, fluxrad, albedo,
    73      &                          naerkind, iaer_dust_doubleq, 
     70     &                          naerkind, iaer_dust_doubleq,
    7471     &                          iaer_stormdust_doubleq, iaer_h2o_ice,
    7572     &                          flux_1AU
    … …  
    7875     &                          dustscaling_mode, dust_rad_adjust,
    7976     &                          freedust, reff_driven_IRtoVIS_scenario
    80       use turb_mod, only: q2, wstar, ustar, sensibFlux, 
     77      use turb_mod, only: q2, wstar, ustar, sensibFlux,
    8178     &                    zmax_th, hfmax_th, turb_resolved
    8279      use planete_h, only: aphelie, periheli, year_day, peri_day,
    8380     &                     obliquit
    8481      use planete_h, only: iniorbit
    85       USE comcstfi_h, only: r, cpp, mugaz, g, rcp, pi, rad 
     82      USE comcstfi_h, only: r, cpp, mugaz, g, rcp, pi, rad
    8683      USE calldrag_noro_mod, ONLY: calldrag_noro
    8784      USE vdifc_mod, ONLY: vdifc
    88       use param_v4_h, only: nreact,n_avog, 
     85      use param_v4_h, only: nreact,n_avog,
    8986     &                      fill_data_thermos, allocate_param_thermos
    9087      use iono_h, only: allocate_param_iono
    … …  
    110107#endif
    111108
    112 #ifdef CPP_XIOS     
     109#ifdef CPP_XIOS
    113110      use xios_output_mod, only: initialize_xios_output,
    114111     &                           update_xios_timestep,
    … …  
    116113      use wxios, only: wxios_context_init, xios_context_finalize
    117114#endif
    118       USE mod_grid_phy_lmdz, ONLY: grid_type, unstructured, 
     115      USE mod_grid_phy_lmdz, ONLY: grid_type, unstructured,
    119116     &                             regular_lonlat
    120117      use ioipsl_getin_p_mod, only: getin_p
    … …  
    132129c   --------
    133130c
    134 c   Organisation of the physical parametrisations of the LMD 
     131c   Organisation of the physical parametrisations of the LMD
    135132c   martian atmospheric general circulation model.
    136133c
    … …  
    164161c           - Saving statistics (if "callstats = .true.")
    165162c           - Dumping eof (if "calleofdump = .true.")
    166 c           - Output any needed variables in "diagfi" 
     163c           - Output any needed variables in "diagfi"
    167164c     11. Diagnostic: mass conservation of tracers
    168 c 
    169 c   author: 
    170 c   ------- 
    171 c           Frederic Hourdin    15/10/93
    172 c           Francois Forget             1994
    173 c           Christophe Hourdin  02/1997
     165c
     166c   author:
     167c   -------
     168c           Frederic Hourdin    15/10/93
     169c           Francois Forget     1994
     170c           Christophe Hourdin  02/1997
    174171c           Subroutine completly rewritten by F.Forget (01/2000)
    175172c           Introduction of the photochemical module: S. Lebonnois (11/2002)
    … …  
    206203c    pt(ngrid,nlayer)      Temperature (K)
    207204c    pq(ngrid,nlayer,nq)   Advected fields
    208 c    pudyn(ngrid,nlayer)    | 
     205c    pudyn(ngrid,nlayer)    |
    209206c    pvdyn(ngrid,nlayer)    | Dynamical temporal derivative for the
    210207c    ptdyn(ngrid,nlayer)    | corresponding variables
    … …  
    250247      REAL,INTENT(in) :: pt(ngrid,nlayer) ! temperature (K)
    251248      REAL,INTENT(in) :: pq(ngrid,nlayer,nq) ! tracers (.../kg_of_air)
    252       REAL,INTENT(in) :: flxw(ngrid,nlayer) ! vertical mass flux (ks/s)
    253                                             ! at lower boundary of layer
     249      REAL,INTENT(in) :: flxw(ngrid,nlayer) ! vertical mass flux (ks/s) at lower boundary of layer
    254250
    255251c   outputs:
    … …  
    262258      REAL,INTENT(out) :: pdpsrf(ngrid) ! surface pressure tendency (Pa/s)
    263259
    264      
    265 
    266260c Local saved variables:
    267261c ----------------------
    … …  
    275269      REAL pq_tmp(ngrid, nlayer, 2) ! To compute tendencies due the dust bomb
    276270#endif
    277            
     271
    278272c     Variables used by the water ice microphysical scheme:
    279273      REAL rice(ngrid,nlayer)    ! Water ice geometric mean radius (m)
    … …  
    284278      real rsedcloudco2(ngrid,nlayer) ! CO2 Cloud sedimentation radius
    285279      real rhocloudco2(ngrid,nlayer)  ! CO2 Cloud density (kg.m-3)
    286       real nuiceco2(ngrid,nlayer) ! Estimated effective variance of the 
     280      real nuiceco2(ngrid,nlayer) ! Estimated effective variance of the
    287281                                  !   size distribution
    288282      REAL inertiesoil_tifeedback(ngrid,nsoilmx,nslope) ! Time varying subsurface
    … …  
    293287      real zdqssed_co2(ngrid)  ! CO2 flux at the surface  (kg.m-2.s-1)
    294288      real zdqssed_ccn(ngrid,nq)  ! CCN flux at the surface  (kg.m-2.s-1)
    295       real zcondicea_co2microp(ngrid,nlayer)
     289      real, dimension(ngrid,nlayer) :: zcondicea_co2microp
    296290c     Variables used by the photochemistry
    297291      REAL surfdust(ngrid,nlayer) ! dust surface area (m2/m3, if photochemistry)
    298292      REAL surfice(ngrid,nlayer)  !  ice surface area (m2/m3, if photochemistry)
    299293c     Variables used by the slope model
    300       REAL sl_ls, sl_lct, sl_lat
     294      REAL sl_lct, sl_lat
    301295      REAL sl_tau, sl_alb, sl_the, sl_psi
    302296      REAL sl_fl0, sl_flu
    … …  
    309303c Local variables :
    310304c -----------------
    311 
    312305      REAL CBRT
    313306      EXTERNAL CBRT
    314307
    315 !      CHARACTER*80 fichier 
    316       INTEGER l,ig,ierr,igout,iq,tapphys,isoil
     308!      CHARACTER*80 fichier
     309      INTEGER l,ig,ierr,igout,iq,isoil
    317310
    318311      REAL fluxsurf_lw(ngrid,nslope)      !incident LW (IR) surface flux (W.m-2)
    … …  
    330323c     rocket dust storm
    331324      REAL totstormfract(ngrid)     ! fraction of the mesh where the dust storm is contained
    332       logical clearatm              ! clearatm used to calculate twice the radiative 
    333                                     ! transfer when rdstorm is active : 
     325      logical clearatm              ! clearatm used to calculate twice the radiative
     326                                    ! transfer when rdstorm is active :
    334327                                    !            - in a mesh with stormdust and background dust (false)
    335328                                    !            - in a mesh with background dust only (true)
    336329c     entrainment by mountain top dust flows
    337       logical nohmons               ! nohmons used to calculate twice the radiative 
    338                                     ! transfer when topflows is active : 
     330      logical nohmons               ! nohmons used to calculate twice the radiative
     331                                    ! transfer when topflows is active :
    339332                                    !            - in a mesh with topdust and background dust (false)
    340333                                    !            - in a mesh with background dust only (true)
    341      
     334
    342335      REAL tau(ngrid,naerkind)     ! Column dust optical depth at each point
    343336                                   ! AS: TBD: this one should be in a module !
    … …  
    398391      REAL tlaymean ! temporary value of mean layer temperature for zzlay
    399392
    400 
    401393c     Variables for the total dust for diagnostics
    402394      REAL qdusttotal(ngrid,nlayer) !it equals to dust + stormdust
    403 
    404       INTEGER iaer
    405395
    406396c local variables only used for diagnostic (output in file "diagfi" or "stats")
    … …  
    419409      real rtopdust(ngrid,nlayer)   ! topdust geometric mean radius (m)
    420410      integer igmin, lmin
    421       logical tdiag
    422 
    423       real co2col(ngrid)        ! CO2 column
     411
    424412      ! pplev and pplay are dynamical inputs and must not be modified in the physics.
    425413      ! instead, use zplay and zplev :
    426       REAL zplev(ngrid,nlayer+1),zplay(ngrid,nlayer) 
     414      REAL zplev(ngrid,nlayer+1),zplay(ngrid,nlayer)
    427415!      REAL zstress(ngrid),cd
    428       real tmean, zlocal(nlayer)
    429416      real rho(ngrid,nlayer)  ! density
    430417      real vmr(ngrid,nlayer)  ! volume mixing ratio
    … …  
    440427      REAL icetotco2(ngrid)     ! Total mass of co2 ice (kg/m2)
    441428      REAL Nccntot(ngrid)       ! Total number of ccn (nbr/m2)
    442       REAL NccnCO2tot(ngrid)    ! Total number of ccnCO2 (nbr/m2)
    443429      REAL Mccntot(ngrid)       ! Total mass of ccn (kg/m2)
    444430      REAL rave(ngrid)          ! Mean water ice effective radius (m)
    … …  
    469455      REAL wspeed(ngrid,nlayer+1) ! vertical velocity stormdust tracer
    470456      REAL wtop(ngrid,nlayer+1) ! vertical velocity topdust tracer
    471        
    472457      REAL dsodust(ngrid,nlayer) ! density scaled opacity for background dust
    473458      REAL dsords(ngrid,nlayer) ! density scaled opacity for stormdust
    … …  
    475460
    476461c Test 1d/3d scavenging
    477       real h2otot(ngrid)
    478       real hdotot(ngrid)
    479462      REAL satu(ngrid,nlayer)  ! satu ratio for output
    480463      REAL zqsat(ngrid,nlayer) ! saturation
    481       REAL satuco2(ngrid,nlayer)  ! co2 satu ratio for output
    482       REAL zqsatco2(ngrid,nlayer) ! saturation co2
    483 
    484464
    485465! Added for new NLTE scheme
    486 
    487466      real co2vmr_gcm(ngrid,nlayer)
    488467      real n2vmr_gcm(ngrid,nlayer)
    … …  
    492471      real*8  varerr
    493472
    494 C  Non-oro GW drag & Calcul of Brunt-Vaisala freq. (BV2)
    495       REAL ztetalev(ngrid,nlayer)
    496       real zdtetalev(ngrid,nlayer), zdzlev(ngrid,nlayer)
    497       REAL bv2(ngrid,nlayer)    ! BV2 at zlev   
    498473c  Non-oro GW tendencies
    499474      REAL d_u_hin(ngrid,nlayer), d_v_hin(ngrid,nlayer)
    … …  
    509484      REAL pdu_th(ngrid,nlayer),pdv_th(ngrid,nlayer)
    510485      REAL pdt_th(ngrid,nlayer),pdq_th(ngrid,nlayer,nq)
    511       INTEGER lmax_th(ngrid),dimout,n_out,n
     486      INTEGER lmax_th(ngrid),n_out,n
    512487      CHARACTER(50) zstring
    513488      REAL dtke_th(ngrid,nlayer+1)
    514       REAL zcdv(ngrid), zcdh(ngrid)
    515489      REAL, ALLOCATABLE, DIMENSION(:,:) :: T_out
    516490      REAL, ALLOCATABLE, DIMENSION(:,:) :: u_out ! Interpolated teta and u at z_out
    … …  
    519493      REAL, ALLOCATABLE, DIMENSION(:) :: z_out     ! height of interpolation between z0 and z1 [meters]
    520494      REAL tstar(ngrid)  ! friction velocity and friction potential temp
    521       REAL L_mo(ngrid),vhf(ngrid),vvv(ngrid)
     495      REAL vhf(ngrid), vvv(ngrid)
    522496      real qdustrds0(ngrid,nlayer),qdustrds1(ngrid,nlayer)
    523       real qstormrds0(ngrid,nlayer),qstormrds1(ngrid,nlayer)   
     497      real qstormrds0(ngrid,nlayer),qstormrds1(ngrid,nlayer)
    524498      real qdusttotal0(ngrid),qdusttotal1(ngrid)
    525499
    … …  
    529503      real zdtswclf(ngrid,nlayer)
    530504      real zdtlwclf(ngrid,nlayer)
    531       real fluxsurf_lwclf(ngrid)     
     505      real fluxsurf_lwclf(ngrid)
    532506      real fluxsurf_dn_swclf(ngrid,2),fluxsurf_up_swclf(ngrid,2)
    533507      real fluxtop_lwclf(ngrid)
    … …  
    551525      logical,save :: check_physics_inputs=.false.
    552526      logical,save :: check_physics_outputs=.false.
    553      
    554 !$OMP THREADPRIVATE(check_physics_inputs,check_physics_outputs)     
     527
     528!$OMP THREADPRIVATE(check_physics_inputs,check_physics_outputs)
    555529
    556530c Sub-grid scale slopes
    … …  
    580554      REAL :: viscoh                              ! kinematic molecular viscosity for heat
    581555
    582 
    583556c=======================================================================
    584557      pdq(:,:,:) = 0.
    … …  
    588561c  1.1   Initialisation only at first call
    589562c  ---------------------------------------
    590 
    591563      IF (firstcall) THEN
    592564
    … …  
    610582#endif
    611583
    612 c        read startfi 
     584c        read startfi
    613585c        ~~~~~~~~~~~~
    614586#ifndef MESOSCALE
    … …  
    623595     &         def_slope,def_slope_mean,subslope_dist)
    624596
    625 !     Sky view: 
     597!     Sky view:
    626598       DO islope=1,nslope
    627599       sky_slope(islope) = (1.+cos(pi*def_slope_mean(islope)/180.))/2.
    … …  
    646618      write(*,*)'check: tsurf ',tsurf(1,:),tsurf(ngrid,:)
    647619      write(*,*)'check: tsoil ',tsoil(1,1,:),tsoil(ngrid,nsoilmx,:)
    648       write(*,*)'check: inert ',inertiedat(1,1),inertiedat(ngrid,nsoilmx)     
     620      write(*,*)'check: inert ',inertiedat(1,1),inertiedat(ngrid,nsoilmx)
    649621      write(*,*)'check: midlayer,layer ', mlayer(:),layer(:)
    650622      write(*,*)'check: tracernames ', noms
    … …  
    675647              !  mlayer(k)=lay1*alpha**(k-1/2)
    676648              lay1=2.e-4
    677                   alpha=2
     649              alpha=2
    678650              do iloop=0,nsoilmx-1
    679                   mlayer(iloop)=lay1*(alpha**(iloop-0.5))
    680               enddo
     651              mlayer(iloop)=lay1*(alpha**(iloop-0.5))
     652          enddo
    681653              lay1=sqrt(mlayer(0)*mlayer(1))
    682654              alpha=mlayer(1)/mlayer(0)
    … …  
    718690         CALL surfini(ngrid,nslope,qsurf)
    719691         CALL iniorbit(aphelie,periheli,year_day,peri_day,obliquit)
    720 c        initialize soil 
     692c        initialize soil
    721693c        ~~~~~~~~~~~~~~~
    722694         IF (callsoil) THEN
    … …  
    744716         ENDIF
    745717         icount=1
    746          
    747          
    748718
    749719#ifndef MESOSCALE
    … …  
    761731c        Initialize rnew cpnew and mmean as constant
    762732c        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    763          call init_r_cp_mu(ngrid,nlayer)     
     733         call init_r_cp_mu(ngrid,nlayer)
    764734
    765735         if(callnlte.and.nltemodel.eq.2) call nlte_setup
    … …  
    768738
    769739        IF (water.AND.(ngrid.NE.1)) THEN
    770           write(*,*)"physiq: water_param Surface water frost albedo:", 
     740          write(*,*)"physiq: water_param Surface water frost albedo:",
    771741     .                  albedo_h2o_frost
    772           write(*,*)"physiq: water_param Surface watercap albedo:", 
     742          write(*,*)"physiq: water_param Surface watercap albedo:",
    773743     .                  albedo_h2o_cap
    774744        ENDIF
    775745
    776 
    777 
    778746#ifndef MESOSCALE
    779          if(ngrid.ne.1) then
    780            ! no need to compute slopes when in 1D; it is an input
    781          if (callslope) call getslopes(ngrid,phisfi)
    782          endif !1D
    783          if (ecritstart.GT.0) then
     747      ! no need to compute slopes when in 1D; it is an input
     748      if (ngrid /= 1 .and. callslope) call getslopes(ngrid,phisfi)
     749      if (ecritstart.GT.0) then
    784750             call physdem0("restartfi.nc",longitude,latitude,
    785751     &                   nsoilmx,ngrid,nlayer,nq,
    … …  
    787753     &                   albedodat,inertiedat,def_slope,
    788754     &                   subslope_dist)
    789           else
     755      else
    790756             call physdem0("restartfi.nc",longitude,latitude,
    791757     &                   nsoilmx,ngrid,nlayer,nq,
    … …  
    793759     &                   albedodat,inertiedat,def_slope,
    794760     &                   subslope_dist)
    795           endif
     761      endif
    796762
    797763c        Initialize mountain mesh fraction for the entrainment by top flows param.
    798764c        ~~~~~~~~~~~~~~~
    799765        if (topflows) call topmons_setup(ngrid)
    800        
     766
    801767c        Parameterization of the ATKE routine
    802768c        ~~~~~~~~~~~~~~~
    … …  
    808774
    809775#endif
    810                  
     776
    811777#ifdef CPP_XIOS
    812778        ! XIOS outputs
    … …  
    826792c ---------------------------------------------------
    827793c
    828 
    829 #ifdef CPP_XIOS     
     794#ifdef CPP_XIOS
    830795      ! update XIOS time/calendar
    831       call update_xios_timestep     
    832 #endif
    833 
    834 
    835 
    836 
     796      call update_xios_timestep
     797#endif
    837798
    838799c     Initialize various variables
    … …  
    851812      dwatercap(:,:)=0
    852813
    853 
    854 
    855814      call compute_meshgridavg(ngrid,nq,albedo,emis,tsurf,qsurf,
    856815     &       albedo_meshavg,emis_meshavg,tsurf_meshavg,qsurf_meshavg)
    857      
     816
    858817      ! Dust scenario conversion coefficient from IRabs to VISext
    859818      IRtoVIScoef(1:ngrid)=2.6 ! initialized with former value from Montabone et al 2015
    … …  
    863822      pq_tmp(:,:,:)=0
    864823#endif
    865       igout=ngrid/2+1 
     824      igout=ngrid/2+1
    866825
    867826
    … …  
    886845      ps(:) = pplev(:,1)
    887846
    888 
    889847#ifndef MESOSCALE
    890848c-----------------------------------------------------------------------
    891849c    1.2.1 Compute mean mass, cp, and R
    892850c    update_r_cp_mu_ak outputs rnew(ngrid,nlayer), cpnew(ngrid,nlayer)
    893 c       , mmean(ngrid,nlayer) and Akknew(ngrid,nlayer)
     851c   , mmean(ngrid,nlayer) and Akknew(ngrid,nlayer)
    894852c    --------------------------------
    895853
    … …  
    904862c     ponderation des altitudes au niveau des couches en dp/p
    905863cc ------------------------------------------
    906         !Calculation zzlev & zzlay for first layer     
     864    !Calculation zzlev & zzlay for first layer
    907865      DO ig=1,ngrid
    908          zzlay(ig,1)=-(log(pplay(ig,1)/ps(ig)))*rnew(ig,1)*pt(ig,1)/g                                     
    909         zzlev(ig,1)=0
     866         zzlay(ig,1)=-(log(pplay(ig,1)/ps(ig)))*rnew(ig,1)*pt(ig,1)/g
     867        zzlev(ig,1)=0
    910868         zzlev(ig,nlayer+1)=1.e7    ! dummy top of last layer above 10000 km...
    911869         gz(ig,1)=g
    912        
    913         DO l=2,nlayer     
     870
     871        DO l=2,nlayer
    914872         ! compute "mean" temperature of the layer
    915873          if(pt(ig,l) .eq. pt(ig,l-1)) then
    … …  
    918876             tlaymean=(pt(ig,l)- pt(ig,l-1))/log(pt(ig,l)/pt(ig,l-1))
    919877          endif
    920                      
     878
    921879          ! compute gravitational acceleration (at altitude zaeroid(nlayer-1))
    922880          gz(ig,l)=g*(rad**2)/(rad+zzlay(ig,l-1)+(phisfi(ig)/g))**2
    923          
    924           zzlay(ig,l)=zzlay(ig,l-1)- 
     881
     882          zzlay(ig,l)=zzlay(ig,l-1)-
    925883     &   (log(pplay(ig,l)/pplay(ig,l-1))*rnew(ig,l)*tlaymean/gz(ig,l))
    926                  
     884
    927885          z1=(zplay(ig,l-1)+zplev(ig,l))/(zplay(ig,l-1)-zplev(ig,l))
    928886          z2=(zplev(ig,l)+zplay(ig,l))/(zplev(ig,l)-zplay(ig,l))
    … …  
    936894c     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    937895      DO l=1,nlayer
    938          DO ig=1,ngrid 
     896         DO ig=1,ngrid
    939897            zpopsk(ig,l)=(zplay(ig,l)/zplev(ig,1))**rcp
    940898            zh(ig,l)=pt(ig,l)/zpopsk(ig,l)
    … …  
    951909      pw(1:ngrid,nlayer)=0.5*flxw(1:ngrid,nlayer) ! since flxw(nlayer+1)=0
    952910      do l=1,nlayer
    953        pw(1:ngrid,l)=(pw(1:ngrid,l)*r*pt(1:ngrid,l)) / 
     911       pw(1:ngrid,l)=(pw(1:ngrid,l)*r*pt(1:ngrid,l)) /
    954912     &               (pplay(1:ngrid,l)*cell_area(1:ngrid))
    955913       ! NB: here we use r and not rnew since this diagnostic comes
    … …  
    1031989
    1032990                 CALL nlte_tcool(ngrid,nlayer,zplay*9.869e-6,
    1033      $                pt,zzlay,co2vmr_gcm, n2vmr_gcm, covmr_gcm, 
     991     $                pt,zzlay,co2vmr_gcm, n2vmr_gcm, covmr_gcm,
    1034992     $                ovmr_gcm,  zdtnlte,ierr_nlte,varerr )
    1035993                 if(ierr_nlte.gt.0) then
    … …  
    10671025       endif
    10681026#endif
    1069          
     1027
    10701028c          Call main radiative transfer scheme
    10711029c          ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    … …  
    11001058           ! case of sub-grid water ice clouds: callradite for the clear case
    11011059            IF (CLFvarying) THEN
    1102                ! ---> PROBLEMS WITH ALLOCATED ARRAYS 
     1060               ! ---> PROBLEMS WITH ALLOCATED ARRAYS
    11031061               ! (temporary solution in callcorrk: do not deallocate
    11041062               ! if
    11051063               ! CLFvarying ...) ?? AP ??
    1106                clearsky=.true. 
     1064               clearsky=.true.
    11071065               CALL callradite(icount,ngrid,nlayer,nq,zday,zls,pq,
    11081066     &              albedo_meshavg,emis_meshavg,mu0,zplev,zplay,pt,
    … …  
    11261084                  ntf_clf=1.-tf_clf
    11271085                  DO islope=1,nslope
    1128                   fluxsurf_lw(ig,islope) = ntf_clf*fluxsurf_lwclf(ig) 
     1086                  fluxsurf_lw(ig,islope) = ntf_clf*fluxsurf_lwclf(ig)
    11291087     &                          + tf_clf*fluxsurf_lw(ig,islope)
    1130                   fluxsurf_dn_sw(ig,1:2,islope) = 
    1131      &                           ntf_clf*fluxsurf_dn_swclf(ig,1:2) 
     1088                  fluxsurf_dn_sw(ig,1:2,islope) =
     1089     &                           ntf_clf*fluxsurf_dn_swclf(ig,1:2)
    11321090     &                          + tf_clf*fluxsurf_dn_sw(ig,1:2,islope)
    11331091                  ENDDO
    1134                   fluxsurf_up_sw(ig,1:2) = 
    1135      &                           ntf_clf*fluxsurf_up_swclf(ig,1:2) 
     1092                  fluxsurf_up_sw(ig,1:2) =
     1093     &                           ntf_clf*fluxsurf_up_swclf(ig,1:2)
    11361094     &                          + tf_clf*fluxsurf_up_sw(ig,1:2)
    1137                   fluxtop_lw(ig)  = ntf_clf*fluxtop_lwclf(ig) 
     1095                  fluxtop_lw(ig)  = ntf_clf*fluxtop_lwclf(ig)
    11381096     &                                      + tf_clf*fluxtop_lw(ig)
    1139                   fluxtop_dn_sw(ig,1:2)=ntf_clf*fluxtop_dn_swclf(ig,1:2) 
     1097                  fluxtop_dn_sw(ig,1:2)=ntf_clf*fluxtop_dn_swclf(ig,1:2)
    11401098     &                                 + tf_clf*fluxtop_dn_sw(ig,1:2)
    1141                   fluxtop_up_sw(ig,1:2)=ntf_clf*fluxtop_up_swclf(ig,1:2) 
     1099                  fluxtop_up_sw(ig,1:2)=ntf_clf*fluxtop_up_swclf(ig,1:2)
    11421100     &                                 + tf_clf*fluxtop_up_sw(ig,1:2)
    1143                   taucloudtes(ig) = ntf_clf*taucloudtesclf(ig) 
     1101                  taucloudtes(ig) = ntf_clf*taucloudtesclf(ig)
    11441102     &                                      + tf_clf*taucloudtes(ig)
    1145                   zdtlw(ig,1:nlayer) = ntf_clf*zdtlwclf(ig,1:nlayer) 
     1103                  zdtlw(ig,1:nlayer) = ntf_clf*zdtlwclf(ig,1:nlayer)
    11461104     &                                      + tf_clf*zdtlw(ig,1:nlayer)
    1147                   zdtsw(ig,1:nlayer) = ntf_clf*zdtswclf(ig,1:nlayer) 
     1105                  zdtsw(ig,1:nlayer) = ntf_clf*zdtswclf(ig,1:nlayer)
    11481106     &                                      + tf_clf*zdtsw(ig,1:nlayer)
    11491107               ENDDO
    11501108
    11511109            ENDIF ! (CLFvarying)
    1152            
     1110
    11531111!============================================================================
    1154            
     1112
    11551113#ifdef DUSTSTORM
    11561114!! specific case: compute the added quantity of dust for perturbation
    11571115       if (firstcall) then
    1158          pdq(1:ngrid,1:nlayer,igcm_dust_mass)= 
    1159      &      pdq(1:ngrid,1:nlayer,igcm_dust_mass)     
     1116         pdq(1:ngrid,1:nlayer,igcm_dust_mass)=
     1117     &      pdq(1:ngrid,1:nlayer,igcm_dust_mass)
    11601118     &      - pq_tmp(1:ngrid,1:nlayer,1)
    11611119     &      + pq(1:ngrid,1:nlayer,igcm_dust_mass)
    … …  
    11891147     &      "physiq","callslope=true but nslope.ne.1",1)
    11901148            endif
    1191             print *, 'Slope scheme is on and computing...'
    1192             DO ig=1,ngrid 
     1149            write(*,*) 'Slope scheme is on and computing...'
     1150            DO ig=1,ngrid
    11931151              sl_the = theta_sl(ig)
    11941152              IF (sl_the .ne. 0.) THEN
    … …  
    12051163                 sl_di0 = 0.
    12061164                 if ((mu0(ig) .gt. 0.).and.(ztim1.gt.0.)) then
    1207                           sl_di0 = mu0(ig)*(exp(-sl_tau/mu0(ig)))
     1165                  sl_di0 = mu0(ig)*(exp(-sl_tau/mu0(ig)))
    12081166                  sl_di0 = sl_di0*flux_1AU/dist_sol/dist_sol
    12091167                  sl_di0 = sl_di0/ztim1
    … …  
    12131171                 if (sl_fl0 .lt. sl_di0) sl_di0=sl_fl0
    12141172                 !!!!!!!!!!!!!!!!!!!!!!!!!!
    1215                  CALL param_slope( mu0(ig), declin, sl_ra, sl_lat, 
     1173                 CALL param_slope( mu0(ig), declin, sl_ra, sl_lat,
    12161174     &                             sl_tau, sl_alb, sl_the, sl_psi,
    12171175     &                             sl_di0, sl_fl0, sl_flu )
    … …  
    12261184          ELSE        ! not calling subslope, nslope might be > 1
    12271185          DO islope = 1,nslope
    1228             sl_the=abs(def_slope_mean(islope)) 
     1186            sl_the=abs(def_slope_mean(islope))
    12291187            IF (sl_the .gt. 1e-6) THEN
    12301188              IF(def_slope_mean(islope).ge.0.) THEN
    … …  
    12331191                psi_sl(:) = 180. !Southward slope
    12341192              ENDIF
    1235               DO ig=1,ngrid 
     1193              DO ig=1,ngrid
    12361194                ztim1=fluxsurf_dn_sw(ig,1,islope)
    12371195     s                +fluxsurf_dn_sw(ig,2,islope)
    … …  
    12461204                 sl_di0 = 0.
    12471205                 if (mu0(ig) .gt. 0.) then
    1248                           sl_di0 = mu0(ig)*(exp(-sl_tau/mu0(ig)))
     1206                  sl_di0 = mu0(ig)*(exp(-sl_tau/mu0(ig)))
    12491207                  sl_di0 = sl_di0*flux_1AU/dist_sol/dist_sol
    12501208                  sl_di0 = sl_di0/ztim1
    … …  
    12541212                 if (sl_fl0 .lt. sl_di0) sl_di0=sl_fl0
    12551213                 !!!!!!!!!!!!!!!!!!!!!!!!!!
    1256                  CALL param_slope( mu0(ig), declin, sl_ra, sl_lat, 
     1214                 CALL param_slope( mu0(ig), declin, sl_ra, sl_lat,
    12571215     &                             sl_tau, sl_alb, sl_the, sl_psi,
    12581216     &                             sl_di0, sl_fl0, sl_flu )
    … …  
    12831241               fluxrad_sky(ig,islope) =
    12841242     $           emis(ig,islope)*fluxsurf_lw(ig,islope)
    1285      $         +fluxsurf_dn_sw(ig,1,islope)*(1.-albedo(ig,1,islope))         
     1243     $         +fluxsurf_dn_sw(ig,1,islope)*(1.-albedo(ig,1,islope))
    12861244     $         +fluxsurf_dn_sw(ig,2,islope)*(1.-albedo(ig,2,islope))
    12871245             ENDDO
    … …  
    13211279                 fluxrad(ig,nslope)=fluxrad(ig,nslope)+
    13221280     $                             (1.-sky)*zplanck(ig)
    1323                ELSE 
     1281               ELSE
    13241282                 fluxrad(ig,islope)=fluxrad(ig,islope) +
    13251283     $              (1.-sky_slope(iflat))*emis(ig,iflat)*
    1326      $              stephan*tsurf(ig,iflat)**4 
     1284     $              stephan*tsurf(ig,iflat)**4
    13271285               ENDIF
    13281286           ENDDO
    … …  
    13791337     &                      pdqrds,wspeed,dsodust,dsords,dsotop,
    13801338     &                      tau_pref_scenario,tau_pref_gcm)
    1381      
     1339
    13821340c      update the tendencies of both dust after vertical transport
    13831341         DO l=1,nlayer
    13841342          DO ig=1,ngrid
    13851343           pdq(ig,l,igcm_stormdust_mass)=
    1386      &              pdq(ig,l,igcm_stormdust_mass)+ 
     1344     &              pdq(ig,l,igcm_stormdust_mass)+
    13871345     &                      pdqrds(ig,l,igcm_stormdust_mass)
    13881346           pdq(ig,l,igcm_stormdust_number)=
    … …  
    13931351     &       pdq(ig,l,igcm_dust_mass)+ pdqrds(ig,l,igcm_dust_mass)
    13941352           pdq(ig,l,igcm_dust_number)=
    1395      &           pdq(ig,l,igcm_dust_number)+ 
     1353     &           pdq(ig,l,igcm_dust_number)+
    13961354     &                  pdqrds(ig,l,igcm_dust_number)
    13971355
    … …  
    14151373!         call write_output('qstormrds1','qstorm after rds',
    14161374!     &           'kg/kg ',qstormrds1(:,:))
    1417 !         
     1375!
    14181376!         call write_output('qdusttotal0','q sum before rds',
    14191377!     &           'kg/m2 ',qdusttotal0(:))
    … …  
    14331391     &                zzlay,zdtsw,zdtlw,
    14341392     &                icount,zday,zls,tsurf(:,iflat),
    1435      &                qsurf_meshavg(:,igcm_co2), 
     1393     &                qsurf_meshavg(:,igcm_co2),
    14361394     &                igout,aerosol,tauscaling,dust_rad_adjust,
    14371395     &                IRtoVIScoef,albedo_meshavg,emis_meshavg,
    … …  
    14411399     &                pdqtop,wtop,dsodust,dsords,dsotop,
    14421400     &                tau_pref_scenario,tau_pref_gcm)
    1443                    
     1401
    14441402c      update the tendencies of both dust after vertical transport
    14451403         DO l=1,nlayer
    … …  
    14661424
    14671425             CALL compute_dtau(ngrid,nlayer,
    1468      &                          zday,pplev,tau_pref_gcm,
    1469      &                          ptimestep,local_time,IRtoVIScoef,
     1426     &                  zday,pplev,tau_pref_gcm,
     1427     &                  ptimestep,local_time,IRtoVIScoef,
    14701428     &                          dustliftday)
    14711429           endif ! end of if (dustinjection.gt.0)
    … …  
    15121470c without using the thermals in gcm and mesoscale can yield problems in
    15131471c weakly unstable situations when winds are near to 0. For those cases, we add
    1514 c a unit subgrid gustiness. Remember that thermals should be used we using the 
     1472c a unit subgrid gustiness. Remember that thermals should be used we using the
    15151473c Richardson based surface layer model.
    1516         IF ( .not.calltherm 
    1517      .       .and. callrichsl 
     1474        IF ( .not.calltherm
     1475     .       .and. callrichsl
    15181476     .       .and. .not.turb_resolved) THEN
    15191477
    … …  
    15251483               wstar(ig)=0.
    15261484               hfmax_th(ig)=0.
    1527              ENDIF     
     1485             ENDIF
    15281486          ENDDO
    15291487        ENDIF
    … …  
    15401498c        Calling vdif (Martian version WITH CO2 condensation)
    15411499         dwatercap_dif(:,:) = 0.
    1542          zcdh(:) = 0.
    1543          zcdv(:) = 0.
    15441500
    15451501         CALL vdifc(ngrid,nlayer,nsoilmx,nq,nqsoil,zpopsk,
    … …  
    15491505     $        zdum1,zdum2,zdh,pdq,zflubid,
    15501506     $        zdudif,zdvdif,zdhdif,zdtsdif,q2,
    1551      &        zdqdif,zdqsdif,wstar,zcdv,zcdh,hfmax_th,
     1507     &        zdqdif,zdqsdif,wstar,hfmax_th,
    15521508     &        zcondicea_co2microp,sensibFlux,
    15531509     &        dustliftday,local_time,watercap,dwatercap_dif)
    … …  
    15921548           !! Specific treatment for lifting in turbulent-resolving mode (AC)
    15931549           IF (lifting .and. doubleq) THEN
    1594              !! lifted dust is injected in the first layer. 
    1595              !! Sedimentation must be called after turbulent mixing, i.e. on next step, after WRF. 
     1550             !! lifted dust is injected in the first layer.
     1551             !! Sedimentation must be called after turbulent mixing, i.e. on next step, after WRF.
    15961552             !! => lifted dust is not incremented before the sedimentation step.
    15971553             zdqdif(1:ngrid,1,1:nq)=0.
    1598              zdqdif(1:ngrid,1,igcm_dust_number) = 
     1554             zdqdif(1:ngrid,1,igcm_dust_number) =
    15991555     .        -zdqsdif_meshavg_tmp(1:ngrid,igcm_dust_number)
    1600              zdqdif(1:ngrid,1,igcm_dust_mass) = 
     1556             zdqdif(1:ngrid,1,igcm_dust_mass) =
    16011557     .        -zdqsdif_meshavg_tmp(1:ngrid,igcm_dust_mass)
    16021558             zdqdif(1:ngrid,2:nlayer,1:nq) = 0.
    … …  
    16121568           ENDIF
    16131569         ENDIF
    1614       ELSE   
     1570      ELSE
    16151571         DO ig=1,ngrid
    16161572           DO islope=1,nslope
    … …  
    16211577
    16221578         IF (turb_resolved) THEN
    1623             write(*,*) 'Turbulent-resolving mode !' 
     1579            write(*,*) 'Turbulent-resolving mode !'
    16241580            write(*,*) 'Please set calldifv to T in callphys.def'
    16251581            call abort_physic("physiq","turbulent-resolving mode",1)
    … …  
    16391595     $ pdu_th,pdv_th,pdt_th,pdq_th,lmax_th,zmax_th,
    16401596     $ dtke_th,zdhdif,hfmax_th,wstar,sensibFlux)
    1641        
     1597
    16421598         DO l=1,nlayer
    16431599           DO ig=1,ngrid
    … …  
    17051661            DO l=1,nlayer
    17061662              DO ig=1,ngrid
    1707                  pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqadj(ig,l,iq) 
     1663                 pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqadj(ig,l,iq)
    17081664              ENDDO
    17091665            ENDDO
    … …  
    17271683         IF (calljliu_gwimix) THEN
    17281684          CALL nonoro_gwd_mix(ngrid,nlayer,ptimestep,
    1729      &               nq,cpnew, rnew, 
    1730      &               zplay, 
    1731      &               zmax_th, 
    1732      &               pt, pu, pv, pq, 
     1685     &               nq,cpnew, rnew,
     1686     &               zplay,
     1687     &               zmax_th,
     1688     &               pt, pu, pv, pq,
    17331689                !loss,  chemical reaction loss rates
    17341690     &               pdt, pdu, pdv, pdq,
    … …  
    17411697     &                         +d_t_hin(1:ngrid,1:nlayer)
    17421698         pdu(1:ngrid,1:nlayer)=pdu(1:ngrid,1:nlayer)
    1743      &                         +d_u_hin(1:ngrid,1:nlayer) 
     1699     &                         +d_u_hin(1:ngrid,1:nlayer)
    17441700         pdv(1:ngrid,1:nlayer)=pdv(1:ngrid,1:nlayer)
    17451701     &                         +d_v_hin(1:ngrid,1:nlayer)
    … …  
    17491705     &                         +d_t_mix(1:ngrid,1:nlayer)
    17501706         pdu(1:ngrid,1:nlayer)=pdu(1:ngrid,1:nlayer)
    1751      &                         +d_u_mix(1:ngrid,1:nlayer) 
     1707     &                         +d_u_mix(1:ngrid,1:nlayer)
    17521708         pdv(1:ngrid,1:nlayer)=pdv(1:ngrid,1:nlayer)
    17531709     &                         +d_v_mix(1:ngrid,1:nlayer)
    … …  
    17601716
    17611717c-----------------------------------------------------------------------
    1762 c   9. Specific parameterizations for tracers 
     1718c   9. Specific parameterizations for tracers
    17631719c:   -----------------------------------------
    17641720
    … …  
    17851741
    17861742! increment water vapour and ice atmospheric tracers tendencies
    1787            pdq(1:ngrid,1:nlayer,igcm_h2o_vap) = 
    1788      &       pdq(1:ngrid,1:nlayer,igcm_h2o_vap) + 
     1743           pdq(1:ngrid,1:nlayer,igcm_h2o_vap) =
     1744     &       pdq(1:ngrid,1:nlayer,igcm_h2o_vap) +
    17891745     &       zdqcloud(1:ngrid,1:nlayer,igcm_h2o_vap)
    1790            pdq(1:ngrid,1:nlayer,igcm_h2o_ice) = 
    1791      &       pdq(1:ngrid,1:nlayer,igcm_h2o_ice) + 
     1746           pdq(1:ngrid,1:nlayer,igcm_h2o_ice) =
     1747     &       pdq(1:ngrid,1:nlayer,igcm_h2o_ice) +
    17921748     &       zdqcloud(1:ngrid,1:nlayer,igcm_h2o_ice)
    17931749
    … …  
    18061762! This is due to single precision rounding problems
    18071763           if (microphys) then
    1808              pdq(1:ngrid,1:nlayer,igcm_ccn_mass) = 
    1809      &         pdq(1:ngrid,1:nlayer,igcm_ccn_mass) + 
     1764             pdq(1:ngrid,1:nlayer,igcm_ccn_mass) =
     1765     &         pdq(1:ngrid,1:nlayer,igcm_ccn_mass) +
    18101766     &         zdqcloud(1:ngrid,1:nlayer,igcm_ccn_mass)
    1811              pdq(1:ngrid,1:nlayer,igcm_ccn_number) = 
    1812      &         pdq(1:ngrid,1:nlayer,igcm_ccn_number) + 
     1767             pdq(1:ngrid,1:nlayer,igcm_ccn_number) =
     1768     &         pdq(1:ngrid,1:nlayer,igcm_ccn_number) +
    18131769     &         zdqcloud(1:ngrid,1:nlayer,igcm_ccn_number)
    1814              where (pq(:,:,igcm_ccn_mass) + 
     1770             where (pq(:,:,igcm_ccn_mass) +
    18151771     &       ptimestep*pdq(:,:,igcm_ccn_mass) < 0.)
    1816                pdq(:,:,igcm_ccn_mass) = 
     1772               pdq(:,:,igcm_ccn_mass) =
    18171773     &         - pq(:,:,igcm_ccn_mass)/ptimestep + 1.e-30
    1818                pdq(:,:,igcm_ccn_number) = 
     1774               pdq(:,:,igcm_ccn_number) =
    18191775     &         - pq(:,:,igcm_ccn_number)/ptimestep + 1.e-30
    18201776             end where
    1821              where (pq(:,:,igcm_ccn_number) + 
     1777             where (pq(:,:,igcm_ccn_number) +
    18221778     &       ptimestep*pdq(:,:,igcm_ccn_number) < 0.)
    1823                pdq(:,:,igcm_ccn_mass) = 
     1779               pdq(:,:,igcm_ccn_mass) =
    18241780     &         - pq(:,:,igcm_ccn_mass)/ptimestep + 1.e-30
    1825                pdq(:,:,igcm_ccn_number) = 
     1781               pdq(:,:,igcm_ccn_number) =
    18261782     &         - pq(:,:,igcm_ccn_number)/ptimestep + 1.e-30
    18271783             end where
    … …  
    18291785
    18301786           if (scavenging) then
    1831              pdq(1:ngrid,1:nlayer,igcm_dust_mass) = 
    1832      &         pdq(1:ngrid,1:nlayer,igcm_dust_mass) + 
     1787             pdq(1:ngrid,1:nlayer,igcm_dust_mass) =
     1788     &         pdq(1:ngrid,1:nlayer,igcm_dust_mass) +
    18331789     &         zdqcloud(1:ngrid,1:nlayer,igcm_dust_mass)
    1834              pdq(1:ngrid,1:nlayer,igcm_dust_number) = 
    1835      &         pdq(1:ngrid,1:nlayer,igcm_dust_number) + 
     1790             pdq(1:ngrid,1:nlayer,igcm_dust_number) =
     1791     &         pdq(1:ngrid,1:nlayer,igcm_dust_number) +
    18361792     &         zdqcloud(1:ngrid,1:nlayer,igcm_dust_number)
    1837              where (pq(:,:,igcm_dust_mass) + 
     1793             where (pq(:,:,igcm_dust_mass) +
    18381794     &       ptimestep*pdq(:,:,igcm_dust_mass) < 0.)
    1839                pdq(:,:,igcm_dust_mass) = 
     1795               pdq(:,:,igcm_dust_mass) =
    18401796     &         - pq(:,:,igcm_dust_mass)/ptimestep + 1.e-30
    1841                pdq(:,:,igcm_dust_number) = 
     1797               pdq(:,:,igcm_dust_number) =
    18421798     &         - pq(:,:,igcm_dust_number)/ptimestep + 1.e-30
    18431799             end where
    1844              where (pq(:,:,igcm_dust_number) + 
     1800             where (pq(:,:,igcm_dust_number) +
    18451801     &       ptimestep*pdq(:,:,igcm_dust_number) < 0.)
    1846                pdq(:,:,igcm_dust_mass) = 
     1802               pdq(:,:,igcm_dust_mass) =
    18471803     &         - pq(:,:,igcm_dust_mass)/ptimestep + 1.e-30
    1848                pdq(:,:,igcm_dust_number) = 
     1804               pdq(:,:,igcm_dust_number) =
    18491805     &         - pq(:,:,igcm_dust_number)/ptimestep + 1.e-30
    18501806             end where
    18511807           endif ! of if scavenging
    1852                      
     1808
    18531809         END IF  ! of IF (water)
    18541810
    … …  
    18591815c flag needed in callphys.def:
    18601816c               co2clouds=.true. is mandatory (default is .false.)
    1861 c               co2useh2o=.true. if you want to allow co2 condensation 
    1862 c                                on water ice particles 
     1817c               co2useh2o=.true. if you want to allow co2 condensation
     1818c                                on water ice particles
    18631819c               meteo_flux=.true. if you want to add a meteoritic
    18641820c                                 supply of CCN
    18651821c               CLFvaryingCO2=.true. if you want to have a sub-grid
    1866 c                                    temperature distribution 
     1822c                                    temperature distribution
    18671823c               spantCO2=integer (i.e. 3) amplitude of the sub-grid T disti
    1868 c               nuiceco2_sed=0.2 variance of the size distribution for the 
     1824c               nuiceco2_sed=0.2 variance of the size distribution for the
    18691825c                                sedimentation
    1870 c               nuiceco2_ref=0.2 variance of the size distribution for the 
     1826c               nuiceco2_ref=0.2 variance of the size distribution for the
    18711827c                                nucleation
    18721828c               imicroco2=50     micro-timestep is 1/50 of physical timestep
    … …  
    18921848c Temperature variation due to latent heat release
    18931849               pdt(1:ngrid,1:nlayer) =
    1894      &              pdt(1:ngrid,1:nlayer) + 
     1850     &              pdt(1:ngrid,1:nlayer) +
    18951851     &              zdtcloudco2(1:ngrid,1:nlayer)
    18961852
    … …  
    19281884!Update water ice clouds values as well
    19291885             if (co2useh2o) then
    1930                 pdq(1:ngrid,1:nlayer,igcm_h2o_ice) = 
    1931      &               pdq(1:ngrid,1:nlayer,igcm_h2o_ice) + 
     1886                pdq(1:ngrid,1:nlayer,igcm_h2o_ice) =
     1887     &               pdq(1:ngrid,1:nlayer,igcm_h2o_ice) +
    19321888     &               zdqcloudco2(1:ngrid,1:nlayer,igcm_h2o_ice)
    1933                 pdq(1:ngrid,1:nlayer,igcm_ccn_mass) = 
    1934      &               pdq(1:ngrid,1:nlayer,igcm_ccn_mass) + 
     1889                pdq(1:ngrid,1:nlayer,igcm_ccn_mass) =
     1890     &               pdq(1:ngrid,1:nlayer,igcm_ccn_mass) +
    19351891     &               zdqcloudco2(1:ngrid,1:nlayer,igcm_ccn_mass)
    1936                 pdq(1:ngrid,1:nlayer,igcm_ccn_number) = 
    1937      &               pdq(1:ngrid,1:nlayer,igcm_ccn_number) + 
     1892                pdq(1:ngrid,1:nlayer,igcm_ccn_number) =
     1893     &               pdq(1:ngrid,1:nlayer,igcm_ccn_number) +
    19381894     &               zdqcloudco2(1:ngrid,1:nlayer,igcm_ccn_number)
    19391895
    … …  
    19531909                where (pq(:,:,igcm_ccn_mass) +
    19541910     &                 ptimestep*pdq(:,:,igcm_ccn_mass) < 0.)
    1955                   pdq(:,:,igcm_ccn_mass) = 
     1911                  pdq(:,:,igcm_ccn_mass) =
    19561912     &               - pq(:,:,igcm_ccn_mass)/ptimestep + 1.e-30
    1957                   pdq(:,:,igcm_ccn_number) = 
     1913                  pdq(:,:,igcm_ccn_number) =
    19581914     &               - pq(:,:,igcm_ccn_number)/ptimestep + 1.e-30
    19591915                end where
    … …  
    19611917                where (pq(:,:,igcm_ccn_number) +
    19621918     &                 ptimestep*pdq(:,:,igcm_ccn_number) < 0.)
    1963                   pdq(:,:,igcm_ccn_mass) = 
     1919                  pdq(:,:,igcm_ccn_mass) =
    19641920     &              - pq(:,:,igcm_ccn_mass)/ptimestep + 1.e-30
    1965                   pdq(:,:,igcm_ccn_number) = 
     1921                  pdq(:,:,igcm_ccn_number) =
    19661922     &              - pq(:,:,igcm_ccn_number)/ptimestep + 1.e-30
    19671923                end where
    … …  
    19971953             endif ! of if (co2useh2o)
    19981954c     Negative values?
    1999              where (pq(:,:,igcm_ccnco2_mass) + 
     1955             where (pq(:,:,igcm_ccnco2_mass) +
    20001956     &              ptimestep*pdq(:,:,igcm_ccnco2_mass) < 0.)
    2001                pdq(:,:,igcm_ccnco2_mass) = 
     1957               pdq(:,:,igcm_ccnco2_mass) =
    20021958     &            - pq(:,:,igcm_ccnco2_mass)/ptimestep + 1.e-30
    2003                pdq(:,:,igcm_ccnco2_number) = 
     1959               pdq(:,:,igcm_ccnco2_number) =
    20041960     &            - pq(:,:,igcm_ccnco2_number)/ptimestep + 1.e-30
    20051961             end where
    2006              where (pq(:,:,igcm_ccnco2_number) + 
     1962             where (pq(:,:,igcm_ccnco2_number) +
    20071963     &              ptimestep*pdq(:,:,igcm_ccnco2_number) < 0.)
    2008                pdq(:,:,igcm_ccnco2_mass) = 
     1964               pdq(:,:,igcm_ccnco2_mass) =
    20091965     &            - pq(:,:,igcm_ccnco2_mass)/ptimestep + 1.e-30
    2010                pdq(:,:,igcm_ccnco2_number) = 
     1966               pdq(:,:,igcm_ccnco2_number) =
    20111967     &            - pq(:,:,igcm_ccnco2_number)/ptimestep + 1.e-30
    20121968             end where
    2013        
     1969
    20141970c     Negative values?
    2015              where (pq(:,:,igcm_dust_mass) + 
     1971             where (pq(:,:,igcm_dust_mass) +
    20161972     &              ptimestep*pdq(:,:,igcm_dust_mass) < 0.)
    2017                pdq(:,:,igcm_dust_mass) = 
     1973               pdq(:,:,igcm_dust_mass) =
    20181974     &           - pq(:,:,igcm_dust_mass)/ptimestep + 1.e-30
    2019                pdq(:,:,igcm_dust_number) = 
     1975               pdq(:,:,igcm_dust_number) =
    20201976     &           - pq(:,:,igcm_dust_number)/ptimestep + 1.e-30
    20211977             end where
    2022              where (pq(:,:,igcm_dust_number) + 
     1978             where (pq(:,:,igcm_dust_number) +
    20231979     &              ptimestep*pdq(:,:,igcm_dust_number) < 0.)
    2024                pdq(:,:,igcm_dust_mass) = 
     1980               pdq(:,:,igcm_dust_mass) =
    20251981     &           - pq(:,:,igcm_dust_mass)/ptimestep + 1.e-30
    2026                pdq(:,:,igcm_dust_number) = 
     1982               pdq(:,:,igcm_dust_number) =
    20271983     &           - pq(:,:,igcm_dust_number)/ptimestep + 1.e-30
    20281984             end where
    … …  
    20502006c        Dust devil :
    20512007c        ----------
    2052          IF(callddevil) then 
     2008         IF(callddevil) then
    20532009           call dustdevil(ngrid,nlayer,nq, zplev,pu,pv,pt, tsurf,q2,
    20542010     &                zdqdev,zdqsdev)
    … …  
    20742030         END IF ! of IF (callddevil)
    20752031
    2076 c        ------------- 
     2032c        -------------
    20772033c        Sedimentation :   acts also on water ice
    20782034c        -------------
    2079          IF (sedimentation) THEN 
     2035         IF (sedimentation) THEN
    20802036           zdqsed(1:ngrid,1:nlayer,1:nq)=0
    20812037           zdqssed(1:ngrid,1:nq)=0
    … …  
    20882044     &                rdust,rstormdust,rtopdust,
    20892045     &                rice,rsedcloud,rhocloud,
    2090      &                pq,pdq,zdqsed,zdqssed,nq, 
     2046     &                pq,pdq,zdqsed,zdqssed,nq,
    20912047     &                tau,tauscaling)
    20922048
    … …  
    20952051           IF (rdstorm) THEN
    20962052c             Storm dust cannot sediment to the surface
    2097               DO ig=1,ngrid 
     2053              DO ig=1,ngrid
    20982054                 zdqsed(ig,1,igcm_stormdust_mass)=
    20992055     &             zdqsed(ig,1,igcm_stormdust_mass)+
    … …  
    21032059     &             zdqsed(ig,1,igcm_stormdust_number)+
    21042060     &             zdqssed(ig,igcm_stormdust_number) /
    2105      &               ((pplev(ig,1)-pplev(ig,2))/g) 
     2061     &               ((pplev(ig,1)-pplev(ig,2))/g)
    21062062                 zdqssed(ig,igcm_stormdust_mass)=0.
    21072063                 zdqssed(ig,igcm_stormdust_number)=0.
    … …  
    21572113!           dust and ice surface area
    21582114            call surfacearea(ngrid, nlayer, naerkind,
    2159      $                       ptimestep, zplay, zzlay, 
    2160      $                       pt, pq, pdq, nq, 
    2161      $                       rdust, rice, tau, tauscaling, 
     2115     $                       ptimestep, zplay, zzlay,
     2116     $                       pt, pq, pdq, nq,
     2117     $                       rdust, rice, tau, tauscaling,
    21622118     $                       surfdust, surfice)
    21632119!           call photochemistry
    … …  
    21832139             ENDDO
    21842140           ENDDO ! of DO iq=1,nq
    2185            
     2141
    21862142           ! add condensation tendency for H2O2
    21872143           if (igcm_h2o2.ne.0) then
    … …  
    22572213         qsurf_tmp(:,igcm_co2) = qsurf_tmp(:,igcm_co2) * 10000.
    22582214      ENDIF
    2259      
     2215
    22602216
    22612217      IF (callcond) THEN
    … …  
    23132269        ENDDO
    23142270        zplev(:,nlayer+1) = 0.
    2315        
    2316         ! Calculation of zzlay and zzlay with udpated pressure and temperature
    2317         DO ig=1,ngrid
    2318         zzlay(ig,1)=-(log(zplay(ig,1)/ps(ig)))*rnew(ig,1)*
    2319      &          (pt(ig,1)+pdt(ig,1)*ptimestep) /g                                     
    2320 
    2321           DO l=2,nlayer
    2322                  
     2271
     2272    ! Calculation of zzlay and zzlay with udpated pressure and temperature
     2273        DO ig=1,ngrid
     2274        zzlay(ig,1)=-(log(zplay(ig,1)/ps(ig)))*rnew(ig,1)*
     2275     &      (pt(ig,1)+pdt(ig,1)*ptimestep) /g
     2276
     2277          DO l=2,nlayer
     2278
    23232279           ! compute "mean" temperature of the layer
    23242280            if((pt(ig,l)+pdt(ig,l)*ptimestep) .eq.
    … …  
    23262282               tlaymean= pt(ig,l)+pdt(ig,l)*ptimestep
    23272283            else
    2328                tlaymean=((pt(ig,l)+pdt(ig,l)*ptimestep)- 
    2329      &            (pt(ig,l-1)+pdt(ig,l-1)*ptimestep))/
    2330      &            log((pt(ig,l)+pdt(ig,l)*ptimestep)/
    2331      &            (pt(ig,l-1)+pdt(ig,l-1)*ptimestep))
     2284               tlaymean=((pt(ig,l)+pdt(ig,l)*ptimestep)-
     2285     &        (pt(ig,l-1)+pdt(ig,l-1)*ptimestep))/
     2286     &        log((pt(ig,l)+pdt(ig,l)*ptimestep)/
     2287     &        (pt(ig,l-1)+pdt(ig,l-1)*ptimestep))
    23322288            endif
    2333                      
     2289
    23342290            ! compute gravitational acceleration (at altitude zaeroid(nlayer-1))
    23352291            gz(ig,l)=g*(rad**2)/(rad+zzlay(ig,l-1)+(phisfi(ig)/g))**2
    2336          
    2337 
    2338             zzlay(ig,l)=zzlay(ig,l-1)- 
     2292
     2293
     2294            zzlay(ig,l)=zzlay(ig,l-1)-
    23392295     &     (log(zplay(ig,l)/zplay(ig,l-1))*rnew(ig,l)*tlaymean/gz(ig,l))
    2340              
    2341        
     2296
     2297
    23422298        !   update layers altitude
    23432299             z1=(zplay(ig,l-1)+zplev(ig,l))/(zplay(ig,l-1)-zplev(ig,l))
    … …  
    23452301             zzlev(ig,l)=(z1*zzlay(ig,l-1)+z2*zzlay(ig,l))/(z1+z2)
    23462302          ENDDO
    2347         ENDDO 
     2303        ENDDO
    23482304#endif
    23492305      ENDIF  ! of IF (callcond)
    … …  
    23512307c-----------------------------------------------------------------------
    23522308c  Updating tracer budget on surface
    2353 c-----------------------------------------------------------------------         
     2309c-----------------------------------------------------------------------
    23542310        DO iq=1, nq
    23552311          DO ig=1,ngrid
    … …  
    23572313            qsurf(ig,iq,islope)=qsurf(ig,iq,islope)+
    23582314     &                ptimestep*dqsurf(ig,iq,islope)
    2359            ENDDO 
     2315           ENDDO
    23602316          ENDDO  ! (ig)
    23612317        ENDDO    ! (iq)
    … …  
    23712327        DO islope = 1,nslope
    23722328            tsurf(ig,islope)=tsurf(ig,islope)+
    2373      &            ptimestep*zdtsurf(ig,islope) 
     2329     &            ptimestep*zdtsurf(ig,islope)
    23742330       ENDDO
    23752331      ENDDO
    … …  
    23812337
    23822338      IF (water) THEN
    2383 !#ifndef MESOSCALE 
     2339!#ifndef MESOSCALE
    23842340!         if (caps.and.(obliquit.lt.27.)) then => now done in co2condens
    23852341           ! NB: Updated surface pressure, at grid point 'ngrid', is
    … …  
    23932349c       Change of surface albedo in case of ground frost
    23942350c       everywhere except on the north permanent cap and in regions
    2395 c       covered by dry ice. 
     2351c       covered by dry ice.
    23962352c              ALWAYS PLACE these lines after co2condens !!!
    23972353c       -------------------------------------------------------------
    … …  
    24372393c     To avoid negative values
    24382394      IF (rdstorm) THEN
    2439            where (pq(:,:,igcm_stormdust_mass) + 
     2395           where (pq(:,:,igcm_stormdust_mass) +
    24402396     &      ptimestep*pdq(:,:,igcm_stormdust_mass) < 0.)
    2441              pdq(:,:,igcm_stormdust_mass) = 
     2397             pdq(:,:,igcm_stormdust_mass) =
    24422398     &       - pq(:,:,igcm_stormdust_mass)/ptimestep + 1.e-30
    2443              pdq(:,:,igcm_stormdust_number) = 
     2399             pdq(:,:,igcm_stormdust_number) =
    24442400     &       - pq(:,:,igcm_stormdust_number)/ptimestep + 1.e-30
    24452401           end where
    2446            where (pq(:,:,igcm_stormdust_number) + 
     2402           where (pq(:,:,igcm_stormdust_number) +
    24472403     &      ptimestep*pdq(:,:,igcm_stormdust_number) < 0.)
    2448              pdq(:,:,igcm_stormdust_mass) = 
     2404             pdq(:,:,igcm_stormdust_mass) =
    24492405     &       - pq(:,:,igcm_stormdust_mass)/ptimestep + 1.e-30
    2450              pdq(:,:,igcm_stormdust_number) = 
     2406             pdq(:,:,igcm_stormdust_number) =
    24512407     &       - pq(:,:,igcm_dust_number)/ptimestep + 1.e-30
    24522408           end where
    24532409
    2454             where (pq(:,:,igcm_dust_mass) + 
     2410            where (pq(:,:,igcm_dust_mass) +
    24552411     &      ptimestep*pdq(:,:,igcm_dust_mass) < 0.)
    2456              pdq(:,:,igcm_dust_mass) = 
     2412             pdq(:,:,igcm_dust_mass) =
    24572413     &       - pq(:,:,igcm_dust_mass)/ptimestep + 1.e-30
    2458              pdq(:,:,igcm_dust_number) = 
     2414             pdq(:,:,igcm_dust_number) =
    24592415     &       - pq(:,:,igcm_dust_number)/ptimestep + 1.e-30
    24602416           end where
    2461            where (pq(:,:,igcm_dust_number) + 
     2417           where (pq(:,:,igcm_dust_number) +
    24622418     &      ptimestep*pdq(:,:,igcm_dust_number) < 0.)
    2463              pdq(:,:,igcm_dust_mass) = 
     2419             pdq(:,:,igcm_dust_mass) =
    24642420     &       - pq(:,:,igcm_dust_mass)/ptimestep + 1.e-30
    2465              pdq(:,:,igcm_dust_number) = 
     2421             pdq(:,:,igcm_dust_number) =
    24662422     &       - pq(:,:,igcm_dust_number)/ptimestep + 1.e-30
    24672423           end where
    2468       ENDIF !(rdstorm)   
     2424      ENDIF !(rdstorm)
    24692425
    24702426c-----------------------------------------------------------------------
    … …  
    25032459c  13. Write output files
    25042460c  ----------------------
    2505 
    25062461      call compute_meshgridavg(ngrid,nq,albedo,emis,tsurf,qsurf,
    25072462     &       albedo_meshavg,emis_meshavg,tsurf_meshavg,qsurf_meshavg)
    … …  
    25682523               ztim1 = pt(ig,l)
    25692524               igmin = ig
    2570                lmin = l 
     2525               lmin = l
    25712526           end if
    25722527        ENDDO
    … …  
    25802535
    25812536c     ---------------------
    2582 c     Outputs to the screen 
     2537c     Outputs to the screen
    25832538c     ---------------------
    25842539
    … …  
    26622617
    26632618         ! default: not writing a restart file at this time step
    2664          write_restart=.false. 
     2619         write_restart=.false.
    26652620         IF (ecritstart.GT.0) THEN
    26662621           ! For when we store multiple time steps in the restart file
    … …  
    26732628           write_restart=.true.
    26742629         ENDIF
    2675          
     2630
    26762631          IF (write_restart) THEN
    26772632           IF (grid_type==unstructured) THEN !IF DYNAMICO
    26782633
    26792634             ! When running Dynamico, no need to add a dynamics time step to ztime_fin
    2680              IF (ptime.LE. 1.E-10) THEN 
     2635             IF (ptime.LE. 1.E-10) THEN
    26812636             ! Residual ptime occurs with Dynamico
    26822637             ztime_fin = pday !+ ptime + ptimestep/(float(iphysiq)*daysec)
    … …  
    26932648
    26942649            if (ecritstart.GT.0) then !IF MULTIPLE RESTARTS nothing change
    2695               ztime_fin = pday - day_ini + ptime 
     2650              ztime_fin = pday - day_ini + ptime
    26962651     &                    + ptimestep/(float(iphysiq)*daysec)
    26972652            else !IF ONE RESTART final time in top of day_end
    2698               ztime_fin = pday - day_ini-(day_end-day_ini)
     2653          ztime_fin = pday - day_ini-(day_end-day_ini)
    26992654     &                    + ptime  + ptimestep/(float(iphysiq)*daysec)
    27002655            endif
    27012656
    27022657           ENDIF ! of IF (grid_type==unstructured)
    2703           write(*,'(A,I7,A,F12.5)') 
     2658          write(*,'(A,I7,A,F12.5)')
    27042659     .         'PHYSIQ: writing in restartfi ; icount=',
    27052660     .          icount,' date=',ztime_fin
    2706            
     2661
    27072662          call physdem1("restartfi.nc",nsoilmx,ngrid,nlayer,nq,nqsoil,
    27082663     .                ptimestep,ztime_fin,
    … …  
    27202675c          diagfi files
    27212676c        -------------------------------------------------------------------
    2722 
    27232677         do ig=1,ngrid
    2724            if(mu0(ig).le.0.01) then
    2725             fluxsurf_dir_dn_sw(ig) = 0.
    2726            else
     2678       if(mu0(ig).le.0.01) then
     2679        fluxsurf_dir_dn_sw(ig) = 0.
     2680       else
    27272681            if (water) then
    27282682             ! both water and dust contribute
    … …  
    27492703
    27502704           if(doubleq) then
    2751               do ig=1,ngrid 
    2752          IF (sedimentation) THEN 
    2753                 dqdustsurf(ig) = 
     2705              do ig=1,ngrid
     2706         IF (sedimentation) THEN
     2707                dqdustsurf(ig) =
    27542708     &                zdqssed(ig,igcm_dust_mass)*tauscaling(ig)
    2755                 dndustsurf(ig) = 
     2709                dndustsurf(ig) =
    27562710     &                zdqssed(ig,igcm_dust_number)*tauscaling(ig)
    27572711         ENDIF
    … …  
    27622716              enddo
    27632717              if (scavenging) then
    2764                 do ig=1,ngrid 
    2765          IF (sedimentation) THEN 
    2766                   dqdustsurf(ig) = dqdustsurf(ig) + 
     2718                do ig=1,ngrid
     2719         IF (sedimentation) THEN
     2720                  dqdustsurf(ig) = dqdustsurf(ig) +
    27672721     &                     zdqssed(ig,igcm_ccn_mass)*tauscaling(ig)
    2768                   dndustsurf(ig) = dndustsurf(ig) + 
     2722                  dndustsurf(ig) = dndustsurf(ig) +
    27692723     &                     zdqssed(ig,igcm_ccn_number)*tauscaling(ig)
    27702724         ENDIF
    … …  
    27812735              mdusttot(:)=0
    27822736              qdusttotal(:,:)=0
    2783               do ig=1,ngrid 
    2784                 rdsdqdustsurf(ig) = 
     2737              do ig=1,ngrid
     2738                rdsdqdustsurf(ig) =
    27852739     &                zdqssed(ig,igcm_stormdust_mass)*tauscaling(ig)
    2786                 rdsdndustsurf(ig) = 
     2740                rdsdndustsurf(ig) =
    27872741     &                zdqssed(ig,igcm_stormdust_number)*tauscaling(ig)
    27882742                rdsndust(ig,:) =
    … …  
    27912745     &                pq(ig,:,igcm_stormdust_mass)*tauscaling(ig)
    27922746                do l=1,nlayer
    2793                     mstormdtot(ig) = mstormdtot(ig) + 
    2794      &                      zq(ig,l,igcm_stormdust_mass) * 
     2747                    mstormdtot(ig) = mstormdtot(ig) +
     2748     &                      zq(ig,l,igcm_stormdust_mass) *
    27952749     &                      (zplev(ig,l) - zplev(ig,l+1)) / g
    2796                     mdusttot(ig) = mdusttot(ig) + 
    2797      &                      zq(ig,l,igcm_dust_mass) * 
     2750                    mdusttot(ig) = mdusttot(ig) +
     2751     &                      zq(ig,l,igcm_dust_mass) *
    27982752     &                      (zplev(ig,l) - zplev(ig,l+1)) / g
    27992753                    qdusttotal(ig,l) = qdust(ig,l)+rdsqdust(ig,l) !calculate total dust
    … …  
    28012755              enddo
    28022756           endif !(rdstorm)
    2803                              
     2757
    28042758           if (water) then
    28052759             mtot(:)=0
    … …  
    28152769             do ig=1,ngrid
    28162770               do l=1,nlayer
    2817                  mtot(ig) = mtot(ig) + 
    2818      &                      zq(ig,l,igcm_h2o_vap) * 
     2771                 mtot(ig) = mtot(ig) +
     2772     &                      zq(ig,l,igcm_h2o_vap) *
    28192773     &                      (zplev(ig,l) - zplev(ig,l+1)) / g
    2820                  icetot(ig) = icetot(ig) + 
    2821      &                        zq(ig,l,igcm_h2o_ice) * 
     2774                 icetot(ig) = icetot(ig) +
     2775     &                        zq(ig,l,igcm_h2o_ice) *
    28222776     &                        (zplev(ig,l) - zplev(ig,l+1)) / g
    28232777                 IF (hdo) then
    … …  
    28352789     &             max(0.4e6*rice(ig,l)*(1.+nuice_ref)-0.05 ,0.),1.2
    28362790     &                        )
    2837                  opTES(ig,l)= 0.75 * Qabsice * 
     2791                 opTES(ig,l)= 0.75 * Qabsice *
    28382792     &             zq(ig,l,igcm_h2o_ice) *
    28392793     &             (zplev(ig,l) - zplev(ig,l+1)) / g
    28402794     &             / (rho_ice * rice(ig,l) * (1.+nuice_ref))
    2841                  tauTES(ig)=tauTES(ig)+ opTES(ig,l) 
     2795                 tauTES(ig)=tauTES(ig)+ opTES(ig,l)
    28422796               enddo
    28432797c              rave(ig)=rave(ig)/max(icetot(ig),1.e-30)       ! mass weight
    … …  
    28512805               Mccntot(:)= 0
    28522806               rave(:)=0
    2853                do ig=1,ngrid 
     2807               do ig=1,ngrid
    28542808                 do l=1,nlayer
    2855                     Nccntot(ig) = Nccntot(ig) + 
     2809                    Nccntot(ig) = Nccntot(ig) +
    28562810     &              zq(ig,l,igcm_ccn_number)*tauscaling(ig)
    28572811     &              *(zplev(ig,l) - zplev(ig,l+1)) / g
    2858                     Mccntot(ig) = Mccntot(ig) + 
     2812                    Mccntot(ig) = Mccntot(ig) +
    28592813     &              zq(ig,l,igcm_ccn_mass)*tauscaling(ig)
    28602814     &              *(zplev(ig,l) - zplev(ig,l+1)) / g
    2861 cccc Column integrated effective ice radius 
    2862 cccc is weighted by total ice surface area (BETTER than total ice mass) 
    2863                     rave(ig) = rave(ig) + 
     2815cccc Column integrated effective ice radius
     2816cccc is weighted by total ice surface area (BETTER than total ice mass)
     2817                    rave(ig) = rave(ig) +
    28642818     &                      tauscaling(ig) *
    28652819     &                      zq(ig,l,igcm_ccn_number) *
    2866      &                      (zplev(ig,l) - zplev(ig,l+1)) / g * 
     2820     &                      (zplev(ig,l) - zplev(ig,l+1)) / g *
    28672821     &                      rice(ig,l) * rice(ig,l)*  (1.+nuice_ref)
    28682822                 enddo
    … …  
    28732827             else ! of if (scavenging)
    28742828               rave(:)=0
    2875                do ig=1,ngrid 
     2829               do ig=1,ngrid
    28762830                 do l=1,nlayer
    2877                  rave(ig) = rave(ig) + 
     2831                 rave(ig) = rave(ig) +
    28782832     &                      zq(ig,l,igcm_h2o_ice) *
    2879      &                      (zplev(ig,l) - zplev(ig,l+1)) / g * 
     2833     &                      (zplev(ig,l) - zplev(ig,l+1)) / g *
    28802834     &                      rice(ig,l) * (1.+nuice_ref)
    28812835                 enddo
    2882                  rave(ig) = max(rave(ig) / 
     2836                 rave(ig) = max(rave(ig) /
    28832837     &             max(icetot(ig),1.e-30),1.e-30) ! mass weight
    28842838               enddo
    … …  
    29042858            do ig=1,ngrid
    29052859              do l=1,nlayer
    2906                 vaptotco2(ig) = vaptotco2(ig) + 
    2907      &                          zq(ig,l,igcm_co2) * 
     2860                vaptotco2(ig) = vaptotco2(ig) +
     2861     &                          zq(ig,l,igcm_co2) *
    29082862     &                          (zplev(ig,l) - zplev(ig,l+1)) / g
    2909                 icetotco2(ig) = icetot(ig) + 
    2910      &                          zq(ig,l,igcm_co2_ice) * 
     2863                icetotco2(ig) = icetot(ig) +
     2864     &                          zq(ig,l,igcm_co2_ice) *
    29112865     &                          (zplev(ig,l) - zplev(ig,l+1)) / g
    29122866              end do
    … …  
    29222876c        which can later be used to make the statistic files of the run:
    29232877c        if flag "callstats" from callphys.def is .true.)
    2924          
     2878
    29252879        call wstats(ngrid,"ps","Surface pressure","Pa",2,ps)
    29262880        call wstats(ngrid,"tsurf","Surface temperature","K",2
    … …  
    29692923          call wstats(ngrid,"fluxsurf_dir_dn_sw",
    29702924     &                "Direct incoming SW flux at surface",
    2971      &                "W.m-2",2,fluxsurf_dir_dn_sw)     
     2925     &                "W.m-2",2,fluxsurf_dir_dn_sw)
    29722926
    29732927          if (calltherm) then
    … …  
    30513005     &                    2,icetotco2)
    30523006             end if
    3053              
    3054              
     3007
     3008
    30553009           if (dustbin.ne.0) then
    3056          
     3010
    30573011             call wstats(ngrid,'tau','taudust','SI',2,tau(1,1))
    3058              
     3012
    30593013             if (doubleq) then
    30603014               call wstats(ngrid,'dqsdust',
    … …  
    30943048     &                        'part/kg',3,nccn)
    30953049             endif ! (scavenging)
    3096          
     3050
    30973051           endif ! (dustbin.ne.0)
    30983052
    … …  
    31403094
    31413095                    do ig = 1,ngrid
    3142                        colden(ig,iq) = 0.                           
     3096                       colden(ig,iq) = 0.
    31433097                    end do
    3144                     do l=1,nlayer                                   
    3145                        do ig=1,ngrid                                 
     3098                    do l=1,nlayer
     3099                       do ig=1,ngrid
    31463100                          colden(ig,iq) = colden(ig,iq) + zq(ig,l,iq)
    3147      $                                   *(zplev(ig,l)-zplev(ig,l+1)) 
    3148      $                                   *6.022e22/(mmol(iq)*g)       
    3149                        end do                                       
    3150                     end do                                         
    3151 
    3152                     call wstats(ngrid,"c_"//trim(noms(iq)),           
    3153      $                          "column","mol cm-2",2,colden(1,iq)) 
     3101     $                                   *(zplev(ig,l)-zplev(ig,l+1))
     3102     $                                   *6.022e22/(mmol(iq)*g)
     3103                       end do
     3104                    end do
     3105
     3106                    call wstats(ngrid,"c_"//trim(noms(iq)),
     3107     $                          "column","mol cm-2",2,colden(1,iq))
    31543108                    call write_output("c_"//trim(noms(iq)),
    31553109     $                          "column","mol cm-2",colden(:,iq))
    31563110
    31573111!                   global mass (g)
    3158                
     3112
    31593113                    call planetwide_sumval(colden(:,iq)/6.022e23
    31603114     $                            *mmol(iq)*1.e4*cell_area(:),mass(iq))
    … …  
    31663120              end do ! of do iq=1,nq
    31673121           end if ! of if (photochem)
    3168 
    31693122
    31703123           IF(lastcall.and.callstats) THEN
    … …  
    31723125             call mkstats(ierr)
    31733126           ENDIF
    3174 
    31753127
    31763128c        (Store EOF for Mars Climate database software)
    … …  
    31823134
    31833135#ifdef MESOSCALE
    3184      
    3185       !! see comm_wrf. 
     3136
     3137      !! see comm_wrf.
    31863138      !! not needed when an array is already in a shared module.
    31873139      !! --> example : hfmax_th, zmax_th
    … …  
    32253177      !state  real  RICE      ikj   misc  1  -  h  "RICE"      "ICE RADIUS"                      "m"
    32263178      comm_RICE(1:ngrid,1:nlayer) = rice(1:ngrid,1:nlayer)
    3227    
     3179
    32283180      !! calculate sensible heat flux in W/m2 for outputs
    32293181      !! -- the one computed in vdifc is not the real one
    … …  
    32333185     .         - capcal(1:ngrid)*zdtsdif(1:ngrid)
    32343186      else
    3235         sensibFlux(1:ngrid) = 
     3187        sensibFlux(1:ngrid) =
    32363188     &   (pplay(1:ngrid,1)/(r*pt(1:ngrid,1)))*cpp
    32373189     &   *sqrt(pu(1:ngrid,1)*pu(1:ngrid,1)+pv(1:ngrid,1)*pv(1:ngrid,1)
    … …  
    32393191     &   *zcdh(1:ngrid)*(tsurf(1:ngrid)-zh(1:ngrid,1))
    32403192      endif
    3241          
     3193
    32423194#else
    32433195#ifndef MESOINI
    … …  
    33793331c        Outputs of the CO2 cycle
    33803332c        ----------------------------------------------------------
    3381 
    33823333      if (igcm_co2.ne.0) then
    33833334        call write_output("co2","co2 mass mixing ratio",
    … …  
    36253576     &                 'Visible dust optical depth at 610Pa in the GCM',
    36263577     &                 'NU',tau_pref_gcm(:))
    3627      
     3578
    36283579      if (reff_driven_IRtoVIS_scenario) then
    36293580             call write_output('IRtoVIScoef',
    … …  
    36493600     &                        'part/kg',ndust(:,:))
    36503601
    3651              select case (trim(dustiropacity))
     3602         select case (trim(dustiropacity))
    36523603              case ("tes")
    36533604               call write_output('dsodust_TES',
    … …  
    37433694     &                        'm',rstormdust(:,:))
    37443695
    3745              select case (trim(dustiropacity))
     3696         select case (trim(dustiropacity))
    37463697              case ("tes")
    37473698               call write_output('dsords_TES',
    … …  
    37633714             call write_output('topdustN','top Dust number',
    37643715     &                        'part/kg',pq(:,:,igcm_topdust_number))
    3765              select case (trim(dustiropacity))
     3716         select case (trim(dustiropacity))
    37663717              case ("tes")
    37673718               call write_output('dsotop_TES',
    … …  
    38413792c        Thermospheric outputs
    38423793c        ----------------------------------------------------------
    3843 
    38443794         if(callthermos) then
    38453795
    … …  
    38953845c        ----------------------------------------------------------
    38963846
    3897 
    38983847c        ----------------------------------------------------------
    38993848c        Output in netcdf file "diagsoil.nc" for subterranean
    39003849c          variables (output every "ecritphy", as for writediagfi)
    39013850c        ----------------------------------------------------------
    3902 
    39033851         ! Write soil temperature
    39043852        call write_output("soiltemp","Soil temperature","K",
    … …  
    39663914
    39673915            DO islope = 1,nslope
    3968             ! Clapeyron law for psat (psat = exp(beta/Th2o+alpha)),following Murphy and Koop 2005
     3916        ! Clapeyron law for psat (psat = exp(beta/Th2o+alpha)),following Murphy and Koop 2005
    39693917             rhowater_surf_sat(ig,islope)  =
    39703918     &         exp(beta_clap_h2o/tsurf(ig,islope)+alpha_clap_h2o)
    39713919     &         / tsurf(ig,islope)
    39723920     &         * mmol(igcm_h2o_vap)/(mugaz*r)
    3973  
     3921
    39743922             if(qsurf(ig,igcm_h2o_ice,islope).gt.(1.e-4)) then
    3975                ! we consider to be at saturation above 1.e-4 kg.m-2
     3923           ! we consider to be at saturation above 1.e-4 kg.m-2
    39763924               rhowater_surf(ig,islope) = rhowater_surf_sat(ig,islope)
    39773925             else
    3978               ! otherwise, use vapor partial pressure
     3926          ! otherwise, use vapor partial pressure
    39793927              rhowater_surf(ig,islope) = pvap_surf(ig)
    39803928     &         / tsurf(ig,islope)
    … …  
    40293977
    40303978c      END IF       ! if(ngrid.ne.1)
    4031 
    40323979
    40333980! test for co2 conservation with co2 microphysics
    … …  
    40574004     &                     'kg', co2conservation)
    40584005! XIOS outputs
    4059 #ifdef CPP_XIOS     
     4006#ifdef CPP_XIOS
    40604007      ! Send fields to XIOS: (NB these fields must also be defined as
    40614008      ! <field id="..." /> in context_lmdz_physics.xml to be correctly used)

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

    r3183 r3203  
    179179if ((theta_s > 90.) .or. (theta_s < 0.)) then
    180180    write(*,*) 'please set theta_s between 0 and 90', theta_s
    181     call abort_physic("param_slope","invalid theta_s",1)
     181    call abort_physic("param_slopes","invalid theta_s",1)
    182182endif
    183183

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

    r3167 r3203  
    1111     $                pdufi,pdvfi,pdhfi,pdqfi,pfluxsrf,
    1212     $                pdudif,pdvdif,pdhdif,pdtsrf,pq2,
    13      $                pdqdif,pdqsdif,wstar,zcdv_true,zcdh_true,
     13     $                pdqdif,pdqsdif,wstar,
    1414     $                hfmax,pcondicea_co2microp,sensibFlux,
    1515     $                dustliftday,local_time,watercap, dwatercap_dif)
    … …  
    114114      REAL zkq(ngrid,nlay+1)
    115115      REAL zcdv(ngrid,nslope),zcdh(ngrid,nslope)
    116       REAL, INTENT(OUT) :: zcdv_true(ngrid,nslope)
    117       REAL, INTENT(OUT) :: zcdh_true(ngrid,nslope)    ! drag coeff are used by the LES to recompute u* and hfx
    118       REAL :: zcdv_tmp(ngrid),zcdh_tmp(ngrid)    ! drag coeffs for the major sub-grid surface
     116      REAL :: zcdv_true(ngrid,nslope)
     117      REAL :: zcdh_true(ngrid,nslope)         ! drag coeff are used by the LES to recompute u* and hfx
     118      REAL :: zcdv_tmp(ngrid),zcdh_tmp(ngrid) ! drag coeffs for the major sub-grid surface
    119119      REAL :: zcdv_true_tmp(ngrid),zcdh_true_tmp(ngrid)    ! drag coeffs (computed with wind gustiness for the major sub-grid surface
    120120      REAL zu(ngrid,nlay),zv(ngrid,nlay)
Note: See TracChangeset for help on using the changeset viewer.