Changeset 1477 for trunk/LMDZ.GENERIC


Ignore:
Timestamp:
Sep 23, 2015, 11:33:47 AM (9 years ago)
Author:
mturbet
Message:

Cleaning of physiq.F90. Condense_cloud becomes Condense_co2

Location:
trunk/LMDZ.GENERIC
Files:
5 edited
1 moved

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.GENERIC/README

    r1470 r1477  
    10831083  aerosol properties files and surface files in datadir) still works.
    10841084
     1085== 21/09/2015 == MT
     1086- Cleanup of physiq.F90 and changed condense_cloud.F90 to condense_co2.F90

  • trunk/LMDZ.GENERIC/libf/phystd/callsedim.F

    r1397 r1477  
    107107      do iq=1,nq
    108108       if((radius(iq).gt.1.e-9).and.(iq.ne.igcm_co2_ice)) then   
    109 !         (no sedim for gases, and co2_ice sedim is done in condense_cloud)     
     109!         (no sedim for gases, and co2_ice sedim is done in condense_co2)     
    110110
    111111! store locally updated tracers

  • trunk/LMDZ.GENERIC/libf/phystd/condense_co2.F90

    r1476 r1477  
    1       subroutine condense_cloud(ngrid,nlayer,nq,ptimestep, &
     1      subroutine condense_co2(ngrid,nlayer,nq,ptimestep, &
    22          pcapcal,pplay,pplev,ptsrf,pt, &
    33          pphi,pdt,pdu,pdv,pdtsrf,pu,pv,pq,pdq, &
    … …  
    164164         enddo
    165165         
    166         write(*,*) "condense_cloud: i_co2ice=",i_co2ice       
     166        write(*,*) "condense_co2: i_co2ice=",i_co2ice       
    167167
    168168        if((i_co2ice.lt.1))then
    … …  
    472472
    473473      return
    474     end subroutine condense_cloud
     474    end subroutine condense_co2
    475475
    476476!-------------------------------------------------------------------------

  • trunk/LMDZ.GENERIC/libf/phystd/physiq.F90

    r1429 r1477  
    5252!   It includes:
    5353!
    54 !      1.  Initialization:
    55 !      1.1 Firstcall initializations
    56 !      1.2 Initialization for every call to physiq
    57 !      1.2.5 Compute mean mass and cp, R and thermal conduction coeff.
    58 !      2. Compute radiative transfer tendencies
    59 !         (longwave and shortwave).
    60 !      4. Vertical diffusion (turbulent mixing):
    61 !      5. Convective adjustment
    62 !      6. Condensation and sublimation of gases (currently just CO2).
    63 !      7. TRACERS
    64 !       7a. water and water ice
    65 !       7c. other schemes for tracer transport (lifting, sedimentation)
    66 !       7d. updates (pressure variations, surface budget)
    67 !      9. Surface and sub-surface temperature calculations
    68 !     10. Write outputs :
    69 !           - "startfi", "histfi" if it's time
    70 !           - Saving statistics if "callstats = .true."
    71 !           - Output any needed variables in "diagfi"
    72 !     10. Diagnostics: mass conservation of tracers, radiative energy balance etc.
     54!      I. Initialization :
     55!         I.1 Firstcall initializations.
     56!         I.2 Initialization for every call to physiq.
     57!
     58!      II. Compute radiative transfer tendencies (longwave and shortwave) :
     59!         II.a Option 1 : Call correlated-k radiative transfer scheme.
     60!         II.b Option 2 : Call Newtonian cooling scheme.
     61!         II.c Option 3 : Atmosphere has no radiative effect.
     62!
     63!      III. Vertical diffusion (turbulent mixing) :
     64!
     65!      IV. Dry Convective adjustment :
     66!
     67!      V. Condensation and sublimation of gases (currently just CO2).
     68!
     69!      VI. Tracers
     70!         VI.1. Water and water ice.
     71!         VI.2. Aerosols and particles.
     72!         VI.3. Updates (pressure variations, surface budget).
     73!         VI.4. Slab Ocean.
     74!         VI.5. Surface Tracer Update.
     75!
     76!      VII. Surface and sub-surface soil temperature.
     77!
     78!      VIII. Perform diagnostics and write output files.
     79!
    7380!
    7481!   arguments
    7582!   ---------
    7683!
    77 !   input
     84!   INPUT
    7885!   -----
    79 !    ecri                  period (in dynamical timestep) to write output
     86!
    8087!    ngrid                 Size of the horizontal grid.
    81 !                          All internal loops are performed on that grid.
    8288!    nlayer                Number of vertical layers.
    83 !    nq                    Number of advected fields
    84 !    firstcall             True at the first call
    85 !    lastcall              True at the last call
    86 !    pday                  Number of days counted from the North. Spring
    87 !                          equinoxe.
    88 !    ptime                 Universal time (0<ptime<1): ptime=0.5 at 12:00 UT
    89 !    ptimestep             timestep (s)
    90 !    pplay(ngrid,nlayer)   Pressure at the middle of the layers (Pa)
    91 !    pplev(ngrid,nlayer+1) intermediate pressure levels (pa)
    92 !    pphi(ngrid,nlayer)    Geopotential at the middle of the layers (m2s-2)
    93 !    pu(ngrid,nlayer)      u component of the wind (ms-1)
    94 !    pv(ngrid,nlayer)      v component of the wind (ms-1)
    95 !    pt(ngrid,nlayer)      Temperature (K)
    96 !    pq(ngrid,nlayer,nq)   Advected fields
     89!    nq                    Number of advected fields.
     90!    nametrac              Name of corresponding advected fields.
     91!
     92!    firstcall             True at the first call.
     93!    lastcall              True at the last call.
     94!
     95!    pday                  Number of days counted from the North. Spring equinoxe.
     96!    ptime                 Universal time (0<ptime<1): ptime=0.5 at 12:00 UT.
     97!    ptimestep             timestep (s).
     98!
     99!    pplay(ngrid,nlayer)   Pressure at the middle of the layers (Pa).
     100!    pplev(ngrid,nlayer+1) Intermediate pressure levels (Pa).
     101!    pphi(ngrid,nlayer)    Geopotential at the middle of the layers (m2.s-2).
     102!
     103!    pu(ngrid,nlayer)      u, zonal component of the wind (ms-1).
     104!    pv(ngrid,nlayer)      v, meridional component of the wind (ms-1).
     105!
     106!    pt(ngrid,nlayer)      Temperature (K).
     107!
     108!    pq(ngrid,nlayer,nq)   Advected fields.
     109!
    97110!    pudyn(ngrid,nlayer)    \
    98111!    pvdyn(ngrid,nlayer)     \ Dynamical temporal derivative for the
    99 !    ptdyn(ngrid,nlayer)     / corresponding variables
     112!    ptdyn(ngrid,nlayer)     / corresponding variables.
    100113!    pqdyn(ngrid,nlayer,nq) /
    101114!    flxw(ngrid,nlayer)      vertical mass flux (kg/s) at layer lower boundary
    102115!
    103 !   output
     116!   OUTPUT
    104117!   ------
    105118!
    … …  
    127140!           Loops converted to F90 matrix format: J. Leconte (2012)
    128141!           No more ngridmx/nqmx, F90 commons and adaptation to parallel: A. Spiga (2012)
     142!           Purge of the code : M. Turbet (2015)
    129143!==================================================================
    130144
    … …  
    138152! -----------
    139153
    140 !   inputs:
     154!   INPUTS:
    141155!   -------
    142       integer,intent(in) :: ngrid ! number of atmospheric columns
    143       integer,intent(in) :: nlayer ! number of atmospheric layers
    144       integer,intent(in) :: nq ! number of tracers
    145       character*20,intent(in) :: nametrac(nq) ! name of the tracer from dynamics
    146       logical,intent(in) :: firstcall ! signals first call to physics
    147       logical,intent(in) :: lastcall ! signals last call to physics
    148       real,intent(in) :: pday ! number of elapsed sols since reference Ls=0
    149       real,intent(in) :: ptime ! "universal time", given as fraction of sol (e.g.: 0.5 for noon)
    150       real,intent(in) :: ptimestep ! physics timestep (s)
    151       real,intent(in) :: pplev(ngrid,nlayer+1) ! inter-layer pressure (Pa)
    152       real,intent(in) :: pplay(ngrid,nlayer) ! mid-layer pressure (Pa)
    153       real,intent(in) :: pphi(ngrid,nlayer) ! geopotential at mid-layer (m2s-2)
    154       real,intent(in) :: pu(ngrid,nlayer) ! zonal wind component (m/s)
    155       real,intent(in) :: pv(ngrid,nlayer) ! meridional wind component (m/s)
    156       real,intent(in) :: pt(ngrid,nlayer) ! temperature (K)
    157       real,intent(in) :: pq(ngrid,nlayer,nq) ! tracers (.../kg_of_air)
    158       real,intent(in) :: flxw(ngrid,nlayer) ! vertical mass flux (ks/s)
    159                                             ! at lower boundary of layer
    160 
    161 
    162 !   outputs:
     156
     157      integer,intent(in) :: ngrid             ! Number of atmospheric columns.
     158      integer,intent(in) :: nlayer            ! Number of atmospheric layers.
     159      integer,intent(in) :: nq                ! Number of tracers.
     160      character*20,intent(in) :: nametrac(nq) ! Names of the tracers taken from dynamics.
     161     
     162      logical,intent(in) :: firstcall ! Signals first call to physics.
     163      logical,intent(in) :: lastcall  ! Signals last call to physics.
     164     
     165      real,intent(in) :: pday                  ! Number of elapsed sols since reference Ls=0.
     166      real,intent(in) :: ptime                 ! "Universal time", given as fraction of sol (e.g.: 0.5 for noon).
     167      real,intent(in) :: ptimestep             ! Physics timestep (s).
     168      real,intent(in) :: pplev(ngrid,nlayer+1) ! Inter-layer pressure (Pa).
     169      real,intent(in) :: pplay(ngrid,nlayer)   ! Mid-layer pressure (Pa).
     170      real,intent(in) :: pphi(ngrid,nlayer)    ! Geopotential at mid-layer (m2s-2).
     171      real,intent(in) :: pu(ngrid,nlayer)      ! Zonal wind component (m/s).
     172      real,intent(in) :: pv(ngrid,nlayer)      ! Meridional wind component (m/s).
     173      real,intent(in) :: pt(ngrid,nlayer)      ! Temperature (K).
     174      real,intent(in) :: pq(ngrid,nlayer,nq)   ! Tracers (kg/kg_of_air).
     175      real,intent(in) :: flxw(ngrid,nlayer)    ! Vertical mass flux (ks/s) at lower boundary of layer
     176
     177!   OUTPUTS:
    163178!   --------
    164 !     physical tendencies
    165       real,intent(out) :: pdu(ngrid,nlayer) ! zonal wind tendency (m/s/s)
    166       real,intent(out) :: pdv(ngrid,nlayer) ! meridional wind tendency (m/s/s)
    167       real,intent(out) :: pdt(ngrid,nlayer) ! temperature tendency (K/s)
    168       real,intent(out) :: pdq(ngrid,nlayer,nq) ! tracer tendencies (../kg/s)
    169       real,intent(out) :: pdpsrf(ngrid) ! surface pressure tendency (Pa/s)
    170       logical,intent(out) :: tracerdyn ! signal to the dynamics to advect tracers or not
     179
     180!     Physical tendencies :
     181
     182      real,intent(out) :: pdu(ngrid,nlayer)    ! Zonal wind tendencies (m/s/s).
     183      real,intent(out) :: pdv(ngrid,nlayer)    ! Meridional wind tendencies (m/s/s).
     184      real,intent(out) :: pdt(ngrid,nlayer)    ! Temperature tendencies (K/s).
     185      real,intent(out) :: pdq(ngrid,nlayer,nq) ! Tracer tendencies (kg/kg_of_air/s).
     186      real,intent(out) :: pdpsrf(ngrid)        ! Surface pressure tendency (Pa/s).
     187      logical,intent(out) :: tracerdyn         ! Signal to the dynamics to advect tracers or not.     
    171188
    172189! Local saved variables:
    173190! ----------------------
    174 !     aerosol (dust or ice) extinction optical depth  at reference wavelength
    175 !     "longrefvis" set in dimradmars.h , for one of the "naerkind"  kind of
    176 !      aerosol optical properties:
    177 !      real aerosol(ngrid,nlayer,naerkind)
    178 !     this is now internal to callcorrk and hence no longer needed here
    179 
    180       integer day_ini                ! Initial date of the run (sol since Ls=0)
    181       integer icount                 ! counter of calls to physiq during the run.
    182       real, dimension(:),allocatable,save ::  tsurf  ! Surface temperature (K)
    183       real, dimension(:,:),allocatable,save ::  tsoil  ! sub-surface temperatures (K)
    184       real, dimension(:),allocatable,save :: albedo ! Surface albedo
     191
     192      integer day_ini                                  ! Initial date of the run (sol since Ls=0).
     193      integer icount                                   ! Counter of calls to physiq during the run.
     194      real, dimension(:),allocatable,save ::  tsurf    ! Surface temperature (K).
     195      real, dimension(:,:),allocatable,save ::  tsoil  ! Sub-surface temperatures (K).
     196      real, dimension(:),allocatable,save :: albedo    ! Surface albedo.
     197     
    185198!$OMP THREADPRIVATE(tsurf,tsoil,albedo)
    186199
    187       real,dimension(:),allocatable,save :: albedo0 ! Surface albedo
    188       real,dimension(:),allocatable,save :: rnat ! added by BC
     200      real,dimension(:),allocatable,save :: albedo0 ! Initial surface albedo.
     201      real,dimension(:),allocatable,save :: rnat    ! Defines the type of the grid (ocean,continent,...). By BC.
     202     
    189203!$OMP THREADPRIVATE(albedo0,rnat)
    190204
    191       real,dimension(:),allocatable,save :: emis ! Thermal IR surface emissivity
    192       real,dimension(:,:),allocatable,save :: dtrad ! Net atm. radiative heating rate (K.s-1)
    193       real,dimension(:),allocatable,save :: fluxrad_sky ! rad. flux from sky absorbed by surface (W.m-2)
    194       real,dimension(:),allocatable,save :: fluxrad ! Net radiative surface flux (W.m-2)
    195       real,dimension(:),allocatable,save :: capcal ! surface heat capacity (J m-2 K-1)
    196       real,dimension(:),allocatable,save :: fluxgrd ! surface conduction flux (W.m-2)
    197       real,dimension(:,:),allocatable,save :: qsurf ! tracer on surface (e.g. kg.m-2)
    198       real,dimension(:,:),allocatable,save :: q2 ! Turbulent Kinetic Energy
     205      real,dimension(:),allocatable,save :: emis        ! Thermal IR surface emissivity.
     206      real,dimension(:,:),allocatable,save :: dtrad     ! Net atmospheric radiative heating rate (K.s-1).
     207      real,dimension(:),allocatable,save :: fluxrad_sky ! Radiative flux from sky absorbed by surface (W.m-2).
     208      real,dimension(:),allocatable,save :: fluxrad     ! Net radiative surface flux (W.m-2).
     209      real,dimension(:),allocatable,save :: capcal      ! Surface heat capacity (J m-2 K-1).
     210      real,dimension(:),allocatable,save :: fluxgrd     ! Surface conduction flux (W.m-2).
     211      real,dimension(:,:),allocatable,save :: qsurf     ! Tracer on surface (e.g. kg.m-2).
     212      real,dimension(:,:),allocatable,save :: q2        ! Turbulent Kinetic Energy.
     213     
    199214!$OMP THREADPRIVATE(emis,dtrad,fluxrad_sky,fluxrad,capcal,fluxgrd,qsurf,q2)
    200215
    201216      save day_ini, icount
     217     
    202218!$OMP THREADPRIVATE(day_ini,icount)
    203219
    … …  
    205221! -----------------
    206222
    207 !     aerosol (dust or ice) extinction optical depth  at reference wavelength
     223!     Aerosol (dust or ice) extinction optical depth  at reference wavelength
    208224!     for the "naerkind" optically active aerosols:
    209       real aerosol(ngrid,nlayer,naerkind)
    210       real zh(ngrid,nlayer)      ! potential temperature (K)
    211       real pw(ngrid,nlayer) ! vertical velocity (m/s) (>0 when downwards)
    212 
    213       character*80 fichier
    214       integer l,ig,ierr,iq,iaer
    215      
    216       !!! this is saved for diagnostic
    217       real,dimension(:),allocatable,save :: fluxsurf_lw ! incident LW (IR) surface flux (W.m-2)
    218       real,dimension(:),allocatable,save :: fluxsurf_sw ! incident SW (stellar) surface flux (W.m-2)
    219       real,dimension(:),allocatable,save :: fluxtop_lw ! Outgoing LW (IR) flux to space (W.m-2)
    220       real,dimension(:),allocatable,save :: fluxabs_sw ! Absorbed SW (stellar) flux (W.m-2)
    221       real,dimension(:),allocatable,save :: fluxtop_dn
    222       real,dimension(:),allocatable,save :: fluxdyn ! horizontal heat transport by dynamics 
    223       real,dimension(:,:),allocatable,save :: OLR_nu ! Outgoing LW radition in each band (Normalized to the band width (W/m2/cm-1)
    224       real,dimension(:,:),allocatable,save :: OSR_nu ! Outgoing SW radition in each band (Normalized to the band width (W/m2/cm-1)
    225       real,dimension(:,:),allocatable,save :: zdtlw ! (K/s)
    226       real,dimension(:,:),allocatable,save :: zdtsw ! (K/s)
    227       real,dimension(:),allocatable,save :: sensibFlux ! turbulent flux given by the atm to the surface
     225
     226      real aerosol(ngrid,nlayer,naerkind) ! Aerosols.
     227      real zh(ngrid,nlayer)               ! Potential temperature (K).
     228      real pw(ngrid,nlayer)               ! Vertical velocity (m/s). (NOTE : >0 WHEN DOWNWARDS !!)
     229
     230      integer l,ig,ierr,iq
     231     
     232      ! FOR DIAGNOSTIC :
     233     
     234      real,dimension(:),allocatable,save :: fluxsurf_lw  ! Incident Long Wave (IR) surface flux (W.m-2).
     235      real,dimension(:),allocatable,save :: fluxsurf_sw  ! incident Short Wave (stellar) surface flux (W.m-2).
     236      real,dimension(:),allocatable,save :: fluxtop_lw   ! Outgoing LW (IR) flux to space (W.m-2).
     237      real,dimension(:),allocatable,save :: fluxabs_sw   ! Absorbed SW (stellar) flux (W.m-2).
     238      real,dimension(:),allocatable,save :: fluxtop_dn   ! Incoming SW (stellar) radiation at the top of the atmosphere (W.m-2).
     239      real,dimension(:),allocatable,save :: fluxdyn      ! Horizontal heat transport by dynamics (W.m-2).
     240      real,dimension(:,:),allocatable,save :: OLR_nu     ! Outgoing LW radiation in each band (Normalized to the band width (W/m2/cm-1)).
     241      real,dimension(:,:),allocatable,save :: OSR_nu     ! Outgoing SW radiation in each band (Normalized to the band width (W/m2/cm-1)).
     242      real,dimension(:,:),allocatable,save :: zdtlw      ! LW heating tendencies (K/s).
     243      real,dimension(:,:),allocatable,save :: zdtsw      ! SW heating tendencies (K/s).
     244      real,dimension(:),allocatable,save :: sensibFlux   ! Turbulent flux given by the atmosphere to the surface (W.m-2).
     245     
    228246!$OMP THREADPRIVATE(fluxsurf_lw,fluxsurf_sw,fluxtop_lw,fluxabs_sw,fluxtop_dn,fluxdyn,OLR_nu,OSR_nu,&
    229         !$OMP zdtlw,zdtsw,sensibFlux)
    230 
    231       real zls                       ! solar longitude (rad)
    232       real zlss                      ! sub solar point longitude (rad)
    233       real zday                      ! date (time since Ls=0, in martian days)
    234       real zzlay(ngrid,nlayer)   ! altitude at the middle of the layers
    235       real zzlev(ngrid,nlayer+1) ! altitude at layer boundaries
    236       real latvl1,lonvl1             ! Viking Lander 1 point (for diagnostic)
    237 
    238 !     Tendencies due to various processes:
    239       real dqsurf(ngrid,nq)
    240       real cldtlw(ngrid,nlayer)                           ! (K/s) LW heating rate for clear areas
    241       real cldtsw(ngrid,nlayer)                           ! (K/s) SW heating rate for clear areas
    242       real zdtsurf(ngrid)                                   ! (K/s)
    243       real dtlscale(ngrid,nlayer)
    244       real zdvdif(ngrid,nlayer),zdudif(ngrid,nlayer)  ! (m.s-2)
    245       real zdhdif(ngrid,nlayer), zdtsdif(ngrid)         ! (K/s)
    246       real zdtdif(ngrid,nlayer)                       ! (K/s)
    247       real zdvadj(ngrid,nlayer),zduadj(ngrid,nlayer)  ! (m.s-2)
    248       real zdhadj(ngrid,nlayer)                           ! (K/s)
    249       real zdtgw(ngrid,nlayer)                            ! (K/s)
    250       real zdtmr(ngrid,nlayer)
    251       real zdugw(ngrid,nlayer),zdvgw(ngrid,nlayer)    ! (m.s-2)
    252       real zdtc(ngrid,nlayer),zdtsurfc(ngrid)
    253       real zdvc(ngrid,nlayer),zduc(ngrid,nlayer)
    254       real zdumr(ngrid,nlayer),zdvmr(ngrid,nlayer)
    255       real zdtsurfmr(ngrid)
    256      
    257       real zdmassmr(ngrid,nlayer),zdpsrfmr(ngrid)
    258       real zdmassmr_col(ngrid)
    259 
    260       real zdqdif(ngrid,nlayer,nq), zdqsdif(ngrid,nq)
    261       real zdqevap(ngrid,nlayer)
    262       real zdqsed(ngrid,nlayer,nq), zdqssed(ngrid,nq)
    263       real zdqdev(ngrid,nlayer,nq), zdqsdev(ngrid,nq)
    264       real zdqadj(ngrid,nlayer,nq)
    265       real zdqc(ngrid,nlayer,nq)
    266       real zdqmr(ngrid,nlayer,nq),zdqsurfmr(ngrid,nq)
    267       real zdqlscale(ngrid,nlayer,nq)
    268       real zdqslscale(ngrid,nq)
    269       real zdqchim(ngrid,nlayer,nq)
    270       real zdqschim(ngrid,nq)
    271 
    272       real zdteuv(ngrid,nlayer)    ! (K/s)
    273       real zdtconduc(ngrid,nlayer) ! (K/s)
    274       real zdumolvis(ngrid,nlayer)
    275       real zdvmolvis(ngrid,nlayer)
    276       real zdqmoldiff(ngrid,nlayer,nq)
    277 
    278 !     Local variables for local calculations:
     247       
     248        !$OMP zdtlw,zdtsw,sensibFlux)
     249
     250      real zls                       ! Solar longitude (radians).
     251      real zlss                      ! Sub solar point longitude (radians).
     252      real zday                      ! Date (time since Ls=0, calculated in sols).
     253      real zzlay(ngrid,nlayer)       ! Altitude at the middle of the atmospheric layers.
     254      real zzlev(ngrid,nlayer+1)     ! Altitude at the atmospheric layer boundaries.
     255
     256! TENDENCIES due to various processes :
     257
     258      ! For Surface Temperature : (K/s)
     259      real zdtsurf(ngrid)     ! Cumulated tendencies.
     260      real zdtsurfmr(ngrid)   ! Mass_redistribution routine.
     261      real zdtsurfc(ngrid)    ! Condense_co2 routine.
     262      real zdtsdif(ngrid)     ! Turbdiff/vdifc routines.
     263      real zdtsurf_hyd(ngrid) ! Hydrol routine.
     264           
     265      ! For Atmospheric Temperatures : (K/s)   
     266      real dtlscale(ngrid,nlayer)                             ! Largescale routine.
     267      real zdtc(ngrid,nlayer)                                 ! Condense_co2 routine.
     268      real zdtdif(ngrid,nlayer)                               ! Turbdiff/vdifc routines.
     269      real zdtmr(ngrid,nlayer)                                ! Mass_redistribution routine.
     270      real zdtrain(ngrid,nlayer)                              ! Rain routine.
     271      real dtmoist(ngrid,nlayer)                              ! Moistadj routine.
     272      real dt_ekman(ngrid,noceanmx), dt_hdiff(ngrid,noceanmx) ! Slab_ocean routine.
     273      real zdtsw1(ngrid,nlayer), zdtlw1(ngrid,nlayer)         ! Callcorrk routine.
     274                             
     275      ! For Surface Tracers : (kg/m2/s)
     276      real dqsurf(ngrid,nq)                 ! Cumulated tendencies.
     277      real zdqslscale(ngrid,nq)             ! Largescale routine.
     278      real zdqsdif(ngrid,nq)                ! Turbdiff/vdifc routines.
     279      real zdqssed(ngrid,nq)                ! Callsedim routine.
     280      real zdqsurfmr(ngrid,nq)              ! Mass_redistribution routine.
     281      real zdqsrain(ngrid), zdqssnow(ngrid) ! Rain routine.
     282      real dqs_hyd(ngrid,nq)                ! Hydrol routine.
     283                 
     284      ! For Tracers : (kg/kg_of_air/s)
     285      real zdqc(ngrid,nlayer,nq)      ! Condense_co2 routine.
     286      real zdqadj(ngrid,nlayer,nq)    ! Convadj routine.
     287      real zdqdif(ngrid,nlayer,nq)    ! Turbdiff/vdifc routines.
     288      real zdqevap(ngrid,nlayer)      ! Turbdiff routine.
     289      real zdqsed(ngrid,nlayer,nq)    ! Callsedim routine.
     290      real zdqmr(ngrid,nlayer,nq)     ! Mass_redistribution routine.
     291      real zdqrain(ngrid,nlayer,nq)   ! Rain routine.
     292      real dqmoist(ngrid,nlayer,nq)   ! Moistadj routine.
     293      real dqvaplscale(ngrid,nlayer)  ! Largescale routine.
     294      real dqcldlscale(ngrid,nlayer)  ! Largescale routine.
     295                       
     296      ! For Winds : (m/s/s)
     297      real zdvc(ngrid,nlayer),zduc(ngrid,nlayer)     ! Condense_co2 routine.
     298      real zdvadj(ngrid,nlayer),zduadj(ngrid,nlayer) ! Convadj routine.
     299      real zdumr(ngrid,nlayer),zdvmr(ngrid,nlayer)   ! Mass_redistribution routine.
     300      real zdvdif(ngrid,nlayer),zdudif(ngrid,nlayer) ! Turbdiff/vdifc routines.
     301      real zdhdif(ngrid,nlayer)                      ! Turbdiff/vdifc routines.
     302      real zdhadj(ngrid,nlayer)                      ! Convadj routine.
     303     
     304      ! For Pressure and Mass :
     305      real zdmassmr(ngrid,nlayer) ! Atmospheric Mass tendency for mass_redistribution (kg_of_air/m2/s).
     306      real zdmassmr_col(ngrid)    ! Atmospheric Column Mass tendency for mass_redistribution (kg_of_air/m2/s).
     307      real zdpsrfmr(ngrid)        ! Pressure tendency for mass_redistribution routine (Pa/s).
     308
     309     
     310   
     311! Local variables for LOCAL CALCULATIONS:
     312! ---------------------------------------
    279313      real zflubid(ngrid)
    280314      real zplanck(ngrid),zpopsk(ngrid,nlayer)
    281       real zdum1(ngrid,nlayer)
    282       real zdum2(ngrid,nlayer)
    283315      real ztim1,ztim2,ztim3, z1,z2
    284316      real ztime_fin
    … …  
    287319      real taux(ngrid),tauy(ngrid)
    288320
    289       integer length
    290       parameter (length=100)
    291 
    292 ! local variables only used for diagnostics (output in file "diagfi" or "stats")
    293 ! ------------------------------------------------------------------------------
    294       real ps(ngrid), zt(ngrid,nlayer)
    295       real zu(ngrid,nlayer),zv(ngrid,nlayer)
    296       real zq(ngrid,nlayer,nq)
    297       character*2 str2
    298       character*5 str5
    299       real zdtadj(ngrid,nlayer)
    300       real zdtdyn(ngrid,nlayer)
    301       real,allocatable,dimension(:,:),save :: ztprevious
     321
     322! local variables for DIAGNOSTICS : (diagfi & stat)
     323! -------------------------------------------------
     324      real ps(ngrid)                                     ! Surface Pressure.
     325      real zt(ngrid,nlayer)                              ! Atmospheric Temperature.
     326      real zu(ngrid,nlayer),zv(ngrid,nlayer)             ! Zonal and Meridional Winds.
     327      real zq(ngrid,nlayer,nq)                           ! Atmospheric Tracers.
     328      real zdtadj(ngrid,nlayer)                          ! Convadj Diagnostic.
     329      real zdtdyn(ngrid,nlayer)                          ! Dynamical Heating (K/s).
     330      real,allocatable,dimension(:,:),save :: ztprevious ! Previous loop Atmospheric Temperature (K)                                                         ! Useful for Dynamical Heating calculation.
    302331!$OMP THREADPRIVATE(ztprevious)
    303       real reff(ngrid,nlayer) ! effective dust radius (used if doubleq=T)
    304       real qtot1,qtot2            ! total aerosol mass
    305       integer igmin, lmin
    306       logical tdiag
    307 
    308       real zplev(ngrid,nlayer+1),zplay(ngrid,nlayer)
    309       real zstress(ngrid), cd
    310       real hco2(nq), tmean, zlocal(nlayer)
    311       real vmr(ngrid,nlayer) ! volume mixing ratio
    312 
     332
     333      real reff(ngrid,nlayer)                       ! Effective dust radius (used if doubleq=T).
     334      real vmr(ngrid,nlayer)                        ! volume mixing ratio
    313335      real time_phys
    314 
    315 !     reinstated by RW for diagnostic
    316       real,allocatable,dimension(:),save :: tau_col
     336      real,allocatable,dimension(:),save :: tau_col ! Total Aerosol Optical Depth.
    317337!$OMP THREADPRIVATE(tau_col)
    318338     
    319 !     included by RW to reduce insanity of code
    320       real ISR,ASR,OLR,GND,DYN,GSR,Ts1,Ts2,Ts3,TsS
    321 
    322 !     included by RW to compute tracer column densities
    323       real qcol(ngrid,nq)
    324 
    325 !     included by RW for H2O precipitation
    326       real zdtrain(ngrid,nlayer)
    327       real zdqrain(ngrid,nlayer,nq)
    328       real zdqsrain(ngrid)
    329       real zdqssnow(ngrid)
    330       real rainout(ngrid)
     339      real ISR,ASR,OLR,GND,DYN,GSR,Ts1,Ts2,Ts3,TsS ! for Diagnostic.
     340     
     341      real qcol(ngrid,nq) ! Tracer Column Mass (kg/m2).
    331342
    332343!     included by RW for H2O Manabe scheme
    333       real dtmoist(ngrid,nlayer)
    334       real dqmoist(ngrid,nlayer,nq)
    335 
    336       real qvap(ngrid,nlayer)
    337       real dqvaplscale(ngrid,nlayer)
    338       real dqcldlscale(ngrid,nlayer)
    339       real rneb_man(ngrid,nlayer)
    340       real rneb_lsc(ngrid,nlayer)
    341 
    342 !     included by RW to account for surface cooling by evaporation
    343       real dtsurfh2olat(ngrid)
     344      real rneb_man(ngrid,nlayer) ! H2O cloud fraction (moistadj).
     345      real rneb_lsc(ngrid,nlayer) ! H2O cloud fraction (large scale).
    344346
    345347
    … …  
    352354      real dEdiffs(ngrid),dEdiff(ngrid)
    353355      real madjdE(ngrid), lscaledE(ngrid),madjdEz(ngrid,nlayer), lscaledEz(ngrid,nlayer)
     356     
    354357!JL12 conservation test for mean flow kinetic energy has been disabled temporarily
    355       real dtmoist_max,dtmoist_min
    356      
     358
     359      real dtmoist_max,dtmoist_min     
    357360      real dItot, dItot_tmp, dVtot, dVtot_tmp
    358361
    359 !     included by BC for evaporation     
    360       real qevap(ngrid,nlayer,nq)
    361       real tevap(ngrid,nlayer)
    362       real dqevap1(ngrid,nlayer)
    363       real dtevap1(ngrid,nlayer)
    364 
    365 !     included by BC for hydrology
    366       real,allocatable,save :: hice(:)
     362      real,allocatable,save :: hice(:) ! Oceanic Ice height. by BC
    367363 !$OMP THREADPRIVATE(hice)     
    368364
    369 !     included by RW to test water conservation (by routine)
    370       real h2otot
     365      real h2otot                      ! Total Amount of water. For diagnostic.
     366      real icesrf,liqsrf,icecol,vapcol ! Total Amounts of water (ice,liq,vap). For diagnostic.
    371367      real dWtot, dWtot_tmp, dWtots, dWtots_tmp
    372       real h2o_surf_all
    373       logical watertest
    374       save watertest
     368      logical,save :: watertest
    375369!$OMP THREADPRIVATE(watertest)
    376370
    377 !     included by RW for RH diagnostic
    378       real qsat(ngrid,nlayer), RH(ngrid,nlayer), H2Omaxcol(ngrid),psat_tmp
    379 
    380 !     included by RW for hydrology
    381       real dqs_hyd(ngrid,nq)
    382       real zdtsurf_hyd(ngrid)
    383 
    384 !     included by RW for water cycle conservation tests
    385       real icesrf,liqsrf,icecol,vapcol
    386 
    387 !     included by BC for double radiative transfer call
    388       logical clearsky
    389       real zdtsw1(ngrid,nlayer)
    390       real zdtlw1(ngrid,nlayer)
    391       real fluxsurf_lw1(ngrid)     
    392       real fluxsurf_sw1(ngrid) 
    393       real fluxtop_lw1(ngrid)   
    394       real fluxabs_sw1(ngrid) 
    395       real tau_col1(ngrid)
    396       real OLR_nu1(ngrid,L_NSPECTI)
    397       real OSR_nu1(ngrid,L_NSPECTV)
     371      real qsat(ngrid,nlayer) ! Water Vapor Volume Mixing Ratio at saturation (kg/kg_of_air).
     372      real RH(ngrid,nlayer)   ! Relative humidity.
     373      real H2Omaxcol(ngrid)   ! Maximum possible H2O column amount (at 100% saturation) (kg/m2).
     374      real psat_tmp
     375     
     376      logical clearsky ! For double radiative transfer call. By BC
     377     
     378      real fluxsurf_lw1(ngrid), fluxsurf_sw1(ngrid), fluxtop_lw1(ngrid), fluxabs_sw1(ngrid) ! For SW/LW flux diagnostics.
     379      real tau_col1(ngrid)                                                                  ! For aerosol optical depth diagnostic.
     380      real OLR_nu1(ngrid,L_NSPECTI), OSR_nu1(ngrid,L_NSPECTV)                               ! For Outgoing Radiation diagnostics.
    398381      real tf, ntf
    399382
    400 !     included by BC for cloud fraction computations
    401       real,allocatable,dimension(:,:),save :: cloudfrac
    402       real,allocatable,dimension(:),save :: totcloudfrac
     383      real,allocatable,dimension(:,:),save :: cloudfrac  ! Fraction of clouds (%).
     384      real,allocatable,dimension(:),save :: totcloudfrac ! Column fraction of clouds (%).
    403385!$OMP THREADPRIVATE(cloudfrac,totcloudfrac)
    404386
    405 !     included by RW for vdifc water conservation test
    406       real nconsMAX
    407       real vdifcncons(ngrid)
    408       real cadjncons(ngrid)
    409 
    410 !      double precision qsurf_hist(ngrid,nq)
     387      real nconsMAX, vdifcncons(ngrid), cadjncons(ngrid) ! Vdfic water conservation test. By RW
     388
    411389      real,allocatable,dimension(:,:),save :: qsurf_hist
    412390!$OMP THREADPRIVATE(qsurf_hist)
    413391
    414 !     included by RW for temp convadj conservation test
    415       real playtest(nlayer)
    416       real plevtest(nlayer)
    417       real ttest(nlayer)
    418       real qtest(nlayer)
    419       integer igtest
    420 
    421 !     included by RW for runway greenhouse 1D study
    422       real muvar(ngrid,nlayer+1)
    423 
    424 !     included by RW for variable H2O particle sizes
    425       real,allocatable,dimension(:,:,:),save :: reffrad ! aerosol effective radius (m)
    426       real,allocatable,dimension(:,:,:),save :: nueffrad ! aerosol effective radius variance
     392      real muvar(ngrid,nlayer+1) ! For Runaway Greenhouse 1D study. By RW
     393
     394      real,allocatable,dimension(:,:,:),save :: reffrad  ! Aerosol effective radius (m). By RW
     395      real,allocatable,dimension(:,:,:),save :: nueffrad ! Aerosol effective radius variance. By RW
    427396!$OMP THREADPRIVATE(reffrad,nueffrad)
    428 !      real :: nueffrad_dummy(ngrid,nlayer,naerkind) !! AS. This is temporary. Check below why.
    429       real :: reffh2oliq(ngrid,nlayer) ! liquid water particles effective radius (m)
    430       real :: reffh2oice(ngrid,nlayer) ! water ice particles effective radius (m)
    431    !   real reffH2O(ngrid,nlayer)
     397
    432398      real reffcol(ngrid,naerkind)
    433399
    434 !     included by RW for sourceevol
     400!     Sourceevol for 'accelerated ice evolution'. By RW
    435401      real,allocatable,dimension(:),save :: ice_initial
    436402      real delta_ice,ice_tot
    437       real,allocatable,dimension(:),save :: ice_min
    438      
     403      real,allocatable,dimension(:),save :: ice_min     
    439404      integer num_run
    440405      logical,save :: ice_update
    441406!$OMP THREADPRIVATE(ice_initial,ice_min,ice_update)
    442407
    443 !     included by BC for slab ocean
     408!     For slab ocean. By BC
    444409      real, dimension(:),allocatable,save ::  pctsrf_sic
    445410      real, dimension(:,:),allocatable,save :: tslab
    … …  
    452417      real :: tsurf2(ngrid)
    453418      real :: flux_o(ngrid),flux_g(ngrid),fluxgrdocean(ngrid)
    454       real :: dt_ekman(ngrid,noceanmx),dt_hdiff(ngrid,noceanmx)
    455419      real :: flux_sens_lat(ngrid)
    456420      real :: qsurfint(ngrid,nq)
    457421
    458422
    459 !=======================================================================
    460 
    461 ! 1. Initialisation
     423!==================================================================================================
     424
    462425! -----------------
    463 
    464 !  1.1   Initialisation only at first call
    465 !  ---------------------------------------
     426! I. INITIALISATION
     427! -----------------
     428
     429! --------------------------------
     430! I.1   First Call Initialisation.
     431! --------------------------------
    466432      if (firstcall) then
    467433
    468         !!! ALLOCATE SAVED ARRAYS
     434        ! Allocate saved arrays.
    469435        ALLOCATE(tsurf(ngrid))
    470436        ALLOCATE(tsoil(ngrid,nsoilmx))   
    … …  
    507473        ALLOCATE(zmasq(ngrid))
    508474        ALLOCATE(knindex(ngrid))
    509 !        ALLOCATE(qsurfint(ngrid,nqmx))
    510 
    511 
    512         !! this is defined in comsaison_h
     475
     476        ! This is defined in comsaison_h
    513477        ALLOCATE(mu0(ngrid))
    514         ALLOCATE(fract(ngrid))
    515 
    516            
    517          !! this is defined in radcommon_h
    518          ALLOCATE(glat(ngrid))
    519 
    520 !        variables set to 0
     478        ALLOCATE(fract(ngrid))           
     479         ! This is defined in radcommon_h
     480        ALLOCATE(glat(ngrid))
     481       
     482
     483!        Variables set to 0
    521484!        ~~~~~~~~~~~~~~~~~~
    522485         dtrad(:,:) = 0.0
    … …  
    526489         zdtlw(:,:) = 0.0
    527490
    528 !        initialize aerosol indexes
    529 !        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    530             call iniaerosol()
    531 
    532      
    533 !        initialize tracer names, indexes and properties
    534 !        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     491
     492!        Initialize aerosol indexes.
     493!        ~~~~~~~~~~~~~~~~~~~~~~~~~~~
     494         call iniaerosol()
     495
     496     
     497!        Initialize tracer names, indexes and properties.
     498!        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    535499         tracerdyn=tracer
    536          IF (.NOT.ALLOCATED(noms)) ALLOCATE(noms(nq)) !! because noms is an argument of physdem1
    537                                                       !! whether or not tracer is on
     500         IF (.NOT.ALLOCATED(noms)) ALLOCATE(noms(nq)) ! (because noms is an argument of physdem1 whether or not tracer is on)
    538501         if (tracer) then
    539502            call initracer(ngrid,nq,nametrac)
    540          endif ! end tracer
    541 
    542 !
    543 
    544 !        read startfi (initial parameters)
     503         endif
     504
     505
     506!        Read 'startfi.nc' file.
    545507!        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    546          call phyetat0(ngrid,nlayer,"startfi.nc",0,0,nsoilmx,nq,   &
    547                day_ini,time_phys,tsurf,tsoil,emis,q2,qsurf,   &
    548                cloudfrac,totcloudfrac,hice,  &
    549                rnat,pctsrf_sic,tslab, tsea_ice,sea_ice)
     508         call phyetat0(ngrid,nlayer,"startfi.nc",0,0,nsoilmx,nq,      &
     509                       day_ini,time_phys,tsurf,tsoil,emis,q2,qsurf,   &
     510                       cloudfrac,totcloudfrac,hice,                   &
     511                       rnat,pctsrf_sic,tslab, tsea_ice,sea_ice)
    550512
    551513         if (pday.ne.day_ini) then
    … …  
    559521         write (*,*) 'In physiq day_ini =', day_ini
    560522
    561 !        Initialize albedo and orbital calculation
     523!        Initialize albedo and orbital calculation.
    562524!        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    563525         call surfini(ngrid,nq,qsurf,albedo0)
    … …  
    571533         endif
    572534
    573 !        Initialize oceanic variables
    574 !        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     535!        Initialize oceanic variables.
     536!        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    575537
    576538         if (ok_slab_ocean)then
    577539
    578540           call ocean_slab_init(ngrid,ptimestep, tslab, &
    579                 sea_ice, pctsrf_sic)
    580 
    581             knindex(:) = 0
    582 
    583             do ig=1,ngrid
     541                                sea_ice, pctsrf_sic)
     542
     543           knindex(:) = 0
     544
     545           do ig=1,ngrid
    584546              zmasq(ig)=1
    585547              knindex(ig) = ig
    … …  
    587549                 zmasq(ig)=0
    588550              endif
    589             enddo
    590 
     551           enddo
    591552
    592553           CALL init_masquv(ngrid,zmasq)
    593554
    594 
    595          endif
    596 
    597 
    598 !        initialize soil
    599 !        ~~~~~~~~~~~~~~~
     555         endif ! end of 'ok_slab_ocean'.
     556
     557
     558!        Initialize soil.
     559!        ~~~~~~~~~~~~~~~~
    600560         if (callsoil) then
     561         
    601562            call soil(ngrid,nsoilmx,firstcall,lastcall,inertiedat, &
    602                 ptimestep,tsurf,tsoil,capcal,fluxgrd)
     563                      ptimestep,tsurf,tsoil,capcal,fluxgrd)
    603564
    604565            if (ok_slab_ocean) then
    605566               do ig=1,ngrid
    606567                  if (nint(rnat(ig)).eq.2) then
    607                     capcal(ig)=capcalocean
    608                     if (pctsrf_sic(ig).gt.0.5) then
    609                       capcal(ig)=capcalseaice
    610                       if (qsurf(ig,igcm_h2o_ice).gt.0.) then
    611                         capcal(ig)=capcalsno
    612                       endif
    613                     endif
     568                     capcal(ig)=capcalocean
     569                     if (pctsrf_sic(ig).gt.0.5) then
     570                        capcal(ig)=capcalseaice
     571                        if (qsurf(ig,igcm_h2o_ice).gt.0.) then
     572                           capcal(ig)=capcalsno
     573                        endif
     574                     endif
    614575                  endif
    615576               enddo
    616             endif
    617 
    618 
    619 
    620 
    621          else
     577            endif ! end of 'ok_slab_ocean'.
     578
     579         else ! else of 'callsoil'.
     580         
    622581            print*,'WARNING! Thermal conduction in the soil turned off'
    623582            capcal(:)=1.e6
    624583            fluxgrd(:)=intheat
    625584            print*,'Flux from ground = ',intheat,' W m^-2'
    626          endif
     585           
     586         endif ! end of 'callsoil'.
     587         
    627588         icount=1
    628589
    629 !        decide whether to update ice at end of run
    630 !        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     590!        Decide whether to update ice at end of run.
     591!        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    631592         ice_update=.false.
    632593         if(sourceevol)then
    … …  
    635596                     status="old",iostat=ierr)
    636597            if (ierr.ne.0) then
    637               write(*,*) "physiq: Error! No num_run file!"
    638               write(*,*) " (which is needed for sourceevol option)"
    639               stop
     598               write(*,*) "physiq: Error! No num_run file!"
     599               write(*,*) " (which is needed for sourceevol option)"
     600               stop
    640601            endif
    641602            read(128,*) num_run
    … …  
    645606       
    646607            if(num_run.ne.0.and.mod(num_run,2).eq.0)then
    647             !if(num_run.ne.0.and.mod(num_run,3).eq.0)then
    648608               print*,'Updating ice at end of this year!'
    649609               ice_update=.true.
    650610               ice_min(:)=1.e4
    651             endif 
    652          endif
    653 
    654 !  In newstart now, will have to be remove (BC 2014)
    655 !        define surface as continent or ocean
    656 !        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     611            endif
     612           
     613         endif ! end of 'sourceevol'.
     614
     615
     616         ! Here is defined the type of the surface : Continent or Ocean.
     617         ! BC2014 : This is now already done in newstart.
     618         ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    657619         if (.not.ok_slab_ocean) then
     620         
    658621           rnat(:)=1.
    659622           do ig=1,ngrid
    660 !             if(iceball.or.oceanball.or.(inertiedat(ig,1).gt.1.E4))then
    661              if(inertiedat(ig,1).gt.1.E4)then
    662                rnat(ig)=0
    663              endif
     623              if(inertiedat(ig,1).gt.1.E4)then
     624                 rnat(ig)=0
     625              endif
    664626           enddo
    665627
    666628           print*,'WARNING! Surface type currently decided by surface inertia'
    667629           print*,'This should be improved e.g. in newstart.F'
    668          endif!(.not.ok_slab_ocean)
    669 
    670 !        initialise surface history variable
     630           
     631         endif ! end of 'ok_slab_ocean'.
     632
     633
     634!        Initialize surface history variable.
    671635!        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    672636         qsurf_hist(:,:)=qsurf(:,:)
    673637
    674 !        initialise variable for dynamical heating diagnostic
     638!        Initialize variable for dynamical heating diagnostic.
    675639!        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    676640         ztprevious(:,:)=pt(:,:)
    … …  
    688652         endif
    689653
    690          if(meanOLR)then
    691             ! to record global radiative balance
    692             call system('rm -f rad_bal.out')
    693             ! to record global mean/max/min temperatures
    694             call system('rm -f tem_bal.out')
    695             ! to record global hydrological balance
    696             call system('rm -f h2o_bal.out')
    697          endif
    698 
    699          watertest=.false.
    700          if(water)then
    701             ! initialise variables for water cycle
    702 
     654         if(meanOLR)then         
     655            call system('rm -f rad_bal.out') ! to record global radiative balance.           
     656            call system('rm -f tem_bal.out') ! to record global mean/max/min temperatures.           
     657            call system('rm -f h2o_bal.out') ! to record global hydrological balance.
     658         endif
     659
     660
     661         watertest=.false.       
     662         if(water)then ! Initialize variables for water cycle
     663           
    703664            if(enertest)then
    704665               watertest = .true.
    … …  
    710671
    711672         endif
     673         
    712674         call su_watercycle ! even if we don't have a water cycle, we might
    713675                            ! need epsi for the wvp definitions in callcorrk.F
    714676
    715          if (ngrid.ne.1) then ! no need to create a restart file in 1d
    716            call physdem0("restartfi.nc",long,lati,nsoilmx,ngrid,nlayer,nq, &
    717                          ptimestep,pday+nday,time_phys,area, &
    718                          albedodat,inertiedat,zmea,zstd,zsig,zgam,zthe)
     677         if (ngrid.ne.1) then ! Note : no need to create a restart file in 1d.
     678            call physdem0("restartfi.nc",long,lati,nsoilmx,ngrid,nlayer,nq, &
     679                          ptimestep,pday+nday,time_phys,area, &
     680                          albedodat,inertiedat,zmea,zstd,zsig,zgam,zthe)
    719681         endif
    720682         
    721       endif        !  (end of "if firstcall")
    722 
    723 ! ---------------------------------------------------
    724 ! 1.2   Initializations done at every physical timestep:
    725 ! ---------------------------------------------------
    726 !     Initialize various variables
    727 !     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    728      
     683      endif ! end of 'firstcall'
     684
     685! ------------------------------------------------------
     686! I.2   Initializations done at every physical timestep:
     687! ------------------------------------------------------
     688
     689      ! Initialize various variables
     690      ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     691     
     692      if ( .not.nearco2cond ) then
     693         pdt(1:ngrid,1:nlayer) = 0.0
     694      endif     
     695      zdtsurf(1:ngrid)      = 0.0
     696      pdq(1:ngrid,1:nlayer,1:nq) = 0.0
     697      dqsurf(1:ngrid,1:nq)= 0.0
    729698      pdu(1:ngrid,1:nlayer) = 0.0
    730699      pdv(1:ngrid,1:nlayer) = 0.0
    731       if ( .not.nearco2cond ) then
    732          pdt(1:ngrid,1:nlayer) = 0.0
    733       endif
    734       pdq(1:ngrid,1:nlayer,1:nq) = 0.0
    735700      pdpsrf(1:ngrid)       = 0.0
    736       zflubid(1:ngrid)      = 0.0
    737       zdtsurf(1:ngrid)      = 0.0
    738       dqsurf(1:ngrid,1:nq)= 0.0
     701      zflubid(1:ngrid)      = 0.0     
    739702      flux_sens_lat(1:ngrid) = 0.0
    740703      taux(1:ngrid) = 0.0
    741704      tauy(1:ngrid) = 0.0
    742705
    743 
    744 
    745 
    746       zday=pday+ptime ! compute time, in sols (and fraction thereof)
    747 
    748 !     Compute Stellar Longitude (Ls), and orbital parameters
    749 !     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     706      zday=pday+ptime ! Compute time, in sols (and fraction thereof).
     707
     708      ! Compute Stellar Longitude (Ls), and orbital parameters.
     709      ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    750710      if (season) then
    751711         call stellarlong(zday,zls)
    … …  
    755715
    756716      call orbite(zls,dist_star,declin,right_ascen)
     717     
    757718      if (tlocked) then
    758719              zlss=Mod(-(2.*pi*(zday/year_day)*nres - right_ascen),2.*pi)
    … …  
    764725
    765726
    766 !    Compute variations of g with latitude (oblate case)
    767 !
    768         if (oblate .eqv. .false.) then
    769            glat(:) = g
    770         else if (flatten .eq. 0.0 .or. J2 .eq. 0.0 .or. Rmean .eq. 0.0 .or. MassPlanet .eq. 0.0) then
    771         print*,'I need values for flatten, J2, Rmean and MassPlanet to compute glat (else set oblate=.false.)'
    772 
    773         call abort
    774         else
    775 
    776         gmplanet = MassPlanet*grav*1e24
    777         do ig=1,ngrid
    778            glat(ig)= gmplanet/(Rmean**2) * (1.D0 + 0.75 *J2 - 2.0*flatten/3. + (2.*flatten - 15./4.* J2) * cos(2. * (pi/2. - lati(ig))))
    779         end do
    780         endif
    781 
    782 !!      write(*,*) 'lati, glat =',lati(1)*180./pi,glat(1), lati(ngrid/2)*180./pi, glat(ngrid/2)
    783 
    784 
    785 
    786 !     Compute geopotential between layers
    787 !     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    788 
     727      ! Compute variations of g with latitude (oblate case).
     728      ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     729      if (oblate .eqv. .false.) then     
     730         glat(:) = g         
     731      else if (flatten .eq. 0.0 .or. J2 .eq. 0.0 .or. Rmean .eq. 0.0 .or. MassPlanet .eq. 0.0) then     
     732         print*,'I need values for flatten, J2, Rmean and MassPlanet to compute glat (else set oblate=.false.)'
     733         call abort
     734      else
     735         gmplanet = MassPlanet*grav*1e24
     736         do ig=1,ngrid
     737            glat(ig)= gmplanet/(Rmean**2) * (1.D0 + 0.75 *J2 - 2.0*flatten/3. + (2.*flatten - 15./4.* J2) * cos(2. * (pi/2. - lati(ig))))
     738         end do
     739      endif
     740
     741      ! Compute geopotential between layers.
     742      ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    789743      zzlay(1:ngrid,1:nlayer)=pphi(1:ngrid,1:nlayer)
    790744      do l=1,nlayer         
    791       zzlay(1:ngrid,l)= zzlay(1:ngrid,l)/glat(1:ngrid)
     745         zzlay(1:ngrid,l)= zzlay(1:ngrid,l)/glat(1:ngrid)
    792746      enddo
    793747
    794748      zzlev(1:ngrid,1)=0.
    795       zzlev(1:ngrid,nlayer+1)=1.e7    ! dummy top of last layer above 10000 km...
     749      zzlev(1:ngrid,nlayer+1)=1.e7 ! Dummy top of last layer above 10000 km...
    796750
    797751      do l=2,nlayer
    … …  
    801755            zzlev(ig,l)=(z1*zzlay(ig,l-1)+z2*zzlay(ig,l))/(z1+z2)
    802756         enddo
    803       enddo
    804 !     Potential temperature calculation may not be the same in physiq and dynamic...
    805 
    806 !     Compute potential temperature
    807 !     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    808 
     757      enddo     
     758
     759      ! Compute potential temperature
     760      ! Note : Potential temperature calculation may not be the same in physiq and dynamic...
     761      ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    809762      do l=1,nlayer         
    810763         do ig=1,ngrid
    … …  
    818771     ! Compute vertical velocity (m/s) from vertical mass flux
    819772     ! w = F / (rho*area) and rho = P/(r*T)
    820      ! but first linearly interpolate mass flux to mid-layers
    821      do l=1,nlayer-1
    822        pw(1:ngrid,l)=0.5*(flxw(1:ngrid,l)+flxw(1:ngrid,l+1))
    823      enddo
    824      pw(1:ngrid,nlayer)=0.5*flxw(1:ngrid,nlayer) ! since flxw(nlayer+1)=0
    825      do l=1,nlayer
    826        pw(1:ngrid,l)=(pw(1:ngrid,l)*r*pt(1:ngrid,l)) /  &
    827                      (pplay(1:ngrid,l)*area(1:ngrid))
    828      enddo
    829 
    830 !-----------------------------------------------------------------------
    831 !    2. Compute radiative tendencies
    832 !-----------------------------------------------------------------------
     773     ! But first linearly interpolate mass flux to mid-layers
     774      do l=1,nlayer-1
     775         pw(1:ngrid,l)=0.5*(flxw(1:ngrid,l)+flxw(1:ngrid,l+1))
     776      enddo
     777      pw(1:ngrid,nlayer)=0.5*flxw(1:ngrid,nlayer) ! since flxw(nlayer+1)=0
     778      do l=1,nlayer
     779         pw(1:ngrid,l)=(pw(1:ngrid,l)*r*pt(1:ngrid,l)) /  &
     780                       (pplay(1:ngrid,l)*area(1:ngrid))
     781      enddo
     782
     783!---------------------------------
     784! II. Compute radiative tendencies
     785!---------------------------------
    833786
    834787      if (callrad) then
    835788         if( mod(icount-1,iradia).eq.0.or.lastcall) then
    836789
    837 !          Compute local stellar zenith angles
    838 !          ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    839 
    840            if (tlocked) then
    841 ! JL14 corrects tidally resonant (and inclined) cases. nres=omega_rot/omega_orb
    842               ztim1=SIN(declin)
    843               ztim2=COS(declin)*COS(zlss)
    844               ztim3=COS(declin)*SIN(zlss)
    845 
    846               call stelang(ngrid,sinlon,coslon,sinlat,coslat,    &
    847                ztim1,ztim2,ztim3,mu0,fract, flatten)
    848 
    849            elseif (diurnal) then
    850               ztim1=SIN(declin)
    851               ztim2=COS(declin)*COS(2.*pi*(zday-.5))
    852               ztim3=-COS(declin)*SIN(2.*pi*(zday-.5))
    853 
    854               call stelang(ngrid,sinlon,coslon,sinlat,coslat,    &
    855                ztim1,ztim2,ztim3,mu0,fract, flatten)
    856            else if(diurnal .eqv. .false.) then
     790          ! Compute local stellar zenith angles
     791          ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     792            if (tlocked) then
     793            ! JL14 corrects tidally resonant (and inclined) cases. nres=omega_rot/omega_orb
     794               ztim1=SIN(declin)
     795               ztim2=COS(declin)*COS(zlss)
     796               ztim3=COS(declin)*SIN(zlss)
     797
     798               call stelang(ngrid,sinlon,coslon,sinlat,coslat,    &
     799                            ztim1,ztim2,ztim3,mu0,fract, flatten)
     800
     801            elseif (diurnal) then
     802               ztim1=SIN(declin)
     803               ztim2=COS(declin)*COS(2.*pi*(zday-.5))
     804               ztim3=-COS(declin)*SIN(2.*pi*(zday-.5))
     805
     806               call stelang(ngrid,sinlon,coslon,sinlat,coslat,    &
     807                            ztim1,ztim2,ztim3,mu0,fract, flatten)
     808            else if(diurnal .eqv. .false.) then
    857809
    858810               call mucorr(ngrid,declin,lati,mu0,fract,10000.,rad,flatten)
    859811               ! WARNING: this function appears not to work in 1D
    860812
    861            endif
     813            endif
    862814           
    863            !! Eclipse incoming sunlight (e.g. Saturn ring shadowing)
    864        
     815            ! Eclipse incoming sunlight (e.g. Saturn ring shadowing).       
    865816            if(rings_shadow) then
    866817                call call_rings(ngrid, ptime, pday, diurnal)
    … …  
    868819
    869820
    870            if (corrk) then
    871 
    872 !          a) Call correlated-k radiative transfer scheme
    873 !          ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    874 
    875               if(kastprof)then
    876                  print*,'kastprof should not = true here'
    877                  call abort
    878               endif
    879               if(water) then
     821            if (corrk) then
     822           
     823! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     824! II.a Call correlated-k radiative transfer scheme
     825! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     826               if(kastprof)then
     827                  print*,'kastprof should not = true here'
     828                  call abort
     829               endif
     830               if(water) then
    880831                  muvar(1:ngrid,1:nlayer)=mugaz/(1.e0+(1.e0/epsi-1.e0)*pq(1:ngrid,1:nlayer,igcm_h2o_vap))
    881832                  muvar(1:ngrid,nlayer+1)=mugaz/(1.e0+(1.e0/epsi-1.e0)*pq(1:ngrid,nlayer,igcm_h2o_vap))
    882833                  ! take into account water effect on mean molecular weight
    883               else
     834               else
    884835                  muvar(1:ngrid,1:nlayer+1)=mugaz
    885               endif
    886      
    887 !             standard callcorrk
    888 
    889               if(ok_slab_ocean) then
    890                 tsurf2(:)=tsurf(:)
    891                 do ig=1,ngrid
    892                   if (nint(rnat(ig))==0) then
    893                     tsurf(ig)=((1.-pctsrf_sic(ig))*tslab(ig,1)**4+pctsrf_sic(ig)*tsea_ice(ig)**4)**0.25
    894                   endif
    895                 enddo
    896               endif!(ok_slab_ocean)
     836               endif
     837     
     838
     839               if(ok_slab_ocean) then
     840                  tsurf2(:)=tsurf(:)
     841                  do ig=1,ngrid
     842                     if (nint(rnat(ig))==0) then
     843                        tsurf(ig)=((1.-pctsrf_sic(ig))*tslab(ig,1)**4+pctsrf_sic(ig)*tsea_ice(ig)**4)**0.25
     844                     endif
     845                  enddo
     846               endif !(ok_slab_ocean)
    897847               
    898 
    899                 clearsky=.false.
    900                 call callcorrk(ngrid,nlayer,pq,nq,qsurf,                   &
    901                       albedo,emis,mu0,pplev,pplay,pt,                    &
    902                       tsurf,fract,dist_star,aerosol,muvar,               &
    903                       zdtlw,zdtsw,fluxsurf_lw,fluxsurf_sw,fluxtop_lw,    &
    904                       fluxabs_sw,fluxtop_dn,OLR_nu,OSR_nu,               &
    905                       tau_col,cloudfrac,totcloudfrac,                    &
    906                       clearsky,firstcall,lastcall)     
    907 
    908 !             Option to call scheme once more for clear regions
    909                 if(CLFvarying)then
    910 
    911                  ! ---> PROBLEMS WITH ALLOCATED ARRAYS
    912                  ! (temporary solution in callcorrk: do not deallocate if CLFvarying ...)
    913                    clearsky=.true.
    914                    call callcorrk(ngrid,nlayer,pq,nq,qsurf,                  &
    915                       albedo,emis,mu0,pplev,pplay,pt,                      &
    916                       tsurf,fract,dist_star,aerosol,muvar,                 &
    917                       zdtlw1,zdtsw1,fluxsurf_lw1,fluxsurf_sw1,fluxtop_lw1, &
    918                       fluxabs_sw1,fluxtop_dn,OLR_nu1,OSR_nu1,              &
    919                       tau_col1,cloudfrac,totcloudfrac,                     &
    920                       clearsky,firstcall,lastcall)
    921                    clearsky = .false.  ! just in case.     
    922 
    923 
    924                  ! Sum the fluxes and heating rates from cloudy/clear cases
    925                    do ig=1,ngrid
    926                       tf=totcloudfrac(ig)
    927                       ntf=1.-tf
     848               ! standard callcorrk
     849               clearsky=.false.
     850               call callcorrk(ngrid,nlayer,pq,nq,qsurf,                          &
     851                              albedo,emis,mu0,pplev,pplay,pt,                    &
     852                              tsurf,fract,dist_star,aerosol,muvar,               &
     853                              zdtlw,zdtsw,fluxsurf_lw,fluxsurf_sw,fluxtop_lw,    &
     854                              fluxabs_sw,fluxtop_dn,OLR_nu,OSR_nu,               &
     855                              tau_col,cloudfrac,totcloudfrac,                    &
     856                              clearsky,firstcall,lastcall)     
     857
     858                ! Option to call scheme once more for clear regions
     859               if(CLFvarying)then
     860
     861                  ! ---> PROBLEMS WITH ALLOCATED ARRAYS : temporary solution in callcorrk: do not deallocate if CLFvarying ...
     862                  clearsky=.true.
     863                  call callcorrk(ngrid,nlayer,pq,nq,qsurf,                            &
     864                                 albedo,emis,mu0,pplev,pplay,pt,                      &
     865                                 tsurf,fract,dist_star,aerosol,muvar,                 &
     866                                 zdtlw1,zdtsw1,fluxsurf_lw1,fluxsurf_sw1,fluxtop_lw1, &
     867                                 fluxabs_sw1,fluxtop_dn,OLR_nu1,OSR_nu1,              &
     868                                 tau_col1,cloudfrac,totcloudfrac,                     &
     869                                 clearsky,firstcall,lastcall)
     870                  clearsky = .false.  ! just in case.     
     871
     872
     873                  ! Sum the fluxes and heating rates from cloudy/clear cases
     874                  do ig=1,ngrid
     875                     tf=totcloudfrac(ig)
     876                     ntf=1.-tf
    928877                   
    929                     fluxsurf_lw(ig) = ntf*fluxsurf_lw1(ig) + tf*fluxsurf_lw(ig)
    930                     fluxsurf_sw(ig) = ntf*fluxsurf_sw1(ig) + tf*fluxsurf_sw(ig)
    931                     fluxtop_lw(ig)  = ntf*fluxtop_lw1(ig)  + tf*fluxtop_lw(ig)
    932                     fluxabs_sw(ig)  = ntf*fluxabs_sw1(ig)  + tf*fluxabs_sw(ig)
    933                     tau_col(ig)     = ntf*tau_col1(ig)     + tf*tau_col(ig)
     878                     fluxsurf_lw(ig) = ntf*fluxsurf_lw1(ig) + tf*fluxsurf_lw(ig)
     879                     fluxsurf_sw(ig) = ntf*fluxsurf_sw1(ig) + tf*fluxsurf_sw(ig)
     880                     fluxtop_lw(ig)  = ntf*fluxtop_lw1(ig)  + tf*fluxtop_lw(ig)
     881                     fluxabs_sw(ig)  = ntf*fluxabs_sw1(ig)  + tf*fluxabs_sw(ig)
     882                     tau_col(ig)     = ntf*tau_col1(ig)     + tf*tau_col(ig)
    934883                   
    935                     zdtlw(ig,1:nlayer) = ntf*zdtlw1(ig,1:nlayer) + tf*zdtlw(ig,1:nlayer)
    936                     zdtsw(ig,1:nlayer) = ntf*zdtsw1(ig,1:nlayer) + tf*zdtsw(ig,1:nlayer)
    937 
    938                     OSR_nu(ig,1:L_NSPECTV) = ntf*OSR_nu1(ig,1:L_NSPECTV) + tf*OSR_nu(ig,1:L_NSPECTV)                   
    939                     OLR_nu(ig,1:L_NSPECTI) = ntf*OLR_nu1(ig,1:L_NSPECTI) + tf*OLR_nu(ig,1:L_NSPECTI)                   
    940 
     884                     zdtlw(ig,1:nlayer) = ntf*zdtlw1(ig,1:nlayer) + tf*zdtlw(ig,1:nlayer)
     885                     zdtsw(ig,1:nlayer) = ntf*zdtsw1(ig,1:nlayer) + tf*zdtsw(ig,1:nlayer)
     886
     887                     OSR_nu(ig,1:L_NSPECTV) = ntf*OSR_nu1(ig,1:L_NSPECTV) + tf*OSR_nu(ig,1:L_NSPECTV)                 
     888                     OLR_nu(ig,1:L_NSPECTI) = ntf*OLR_nu1(ig,1:L_NSPECTI) + tf*OLR_nu(ig,1:L_NSPECTI)                 
    941889                 enddo
    942890
    943               endif !CLFvarying
    944 
    945               if(ok_slab_ocean) then
    946                 tsurf(:)=tsurf2(:)
    947               endif!(ok_slab_ocean)
    948 
    949 
    950 !             Radiative flux from the sky absorbed by the surface (W.m-2)
    951               GSR=0.0
    952               fluxrad_sky(1:ngrid)=emis(1:ngrid)*fluxsurf_lw(1:ngrid)+fluxsurf_sw(1:ngrid)*(1.-albedo(1:ngrid))
    953 
    954               !if(noradsurf)then ! no lower surface; SW flux just disappears
    955               !   GSR = SUM(fluxsurf_sw(1:ngrid)*area(1:ngrid))/totarea
    956               !   fluxrad_sky(1:ngrid)=emis(1:ngrid)*fluxsurf_lw(1:ngrid)
    957               !   print*,'SW lost in deep atmosphere = ',GSR,' W m^-2'
    958               !endif
    959 
    960 !             Net atmospheric radiative heating rate (K.s-1)
    961               dtrad(1:ngrid,1:nlayer)=zdtsw(1:ngrid,1:nlayer)+zdtlw(1:ngrid,1:nlayer)
    962 
    963            elseif(newtonian)then
    964 
    965 !          b) Call Newtonian cooling scheme
    966 !          ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    967               call newtrelax(ngrid,nlayer,mu0,sinlat,zpopsk,pt,pplay,pplev,dtrad,firstcall)
    968 
    969               zdtsurf(1:ngrid) = +(pt(1:ngrid,1)-tsurf(1:ngrid))/ptimestep
    970               ! e.g. surface becomes proxy for 1st atmospheric layer ?
    971 
    972            else
    973 
    974 !          c) Atmosphere has no radiative effect
    975 !          ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    976               fluxtop_dn(1:ngrid)  = fract(1:ngrid)*mu0(1:ngrid)*Fat1AU/dist_star**2
    977               if(ngrid.eq.1)then ! / by 4 globally in 1D case...
    978                  fluxtop_dn(1)  = fract(1)*Fat1AU/dist_star**2/2.0
    979               endif
    980               fluxsurf_sw(1:ngrid) = fluxtop_dn(1:ngrid)
    981               fluxrad_sky(1:ngrid) = fluxtop_dn(1:ngrid)*(1.-albedo(1:ngrid))
    982               fluxtop_lw(1:ngrid)  = emis(1:ngrid)*sigma*tsurf(1:ngrid)**4
    983               ! radiation skips the atmosphere entirely
    984 
    985 
    986               dtrad(1:ngrid,1:nlayer)=0.0
    987               ! hence no atmospheric radiative heating
    988 
    989            endif                ! if corrk
    990 
    991         endif ! of if(mod(icount-1,iradia).eq.0)
     891               endif ! end of CLFvarying.
     892
     893               if(ok_slab_ocean) then
     894                  tsurf(:)=tsurf2(:)
     895               endif
     896
     897
     898               ! Radiative flux from the sky absorbed by the surface (W.m-2)
     899               GSR=0.0
     900               fluxrad_sky(1:ngrid)=emis(1:ngrid)*fluxsurf_lw(1:ngrid)+fluxsurf_sw(1:ngrid)*(1.-albedo(1:ngrid))
     901
     902                            !if(noradsurf)then ! no lower surface; SW flux just disappears
     903                            !   GSR = SUM(fluxsurf_sw(1:ngrid)*area(1:ngrid))/totarea
     904                            !   fluxrad_sky(1:ngrid)=emis(1:ngrid)*fluxsurf_lw(1:ngrid)
     905                            !   print*,'SW lost in deep atmosphere = ',GSR,' W m^-2'
     906                            !endif
     907
     908               ! Net atmospheric radiative heating rate (K.s-1)
     909               dtrad(1:ngrid,1:nlayer)=zdtsw(1:ngrid,1:nlayer)+zdtlw(1:ngrid,1:nlayer)
     910
     911            elseif(newtonian)then
     912
     913               ! II.b Call Newtonian cooling scheme
     914               ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     915               call newtrelax(ngrid,nlayer,mu0,sinlat,zpopsk,pt,pplay,pplev,dtrad,firstcall)
     916
     917               zdtsurf(1:ngrid) = +(pt(1:ngrid,1)-tsurf(1:ngrid))/ptimestep
     918               ! e.g. surface becomes proxy for 1st atmospheric layer ?
     919
     920            else
     921           
     922! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     923! II.c Atmosphere has no radiative effect
     924! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     925               fluxtop_dn(1:ngrid)  = fract(1:ngrid)*mu0(1:ngrid)*Fat1AU/dist_star**2
     926               if(ngrid.eq.1)then ! / by 4 globally in 1D case...
     927                  fluxtop_dn(1)  = fract(1)*Fat1AU/dist_star**2/2.0
     928               endif
     929               fluxsurf_sw(1:ngrid) = fluxtop_dn(1:ngrid)
     930               fluxrad_sky(1:ngrid) = fluxtop_dn(1:ngrid)*(1.-albedo(1:ngrid))
     931               fluxtop_lw(1:ngrid)  = emis(1:ngrid)*sigma*tsurf(1:ngrid)**4
     932
     933               dtrad(1:ngrid,1:nlayer)=0.0 ! no atmospheric radiative heating
     934
     935            endif ! end of corrk
     936
     937         endif ! of if(mod(icount-1,iradia).eq.0)
    992938       
    993939
    994 !    Transformation of the radiative tendencies
    995 !    ------------------------------------------
    996 
    997         zplanck(1:ngrid)=tsurf(1:ngrid)*tsurf(1:ngrid)
    998         zplanck(1:ngrid)=emis(1:ngrid)*sigma*zplanck(1:ngrid)*zplanck(1:ngrid)
    999         fluxrad(1:ngrid)=fluxrad_sky(1:ngrid)-zplanck(1:ngrid)
    1000         pdt(1:ngrid,1:nlayer)=pdt(1:ngrid,1:nlayer)+dtrad(1:ngrid,1:nlayer)
    1001 
    1002 !-------------------------
    1003 ! test energy conservation
     940         ! Transformation of the radiative tendencies
     941         ! ------------------------------------------
     942         zplanck(1:ngrid)=tsurf(1:ngrid)*tsurf(1:ngrid)
     943         zplanck(1:ngrid)=emis(1:ngrid)*sigma*zplanck(1:ngrid)*zplanck(1:ngrid)
     944         fluxrad(1:ngrid)=fluxrad_sky(1:ngrid)-zplanck(1:ngrid)
     945         pdt(1:ngrid,1:nlayer)=pdt(1:ngrid,1:nlayer)+dtrad(1:ngrid,1:nlayer)
     946         
     947         ! Test of energy conservation
     948         !----------------------------
    1004949         if(enertest)then
    1005950            call planetwide_sumval(cpp*massarea(:,:)*zdtsw(:,:)/totarea_planet,dEtotSW)
    … …  
    1010955            dEzRadlw(:,:)=cpp*mass(:,:)*zdtlw(:,:)
    1011956            if (is_master) then
    1012                 print*,'---------------------------------------------------------------'
    1013                 print*,'In corrk SW atmospheric heating       =',dEtotSW,' W m-2'
    1014                 print*,'In corrk LW atmospheric heating       =',dEtotLW,' W m-2'
    1015                 print*,'atmospheric net rad heating (SW+LW)   =',dEtotLW+dEtotSW,' W m-2'
    1016                 print*,'In corrk SW surface heating           =',dEtotsSW,' W m-2'
    1017                 print*,'In corrk LW surface heating           =',dEtotsLW,' W m-2'
    1018                 print*,'surface net rad heating (SW+LW)       =',dEtotsLW+dEtotsSW,' W m-2'
     957               print*,'---------------------------------------------------------------'
     958               print*,'In corrk SW atmospheric heating       =',dEtotSW,' W m-2'
     959               print*,'In corrk LW atmospheric heating       =',dEtotLW,' W m-2'
     960               print*,'atmospheric net rad heating (SW+LW)   =',dEtotLW+dEtotSW,' W m-2'
     961               print*,'In corrk SW surface heating           =',dEtotsSW,' W m-2'
     962               print*,'In corrk LW surface heating           =',dEtotsLW,' W m-2'
     963               print*,'surface net rad heating (SW+LW)       =',dEtotsLW+dEtotsSW,' W m-2'
    1019964            endif
    1020          endif
    1021 !-------------------------
     965         endif ! end of 'enertest'
    1022966
    1023967      endif ! of if (callrad)
    1024968
    1025 !-----------------------------------------------------------------------
    1026 !    4. Vertical diffusion (turbulent mixing):
    1027 !    -----------------------------------------
     969
     970
     971!  --------------------------------------------
     972!  III. Vertical diffusion (turbulent mixing) :
     973!  --------------------------------------------
    1028974
    1029975      if (calldifv) then
    … …  
    1031977         zflubid(1:ngrid)=fluxrad(1:ngrid)+fluxgrd(1:ngrid)
    1032978
    1033          zdum1(1:ngrid,1:nlayer)=0.0
    1034          zdum2(1:ngrid,1:nlayer)=0.0
    1035 
    1036 
    1037 !JL12 the following if test is temporarily there to allow us to compare the old vdifc with turbdiff     
     979         ! JL12 the following if test is temporarily there to allow us to compare the old vdifc with turbdiff.
    1038980         if (UseTurbDiff) then
    1039981         
    1040            call turbdiff(ngrid,nlayer,nq,rnat,       &
    1041               ptimestep,capcal,lwrite,               &
    1042               pplay,pplev,zzlay,zzlev,z0,            &
    1043               pu,pv,pt,zpopsk,pq,tsurf,emis,qsurf,   &
    1044               zdum1,zdum2,pdt,pdq,zflubid,           &
    1045               zdudif,zdvdif,zdtdif,zdtsdif,          &
    1046               sensibFlux,q2,zdqdif,zdqevap,zdqsdif,  &
    1047               taux,tauy,lastcall)
     982            call turbdiff(ngrid,nlayer,nq,rnat,                  &
     983                          ptimestep,capcal,lwrite,               &
     984                          pplay,pplev,zzlay,zzlev,z0,            &
     985                          pu,pv,pt,zpopsk,pq,tsurf,emis,qsurf,   &
     986                          pdt,pdq,zflubid,                       &
     987                          zdudif,zdvdif,zdtdif,zdtsdif,          &
     988                          sensibFlux,q2,zdqdif,zdqevap,zdqsdif,  &
     989                          taux,tauy,lastcall)
    1048990
    1049991         else
    1050992         
    1051            zdh(1:ngrid,1:nlayer)=pdt(1:ngrid,1:nlayer)/zpopsk(1:ngrid,1:nlayer)
     993            zdh(1:ngrid,1:nlayer)=pdt(1:ngrid,1:nlayer)/zpopsk(1:ngrid,1:nlayer)
    1052994 
    1053            call vdifc(ngrid,nlayer,nq,rnat,zpopsk,   &
    1054               ptimestep,capcal,lwrite,               &
    1055               pplay,pplev,zzlay,zzlev,z0,            &
    1056               pu,pv,zh,pq,tsurf,emis,qsurf,          &
    1057               zdum1,zdum2,zdh,pdq,zflubid,           &
    1058               zdudif,zdvdif,zdhdif,zdtsdif,          &
    1059               sensibFlux,q2,zdqdif,zdqsdif,          &
    1060               taux,tauy,lastcall)
    1061 
    1062            zdtdif(1:ngrid,1:nlayer)=zdhdif(1:ngrid,1:nlayer)*zpopsk(1:ngrid,1:nlayer) ! for diagnostic only
    1063            zdqevap(1:ngrid,1:nlayer)=0.
    1064 
    1065          end if !turbdiff
     995            call vdifc(ngrid,nlayer,nq,rnat,zpopsk,           &
     996                       ptimestep,capcal,lwrite,               &
     997                       pplay,pplev,zzlay,zzlev,z0,            &
     998                       pu,pv,zh,pq,tsurf,emis,qsurf,          &
     999                       zdh,pdq,zflubid,                       &
     1000                       zdudif,zdvdif,zdhdif,zdtsdif,          &
     1001                       sensibFlux,q2,zdqdif,zdqsdif,          &
     1002                       taux,tauy,lastcall)
     1003
     1004            zdtdif(1:ngrid,1:nlayer)=zdhdif(1:ngrid,1:nlayer)*zpopsk(1:ngrid,1:nlayer) ! for diagnostic only
     1005            zdqevap(1:ngrid,1:nlayer)=0.
     1006
     1007         end if !end of 'UseTurbDiff'
     1008
    10661009
    10671010         pdv(1:ngrid,1:nlayer)=pdv(1:ngrid,1:nlayer)+zdvdif(1:ngrid,1:nlayer)
    … …  
    10701013         zdtsurf(1:ngrid)=zdtsurf(1:ngrid)+zdtsdif(1:ngrid)
    10711014
     1015
    10721016         if(ok_slab_ocean)then
    10731017            flux_sens_lat(1:ngrid)=(zdtsdif(1:ngrid)*capcal(1:ngrid)-fluxrad(1:ngrid))
    10741018         endif
    1075 
    10761019
    10771020
    … …  
    10811024         end if ! of if (tracer)
    10821025
     1026
     1027         ! test energy conservation
    10831028         !-------------------------
    1084          ! test energy conservation
    10851029         if(enertest)then
     1030         
    10861031            dEzdiff(:,:)=cpp*mass(:,:)*zdtdif(:,:)
    10871032            do ig = 1, ngrid
    … …  
    10891034               dEzdiff(ig,1)= dEzdiff(ig,1)+ sensibFlux(ig)! subtract flux to the ground
    10901035            enddo
     1036           
    10911037            call planetwide_sumval(dEdiff(:)*area(:)/totarea_planet,dEtot)
    10921038            dEdiffs(:)=capcal(:)*zdtsdif(:)-zflubid(:)-sensibFlux(:)
    10931039            call planetwide_sumval(dEdiffs(:)*area(:)/totarea_planet,dEtots)
    10941040            call planetwide_sumval(sensibFlux(:)*area(:)/totarea_planet,AtmToSurf_TurbFlux)
     1041           
    10951042            if (is_master) then
    1096              if (UseTurbDiff) then
    1097                 print*,'In TurbDiff sensible flux (atm=>surf) =',AtmToSurf_TurbFlux,' W m-2'
    1098                 print*,'In TurbDiff non-cons atm nrj change   =',dEtot,' W m-2'
    1099                 print*,'In TurbDiff (correc rad+latent heat) surf nrj change =',dEtots,' W m-2'
    1100              else
    1101                 print*,'In vdifc sensible flux (atm=>surf)    =',AtmToSurf_TurbFlux,' W m-2'
    1102                 print*,'In vdifc non-cons atm nrj change      =',dEtot,' W m-2'
    1103                 print*,'In vdifc (correc rad+latent heat) surf nrj change =',dEtots,' W m-2'
    1104              end if
    1105             endif ! of if (is_master)
    1106 ! JL12 note that the black body radiative flux emitted by the surface has been updated by the implicit scheme
    1107 !    but not given back elsewhere
    1108          endif
    1109          !-------------------------
    1110 
    1111          !-------------------------
    1112          ! test water conservation
     1043           
     1044               if (UseTurbDiff) then
     1045                  print*,'In TurbDiff sensible flux (atm=>surf) =',AtmToSurf_TurbFlux,' W m-2'
     1046                  print*,'In TurbDiff non-cons atm nrj change   =',dEtot,' W m-2'
     1047                  print*,'In TurbDiff (correc rad+latent heat) surf nrj change =',dEtots,' W m-2'
     1048               else
     1049                  print*,'In vdifc sensible flux (atm=>surf)    =',AtmToSurf_TurbFlux,' W m-2'
     1050                  print*,'In vdifc non-cons atm nrj change      =',dEtot,' W m-2'
     1051                  print*,'In vdifc (correc rad+latent heat) surf nrj change =',dEtots,' W m-2'
     1052               end if
     1053            endif ! end of 'is_master'
     1054           
     1055         ! JL12 : note that the black body radiative flux emitted by the surface has been updated by the implicit scheme but not given back elsewhere.
     1056         endif ! end of 'enertest'
     1057
     1058
     1059         ! Test water conservation.
    11131060         if(watertest.and.water)then
     1061         
    11141062            call planetwide_sumval(massarea(:,:)*zdqdif(:,:,igcm_h2o_vap)*ptimestep/totarea_planet,dWtot_tmp)
    11151063            call planetwide_sumval(zdqsdif(:,igcm_h2o_vap)*area(:)*ptimestep/totarea_planet,dWtots_tmp)
    11161064            do ig = 1, ngrid
    11171065               vdifcncons(ig)=SUM(mass(ig,:)*zdqdif(ig,:,igcm_h2o_vap))
    1118             Enddo
     1066            enddo
    11191067            call planetwide_sumval(massarea(:,:)*zdqdif(:,:,igcm_h2o_ice)*ptimestep/totarea_planet,dWtot)
    11201068            call planetwide_sumval(zdqsdif(:,igcm_h2o_ice)*area(:)*ptimestep/totarea_planet,dWtots)
    … …  
    11231071            do ig = 1, ngrid
    11241072               vdifcncons(ig)=vdifcncons(ig) + SUM(mass(ig,:)*zdqdif(ig,:,igcm_h2o_ice))
    1125             Enddo           
     1073            enddo           
    11261074            call planetwide_maxval(vdifcncons(:),nconsMAX)
    11271075
    11281076            if (is_master) then
    1129                 print*,'---------------------------------------------------------------'
    1130                 print*,'In difv atmospheric water change        =',dWtot,' kg m-2'
    1131                 print*,'In difv surface water change            =',dWtots,' kg m-2'
    1132                 print*,'In difv non-cons factor                 =',dWtot+dWtots,' kg m-2'
    1133                 print*,'In difv MAX non-cons factor             =',nconsMAX,' kg m-2 s-1'
     1077               print*,'---------------------------------------------------------------'
     1078               print*,'In difv atmospheric water change        =',dWtot,' kg m-2'
     1079               print*,'In difv surface water change            =',dWtots,' kg m-2'
     1080               print*,'In difv non-cons factor                 =',dWtot+dWtots,' kg m-2'
     1081               print*,'In difv MAX non-cons factor             =',nconsMAX,' kg m-2 s-1'
    11341082            endif
    11351083
    1136          endif
     1084         endif ! end of 'watertest'
    11371085         !-------------------------
    11381086
    1139       else    
     1087      else ! calldifv
    11401088
    11411089         if(.not.newtonian)then
    … …  
    11451093         endif
    11461094
    1147       endif ! of if (calldifv)
    1148 
    1149 
    1150 !-----------------------------------------------------------------------
    1151 !   5. Dry convective adjustment:
    1152 !   -----------------------------
     1095      endif ! end of 'calldifv'
     1096
     1097
     1098!----------------------------------
     1099!   IV. Dry convective adjustment :
     1100!----------------------------------
    11531101
    11541102      if(calladj) then
    … …  
    11601108
    11611109
    1162          call convadj(ngrid,nlayer,nq,ptimestep,    &
    1163               pplay,pplev,zpopsk,                   &
    1164               pu,pv,zh,pq,                          &
    1165               pdu,pdv,zdh,pdq,                      &
    1166               zduadj,zdvadj,zdhadj,                 &
    1167               zdqadj)
     1110         call convadj(ngrid,nlayer,nq,ptimestep,            &
     1111                      pplay,pplev,zpopsk,                   &
     1112                      pu,pv,zh,pq,                          &
     1113                      pdu,pdv,zdh,pdq,                      &
     1114                      zduadj,zdvadj,zdhadj,                 &
     1115                      zdqadj)
    11681116
    11691117         pdu(1:ngrid,1:nlayer) = pdu(1:ngrid,1:nlayer) + zduadj(1:ngrid,1:nlayer)
    … …  
    11761124         end if
    11771125
    1178          !-------------------------
    1179          ! test energy conservation
     1126         ! Test energy conservation
    11801127         if(enertest)then
    11811128            call planetwide_sumval(cpp*massarea(:,:)*zdtadj(:,:)/totarea_planet,dEtot)
    11821129            if (is_master) print*,'In convadj atmospheric energy change  =',dEtot,' W m-2'
    11831130         endif
    1184          !-------------------------
    1185 
    1186          !-------------------------
    1187          ! test water conservation
     1131
     1132         ! Test water conservation
    11881133         if(watertest)then
    11891134            call planetwide_sumval(massarea(:,:)*zdqadj(:,:,igcm_h2o_vap)*ptimestep/totarea_planet,dWtot_tmp)
    11901135            do ig = 1, ngrid
    11911136               cadjncons(ig)=SUM(mass(ig,:)*zdqadj(ig,:,igcm_h2o_vap))
    1192             Enddo
     1137            enddo
    11931138            call planetwide_sumval(massarea(:,:)*zdqadj(:,:,igcm_h2o_ice)*ptimestep/totarea_planet,dWtot)
    11941139            dWtot = dWtot + dWtot_tmp
    11951140            do ig = 1, ngrid
    11961141               cadjncons(ig)=cadjncons(ig) + SUM(mass(ig,:)*zdqadj(ig,:,igcm_h2o_ice))
    1197             Enddo           
     1142            enddo           
    11981143            call planetwide_maxval(cadjncons(:),nconsMAX)
    11991144
    12001145            if (is_master) then
    1201                 print*,'In convadj atmospheric water change     =',dWtot,' kg m-2'
    1202                 print*,'In convadj MAX non-cons factor          =',nconsMAX,' kg m-2 s-1'
     1146               print*,'In convadj atmospheric water change     =',dWtot,' kg m-2'
     1147               print*,'In convadj MAX non-cons factor          =',nconsMAX,' kg m-2 s-1'
    12031148            endif
    1204          endif
    1205          !-------------------------
     1149            
     1150         endif ! end of 'watertest'
    12061151         
    1207       endif ! of if(calladj)
    1208 
    1209 !-----------------------------------------------------------------------
    1210 !   6. Carbon dioxide condensation-sublimation:
    1211 !   -------------------------------------------
     1152      endif ! end of 'calladj'
     1153
     1154!-----------------------------------------------
     1155!   V. Carbon dioxide condensation-sublimation :
     1156!-----------------------------------------------
    12121157
    12131158      if (co2cond) then
     1159     
    12141160         if (.not.tracer) then
    12151161            print*,'We need a CO2 ice tracer to condense CO2'
    12161162            call abort
    12171163         endif
    1218          call condense_cloud(ngrid,nlayer,nq,ptimestep,   &
    1219               capcal,pplay,pplev,tsurf,pt,                  &
    1220               pphi,pdt,pdu,pdv,zdtsurf,pu,pv,pq,pdq,        &
    1221               qsurf(1,igcm_co2_ice),albedo,emis,            &
    1222               zdtc,zdtsurfc,pdpsrf,zduc,zdvc,               &
    1223               zdqc)
     1164         call condense_co2(ngrid,nlayer,nq,ptimestep,                    &
     1165                           capcal,pplay,pplev,tsurf,pt,                  &
     1166                           pphi,pdt,pdu,pdv,zdtsurf,pu,pv,pq,pdq,        &
     1167                           qsurf(1,igcm_co2_ice),albedo,emis,            &
     1168                           zdtc,zdtsurfc,pdpsrf,zduc,zdvc,               &
     1169                           zdqc)
    12241170
    12251171         pdt(1:ngrid,1:nlayer)=pdt(1:ngrid,1:nlayer)+zdtc(1:ngrid,1:nlayer)
    … …  
    12291175
    12301176         pdq(1:ngrid,1:nlayer,1:nq)=pdq(1:ngrid,1:nlayer,1:nq)+ zdqc(1:ngrid,1:nlayer,1:nq)
    1231          ! Note: we do not add surface co2ice tendency
    1232          ! because qsurf(:,igcm_co2_ice) is updated in condens_co2cloud
    1233 
    1234          !-------------------------
     1177
    12351178         ! test energy conservation
    12361179         if(enertest)then
    … …  
    12381181            call planetwide_sumval(capcal(:)*zdtsurfc(:)*area(:)/totarea_planet,dEtots)
    12391182            if (is_master) then
    1240                 print*,'In co2cloud atmospheric energy change   =',dEtot,' W m-2'
    1241                 print*,'In co2cloud surface energy change       =',dEtots,' W m-2'
     1183               print*,'In co2cloud atmospheric energy change   =',dEtot,' W m-2'
     1184               print*,'In co2cloud surface energy change       =',dEtots,' W m-2'
    12421185            endif
    12431186         endif
    1244          !-------------------------
    1245 
    1246       endif  ! of if (co2cond)
    1247 
    1248 
    1249 !-----------------------------------------------------------------------
    1250 !   7. Specific parameterizations for tracers
    1251 !   -----------------------------------------
    1252 
    1253       if (tracer) then
    1254 
    1255 !   7a. Water and ice
    1256 !     ---------------
     1187
     1188      endif  ! end of 'co2cond'
     1189
     1190
     1191!---------------------------------------------
     1192!   VI. Specific parameterizations for tracers
     1193!---------------------------------------------
     1194
     1195      if (tracer) then
     1196     
     1197  ! ---------------------
     1198  !   VI.1. Water and ice
     1199  ! ---------------------
    12571200         if (water) then
    12581201
    1259 !        ----------------------------------------
    1260 !        Water ice condensation in the atmosphere
    1261 !        ----------------------------------------
     1202            ! Water ice condensation in the atmosphere
    12621203            if(watercond.and.(RLVTT.gt.1.e-8))then
    1263 
    1264 !             ----------------
    1265 !             Moist convection
    1266 !             ----------------
    12671204
    12681205               dqmoist(1:ngrid,1:nlayer,1:nq)=0.
    12691206               dtmoist(1:ngrid,1:nlayer)=0.
    1270 
     1207               
     1208               ! Moist Convective Adjustment.
     1209               ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    12711210               call moistadj(ngrid,nlayer,nq,pt,pq,pdq,pplev,pplay,dtmoist,dqmoist,ptimestep,rneb_man)
    12721211
    1273                pdq(1:ngrid,1:nlayer,igcm_h2o_vap) = pdq(1:ngrid,1:nlayer,igcm_h2o_vap)   &
    1274                           +dqmoist(1:ngrid,1:nlayer,igcm_h2o_vap)
    1275                pdq(1:ngrid,1:nlayer,igcm_h2o_ice) =pdq(1:ngrid,1:nlayer,igcm_h2o_ice)     &
    1276                           +dqmoist(1:ngrid,1:nlayer,igcm_h2o_ice)
     1212               pdq(1:ngrid,1:nlayer,igcm_h2o_vap) = pdq(1:ngrid,1:nlayer,igcm_h2o_vap)     &
     1213                                                  + dqmoist(1:ngrid,1:nlayer,igcm_h2o_vap)
     1214               pdq(1:ngrid,1:nlayer,igcm_h2o_ice) = pdq(1:ngrid,1:nlayer,igcm_h2o_ice)     &
     1215                                                  + dqmoist(1:ngrid,1:nlayer,igcm_h2o_ice)
    12771216               pdt(1:ngrid,1:nlayer) = pdt(1:ngrid,1:nlayer)+dtmoist(1:ngrid,1:nlayer)
    12781217
    1279                !-------------------------
    1280                ! test energy conservation
     1218               ! Test energy conservation.
    12811219               if(enertest)then
     1220               
    12821221                  call planetwide_sumval(cpp*massarea(:,:)*dtmoist(:,:)/totarea_planet,dEtot)
    12831222                  call planetwide_maxval(dtmoist(:,:),dtmoist_max)
    … …  
    12871226                     madjdE(ig) = cpp*SUM(mass(:,:)*dtmoist(:,:))
    12881227                  enddo
     1228                 
    12891229                  if (is_master) then
    1290                         print*,'In moistadj atmospheric energy change   =',dEtot,' W m-2'
    1291                         print*,'In moistadj MAX atmospheric energy change   =',dtmoist_max*ptimestep,'K/step'
    1292                         print*,'In moistadj MIN atmospheric energy change   =',dtmoist_min*ptimestep,'K/step'
     1230                     print*,'In moistadj atmospheric energy change   =',dEtot,' W m-2'
     1231                     print*,'In moistadj MAX atmospheric energy change   =',dtmoist_max*ptimestep,'K/step'
     1232                     print*,'In moistadj MIN atmospheric energy change   =',dtmoist_min*ptimestep,'K/step'
    12931233                  endif
    12941234                 
    1295                 ! test energy conservation
    12961235                  call planetwide_sumval(massarea(:,:)*dqmoist(:,:,igcm_h2o_vap)*ptimestep/totarea_planet+      &
    1297                         massarea(:,:)*dqmoist(:,:,igcm_h2o_ice)*ptimestep/totarea_planet,dWtot)
     1236                                         massarea(:,:)*dqmoist(:,:,igcm_h2o_ice)*ptimestep/totarea_planet,dWtot)
    12981237                  if (is_master) print*,'In moistadj atmospheric water change    =',dWtot,' kg m-2'
    1299                endif
    1300                !-------------------------
     1238                  
     1239               endif ! end of 'enertest'
    13011240               
    13021241
    1303 !        --------------------------------
    1304 !        Large scale condensation/evaporation
    1305 !        --------------------------------
     1242               ! Large scale condensation/evaporation.
     1243               ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    13061244               call largescale(ngrid,nlayer,nq,ptimestep,pplev,pplay,pt,pq,pdt,pdq,dtlscale,dqvaplscale,dqcldlscale,rneb_lsc)
    13071245
    … …  
    13101248               pdq(1:ngrid,1:nlayer,igcm_h2o_ice) = pdq(1:ngrid,1:nlayer,igcm_h2o_ice)+dqcldlscale(1:ngrid,1:nlayer)
    13111249
    1312                !-------------------------
    1313                ! test energy conservation
     1250               ! Test energy conservation.
    13141251               if(enertest)then
    13151252                  lscaledEz(:,:) = cpp*mass(:,:)*dtlscale(:,:)
    … …  
    13181255                  enddo
    13191256                  call planetwide_sumval(cpp*massarea(:,:)*dtlscale(:,:)/totarea_planet,dEtot)
    1320 !                 if(isnan(dEtot)) then ! NB: isnan() is not a standard function...
    1321 !                    print*,'Nan in largescale, abort'
    1322 !                     STOP
    1323 !                 endif
     1257
    13241258                  if (is_master) print*,'In largescale atmospheric energy change =',dEtot,' W m-2'
    13251259
    1326                ! test water conservation
     1260                  ! Test water conservation.
    13271261                  call planetwide_sumval(massarea(:,:)*dqvaplscale(:,:)*ptimestep/totarea_planet+       &
    1328                         SUM(massarea(:,:)*dqcldlscale(:,:))*ptimestep/totarea_planet,dWtot)
     1262                                         SUM(massarea(:,:)*dqcldlscale(:,:))*ptimestep/totarea_planet,dWtot)
     1263                       
    13291264                  if (is_master) print*,'In largescale atmospheric water change  =',dWtot,' kg m-2'
    1330                endif
    1331                !-------------------------
    1332 
    1333                ! compute cloud fraction
     1265               endif ! end of 'enertest'
     1266
     1267               ! Compute cloud fraction.
    13341268               do l = 1, nlayer
    13351269                  do ig=1,ngrid
    … …  
    13381272               enddo
    13391273
    1340             endif                ! of if (watercondense)
     1274            endif ! end of 'watercond'
    13411275           
    1342 
    1343 !        --------------------------------
    1344 !        Water ice / liquid precipitation
    1345 !        --------------------------------
     1276            ! Water ice / liquid precipitation.
     1277            ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    13461278            if(waterrain)then
    13471279
    … …  
    13511283
    13521284               call rain(ngrid,nlayer,nq,ptimestep,pplev,pplay,pt,pdt,pq,pdq,            &
    1353                    zdtrain,zdqrain,zdqsrain,zdqssnow,cloudfrac)
     1285                         zdtrain,zdqrain,zdqsrain,zdqssnow,cloudfrac)
    13541286
    13551287               pdq(1:ngrid,1:nlayer,igcm_h2o_vap) = pdq(1:ngrid,1:nlayer,igcm_h2o_vap) &
    1356                      +zdqrain(1:ngrid,1:nlayer,igcm_h2o_vap)
     1288                                                  + zdqrain(1:ngrid,1:nlayer,igcm_h2o_vap)
    13571289               pdq(1:ngrid,1:nlayer,igcm_h2o_ice) = pdq(1:ngrid,1:nlayer,igcm_h2o_ice) &
    1358                      +zdqrain(1:ngrid,1:nlayer,igcm_h2o_ice)
     1290                                                  + zdqrain(1:ngrid,1:nlayer,igcm_h2o_ice)
    13591291               pdt(1:ngrid,1:nlayer) = pdt(1:ngrid,1:nlayer)+zdtrain(1:ngrid,1:nlayer)
    1360                dqsurf(1:ngrid,igcm_h2o_vap) = dqsurf(1:ngrid,igcm_h2o_vap)+zdqsrain(1:ngrid) ! a bug was here
    1361                dqsurf(1:ngrid,igcm_h2o_ice) = dqsurf(1:ngrid,igcm_h2o_ice)+zdqssnow(1:ngrid)
    1362                rainout(1:ngrid)             = zdqsrain(1:ngrid)+zdqssnow(1:ngrid) ! diagnostic   
    1363 
    1364                !-------------------------
    1365                ! test energy conservation
     1292               
     1293               dqsurf(1:ngrid,igcm_h2o_vap) = dqsurf(1:ngrid,igcm_h2o_vap)+zdqsrain(1:ngrid)
     1294               dqsurf(1:ngrid,igcm_h2o_ice) = dqsurf(1:ngrid,igcm_h2o_ice)+zdqssnow(1:ngrid)
     1295
     1296               ! Test energy conservation.
    13661297               if(enertest)then
     1298               
    13671299                  call planetwide_sumval(cpp*massarea(:,:)*zdtrain(:,:)/totarea_planet,dEtot)
    13681300                  if (is_master) print*,'In rain atmospheric T energy change       =',dEtot,' W m-2'
    … …  
    13741306                  dVtot = dVtot + dVtot_tmp
    13751307                  dEtot = dItot + dVtot
     1308                 
    13761309                  if (is_master) then
    1377                         print*,'In rain dItot =',dItot,' W m-2'
    1378                         print*,'In rain dVtot =',dVtot,' W m-2'
    1379                         print*,'In rain atmospheric L energy change       =',dEtot,' W m-2'
     1310                     print*,'In rain dItot =',dItot,' W m-2'
     1311                     print*,'In rain dVtot =',dVtot,' W m-2'
     1312                     print*,'In rain atmospheric L energy change       =',dEtot,' W m-2'
    13801313                  endif
    13811314
    1382                ! test water conservation
     1315                  ! Test water conservation
    13831316                  call planetwide_sumval(massarea(:,:)*zdqrain(:,:,igcm_h2o_vap)*ptimestep/totarea_planet+      &
    13841317                        massarea(:,:)*zdqrain(:,:,igcm_h2o_ice)*ptimestep/totarea_planet,dWtot)
    13851318                  call planetwide_sumval((zdqsrain(:)+zdqssnow(:))*area(:)*ptimestep/totarea_planet,dWtots)
     1319                 
    13861320                  if (is_master) then
    13871321                        print*,'In rain atmospheric water change        =',dWtot,' kg m-2'
    … …  
    13891323                        print*,'In rain non-cons factor                 =',dWtot+dWtots,' kg m-2'
    13901324                  endif
    1391               endif
    1392               !-------------------------
    1393 
    1394             end if                 ! of if (waterrain)
    1395          end if                    ! of if (water)
    1396 
    1397 
    1398 !   7c. Aerosol particles
    1399 !     -------------------
    1400 !        -------------
    1401 !        Sedimentation
    1402 !        -------------
    1403         if (sedimentation) then
    1404            zdqsed(1:ngrid,1:nlayer,1:nq) = 0.0
    1405            zdqssed(1:ngrid,1:nq)  = 0.0
    1406 
    1407 
    1408            !-------------------------
    1409            ! find qtot
    1410            if(watertest)then
    1411               iq=igcm_h2o_ice
    1412               call planetwide_sumval(massarea(:,:)*pq(:,:,iq)*ptimestep/totarea_planet,dWtot)
    1413               call planetwide_sumval(massarea(:,:)*pdq(:,:,iq)*ptimestep/totarea_planet,dWtots)
    1414               if (is_master) then
    1415                 print*,'Before sedim pq  =',dWtot,' kg m-2'
    1416                 print*,'Before sedim pdq =',dWtots,' kg m-2'
    1417               endif
    1418            endif
    1419            !-------------------------
    1420 
    1421            call callsedim(ngrid,nlayer,ptimestep,           &
    1422                 pplev,zzlev,pt,pdt,pq,pdq,zdqsed,zdqssed,nq)
    1423 
    1424            !-------------------------
    1425            ! find qtot
    1426            if(watertest)then
    1427               iq=igcm_h2o_ice
    1428               call planetwide_sumval(massarea(:,:)*pq(:,:,iq)*ptimestep/totarea_planet,dWtot)
    1429               call planetwide_sumval(massarea(:,:)*pdq(:,:,iq)*ptimestep/totarea_planet,dWtots)
    1430               if (is_master) then
    1431                 print*,'After sedim pq  =',dWtot,' kg m-2'
    1432                 print*,'After sedim pdq =',dWtots,' kg m-2'
    1433               endif
    1434            endif
    1435            !-------------------------
    1436 
    1437            ! for now, we only allow H2O ice to sediment
    1438               ! and as in rain.F90, whether it falls as rain or snow depends
    1439               ! only on the surface temperature
    1440            pdq(1:ngrid,1:nlayer,1:nq) = pdq(1:ngrid,1:nlayer,1:nq) + zdqsed(1:ngrid,1:nlayer,1:nq)
    1441            dqsurf(1:ngrid,1:nq) = dqsurf(1:ngrid,1:nq) + zdqssed(1:ngrid,1:nq)
    1442 
    1443            !-------------------------
    1444            ! test water conservation
    1445            if(watertest)then
    1446               call planetwide_sumval(massarea(:,:)*(zdqsed(:,:,igcm_h2o_vap)+zdqsed(:,:,igcm_h2o_ice))*ptimestep/totarea_planet,dWtot)
    1447               call planetwide_sumval((zdqssed(:,igcm_h2o_vap)+zdqssed(:,igcm_h2o_ice))*area(:)*ptimestep/totarea_planet,dWtots)
    1448               if (is_master) then
    1449                 print*,'In sedim atmospheric ice change         =',dWtot,' kg m-2'
    1450                 print*,'In sedim surface ice change             =',dWtots,' kg m-2'
    1451                 print*,'In sedim non-cons factor                =',dWtot+dWtots,' kg m-2'
    1452               endif
    1453            endif
    1454            !-------------------------
    1455 
    1456         end if   ! of if (sedimentation)
    1457 
    1458 
    1459 !   7d. Updates
    1460 !     ---------
    1461 
    1462 !       -----------------------------------
    1463 !       Updating atm mass and tracer budget
    1464 !       -----------------------------------
    1465 
     1325                 
     1326               endif ! end of 'enertest'
     1327
     1328            end if ! enf of 'waterrain'
     1329           
     1330         end if ! end of 'water'
     1331
     1332  ! -------------------------
     1333  !   VI.2. Aerosol particles
     1334  ! -------------------------
     1335
     1336         ! Sedimentation.
     1337         if (sedimentation) then
     1338       
     1339            zdqsed(1:ngrid,1:nlayer,1:nq) = 0.0
     1340            zdqssed(1:ngrid,1:nq)  = 0.0
     1341
     1342            if(watertest)then
     1343           
     1344               iq=igcm_h2o_ice
     1345               call planetwide_sumval(massarea(:,:)*pq(:,:,iq)*ptimestep/totarea_planet,dWtot)
     1346               call planetwide_sumval(massarea(:,:)*pdq(:,:,iq)*ptimestep/totarea_planet,dWtots)
     1347               if (is_master) then
     1348                  print*,'Before sedim pq  =',dWtot,' kg m-2'
     1349                  print*,'Before sedim pdq =',dWtots,' kg m-2'
     1350               endif
     1351            endif
     1352           
     1353            call callsedim(ngrid,nlayer,ptimestep,           &
     1354                          pplev,zzlev,pt,pdt,pq,pdq,        &
     1355                          zdqsed,zdqssed,nq)
     1356
     1357            if(watertest)then
     1358               iq=igcm_h2o_ice
     1359               call planetwide_sumval(massarea(:,:)*pq(:,:,iq)*ptimestep/totarea_planet,dWtot)
     1360               call planetwide_sumval(massarea(:,:)*pdq(:,:,iq)*ptimestep/totarea_planet,dWtots)
     1361               if (is_master) then
     1362                  print*,'After sedim pq  =',dWtot,' kg m-2'
     1363                  print*,'After sedim pdq =',dWtots,' kg m-2'
     1364               endif
     1365            endif
     1366
     1367            ! Whether it falls as rain or snow depends only on the surface temperature
     1368            pdq(1:ngrid,1:nlayer,1:nq) = pdq(1:ngrid,1:nlayer,1:nq) + zdqsed(1:ngrid,1:nlayer,1:nq)
     1369            dqsurf(1:ngrid,1:nq) = dqsurf(1:ngrid,1:nq) + zdqssed(1:ngrid,1:nq)
     1370
     1371            ! Test water conservation
     1372            if(watertest)then
     1373               call planetwide_sumval(massarea(:,:)*(zdqsed(:,:,igcm_h2o_vap)+zdqsed(:,:,igcm_h2o_ice))*ptimestep/totarea_planet,dWtot)
     1374               call planetwide_sumval((zdqssed(:,igcm_h2o_vap)+zdqssed(:,igcm_h2o_ice))*area(:)*ptimestep/totarea_planet,dWtots)
     1375               if (is_master) then
     1376                  print*,'In sedim atmospheric ice change         =',dWtot,' kg m-2'
     1377                  print*,'In sedim surface ice change             =',dWtots,' kg m-2'
     1378                  print*,'In sedim non-cons factor                =',dWtot+dWtots,' kg m-2'
     1379               endif
     1380            endif
     1381
     1382         end if ! end of 'sedimentation'
     1383
     1384
     1385  ! ---------------
     1386  !   VI.3. Updates
     1387  ! ---------------
     1388
     1389         ! Updating Atmospheric Mass and Tracers budgets.
    14661390         if(mass_redistrib) then
    14671391
    1468             zdmassmr(1:ngrid,1:nlayer) = mass(1:ngrid,1:nlayer) * &
     1392            zdmassmr(1:ngrid,1:nlayer) = mass(1:ngrid,1:nlayer) *     &
    14691393               (   zdqevap(1:ngrid,1:nlayer)                          &
    14701394                 + zdqrain(1:ngrid,1:nlayer,igcm_h2o_vap)             &
    … …  
    14781402            call writediagfi(ngrid,"mass_evap","mass gain"," ",3,zdmassmr)
    14791403            call writediagfi(ngrid,"mass_evap_col","mass gain col"," ",2,zdmassmr_col)
    1480             call writediagfi(ngrid,"mass","mass"," ",3,mass)
    1481 
    1482             call mass_redistribution(ngrid,nlayer,nq,ptimestep,   &
    1483                        rnat,capcal,pplay,pplev,pt,tsurf,pq,qsurf,     &
    1484                        pu,pv,pdt,zdtsurf,pdq,pdu,pdv,zdmassmr,  &
    1485                        zdtmr,zdtsurfmr,zdpsrfmr,zdumr,zdvmr,zdqmr,zdqsurfmr)
     1404            call writediagfi(ngrid,"mass","mass","kg/m2",3,mass)
     1405
     1406            call mass_redistribution(ngrid,nlayer,nq,ptimestep,                     &
     1407                                     rnat,capcal,pplay,pplev,pt,tsurf,pq,qsurf,     &
     1408                                     pu,pv,pdt,zdtsurf,pdq,pdu,pdv,zdmassmr,        &
     1409                                     zdtmr,zdtsurfmr,zdpsrfmr,zdumr,zdvmr,zdqmr,zdqsurfmr)
    14861410         
    1487 
    14881411            pdq(1:ngrid,1:nlayer,1:nq) = pdq(1:ngrid,1:nlayer,1:nq) + zdqmr(1:ngrid,1:nlayer,1:nq)
    1489             dqsurf(1:ngrid,1:nq)         = dqsurf(1:ngrid,1:nq) + zdqsurfmr(1:ngrid,1:nq)
    1490             pdt(1:ngrid,1:nlayer)        = pdt(1:ngrid,1:nlayer) + zdtmr(1:ngrid,1:nlayer)
    1491             pdu(1:ngrid,1:nlayer)        = pdu(1:ngrid,1:nlayer) + zdumr(1:ngrid,1:nlayer)
    1492             pdv(1:ngrid,1:nlayer)        = pdv(1:ngrid,1:nlayer) + zdvmr(1:ngrid,1:nlayer)
    1493             pdpsrf(1:ngrid)                = pdpsrf(1:ngrid) + zdpsrfmr(1:ngrid)
    1494             zdtsurf(1:ngrid)               = zdtsurf(1:ngrid) + zdtsurfmr(1:ngrid)
     1412            dqsurf(1:ngrid,1:nq)       = dqsurf(1:ngrid,1:nq) + zdqsurfmr(1:ngrid,1:nq)
     1413            pdt(1:ngrid,1:nlayer)      = pdt(1:ngrid,1:nlayer) + zdtmr(1:ngrid,1:nlayer)
     1414            pdu(1:ngrid,1:nlayer)      = pdu(1:ngrid,1:nlayer) + zdumr(1:ngrid,1:nlayer)
     1415            pdv(1:ngrid,1:nlayer)      = pdv(1:ngrid,1:nlayer) + zdvmr(1:ngrid,1:nlayer)
     1416            pdpsrf(1:ngrid)            = pdpsrf(1:ngrid) + zdpsrfmr(1:ngrid)
     1417            zdtsurf(1:ngrid)           = zdtsurf(1:ngrid) + zdtsurfmr(1:ngrid)
    14951418           
    1496 !           print*,'after mass redistrib, q=',pq(211,1:nlayer,igcm_h2o_vap)+ptimestep*pdq(211,1:nlayer,igcm_h2o_vap)
    14971419         endif
    14981420
    1499 
    1500 !   7e. Ocean
    1501 !     ---------
    1502 
    1503 !       ---------------------------------
    1504 !       Slab_ocean
    1505 !       ---------------------------------
    1506       if (ok_slab_ocean)then
    1507 
    1508          do ig=1,ngrid
    1509             qsurfint(:,igcm_h2o_ice)=qsurf(:,igcm_h2o_ice)
    1510          enddo
    1511 
    1512          call ocean_slab_ice(ptimestep, &
    1513                ngrid, knindex, tsea_ice, fluxrad, &
    1514                zdqssnow, qsurf(:,igcm_h2o_ice), &
    1515                -zdqsdif(:,igcm_h2o_vap), &
    1516                flux_sens_lat,tsea_ice, pctsrf_sic, &
    1517                taux,tauy,icount)
    1518 
    1519 
    1520          call ocean_slab_get_vars(ngrid,tslab, &
    1521                sea_ice, flux_o, &
    1522                flux_g, dt_hdiff, &
    1523                dt_ekman)
    1524 
     1421  ! ------------------
     1422  !   VI.4. Slab Ocean
     1423  ! ------------------
     1424
     1425         if (ok_slab_ocean)then
     1426
     1427            do ig=1,ngrid
     1428               qsurfint(:,igcm_h2o_ice)=qsurf(:,igcm_h2o_ice)
     1429            enddo
     1430
     1431            call ocean_slab_ice(ptimestep,                          &
     1432                                ngrid, knindex, tsea_ice, fluxrad,  &
     1433                                zdqssnow, qsurf(:,igcm_h2o_ice),    &
     1434                              - zdqsdif(:,igcm_h2o_vap),            &
     1435                                flux_sens_lat,tsea_ice, pctsrf_sic, &
     1436                                taux,tauy,icount)
     1437
     1438
     1439            call ocean_slab_get_vars(ngrid,tslab,      &
     1440                                     sea_ice, flux_o,  &
     1441                                     flux_g, dt_hdiff, &
     1442                                     dt_ekman)
     1443   
    15251444            do ig=1,ngrid
    15261445               if (nint(rnat(ig)).eq.1)then
    1527                tslab(ig,1) = 0.
    1528                tslab(ig,2) = 0.
    1529                tsea_ice(ig) = 0.
    1530                sea_ice(ig) = 0.
    1531                pctsrf_sic(ig) = 0.
    1532                qsurf(ig,igcm_h2o_ice)=qsurfint(ig,igcm_h2o_ice)
     1446                  tslab(ig,1) = 0.
     1447                  tslab(ig,2) = 0.
     1448                  tsea_ice(ig) = 0.
     1449                  sea_ice(ig) = 0.
     1450                  pctsrf_sic(ig) = 0.
     1451                  qsurf(ig,igcm_h2o_ice) = qsurfint(ig,igcm_h2o_ice)
    15331452               endif
    15341453            enddo
    15351454
    1536 
    1537       endif! (ok_slab_ocean)
    1538 
    1539 !       ---------------------------------
    1540 !       Updating tracer budget on surface
    1541 !       ---------------------------------
     1455         endif ! end of 'ok_slab_ocean'
     1456
     1457  ! -----------------------------
     1458  !   VI.5. Surface Tracer Update
     1459  ! -----------------------------
    15421460
    15431461         if(hydrology)then
    15441462           
    1545             call hydrol(ngrid,nq,ptimestep,rnat,tsurf,qsurf,dqsurf,dqs_hyd,    &
    1546                capcal,albedo0,albedo,mu0,zdtsurf,zdtsurf_hyd,hice,pctsrf_sic,  &
    1547                sea_ice)
    1548             ! note: for now, also changes albedo in the subroutine
    1549 
    1550             zdtsurf(1:ngrid) = zdtsurf(1:ngrid) + zdtsurf_hyd(1:ngrid)
     1463            call hydrol(ngrid,nq,ptimestep,rnat,tsurf,qsurf,dqsurf,dqs_hyd,             &
     1464                        capcal,albedo0,albedo,mu0,zdtsurf,zdtsurf_hyd,hice,pctsrf_sic,  &
     1465                        sea_ice)
     1466
     1467            zdtsurf(1:ngrid)    = zdtsurf(1:ngrid) + zdtsurf_hyd(1:ngrid)
    15511468            qsurf(1:ngrid,1:nq) = qsurf(1:ngrid,1:nq) + ptimestep*dqs_hyd(1:ngrid,1:nq)
    1552             ! when hydrology is used, other dqsurf tendencies are all added to dqs_hyd inside
    1553 
    1554             !-------------------------
    1555             ! test energy conservation
     1469
     1470            ! Test energy conservation
    15561471            if(enertest)then
    15571472               call planetwide_sumval(area(:)*capcal(:)*zdtsurf_hyd(:)/totarea_planet,dEtots)
    15581473               if (is_master) print*,'In hydrol surface energy change     =',dEtots,' W m-2'
    15591474            endif
    1560             !-------------------------
    1561 
    1562             !-------------------------
     1475
    15631476            ! test water conservation
    15641477            if(watertest)then
    15651478               call planetwide_sumval(dqs_hyd(:,igcm_h2o_ice)*area(:)*ptimestep/totarea_planet,dWtots)
    15661479               if (is_master) print*,'In hydrol surface ice change            =',dWtots,' kg m-2'
    1567                call planetwide_sumval(dqs_hyd(:,igcm_h2o_vap)*area(:)*ptimestep/totarea_planet,dWtots)
     1480                  call planetwide_sumval(dqs_hyd(:,igcm_h2o_vap)*area(:)*ptimestep/totarea_planet,dWtots)
    15681481               if (is_master) then
    1569                 print*,'In hydrol surface water change          =',dWtots,' kg m-2'
    1570                 print*,'---------------------------------------------------------------'
     1482                  print*,'In hydrol surface water change          =',dWtots,' kg m-2'
     1483                  print*,'---------------------------------------------------------------'
    15711484               endif
    15721485            endif
    1573             !-------------------------
    1574 
    1575          ELSE     ! of if (hydrology)
     1486
     1487         else ! of if (hydrology)
    15761488
    15771489            qsurf(1:ngrid,1:nq)=qsurf(1:ngrid,1:nq)+ptimestep*dqsurf(1:ngrid,1:nq)
    15781490
    1579          END IF   ! of if (hydrology)
    1580 
    1581          ! Add qsurf to qsurf_hist, which is what we save in
    1582          ! diagfi.nc etc. At the same time, we set the water
    1583          ! content of ocean gridpoints back to zero, in order
    1584          ! to avoid rounding errors in vdifc, rain
     1491         end if ! of if (hydrology)
     1492
     1493         ! Add qsurf to qsurf_hist, which is what we save in diagfi.nc. At the same time, we set the water
     1494         ! content of ocean gridpoints back to zero, in order to avoid rounding errors in vdifc, rain.
    15851495         qsurf_hist(:,:) = qsurf(:,:)
    15861496
    … …  
    15891499         endif
    15901500
    1591       endif   !  of if (tracer)
    1592 
    1593 !-----------------------------------------------------------------------
    1594 !   9. Surface and sub-surface soil temperature
    1595 !-----------------------------------------------------------------------
    1596 
    1597 !     Increment surface temperature
    1598 
     1501      endif! end of if 'tracer'
     1502
     1503
     1504!------------------------------------------------
     1505!   VII. Surface and sub-surface soil temperature
     1506!------------------------------------------------
     1507
     1508
     1509      ! Increment surface temperature
    15991510      if(ok_slab_ocean)then
    16001511         do ig=1,ngrid
    … …  
    16081519      else
    16091520        tsurf(1:ngrid)=tsurf(1:ngrid)+ptimestep*zdtsurf(1:ngrid)
    1610       endif!(ok_slab_ocean)
    1611 
    1612 !     Compute soil temperatures and subsurface heat flux
     1521      endif ! end of 'ok_slab_ocean'
     1522
     1523
     1524      ! Compute soil temperatures and subsurface heat flux.
    16131525      if (callsoil) then
    16141526         call soil(ngrid,nsoilmx,.false.,lastcall,inertiedat,   &
    1615                 ptimestep,tsurf,tsoil,capcal,fluxgrd)     
     1527                   ptimestep,tsurf,tsoil,capcal,fluxgrd)     
    16161528      endif
    16171529
    16181530
    16191531      if (ok_slab_ocean) then
    1620             do ig=1,ngrid
    1621                fluxgrdocean(ig)=fluxgrd(ig)
    1622                if (nint(rnat(ig)).eq.0) then
     1532     
     1533         do ig=1,ngrid
     1534         
     1535            fluxgrdocean(ig)=fluxgrd(ig)
     1536            if (nint(rnat(ig)).eq.0) then
    16231537               capcal(ig)=capcalocean
    16241538               fluxgrd(ig)=0.
    16251539               fluxgrdocean(ig)=pctsrf_sic(ig)*flux_g(ig)+(1-pctsrf_sic(ig))*(dt_hdiff(ig,1)+dt_ekman(ig,1))
    1626                  do iq=1,nsoilmx
    1627                     tsoil(ig,iq)=tsurf(ig)
    1628                  enddo
    1629                  if (pctsrf_sic(ig).gt.0.5) then
    1630                    capcal(ig)=capcalseaice
    1631                    if (qsurf(ig,igcm_h2o_ice).gt.0.) then
    1632                    capcal(ig)=capcalsno
    1633                    endif
    1634                  endif
    1635                endif
    1636             enddo
    1637       endif !(ok_slab_ocean)
    1638 
    1639 !-------------------------
    1640 ! test energy conservation
     1540               do iq=1,nsoilmx
     1541                  tsoil(ig,iq)=tsurf(ig)
     1542               enddo
     1543               if (pctsrf_sic(ig).gt.0.5) then
     1544                  capcal(ig)=capcalseaice
     1545                  if (qsurf(ig,igcm_h2o_ice).gt.0.) then
     1546                     capcal(ig)=capcalsno
     1547                  endif
     1548               endif               
     1549            endif
     1550           
     1551         enddo
     1552         
     1553      endif !end of 'ok_slab_ocean'
     1554
     1555
     1556      ! Test energy conservation
    16411557      if(enertest)then
    16421558         call planetwide_sumval(area(:)*capcal(:)*zdtsurf(:)/totarea_planet,dEtots)     
    16431559         if (is_master) print*,'Surface energy change                 =',dEtots,' W m-2'
    16441560      endif
    1645 !-------------------------
    1646 
    1647 !-----------------------------------------------------------------------
    1648 !  10. Perform diagnostics and write output files
    1649 !-----------------------------------------------------------------------
    1650 
    1651 !    -------------------------------
    1652 !    Dynamical fields incrementation
    1653 !    -------------------------------
    1654 !    For output only: the actual model integration is performed in the dynamics
     1561
     1562
     1563!---------------------------------------------------
     1564!   VIII. Perform diagnostics and write output files
     1565!---------------------------------------------------
     1566
     1567      ! Note : For output only: the actual model integration is performed in the dynamics.
     1568
     1569
    16551570 
    1656       ! temperature, zonal and meridional wind
     1571      ! Temperature, zonal and meridional winds.
    16571572      zt(1:ngrid,1:nlayer) = pt(1:ngrid,1:nlayer) + pdt(1:ngrid,1:nlayer)*ptimestep
    16581573      zu(1:ngrid,1:nlayer) = pu(1:ngrid,1:nlayer) + pdu(1:ngrid,1:nlayer)*ptimestep
    16591574      zv(1:ngrid,1:nlayer) = pv(1:ngrid,1:nlayer) + pdv(1:ngrid,1:nlayer)*ptimestep
    16601575
    1661       ! diagnostic
     1576      ! Diagnostic.
    16621577      zdtdyn(1:ngrid,1:nlayer)     = pt(1:ngrid,1:nlayer)-ztprevious(1:ngrid,1:nlayer)
    16631578      ztprevious(1:ngrid,1:nlayer) = zt(1:ngrid,1:nlayer)
    … …  
    16671582      endif
    16681583
    1669       ! dynamical heating diagnostic
     1584      ! Dynamical heating diagnostic.
    16701585      do ig=1,ngrid
    16711586         fluxdyn(ig)= SUM(zdtdyn(ig,:) *mass(ig,:))*cpp/ptimestep
    16721587      enddo
    16731588
    1674       ! tracers
     1589      ! Tracers.
    16751590      zq(1:ngrid,1:nlayer,1:nq) = pq(1:ngrid,1:nlayer,1:nq) + pdq(1:ngrid,1:nlayer,1:nq)*ptimestep
    16761591
    1677       ! surface pressure
     1592      ! Surface pressure.
    16781593      ps(1:ngrid) = pplev(1:ngrid,1) + pdpsrf(1:ngrid)*ptimestep
    16791594
    1680       ! pressure
    1681       do l=1,nlayer
    1682           zplev(1:ngrid,l) = pplev(1:ngrid,l)/pplev(1:ngrid,1)*ps(:)
    1683           zplay(1:ngrid,l) = pplay(1:ngrid,l)/pplev(1:ngrid,1)*ps(1:ngrid)
    1684       enddo
    1685 
    1686 !     ---------------------------------------------------------
    1687 !     Surface and soil temperature information
    1688 !     ---------------------------------------------------------
    1689 
     1595
     1596
     1597      ! Surface and soil temperature information
    16901598      call planetwide_sumval(area(:)*tsurf(:)/totarea_planet,Ts1)
    16911599      call planetwide_minval(tsurf(:),Ts2)
    … …  
    16931601      if(callsoil)then
    16941602         TsS = SUM(area(:)*tsoil(:,nsoilmx))/totarea        ! mean temperature at bottom soil layer
    1695            print*,'  ave[Tsurf]     min[Tsurf]     max[Tsurf]     ave[Tdeep]'
    1696            print*,Ts1,Ts2,Ts3,TsS
     1603         print*,'  ave[Tsurf]     min[Tsurf]     max[Tsurf]     ave[Tdeep]'
     1604         print*,Ts1,Ts2,Ts3,TsS
    16971605      else
    1698         if (is_master) then
    1699            print*,'  ave[Tsurf]     min[Tsurf]     max[Tsurf]'
    1700            print*,Ts1,Ts2,Ts3
    1701         endif
     1606         if (is_master) then
     1607            print*,'  ave[Tsurf]     min[Tsurf]     max[Tsurf]'
     1608            print*,Ts1,Ts2,Ts3
     1609         endif
    17021610      end if
    17031611
    1704 !     ---------------------------------------------------------
    1705 !     Check the energy balance of the simulation during the run
    1706 !     ---------------------------------------------------------
    1707 
     1612
     1613      ! Check the energy balance of the simulation during the run
    17081614      if(corrk)then
    17091615
    … …  
    17301636         
    17311637         if (is_master) then
    1732                 print*,'  ISR            ASR            OLR            GND            DYN [W m^-2]'
    1733                 print*, ISR,ASR,OLR,GND,DYN
     1638            print*,'  ISR            ASR            OLR            GND            DYN [W m^-2]'
     1639            print*, ISR,ASR,OLR,GND,DYN
    17341640         endif
    17351641
    … …  
    17481654               open(93,file="tem_bal.out",form='formatted',position='append')
    17491655               if(callsoil)then
    1750                 write(93,*) zday,Ts1,Ts2,Ts3,TsS
     1656                  write(93,*) zday,Ts1,Ts2,Ts3,TsS
    17511657               else
    1752                 write(93,*) zday,Ts1,Ts2,Ts3
     1658                  write(93,*) zday,Ts1,Ts2,Ts3
    17531659               endif
    17541660               close(93)
    … …  
    17561662         endif
    17571663
    1758       endif
    1759 
    1760 
    1761 !     ------------------------------------------------------------------
    1762 !     Diagnostic to test radiative-convective timescales in code
    1763 !     ------------------------------------------------------------------
     1664      endif ! end of 'corrk'
     1665
     1666
     1667      ! Diagnostic to test radiative-convective timescales in code.
    17641668      if(testradtimes)then
    17651669         open(38,file="tau_phys.out",form='formatted',position='append')
    … …  
    17741678      endif
    17751679
    1776 !     ---------------------------------------------------------
    1777 !     Compute column amounts (kg m-2) if tracers are enabled
    1778 !     ---------------------------------------------------------
     1680
     1681      ! Compute column amounts (kg m-2) if tracers are enabled.
    17791682      if(tracer)then
    17801683         qcol(1:ngrid,1:nq)=0.0
    17811684         do iq=1,nq
    1782            do ig=1,ngrid
    1783               qcol(ig,iq) = SUM( zq(ig,1:nlayer,iq) * mass(ig,1:nlayer))
    1784            enddo
     1685            do ig=1,ngrid
     1686               qcol(ig,iq) = SUM( zq(ig,1:nlayer,iq) * mass(ig,1:nlayer))
     1687            enddo
    17851688         enddo
    17861689
    1787          ! Generalised for arbitrary aerosols now. (LK)
     1690         ! Generalised for arbitrary aerosols now. By LK
    17881691         reffcol(1:ngrid,1:naerkind)=0.0
    17891692         if(co2cond.and.(iaero_co2.ne.0))then
    … …  
    18011704         endif
    18021705
    1803       endif
    1804 
    1805 !     ---------------------------------------------------------
    1806 !     Test for water conservation if water is enabled
    1807 !     ---------------------------------------------------------
    1808 
     1706      endif ! end of 'tracer'
     1707
     1708
     1709      ! Test for water conservation.
    18091710      if(water)then
    18101711
    … …  
    18171718         
    18181719         if (is_master) then
    1819                 print*,' Total water amount [kg m^-2]: ',h2otot
    1820                 print*,' Surface ice    Surface liq.   Atmos. con.     Atmos. vap. [kg m^-2] '
    1821                 print*, icesrf,liqsrf,icecol,vapcol
     1720            print*,' Total water amount [kg m^-2]: ',h2otot
     1721            print*,' Surface ice    Surface liq.   Atmos. con.     Atmos. vap. [kg m^-2] '
     1722            print*, icesrf,liqsrf,icecol,vapcol
    18221723         endif
    18231724
    … …  
    18331734      endif
    18341735
    1835 !     ---------------------------------------------------------
    1836 !     Calculate RH for diagnostic if water is enabled
    1837 !     ---------------------------------------------------------
    1838 
     1736
     1737      ! Calculate RH (Relative Humidity) for diagnostic.
    18391738      if(water)then
     1739     
    18401740         do l = 1, nlayer
    18411741            do ig=1,ngrid
    … …  
    18451745         enddo
    18461746
    1847          ! compute maximum possible H2O column amount (100% saturation)
     1747         ! Compute maximum possible H2O column amount (100% saturation).
    18481748         do ig=1,ngrid
    1849                H2Omaxcol(ig) = SUM( qsat(ig,:) * mass(ig,:))
     1749            H2Omaxcol(ig) = SUM( qsat(ig,:) * mass(ig,:))
    18501750         enddo
    18511751
    1852       endif
    1853 
    1854 
    1855            print*,'--> Ls =',zls*180./pi
    1856 !        -------------------------------------------------------------------
     1752      endif ! end of 'water'
     1753
     1754
     1755      print*,'--> Ls =',zls*180./pi
     1756     
     1757     
     1758!----------------------------------------------------------------------
    18571759!        Writing NetCDF file  "RESTARTFI" at the end of the run
    1858 !        -------------------------------------------------------------------
     1760!----------------------------------------------------------------------
     1761
    18591762!        Note: 'restartfi' is stored just before dynamics are stored
    18601763!              in 'restart'. Between now and the writting of 'restart',
    … …  
    18631766!              thus we store for time=time+dtvr
    18641767
    1865          if(lastcall) then
    1866             ztime_fin = ptime + ptimestep/(float(iphysiq)*daysec)
    1867 
    1868 
    1869 !           Update surface ice distribution to iterate to steady state if requested
    1870             if(ice_update)then
    1871 
    1872                do ig=1,ngrid
    1873 
    1874                   delta_ice = (qsurf(ig,igcm_h2o_ice)-ice_initial(ig))
    1875 
    1876                   ! add multiple years of evolution
    1877                   qsurf_hist(ig,igcm_h2o_ice) = qsurf_hist(ig,igcm_h2o_ice) + delta_ice*icetstep
    1878 
    1879                   ! if ice has gone -ve, set to zero
    1880                   if(qsurf_hist(ig,igcm_h2o_ice).lt.0.0)then
    1881                      qsurf_hist(ig,igcm_h2o_ice) = 0.0
    1882                   endif
    1883 
    1884                   ! if ice is seasonal, set to zero (NEW)
    1885                   if(ice_min(ig).lt.0.01)then
    1886                      qsurf_hist(ig,igcm_h2o_ice) = 0.0
    1887                   endif
    1888 
    1889                enddo
    1890 
    1891                ! enforce ice conservation
    1892                ice_tot= SUM(qsurf_hist(:,igcm_h2o_ice)*area(:) )
    1893                qsurf_hist(:,igcm_h2o_ice) = qsurf_hist(:,igcm_h2o_ice)*(icesrf/ice_tot)
    1894 
    1895             endif
    1896 
    1897             if (ngrid.ne.1) then
    1898               write(*,*)'PHYSIQ: for physdem ztime_fin =',ztime_fin
    1899 !#ifdef CPP_PARA
    1900 !! for now in parallel we use a different routine to write restartfi.nc
    1901 !            call phyredem(ngrid,"restartfi.nc",           &
    1902 !                    ptimestep,pday,ztime_fin,tsurf,tsoil,emis,q2,qsurf_hist, &
    1903 !                    cloudfrac,totcloudfrac,hice)
    1904 !#else
    1905 !            call physdem1(ngrid,"restartfi.nc",long,lati,nsoilmx,nq,            &
    1906 !                    ptimestep,pday,ztime_fin,tsurf,tsoil,emis,q2,qsurf_hist, &
    1907 !                    area,albedodat,inertiedat,zmea,zstd,zsig,zgam,zthe,      &
    1908 !                    cloudfrac,totcloudfrac,hice,noms)
    1909 !#endif
    1910               call physdem1("restartfi.nc",nsoilmx,ngrid,nlayer,nq, &
    1911                       ptimestep,ztime_fin, &
    1912                       tsurf,tsoil,emis,q2,qsurf_hist, &
    1913                       cloudfrac,totcloudfrac,hice, &
    1914                       rnat,pctsrf_sic,tslab,tsea_ice,sea_ice)
    1915             endif
    1916 
    1917             if(ok_slab_ocean) then
    1918               call ocean_slab_final!(tslab, seaice)
    1919             end if
     1768
     1769
     1770      if(lastcall) then
     1771         ztime_fin = ptime + ptimestep/(float(iphysiq)*daysec)
     1772
     1773         ! Update surface ice distribution to iterate to steady state if requested
     1774         if(ice_update)then
     1775
     1776            do ig=1,ngrid
     1777
     1778               delta_ice = (qsurf(ig,igcm_h2o_ice)-ice_initial(ig))
     1779               
     1780               ! add multiple years of evolution
     1781               qsurf_hist(ig,igcm_h2o_ice) = qsurf_hist(ig,igcm_h2o_ice) + delta_ice*icetstep
     1782
     1783               ! if ice has gone -ve, set to zero
     1784               if(qsurf_hist(ig,igcm_h2o_ice).lt.0.0)then
     1785                  qsurf_hist(ig,igcm_h2o_ice) = 0.0
     1786               endif
     1787
     1788               ! if ice is seasonal, set to zero (NEW)
     1789               if(ice_min(ig).lt.0.01)then
     1790                  qsurf_hist(ig,igcm_h2o_ice) = 0.0
     1791               endif
     1792
     1793            enddo
     1794
     1795            ! enforce ice conservation
     1796            ice_tot= SUM(qsurf_hist(:,igcm_h2o_ice)*area(:) )
     1797            qsurf_hist(:,igcm_h2o_ice) = qsurf_hist(:,igcm_h2o_ice)*(icesrf/ice_tot)
     1798
     1799         endif
     1800
     1801         if (ngrid.ne.1) then
     1802            write(*,*)'PHYSIQ: for physdem ztime_fin =',ztime_fin
     1803           
     1804            call physdem1("restartfi.nc",nsoilmx,ngrid,nlayer,nq, &
     1805                          ptimestep,ztime_fin,                    &
     1806                          tsurf,tsoil,emis,q2,qsurf_hist,         &
     1807                          cloudfrac,totcloudfrac,hice,            &
     1808                          rnat,pctsrf_sic,tslab,tsea_ice,sea_ice)
     1809         endif
     1810
     1811         if(ok_slab_ocean) then
     1812            call ocean_slab_final!(tslab, seaice)
     1813         end if
    19201814
    19211815         
    1922          endif ! of if (lastcall)
    1923 
    1924 
    1925 !        -----------------------------------------------------------------
     1816      endif ! end of 'lastcall'
     1817
     1818
     1819!-----------------------------------
    19261820!        Saving statistics :
    1927 !        -----------------------------------------------------------------
    1928 !        ("stats" stores and accumulates 8 key variables in file "stats.nc"
    1929 !        which can later be used to make the statistic files of the run:
    1930 !        "stats")          only possible in 3D runs !
     1821!-----------------------------------
     1822
     1823!    Note :("stats" stores and accumulates 8 key variables in file "stats.nc"
     1824!           which can later be used to make the statistic files of the run:
     1825!           "stats")          only possible in 3D runs !!!
    19311826
    19321827         
    1933          if (callstats) then
    1934 
    1935             call wstats(ngrid,"ps","Surface pressure","Pa",2,ps)
    1936             call wstats(ngrid,"tsurf","Surface temperature","K",2,tsurf)
    1937             call wstats(ngrid,"fluxsurf_lw",                               &
    1938                         "Thermal IR radiative flux to surface","W.m-2",2,  &
    1939                          fluxsurf_lw)
     1828      if (callstats) then
     1829
     1830         call wstats(ngrid,"ps","Surface pressure","Pa",2,ps)
     1831         call wstats(ngrid,"tsurf","Surface temperature","K",2,tsurf)
     1832         call wstats(ngrid,"fluxsurf_lw",                               &
     1833                     "Thermal IR radiative flux to surface","W.m-2",2,  &
     1834                     fluxsurf_lw)
     1835         call wstats(ngrid,"fluxtop_lw",                                &
     1836                     "Thermal IR radiative flux to space","W.m-2",2,    &
     1837                     fluxtop_lw)
     1838                     
    19401839!            call wstats(ngrid,"fluxsurf_sw",                               &
    19411840!                        "Solar radiative flux to surface","W.m-2",2,       &
    1942 !                         fluxsurf_sw_tot)
    1943             call wstats(ngrid,"fluxtop_lw",                                &
    1944                         "Thermal IR radiative flux to space","W.m-2",2,    &
    1945                         fluxtop_lw)
     1841!                         fluxsurf_sw_tot)                     
    19461842!            call wstats(ngrid,"fluxtop_sw",                                &
    19471843!                        "Solar radiative flux to space","W.m-2",2,         &
    19481844!                        fluxtop_sw_tot)
    19491845
    1950             call wstats(ngrid,"ISR","incoming stellar rad.","W m-2",2,fluxtop_dn)
    1951             call wstats(ngrid,"ASR","absorbed stellar rad.","W m-2",2,fluxabs_sw)
    1952             call wstats(ngrid,"OLR","outgoing longwave rad.","W m-2",2,fluxtop_lw)
    1953             call wstats(ngrid,"ALB","Surface albedo"," ",2,albedo)
    1954             call wstats(ngrid,"p","Pressure","Pa",3,pplay)
    1955             call wstats(ngrid,"temp","Atmospheric temperature","K",3,zt)
    1956             call wstats(ngrid,"u","Zonal (East-West) wind","m.s-1",3,zu)
    1957             call wstats(ngrid,"v","Meridional (North-South) wind","m.s-1",3,zv)
    1958             call wstats(ngrid,"w","Vertical (down-up) wind","m.s-1",3,pw)
    1959             call wstats(ngrid,"q2","Boundary layer eddy kinetic energy","m2.s-2",3,q2)
    1960 
    1961            if (tracer) then
    1962              do iq=1,nq
    1963                 call wstats(ngrid,noms(iq),noms(iq),'kg/kg',3,zq(1,1,iq))
    1964                 call wstats(ngrid,trim(noms(iq))//'_surf',trim(noms(iq))//'_surf',  &
    1965                           'kg m^-2',2,qsurf(1,iq) )
    1966 
    1967                 call wstats(ngrid,trim(noms(iq))//'_col',trim(noms(iq))//'_col',    &
     1846
     1847         call wstats(ngrid,"ISR","incoming stellar rad.","W m-2",2,fluxtop_dn)
     1848         call wstats(ngrid,"ASR","absorbed stellar rad.","W m-2",2,fluxabs_sw)
     1849         call wstats(ngrid,"OLR","outgoing longwave rad.","W m-2",2,fluxtop_lw)
     1850         call wstats(ngrid,"ALB","Surface albedo"," ",2,albedo)
     1851         call wstats(ngrid,"p","Pressure","Pa",3,pplay)
     1852         call wstats(ngrid,"temp","Atmospheric temperature","K",3,zt)
     1853         call wstats(ngrid,"u","Zonal (East-West) wind","m.s-1",3,zu)
     1854         call wstats(ngrid,"v","Meridional (North-South) wind","m.s-1",3,zv)
     1855         call wstats(ngrid,"w","Vertical (down-up) wind","m.s-1",3,pw)
     1856         call wstats(ngrid,"q2","Boundary layer eddy kinetic energy","m2.s-2",3,q2)
     1857
     1858         if (tracer) then
     1859            do iq=1,nq
     1860               call wstats(ngrid,noms(iq),noms(iq),'kg/kg',3,zq(1,1,iq))
     1861               call wstats(ngrid,trim(noms(iq))//'_surf',trim(noms(iq))//'_surf',  &
     1862                           'kg m^-2',2,qsurf(1,iq) )
     1863               call wstats(ngrid,trim(noms(iq))//'_col',trim(noms(iq))//'_col',    &
    19681864                          'kg m^-2',2,qcol(1,iq) )
    1969 !                call wstats(ngrid,trim(noms(iq))//'_reff',                          &
     1865                         
     1866!              call wstats(ngrid,trim(noms(iq))//'_reff',                          &
    19701867!                          trim(noms(iq))//'_reff',                                   &
    19711868!                          'm',3,reffrad(1,1,iq))
    1972               end do
     1869
     1870            end do
     1871           
    19731872            if (water) then
    19741873               vmr=zq(1:ngrid,1:nlayer,igcm_h2o_vap)*mugaz/mmol(igcm_h2o_vap)
    1975                call wstats(ngrid,"vmr_h2ovapor",                           &
    1976                           "H2O vapour volume mixing ratio","mol/mol",       &
    1977                           3,vmr)
    1978             endif ! of if (water)
    1979 
    1980            endif !tracer
    1981 
    1982           if(watercond.and.CLFvarying)then
    1983              call wstats(ngrid,"rneb_man","H2O cloud fraction (conv)"," ",3,rneb_man)
    1984              call wstats(ngrid,"rneb_lsc","H2O cloud fraction (large scale)"," ",3,rneb_lsc)
    1985              call wstats(ngrid,"CLF","H2O cloud fraction"," ",3,cloudfrac)
    1986              call wstats(ngrid,"CLFt","H2O column cloud fraction"," ",2,totcloudfrac)
    1987              call wstats(ngrid,"RH","relative humidity"," ",3,RH)
    1988           endif
    1989 
    1990 
    1991 
    1992            if (ok_slab_ocean) then
     1874               call wstats(ngrid,"vmr_h2ovapor",                             &
     1875                           "H2O vapour volume mixing ratio","mol/mol",       &
     1876                           3,vmr)
     1877            endif
     1878
     1879         endif ! end of 'tracer'
     1880
     1881         if(watercond.and.CLFvarying)then
     1882            call wstats(ngrid,"rneb_man","H2O cloud fraction (conv)"," ",3,rneb_man)
     1883            call wstats(ngrid,"rneb_lsc","H2O cloud fraction (large scale)"," ",3,rneb_lsc)
     1884            call wstats(ngrid,"CLF","H2O cloud fraction"," ",3,cloudfrac)
     1885            call wstats(ngrid,"CLFt","H2O column cloud fraction"," ",2,totcloudfrac)
     1886            call wstats(ngrid,"RH","relative humidity"," ",3,RH)
     1887         endif
     1888
     1889         if (ok_slab_ocean) then
    19931890            call wstats(ngrid,"dt_hdiff1","dt_hdiff1","K/s",2,dt_hdiff(:,1))
    19941891            call wstats(ngrid,"dt_hdiff2","dt_hdiff2","K/s",2,dt_hdiff(:,2))
    … …  
    20011898            call wstats(ngrid,"sea_ice","sea ice","kg/m2",2,sea_ice)
    20021899            call wstats(ngrid,"rnat","nature of the surface","",2,rnat)
    2003            endif! (ok_slab_ocean)
    2004 
    2005             if(lastcall) then
    2006               write (*,*) "Writing stats..."
    2007               call mkstats(ierr)
    2008             endif
    2009           endif !if callstats
    2010 
    2011 
    2012 !        ----------------------------------------------------------------------
    2013 !        output in netcdf file "DIAGFI", containing any variable for diagnostic
    2014 !        (output with  period "ecritphy", set in "run.def")
    2015 !        ----------------------------------------------------------------------
    2016 !        writediagfi can also be called from any other subroutine for any variable.
    2017 !        but its preferable to keep all the calls in one place...
    2018 
    2019         call writediagfi(ngrid,"Ls","solar longitude","deg",0,zls*180./pi)
    2020         call writediagfi(ngrid,"Lss","sub solar longitude","deg",0,zlss*180./pi)
    2021         call writediagfi(ngrid,"RA","right ascension","deg",0,right_ascen*180./pi)
    2022         call writediagfi(ngrid,"Declin","solar declination","deg",0,declin*180./pi)
    2023         call writediagfi(ngrid,"tsurf","Surface temperature","K",2,tsurf)
    2024         call writediagfi(ngrid,"ps","Surface pressure","Pa",2,ps)
    2025         call writediagfi(ngrid,"temp","temperature","K",3,zt)
    2026         call writediagfi(ngrid,"teta","potential temperature","K",3,zh)
    2027         call writediagfi(ngrid,"u","Zonal wind","m.s-1",3,zu)
    2028         call writediagfi(ngrid,"v","Meridional wind","m.s-1",3,zv)
    2029         call writediagfi(ngrid,"w","Vertical wind","m.s-1",3,pw)
    2030         call writediagfi(ngrid,"p","Pressure","Pa",3,pplay)
     1900         endif
     1901
     1902         if(lastcall) then
     1903            write (*,*) "Writing stats..."
     1904            call mkstats(ierr)
     1905         endif
     1906         
     1907      endif ! end of 'callstats'
     1908
     1909
     1910!-----------------------------------------------------------------------------------------------------
     1911!           OUTPUT in netcdf file "DIAGFI.NC", containing any variable for diagnostic
     1912!
     1913!             Note 1 : output with  period "ecritphy", set in "run.def"
     1914!
     1915!             Note 2 : writediagfi can also be called from any other subroutine for any variable,
     1916!                      but its preferable to keep all the calls in one place ...
     1917!-----------------------------------------------------------------------------------------------------
     1918
     1919
     1920      call writediagfi(ngrid,"Ls","solar longitude","deg",0,zls*180./pi)
     1921      call writediagfi(ngrid,"Lss","sub solar longitude","deg",0,zlss*180./pi)
     1922      call writediagfi(ngrid,"RA","right ascension","deg",0,right_ascen*180./pi)
     1923      call writediagfi(ngrid,"Declin","solar declination","deg",0,declin*180./pi)
     1924      call writediagfi(ngrid,"tsurf","Surface temperature","K",2,tsurf)
     1925      call writediagfi(ngrid,"ps","Surface pressure","Pa",2,ps)
     1926      call writediagfi(ngrid,"temp","temperature","K",3,zt)
     1927      call writediagfi(ngrid,"teta","potential temperature","K",3,zh)
     1928      call writediagfi(ngrid,"u","Zonal wind","m.s-1",3,zu)
     1929      call writediagfi(ngrid,"v","Meridional wind","m.s-1",3,zv)
     1930      call writediagfi(ngrid,"w","Vertical wind","m.s-1",3,pw)
     1931      call writediagfi(ngrid,"p","Pressure","Pa",3,pplay)
    20311932
    20321933!     Subsurface temperatures
    … …  
    20341935!        call writediagsoil(ngrid,"temp","temperature","K",3,tsoil)
    20351936
    2036 !     Oceanic layers
    2037         if(ok_slab_ocean) then
    2038           call writediagfi(ngrid,"pctsrf_sic","pct ice/sea","",2,pctsrf_sic)
    2039           call writediagfi(ngrid,"tsea_ice","tsea_ice","K",2,tsea_ice)
    2040           call writediagfi(ngrid,"sea_ice","sea ice","kg/m2",2,sea_ice)
    2041           call writediagfi(ngrid,"tslab1","tslab1","K",2,tslab(:,1))
    2042           call writediagfi(ngrid,"tslab2","tslab2","K",2,tslab(:,2))
    2043           call writediagfi(ngrid,"dt_hdiff1","dt_hdiff1","K/s",2,dt_hdiff(:,1))
    2044           call writediagfi(ngrid,"dt_hdiff2","dt_hdiff2","K/s",2,dt_hdiff(:,2))
    2045           call writediagfi(ngrid,"dt_ekman1","dt_ekman1","K/s",2,dt_ekman(:,1))
    2046           call writediagfi(ngrid,"dt_ekman2","dt_ekman2","K/s",2,dt_ekman(:,2))
    2047           call writediagfi(ngrid,"rnat","nature of the surface","",2,rnat)
    2048           call writediagfi(ngrid,"sensibFlux","sensible heat flux","w.m^-2",2,sensibFlux)
    2049           call writediagfi(ngrid,"latentFlux","latent heat flux","w.m^-2",2,zdqsdif(:,igcm_h2o_vap)*RLVTT)
    2050        endif! (ok_slab_ocean)
    2051 
    2052 !     Total energy balance diagnostics
    2053         if(callrad.and.(.not.newtonian))then
    2054            call writediagfi(ngrid,"ALB","Surface albedo"," ",2,albedo)
    2055            call writediagfi(ngrid,"ISR","incoming stellar rad.","W m-2",2,fluxtop_dn)
    2056            call writediagfi(ngrid,"ASR","absorbed stellar rad.","W m-2",2,fluxabs_sw)
    2057            call writediagfi(ngrid,"OLR","outgoing longwave rad.","W m-2",2,fluxtop_lw)
    2058            call writediagfi(ngrid,"shad","rings"," ", 2, fract)
     1937      ! Oceanic layers
     1938      if(ok_slab_ocean) then
     1939         call writediagfi(ngrid,"pctsrf_sic","pct ice/sea","",2,pctsrf_sic)
     1940         call writediagfi(ngrid,"tsea_ice","tsea_ice","K",2,tsea_ice)
     1941         call writediagfi(ngrid,"sea_ice","sea ice","kg/m2",2,sea_ice)
     1942         call writediagfi(ngrid,"tslab1","tslab1","K",2,tslab(:,1))
     1943         call writediagfi(ngrid,"tslab2","tslab2","K",2,tslab(:,2))
     1944         call writediagfi(ngrid,"dt_hdiff1","dt_hdiff1","K/s",2,dt_hdiff(:,1))
     1945         call writediagfi(ngrid,"dt_hdiff2","dt_hdiff2","K/s",2,dt_hdiff(:,2))
     1946         call writediagfi(ngrid,"dt_ekman1","dt_ekman1","K/s",2,dt_ekman(:,1))
     1947         call writediagfi(ngrid,"dt_ekman2","dt_ekman2","K/s",2,dt_ekman(:,2))
     1948         call writediagfi(ngrid,"rnat","nature of the surface","",2,rnat)
     1949         call writediagfi(ngrid,"sensibFlux","sensible heat flux","w.m^-2",2,sensibFlux)
     1950         call writediagfi(ngrid,"latentFlux","latent heat flux","w.m^-2",2,zdqsdif(:,igcm_h2o_vap)*RLVTT)
     1951      endif
     1952     
     1953
     1954      ! Total energy balance diagnostics
     1955      if(callrad.and.(.not.newtonian))then
     1956     
     1957         call writediagfi(ngrid,"ALB","Surface albedo"," ",2,albedo)
     1958         call writediagfi(ngrid,"ISR","incoming stellar rad.","W m-2",2,fluxtop_dn)
     1959         call writediagfi(ngrid,"ASR","absorbed stellar rad.","W m-2",2,fluxabs_sw)
     1960         call writediagfi(ngrid,"OLR","outgoing longwave rad.","W m-2",2,fluxtop_lw)
     1961         call writediagfi(ngrid,"shad","rings"," ", 2, fract)
     1962         
    20591963!           call writediagfi(ngrid,"ASRcs","absorbed stellar rad (cs).","W m-2",2,fluxabs_sw1)
    20601964!           call writediagfi(ngrid,"OLRcs","outgoing longwave rad (cs).","W m-2",2,fluxtop_lw1)
    … …  
    20631967!           call writediagfi(ngrid,"fluxsurfswcs","sw surface flux (cs).","W m-2",2,fluxsurf_sw1)
    20641968!           call writediagfi(ngrid,"fluxsurflwcs","lw back radiation (cs).","W m-2",2,fluxsurf_lw1)
    2065            if(ok_slab_ocean) then
    2066              call writediagfi(ngrid,"GND","heat flux from ground","W m-2",2,fluxgrdocean)
    2067            else
    2068              call writediagfi(ngrid,"GND","heat flux from ground","W m-2",2,fluxgrd)
    2069            endif!(ok_slab_ocean)
    2070            call writediagfi(ngrid,"DYN","dynamical heat input","W m-2",2,fluxdyn)
    2071          endif
     1969
     1970         if(ok_slab_ocean) then
     1971            call writediagfi(ngrid,"GND","heat flux from ground","W m-2",2,fluxgrdocean)
     1972         else
     1973            call writediagfi(ngrid,"GND","heat flux from ground","W m-2",2,fluxgrd)
     1974         endif
     1975         
     1976         call writediagfi(ngrid,"DYN","dynamical heat input","W m-2",2,fluxdyn)
     1977         
     1978      endif ! end of 'callrad'
    20721979       
    2073         if(enertest) then
    2074           if (calldifv) then
    2075              call writediagfi(ngrid,"q2","turbulent kinetic energy","J.kg^-1",3,q2)
     1980      if(enertest) then
     1981     
     1982         if (calldifv) then
     1983         
     1984            call writediagfi(ngrid,"q2","turbulent kinetic energy","J.kg^-1",3,q2)
     1985            call writediagfi(ngrid,"sensibFlux","sensible heat flux","w.m^-2",2,sensibFlux)
     1986           
    20761987!            call writediagfi(ngrid,"dEzdiff","turbulent diffusion heating (-sensible flux)","w.m^-2",3,dEzdiff)
    20771988!            call writediagfi(ngrid,"dEdiff","integrated turbulent diffusion heating (-sensible flux)","w.m^-2",2,dEdiff)
    20781989!            call writediagfi(ngrid,"dEdiffs","In TurbDiff (correc rad+latent heat) surf nrj change","w.m^-2",2,dEdiffs)
    2079              call writediagfi(ngrid,"sensibFlux","sensible heat flux","w.m^-2",2,sensibFlux)
    2080           endif
    2081           if (corrk) then
    2082              call writediagfi(ngrid,"dEzradsw","radiative heating","w.m^-2",3,dEzradsw)
    2083              call writediagfi(ngrid,"dEzradlw","radiative heating","w.m^-2",3,dEzradlw)
    2084           endif
    2085           if(watercond) then
     1990             
     1991         endif
     1992         
     1993         if (corrk) then
     1994            call writediagfi(ngrid,"dEzradsw","radiative heating","w.m^-2",3,dEzradsw)
     1995            call writediagfi(ngrid,"dEzradlw","radiative heating","w.m^-2",3,dEzradlw)
     1996         endif
     1997         
     1998         if(watercond) then
     1999
     2000            call writediagfi(ngrid,"lscaledE","heat from largescale","W m-2",2,lscaledE)
     2001            call writediagfi(ngrid,"madjdE","heat from moistadj","W m-2",2,madjdE)
     2002            call writediagfi(ngrid,"qsatatm","atm qsat"," ",3,qsat)
     2003             
    20862004!            call writediagfi(ngrid,"lscaledEz","heat from largescale","W m-2",3,lscaledEz)
    2087 !            call writediagfi(ngrid,"madjdEz","heat from moistadj","W m-2",3,madjdEz)
    2088              call writediagfi(ngrid,"lscaledE","heat from largescale","W m-2",2,lscaledE)
    2089              call writediagfi(ngrid,"madjdE","heat from moistadj","W m-2",2,madjdE)
    2090              call writediagfi(ngrid,"qsatatm","atm qsat"," ",3,qsat)
     2005!            call writediagfi(ngrid,"madjdEz","heat from moistadj","W m-2",3,madjdEz)             
    20912006!            call writediagfi(ngrid,"h2o_max_col","maximum H2O column amount","kg.m^-2",2,H2Omaxcol)
    2092           endif
    2093         endif
    2094 
    2095 !     Temporary inclusions for heating diagnostics
    2096 !        call writediagfi(ngrid,"zdtdyn","Dyn. heating","T s-1",3,zdtdyn)
    2097         call writediagfi(ngrid,"zdtsw","SW heating","T s-1",3,zdtsw)
    2098         call writediagfi(ngrid,"zdtlw","LW heating","T s-1",3,zdtlw)
    2099         call writediagfi(ngrid,"dtrad","radiative heating","K s-1",3,dtrad)
    2100 
    2101         ! debugging
     2007
     2008         endif
     2009         
     2010      endif ! end of 'enertest'
     2011
     2012
     2013        ! Temporary inclusions for heating diagnostics.
     2014        !call writediagfi(ngrid,"zdtsw","SW heating","T s-1",3,zdtsw)
     2015        !call writediagfi(ngrid,"zdtlw","LW heating","T s-1",3,zdtlw)
     2016        !call writediagfi(ngrid,"dtrad","radiative heating","K s-1",3,dtrad)
     2017        ! call writediagfi(ngrid,"zdtdyn","Dyn. heating","T s-1",3,zdtdyn)
     2018       
     2019        ! For Debugging.
    21022020        !call writediagfi(ngrid,'rnat','Terrain type',' ',2,real(rnat))
    21032021        !call writediagfi(ngrid,'pphi','Geopotential',' ',3,pphi)
    2104 
    2105 !     Output aerosols
    2106         if (igcm_co2_ice.ne.0.and.iaero_co2.ne.0) &
    2107           call writediagfi(ngrid,'CO2ice_reff','CO2ice_reff','m',3,reffrad(1,1,iaero_co2))
    2108         if (igcm_h2o_ice.ne.0.and.iaero_h2o.ne.0) &
    2109           call writediagfi(ngrid,'H2Oice_reff','H2Oice_reff','m',3,reffrad(:,:,iaero_h2o))
    2110         if (igcm_co2_ice.ne.0.and.iaero_co2.ne.0) &
    2111           call writediagfi(ngrid,'CO2ice_reffcol','CO2ice_reffcol','um kg m^-2',2,reffcol(1,iaero_co2))
    2112         if (igcm_h2o_ice.ne.0.and.iaero_h2o.ne.0) &
    2113           call writediagfi(ngrid,'H2Oice_reffcol','H2Oice_reffcol','um kg m^-2',2,reffcol(1,iaero_h2o))
    2114 
    2115 !     Output tracers
    2116        if (tracer) then
    2117         do iq=1,nq
    2118           call writediagfi(ngrid,noms(iq),noms(iq),'kg/kg',3,zq(1,1,iq))
     2022       
     2023
     2024      ! Output aerosols.
     2025      if (igcm_co2_ice.ne.0.and.iaero_co2.ne.0) &
     2026         call writediagfi(ngrid,'CO2ice_reff','CO2ice_reff','m',3,reffrad(1,1,iaero_co2))
     2027      if (igcm_h2o_ice.ne.0.and.iaero_h2o.ne.0) &
     2028         call writediagfi(ngrid,'H2Oice_reff','H2Oice_reff','m',3,reffrad(:,:,iaero_h2o))
     2029      if (igcm_co2_ice.ne.0.and.iaero_co2.ne.0) &
     2030         call writediagfi(ngrid,'CO2ice_reffcol','CO2ice_reffcol','um kg m^-2',2,reffcol(1,iaero_co2))
     2031      if (igcm_h2o_ice.ne.0.and.iaero_h2o.ne.0) &
     2032         call writediagfi(ngrid,'H2Oice_reffcol','H2Oice_reffcol','um kg m^-2',2,reffcol(1,iaero_h2o))
     2033
     2034      ! Output tracers.
     2035      if (tracer) then
     2036     
     2037         do iq=1,nq
     2038            call writediagfi(ngrid,noms(iq),noms(iq),'kg/kg',3,zq(1,1,iq))
     2039            call writediagfi(ngrid,trim(noms(iq))//'_surf',trim(noms(iq))//'_surf',  &
     2040                             'kg m^-2',2,qsurf_hist(1,iq) )
     2041            call writediagfi(ngrid,trim(noms(iq))//'_col',trim(noms(iq))//'_col',    &
     2042                            'kg m^-2',2,qcol(1,iq) )
     2043                         
    21192044!          call writediagfi(ngrid,trim(noms(iq))//'_surf',trim(noms(iq))//'_surf',  &
    2120 !                          'kg m^-2',2,qsurf(1,iq) )
    2121           call writediagfi(ngrid,trim(noms(iq))//'_surf',trim(noms(iq))//'_surf',  &
    2122                           'kg m^-2',2,qsurf_hist(1,iq) )
    2123           call writediagfi(ngrid,trim(noms(iq))//'_col',trim(noms(iq))//'_col',    &
    2124                           'kg m^-2',2,qcol(1,iq) )
    2125 
    2126           if(watercond.or.CLFvarying)then
    2127              call writediagfi(ngrid,"rneb_man","H2O cloud fraction (conv)"," ",3,rneb_man)
    2128              call writediagfi(ngrid,"rneb_lsc","H2O cloud fraction (large scale)"," ",3,rneb_lsc)
    2129              call writediagfi(ngrid,"CLF","H2O cloud fraction"," ",3,cloudfrac)
    2130              call writediagfi(ngrid,"CLFt","H2O column cloud fraction"," ",2,totcloudfrac)
    2131              call writediagfi(ngrid,"RH","relative humidity"," ",3,RH)
    2132           endif
    2133 
    2134           if(waterrain)then
    2135              call writediagfi(ngrid,"rain","rainfall","kg m-2 s-1",2,zdqsrain)
    2136              call writediagfi(ngrid,"snow","snowfall","kg m-2 s-1",2,zdqssnow)
    2137           endif
    2138 
    2139           if((hydrology).and.(.not.ok_slab_ocean))then
    2140              call writediagfi(ngrid,"hice","oceanic ice height","m",2,hice)
    2141           endif
    2142 
    2143           if(ice_update)then
    2144              call writediagfi(ngrid,"ice_min","min annual ice","m",2,ice_min)
    2145              call writediagfi(ngrid,"ice_ini","initial annual ice","m",2,ice_initial)
    2146           endif
    2147 
    2148           do ig=1,ngrid
    2149              if(tau_col(ig).gt.1.e3)then
    2150 !                print*,'WARNING: tau_col=',tau_col(ig)
    2151 !                print*,'at ig=',ig,'in PHYSIQ'
    2152              endif
    2153           end do
    2154 
    2155           call writediagfi(ngrid,"tau_col","Total aerosol optical depth","[]",2,tau_col)
    2156 
    2157         enddo
    2158        endif
    2159 
    2160 !      output spectrum
     2045!                          'kg m^-2',2,qsurf(1,iq) )                         
     2046
     2047            if(watercond.or.CLFvarying)then
     2048               call writediagfi(ngrid,"rneb_man","H2O cloud fraction (conv)"," ",3,rneb_man)
     2049               call writediagfi(ngrid,"rneb_lsc","H2O cloud fraction (large scale)"," ",3,rneb_lsc)
     2050               call writediagfi(ngrid,"CLF","H2O cloud fraction"," ",3,cloudfrac)
     2051               call writediagfi(ngrid,"CLFt","H2O column cloud fraction"," ",2,totcloudfrac)
     2052               call writediagfi(ngrid,"RH","relative humidity"," ",3,RH)
     2053            endif
     2054
     2055            if(waterrain)then
     2056               call writediagfi(ngrid,"rain","rainfall","kg m-2 s-1",2,zdqsrain)
     2057               call writediagfi(ngrid,"snow","snowfall","kg m-2 s-1",2,zdqssnow)
     2058            endif
     2059
     2060            if((hydrology).and.(.not.ok_slab_ocean))then
     2061               call writediagfi(ngrid,"hice","oceanic ice height","m",2,hice)
     2062            endif
     2063
     2064            if(ice_update)then
     2065               call writediagfi(ngrid,"ice_min","min annual ice","m",2,ice_min)
     2066               call writediagfi(ngrid,"ice_ini","initial annual ice","m",2,ice_initial)
     2067            endif
     2068
     2069            do ig=1,ngrid
     2070               if(tau_col(ig).gt.1.e3)then
     2071                  print*,'WARNING: tau_col=',tau_col(ig)
     2072                  print*,'at ig=',ig,'in PHYSIQ'
     2073               endif         
     2074            end do
     2075
     2076            call writediagfi(ngrid,"tau_col","Total aerosol optical depth","[]",2,tau_col)
     2077
     2078         enddo ! end of 'nq' loop
     2079         
     2080       endif ! end of 'tracer'
     2081
     2082
     2083      ! Output spectrum.
    21612084      if(specOLR.and.corrk)then     
    21622085         call writediagspecIR(ngrid,"OLR3D","OLR(lon,lat,band)","W/m^2/cm^-1",3,OLR_nu)
    … …  
    21662089      icount=icount+1
    21672090
    2168       if (lastcall) then
    2169 
    2170         ! deallocate gas variables
    2171 !$OMP BARRIER
    2172 !$OMP MASTER
    2173         IF ( ALLOCATED( gnom ) ) DEALLOCATE( gnom )
    2174         IF ( ALLOCATED( gfrac ) ) DEALLOCATE( gfrac )  ! both allocated in su_gases.F90
    2175 !$OMP END MASTER
    2176 !$OMP BARRIER
    2177 
    2178         ! deallocate saved arrays
    2179         ! this is probably not that necessary
    2180         ! ... but probably a good idea to clean the place before we leave
    2181         IF ( ALLOCATED(tsurf)) DEALLOCATE(tsurf)
    2182         IF ( ALLOCATED(tsoil)) DEALLOCATE(tsoil)
    2183         IF ( ALLOCATED(albedo)) DEALLOCATE(albedo)
    2184         IF ( ALLOCATED(albedo0)) DEALLOCATE(albedo0)
    2185         IF ( ALLOCATED(rnat)) DEALLOCATE(rnat)
    2186         IF ( ALLOCATED(emis)) DEALLOCATE(emis)
    2187         IF ( ALLOCATED(dtrad)) DEALLOCATE(dtrad)
    2188         IF ( ALLOCATED(fluxrad_sky)) DEALLOCATE(fluxrad_sky)
    2189         IF ( ALLOCATED(fluxrad)) DEALLOCATE(fluxrad)
    2190         IF ( ALLOCATED(capcal)) DEALLOCATE(capcal)
    2191         IF ( ALLOCATED(fluxgrd)) DEALLOCATE(fluxgrd)
    2192         IF ( ALLOCATED(qsurf)) DEALLOCATE(qsurf)
    2193         IF ( ALLOCATED(q2)) DEALLOCATE(q2)
    2194         IF ( ALLOCATED(ztprevious)) DEALLOCATE(ztprevious)
    2195         IF ( ALLOCATED(hice)) DEALLOCATE(hice)
    2196         IF ( ALLOCATED(cloudfrac)) DEALLOCATE(cloudfrac)
    2197         IF ( ALLOCATED(totcloudfrac)) DEALLOCATE(totcloudfrac)
    2198         IF ( ALLOCATED(qsurf_hist)) DEALLOCATE(qsurf_hist)
    2199         IF ( ALLOCATED(reffrad)) DEALLOCATE(reffrad)
    2200         IF ( ALLOCATED(nueffrad)) DEALLOCATE(nueffrad)
    2201         IF ( ALLOCATED(ice_initial)) DEALLOCATE(ice_initial)
    2202         IF ( ALLOCATED(ice_min)) DEALLOCATE(ice_min)
    2203 
    2204         IF ( ALLOCATED(fluxsurf_lw)) DEALLOCATE(fluxsurf_lw)
    2205         IF ( ALLOCATED(fluxsurf_sw)) DEALLOCATE(fluxsurf_sw)
    2206         IF ( ALLOCATED(fluxtop_lw)) DEALLOCATE(fluxtop_lw)
    2207         IF ( ALLOCATED(fluxabs_sw)) DEALLOCATE(fluxabs_sw)
    2208         IF ( ALLOCATED(fluxtop_dn)) DEALLOCATE(fluxtop_dn)
    2209         IF ( ALLOCATED(fluxdyn)) DEALLOCATE(fluxdyn)
    2210         IF ( ALLOCATED(OLR_nu)) DEALLOCATE(OLR_nu)
    2211         IF ( ALLOCATED(OSR_nu)) DEALLOCATE(OSR_nu)
    2212         IF ( ALLOCATED(sensibFlux)) DEALLOCATE(sensibFlux)
    2213         IF ( ALLOCATED(zdtlw)) DEALLOCATE(zdtlw)
    2214         IF ( ALLOCATED(zdtsw)) DEALLOCATE(zdtsw)
    2215         IF ( ALLOCATED(tau_col)) DEALLOCATE(tau_col)
    2216         IF ( ALLOCATED(pctsrf_sic)) DEALLOCATE(pctsrf_sic)
    2217         IF ( ALLOCATED(tslab)) DEALLOCATE(tslab)
    2218         IF ( ALLOCATED(tsea_ice)) DEALLOCATE(tsea_ice)
    2219         IF ( ALLOCATED(sea_ice)) DEALLOCATE(sea_ice)
    2220         IF ( ALLOCATED(zmasq)) DEALLOCATE(zmasq)
    2221         IF ( ALLOCATED(knindex)) DEALLOCATE(knindex)
    2222 
    2223         !! this is defined in comsaison_h
    2224         IF ( ALLOCATED(mu0)) DEALLOCATE(mu0)
    2225         IF ( ALLOCATED(fract)) DEALLOCATE(fract)
    2226 
    2227         !! this is defined in comsoil_h
    2228         IF ( ALLOCATED(layer)) DEALLOCATE(layer)
    2229         IF ( ALLOCATED(mlayer)) DEALLOCATE(mlayer)
    2230         IF ( ALLOCATED(inertiedat)) DEALLOCATE(inertiedat)
    2231 
    2232         !! this is defined in surfdat_h
    2233         IF ( ALLOCATED(albedodat)) DEALLOCATE(albedodat)
    2234         IF ( ALLOCATED(phisfi)) DEALLOCATE(phisfi)
    2235         IF ( ALLOCATED(zmea)) DEALLOCATE(zmea)
    2236         IF ( ALLOCATED(zstd)) DEALLOCATE(zstd)
    2237         IF ( ALLOCATED(zsig)) DEALLOCATE(zsig)
    2238         IF ( ALLOCATED(zgam)) DEALLOCATE(zgam)
    2239         IF ( ALLOCATED(zthe)) DEALLOCATE(zthe)
    2240         IF ( ALLOCATED(dryness)) DEALLOCATE(dryness)
    2241         IF ( ALLOCATED(watercaptag)) DEALLOCATE(watercaptag)
    2242  
    2243         !! this is defined in tracer_h
    2244         IF ( ALLOCATED(noms)) DEALLOCATE(noms)
    2245         IF ( ALLOCATED(mmol)) DEALLOCATE(mmol)
    2246         IF ( ALLOCATED(radius)) DEALLOCATE(radius)
    2247         IF ( ALLOCATED(rho_q)) DEALLOCATE(rho_q)
    2248         IF ( ALLOCATED(qext)) DEALLOCATE(qext)
    2249         IF ( ALLOCATED(alpha_lift)) DEALLOCATE(alpha_lift)
    2250         IF ( ALLOCATED(alpha_devil)) DEALLOCATE(alpha_devil)
    2251         IF ( ALLOCATED(qextrhor)) DEALLOCATE(qextrhor)
    2252         IF ( ALLOCATED(igcm_dustbin)) DEALLOCATE(igcm_dustbin)
    2253 
    2254         !! this is defined in comgeomfi_h
    2255         IF ( ALLOCATED(lati)) DEALLOCATE(lati)
    2256         IF ( ALLOCATED(long)) DEALLOCATE(long)
    2257         IF ( ALLOCATED(area)) DEALLOCATE(area)
    2258         IF ( ALLOCATED(sinlat)) DEALLOCATE(sinlat)
    2259         IF ( ALLOCATED(coslat)) DEALLOCATE(coslat)
    2260         IF ( ALLOCATED(sinlon)) DEALLOCATE(sinlon)
    2261         IF ( ALLOCATED(coslon)) DEALLOCATE(coslon)
    2262       endif
    2263 
    2264       return
    22652091    end subroutine physiq

  • trunk/LMDZ.GENERIC/libf/phystd/turbdiff.F90

    r1397 r1477  
    33          pplay,pplev,pzlay,pzlev,pz0,                 &
    44          pu,pv,pt,ppopsk,pq,ptsrf,pemis,pqsurf,       &
    5           pdufi,pdvfi,pdtfi,pdqfi,pfluxsrf,            &
     5          pdtfi,pdqfi,pfluxsrf,            &
    66          Pdudif,pdvdif,pdtdif,pdtsrf,sensibFlux,pq2,  &
    77          pdqdif,pdqevap,pdqsdif,flux_u,flux_v,lastcall)
    … …  
    5454      REAL,INTENT(IN) :: ptsrf(ngrid) ! surface temperature (K)
    5555      REAL,INTENT(IN) :: pemis(ngrid)
    56       REAL,INTENT(IN) :: pdufi(ngrid,nlay),pdvfi(ngrid,nlay)
    5756      REAL,INTENT(IN) :: pdtfi(ngrid,nlay)
    5857      REAL,INTENT(IN) :: pfluxsrf(ngrid)
    … …  
    189188      DO ilev=1,nlay
    190189         DO ig=1,ngrid
    191             zu(ig,ilev)=pu(ig,ilev)+pdufi(ig,ilev)*ptimestep
    192             zv(ig,ilev)=pv(ig,ilev)+pdvfi(ig,ilev)*ptimestep
     190            zu(ig,ilev)=pu(ig,ilev)
     191            zv(ig,ilev)=pv(ig,ilev)
    193192            zt(ig,ilev)=pt(ig,ilev)+pdtfi(ig,ilev)*ptimestep
    194193            zh(ig,ilev)=pt(ig,ilev)*zovExner(ig,ilev) !for call vdif_kc, but could be moved and computed there
    … …  
    689688      do ilev = 1, nlay
    690689         do ig=1,ngrid
    691             pdudif(ig,ilev)=(zu(ig,ilev)-(pu(ig,ilev)+pdufi(ig,ilev)*ptimestep))/ptimestep
    692             pdvdif(ig,ilev)=(zv(ig,ilev)-(pv(ig,ilev)+pdvfi(ig,ilev)*ptimestep))/ptimestep
     690            pdudif(ig,ilev)=(zu(ig,ilev)-(pu(ig,ilev)))/ptimestep
     691            pdvdif(ig,ilev)=(zv(ig,ilev)-(pv(ig,ilev)))/ptimestep
    693692            pdtdif(ig,ilev)=( zt(ig,ilev)- pt(ig,ilev))/ptimestep-pdtfi(ig,ilev)
    694693         enddo

  • trunk/LMDZ.GENERIC/libf/phystd/vdifc.F

    r1397 r1477  
    33     &     pplay,pplev,pzlay,pzlev,pz0,
    44     &     pu,pv,ph,pq,ptsrf,pemis,pqsurf,
    5      &     pdufi,pdvfi,pdhfi,pdqfi,pfluxsrf,
     5     &     pdhfi,pdqfi,pfluxsrf,
    66     &     pdudif,pdvdif,pdhdif,pdtsrf,sensibFlux,pq2,
    77     &     pdqdif,pdqsdif,lastcall)
    … …  
    4949      REAL pu(ngrid,nlay),pv(ngrid,nlay),ph(ngrid,nlay)
    5050      REAL ptsrf(ngrid),pemis(ngrid)
    51       REAL pdufi(ngrid,nlay),pdvfi(ngrid,nlay),pdhfi(ngrid,nlay)
     51      REAL pdhfi(ngrid,nlay)
    5252      REAL pfluxsrf(ngrid)
    5353      REAL pdudif(ngrid,nlay),pdvdif(ngrid,nlay),pdhdif(ngrid,nlay)
    … …  
    205205      DO ilev=1,nlay
    206206         DO ig=1,ngrid
    207             zu(ig,ilev)=pu(ig,ilev)+pdufi(ig,ilev)*ptimestep
    208             zv(ig,ilev)=pv(ig,ilev)+pdvfi(ig,ilev)*ptimestep
     207            zu(ig,ilev)=pu(ig,ilev)
     208            zv(ig,ilev)=pv(ig,ilev)
    209209            zh(ig,ilev)=ph(ig,ilev)+pdhfi(ig,ilev)*ptimestep
    210210         ENDDO
    … …  
    662662         do ig=1,ngrid
    663663            pdudif(ig,ilev)=(zu(ig,ilev)-
    664      &           (pu(ig,ilev)+pdufi(ig,ilev)*ptimestep))/ptimestep
     664     &           (pu(ig,ilev)))/ptimestep
    665665            pdvdif(ig,ilev)=(zv(ig,ilev)-
    666      &           (pv(ig,ilev)+pdvfi(ig,ilev)*ptimestep))/ptimestep
     666     &           (pv(ig,ilev)))/ptimestep
    667667            hh = ph(ig,ilev)+pdhfi(ig,ilev)*ptimestep
    668668
Note: See TracChangeset for help on using the changeset viewer.