Changeset 3232 for trunk/LMDZ.GENERIC

Timestamp:

Feb 21, 2024, 5:45:11 PM (10 months ago)

Author:

afalco

Message:

Pluto PCM:
1D functional
AF

File:

• trunk/LMDZ.GENERIC/libf/phystd/physiq_mod.F90 (modified) (140 diffs)

trunk/LMDZ.GENERIC/libf/phystd/physiq_mod.F90

@@ -1 +1 @@
       module physiq_mod
-      
+
       implicit none
-      
+
       contains
-      
+
       subroutine physiq(ngrid,nlayer,nq,   &
                   firstcall,lastcall,      &
@@ -14 +14 @@
 
 !!
-      use write_field_phy, only: Writefield_phy  
+      use write_field_phy, only: Writefield_phy
 !!
-      use ioipsl_getin_p_mod, only: getin_p 
+      use ioipsl_getin_p_mod, only: getin_p
       use radinc_h, only : L_NSPECTI,L_NSPECTV,naerkind, corrkdir, banddir
       use watercommon_h, only : RLVTT, Psat_water,epsi,su_watercycle, RV, T_h2o_ice_liq
@@ -80 +80 @@
       use turb_mod, only : q2,sensibFlux,turb_resolved
       use mass_redistribution_mod, only: mass_redistribution
-      use condensation_generic_mod, only: condensation_generic 
+      use condensation_generic_mod, only: condensation_generic
       use datafile_mod, only: datadir
 #ifndef MESOSCALE
@@ -98 +98 @@
 #endif
 
-#ifdef CPP_XIOS      
+#ifdef CPP_XIOS
       use xios_output_mod, only: initialize_xios_output, &
                                  update_xios_timestep, &
@@ -109 +109 @@
 
 !==================================================================
-!     
+!
 !     Purpose
 !     -------
@@ -199 +199 @@
 !           Frederic Hourdin        15/10/93
 !           Francois Forget        1994
-!           Christophe Hourdin        02/1997 
+!           Christophe Hourdin        02/1997
 !           Subroutine completely rewritten by F. Forget (01/2000)
 !           Water ice clouds: Franck Montmessin (update 06/2003)
 !           Radiatively active tracers: J.-B. Madeleine (10/2008-06/2009)
 !           New correlated-k radiative scheme: R. Wordsworth (2009)
-!           Many specifically Martian subroutines removed: R. Wordsworth (2009)       
+!           Many specifically Martian subroutines removed: R. Wordsworth (2009)
 !           Improved water cycle: R. Wordsworth / B. Charnay (2010)
 !           To F90: R. Wordsworth (2010)
@@ -229 +229 @@
       integer,intent(in) :: nlayer            ! Number of atmospheric layers.
       integer,intent(in) :: nq                ! Number of tracers.
-      
+
       logical,intent(in) :: firstcall ! Signals first call to physics.
       logical,intent(in) :: lastcall  ! Signals last call to physics.
-      
+
       real,intent(in) :: pday                  ! Number of elapsed sols since reference Ls=0.
       real,intent(in) :: ptime                 ! "Universal time", given as fraction of sol (e.g.: 0.5 for noon).
@@ -265 +265 @@
 ! -----------------
 
-!     Aerosol (dust or ice) extinction optical depth  at reference wavelength 
+!     Aerosol (dust or ice) extinction optical depth  at reference wavelength
 !     for the "naerkind" optically active aerosols:
 
@@ -276 +276 @@
       integer l,ig,ierr,iq,nw,isoil,iesp, igcm_generic_vap, igcm_generic_ice
       logical call_ice_vap_generic ! to call only one time the ice/vap pair of a tracer
-      
+
       real zls                       ! Solar longitude (radians).
       real zlss                      ! Sub solar point longitude (radians).
@@ -284 +284 @@
 
 ! VARIABLES for the thermal plume model
-      
+
       real f(ngrid)                       ! Mass flux norm
       real fm(ngrid,nlayer+1)             ! Mass flux
@@ -297 +297 @@
       real zw2(ngrid,nlayer+1)            ! Vertical speed
       real zw2_bis(ngrid,nlayer)          ! Recasted zw2
-      
+
 ! TENDENCIES due to various processes :
 
@@ -306 +306 @@
       real zdtsdif(ngrid)     ! Turbdiff/vdifc routines.
       real zdtsurf_hyd(ngrid) ! Hydrol routine.
-            
-      ! For Atmospheric Temperatures : (K/s)    
+
+      ! For Atmospheric Temperatures : (K/s)
       real dtlscale(ngrid,nlayer)                             ! Largescale routine.
-      real dt_generic_condensation(ngrid,nlayer)              ! condensation_generic routine. 
+      real dt_generic_condensation(ngrid,nlayer)              ! condensation_generic routine.
       real zdtc(ngrid,nlayer)                                 ! Condense_co2 routine.
       real zdtdif(ngrid,nlayer)                               ! Turbdiff/vdifc routines.
@@ -320 +320 @@
       real zdtsw1(ngrid,nlayer), zdtlw1(ngrid,nlayer)         ! Callcorrk routine.
       real zdtchim(ngrid,nlayer)                              ! Calchim routine.
-                              
+
       ! For Surface Tracers : (kg/m2/s)
       real dqsurf(ngrid,nq)                 ! Cumulated tendencies.
@@ -332 +332 @@
       real reevap_precip(ngrid)             ! re-evaporation flux of precipitation (integrated over the atmospheric column)
       real reevap_precip_generic(ngrid,nq)
-                  
+
       ! For Tracers : (kg/kg_of_air/s)
       real zdqc(ngrid,nlayer,nq)      ! Condense_co2 routine.
@@ -355 +355 @@
       REAL array_zero1(ngrid)
       REAL array_zero2(ngrid,nlayer)
-                        
+
       ! For Winds : (m/s/s)
       real zdvadj(ngrid,nlayer), zduadj(ngrid,nlayer)       ! Convadj routine.
@@ -363 +363 @@
       real zdhdif(ngrid,nlayer)                             ! Turbdiff/vdifc routines.
       real zdhadj(ngrid,nlayer)                             ! Convadj routine.
-     
+
       ! For Pressure and Mass :
       real zdmassmr(ngrid,nlayer) ! Atmospheric Mass tendency for mass_redistribution (kg_of_air/m2/s).
@@ -369 +369 @@
       real zdpsrfmr(ngrid)        ! Pressure tendency for mass_redistribution routine (Pa/s).
 
-      
-    
+
+
 ! Local variables for LOCAL CALCULATIONS:
 ! ---------------------------------------
@@ -395 +395 @@
       real vmr(ngrid,nlayer)                        ! volume mixing ratio
       real time_phys
-      
+
       real ISR,ASR,OLR,GND,DYN,GSR,Ts1,Ts2,Ts3,TsS ! for Diagnostic.
-      
+
       real qcol(ngrid,nq) ! Tracer Column Mass (kg/m2).
 
@@ -410 +410 @@
       real,save :: dEtotSW, dEtotsSW, dEtotLW, dEtotsLW
 !$OMP THREADPRIVATE(dEtotSW, dEtotsSW, dEtotLW, dEtotsLW)
-      
+
 !JL12 conservation test for mean flow kinetic energy has been disabled temporarily
 
-      real dtmoist_max,dtmoist_min     
+      real dtmoist_max,dtmoist_min
       real dItot, dItot_tmp, dVtot, dVtot_tmp
 
@@ -437 +437 @@
 
       logical clearsky ! For double radiative transfer call. By BC
-      
+
       ! For Clear Sky Case.
       real fluxsurfabs_sw1(ngrid)  ! For SW/LW flux.
@@ -466 +466 @@
       real :: flux_sens_lat(ngrid)
       real :: qsurfint(ngrid,nq)
-#ifdef MESOSCALE 
+#ifdef MESOSCALE
       REAL :: lsf_dt(nlayer)
       REAL :: lsf_dq(nlayer)
@@ -474 +474 @@
       logical, save ::  check_physics_inputs=.false.
       logical, save ::  check_physics_outputs=.false.
-!$OPM THREADPRIVATE(check_physics_inputs,check_physics_outputs)     
+!$OPM THREADPRIVATE(check_physics_inputs,check_physics_outputs)
 
       ! Misc
       character*2 :: str2
-      character(len=10) :: tmp1 
+      character(len=10) :: tmp1
       character(len=10) :: tmp2
 !==================================================================================================
@@ -534 +534 @@
 !        Read 'startfi.nc' file.
 !        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-#ifndef MESOSCALE 
+#ifndef MESOSCALE
          call phyetat0(startphy_file,                                 &
                        ngrid,nlayer,"startfi.nc",0,0,nsoilmx,nq,      &
@@ -541 +541 @@
                        rnat,pctsrf_sic,tslab, tsea_ice,sea_ice)
 
-#else   
+#else
 
          day_ini = pday
@@ -555 +555 @@
              tsoil(1:ngrid,isoil)=tsurf(1:ngrid)
            enddo
-           if (is_master) write(*,*) "Physiq: initializing day_ini to pdat !"
+           if (is_master) write(*,*) "Physiq: initializing day_ini to pday !"
            day_ini=pday
          endif
@@ -569 +569 @@
          write (*,*) 'In physiq day_ini =', day_ini
 
-!        Initialize albedo calculation. 
+!        Initialize albedo calculation.
 !        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
          albedo(:,:)=0.0
@@ -575 +575 @@
          albedo_snow_SPECTV(:)=0.0
          albedo_co2_ice_SPECTV(:)=0.0
-         call surfini(ngrid,nq,qsurf,albedo,albedo_bareground,albedo_snow_SPECTV,albedo_co2_ice_SPECTV) 
-
-!        Initialize orbital calculation. 
+         call surfini(ngrid,nq,qsurf,albedo,albedo_bareground,albedo_snow_SPECTV,albedo_co2_ice_SPECTV)
+
+!        Initialize orbital calculation.
 !        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
          call iniorbit(apoastr,periastr,year_day,peri_day,obliquit)
@@ -596 +596 @@
 
 !!           call ini_surf_heat_transp_mod()
-           
+
            knindex(:) = 0
            knon = 0
@@ -620 +620 @@
 !        ~~~~~~~~~~~~~~~~
          if (callsoil) then
-         
+
             call soil(ngrid,nsoilmx,firstcall,lastcall,inertiedat, &
                       ptimestep,tsurf,tsoil,capcal,fluxgrd)
@@ -639 +639 @@
 
          else ! else of 'callsoil'.
-         
+
             print*,'WARNING! Thermal conduction in the soil turned off'
             capcal(:)=1.e6
             fluxgrd(:)=intheat
             print*,'Flux from ground = ',intheat,' W m^-2'
-            
+
          endif ! end of 'callsoil'.
-         
+
          icount=1
 
@@ -653 +653 @@
          ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
          if (.not.ok_slab_ocean) then
-         
+
            rnat(:)=1.
            if (.not.nosurf) then ! inertiedat only defined if there is a surface
@@ -691 +691 @@
          endif
 
-         if(meanOLR)then          
-            call system('rm -f rad_bal.out') ! to record global radiative balance.           
-            call system('rm -f tem_bal.out') ! to record global mean/max/min temperatures.            
+         if(meanOLR)then
+            call system('rm -f rad_bal.out') ! to record global radiative balance.
+            call system('rm -f tem_bal.out') ! to record global mean/max/min temperatures.
             call system('rm -f h2o_bal.out') ! to record global hydrological balance.
          endif
 
 
-         watertest=.false.        
+         watertest=.false.
          if(water)then ! Initialize variables for water cycle
-            
+
             if(enertest)then
                watertest = .true.
@@ -711 +711 @@
          metallicity=0.0 ! default value --- is not used here but necessary to call function Psat_generic
          call getin_p("metallicity",metallicity) ! --- is not used here but necessary to call function Psat_generic
-         
+
 !        Set some parameters for the thermal plume model
 !        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -731 +731 @@
 
 !!         call WriteField_phy("post_physdem_qsurf",qsurf(1:ngrid,igcm_h2o_vap),1)
-#endif         
-         if (corrk) then 
-            ! We initialise the spectral grid here instead of 
-            ! at firstcall of callcorrk so we can output XspecIR, XspecVI 
+#endif
+         if (corrk) then
+            ! We initialise the spectral grid here instead of
+            ! at firstcall of callcorrk so we can output XspecIR, XspecVI
             ! when using Dynamico
             print*, "physiq_mod: Correlated-k data base folder:",trim(datadir)
@@ -757 +757 @@
                                      year_day,presnivs,pseudoalt,WNOI,WNOV)
 #endif
-         
+
 !!         call WriteField_phy("post_xios_qsurf",qsurf(1:ngrid,igcm_h2o_vap),1)
 
@@ -764 +764 @@
 
 !!      call WriteField_phy("post_firstcall_qsurf",qsurf(1:ngrid,igcm_h2o_vap),1)
-!!      call writediagfi(ngrid,"firstcall_post_qsurf"," "," ",2,qsurf(1:ngrid,igcm_h2o_vap))  
+!!      call writediagfi(ngrid,"firstcall_post_qsurf"," "," ",2,qsurf(1:ngrid,igcm_h2o_vap))
 
       if (check_physics_inputs) then
@@ -778 +778 @@
 ! ------------------------------------------------------
 
-#ifdef CPP_XIOS      
+#ifdef CPP_XIOS
       ! update XIOS time/calendar
       call update_xios_timestep
-#endif      
+#endif
 
       ! Initialize various variables
       ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-      
+
       if ( .not.nearco2cond ) then
          pdt(1:ngrid,1:nlayer) = 0.0
-      endif      
+      endif
       zdtsurf(1:ngrid)      = 0.0
       pdq(1:ngrid,1:nlayer,1:nq) = 0.0
@@ -795 +795 @@
       pdv(1:ngrid,1:nlayer) = 0.0
       pdpsrf(1:ngrid)       = 0.0
-      zflubid(1:ngrid)      = 0.0      
+      zflubid(1:ngrid)      = 0.0
       flux_sens_lat(1:ngrid) = 0.0
       taux(1:ngrid) = 0.0
@@ -811 +811 @@
 
       call orbite(zls,dist_star,declin,right_ascen)
-      
+
       if (tlocked) then
               zlss=Mod(-(2.*pi*(zday/year_day)*nres - right_ascen),2.*pi)
@@ -818 +818 @@
       else if(diurnal .eqv. .false.) then
               zlss=9999.
-      endif 
+      endif
 
 
       ! Compute variations of g with latitude (oblate case).
       ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-      if (oblate .eqv. .false.) then      
-         glat(:) = g         
-      else if (flatten .eq. 0.0 .or. J2 .eq. 0.0 .or. Rmean .eq. 0.0 .or. MassPlanet .eq. 0.0) then      
+      if (oblate .eqv. .false.) then
+         glat(:) = g
+      else if (flatten .eq. 0.0 .or. J2 .eq. 0.0 .or. Rmean .eq. 0.0 .or. MassPlanet .eq. 0.0) then
          print*,'I need values for flatten, J2, Rmean and MassPlanet to compute glat (else set oblate=.false.)'
          call abort
@@ -831 +831 @@
          gmplanet = MassPlanet*grav*1e24
          do ig=1,ngrid
-            glat(ig)= gmplanet/(Rmean**2) * (1.D0 + 0.75 *J2 - 2.0*flatten/3. + (2.*flatten - 15./4.* J2) * cos(2. * (pi/2. - latitude(ig)))) 
+            glat(ig)= gmplanet/(Rmean**2) * (1.D0 + 0.75 *J2 - 2.0*flatten/3. + (2.*flatten - 15./4.* J2) * cos(2. * (pi/2. - latitude(ig))))
          end do
       endif
@@ -838 +838 @@
       ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
       zzlay(1:ngrid,1:nlayer)=pphi(1:ngrid,1:nlayer)
-      do l=1,nlayer         
+      do l=1,nlayer
          zzlay(1:ngrid,l)= zzlay(1:ngrid,l)/glat(1:ngrid)
       enddo
@@ -850 +850 @@
             zzlev(ig,l)=(z1*zzlay(ig,l-1)+z2*zzlay(ig,l))/(z1+z2)
          enddo
-      enddo     
+      enddo
 
       !Altitude of top interface (nlayer+1), using the thicknesss of the level below the top one. LT22
-      
+
       zzlev(1:ngrid,nlayer+1) = 2*zzlev(1:ngrid,nlayer)-zzlev(1:ngrid,nlayer-1)
 
@@ -859 +859 @@
       ! Note : Potential temperature calculation may not be the same in physiq and dynamic...
       ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-      do l=1,nlayer         
+      do l=1,nlayer
          do ig=1,ngrid
             zpopsk(ig,l)=(pplay(ig,l)/pplev(ig,1))**rcp
@@ -894 +894 @@
       if(photochem) then
          call concentrations(ngrid,nlayer,nq,pplay,pt,pdt,pq,pdq,ptimestep)
-      endif 
+      endif
 #endif
 
@@ -928 +928 @@
          endif
 
-      endif 
+      endif
 
       if (callrad) then
          if( mod(icount-1,iradia).eq.0.or.lastcall) then
-           
-            ! Eclipse incoming sunlight (e.g. Saturn ring shadowing).       
+
+            ! Eclipse incoming sunlight (e.g. Saturn ring shadowing).
             if(rings_shadow) then
                 call call_rings(ngrid, ptime, pday, diurnal)
-            endif   
+            endif
 
 !!            call writediagfi(ngrid,"corrk_pre_dqsurf"," "," ",2,dqsurf(1:ngrid,igcm_h2o_vap))
-!!            call writediagfi(ngrid,"corrk_pre_qsurf"," "," ",2,qsurf(1:ngrid,igcm_h2o_vap)) 
+!!            call writediagfi(ngrid,"corrk_pre_qsurf"," "," ",2,qsurf(1:ngrid,igcm_h2o_vap))
 
 
             if (corrk) then
-            
+
 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 ! II.a Call correlated-k radiative transfer scheme
@@ -952 +952 @@
                endif
                if(water) then
-                  muvar(1:ngrid,1:nlayer)=mugaz/(1.e0+(1.e0/epsi-1.e0)*pq(1:ngrid,1:nlayer,igcm_h2o_vap)) 
-                  muvar(1:ngrid,nlayer+1)=mugaz/(1.e0+(1.e0/epsi-1.e0)*pq(1:ngrid,nlayer,igcm_h2o_vap)) 
+                  muvar(1:ngrid,1:nlayer)=mugaz/(1.e0+(1.e0/epsi-1.e0)*pq(1:ngrid,1:nlayer,igcm_h2o_vap))
+                  muvar(1:ngrid,nlayer+1)=mugaz/(1.e0+(1.e0/epsi-1.e0)*pq(1:ngrid,nlayer,igcm_h2o_vap))
                   ! take into account water effect on mean molecular weight
                elseif(generic_condensation) then
@@ -959 +959 @@
 
                      call generic_tracer_index(nq,iq,igcm_generic_vap,igcm_generic_ice,call_ice_vap_generic)
-                     
+
                      if (call_ice_vap_generic) then ! to call only one time the ice/vap pair of a tracer
 
                         epsi_generic=constants_epsi_generic(iq)
 
-                        muvar(1:ngrid,1:nlayer)=mugaz/(1.e0+(1.e0/epsi_generic-1.e0)*pq(1:ngrid,1:nlayer,igcm_generic_vap)) 
-                        muvar(1:ngrid,nlayer+1)=mugaz/(1.e0+(1.e0/epsi_generic-1.e0)*pq(1:ngrid,nlayer,igcm_generic_vap))  
+                        muvar(1:ngrid,1:nlayer)=mugaz/(1.e0+(1.e0/epsi_generic-1.e0)*pq(1:ngrid,1:nlayer,igcm_generic_vap))
+                        muvar(1:ngrid,nlayer+1)=mugaz/(1.e0+(1.e0/epsi_generic-1.e0)*pq(1:ngrid,nlayer,igcm_generic_vap))
 
                      endif
-                  end do ! do iq=1,nq loop on tracers 
+                  end do ! do iq=1,nq loop on tracers
                ! take into account generic condensable specie (GCS) effect on mean molecular weight
-               
+
                else
                   muvar(1:ngrid,1:nlayer+1)=mugaz
-               endif 
-      
+               endif
+
 
                if(ok_slab_ocean) then
-                  tsurf2(:)=tsurf(:) 
+                  tsurf2(:)=tsurf(:)
                   do ig=1,ngrid
                      if (nint(rnat(ig))==0) then
@@ -984 +984 @@
                   enddo
                endif !(ok_slab_ocean)
-                
+
                ! standard callcorrk
                clearsky=.false.
@@ -995 +995 @@
                               int_dtaui,int_dtauv,                                &
                               tau_col,cloudfrac,totcloudfrac,                     &
-                              clearsky,firstcall,lastcall)     
+                              clearsky,firstcall,lastcall)
 
                 !GG (feb2021): Option to "artificially" decrease the raditive time scale in
                 !the deep atmosphere  press > 0.1 bar. Suggested by MT.
-                !! COEFF_RAD to be "tuned" to facilitate convergence of tendency 
-        
+                !! COEFF_RAD to be "tuned" to facilitate convergence of tendency
+
                 !coeff_rad=0.   ! 0 values, it doesn't accelerate the convergence
                 !coeff_rad=0.5
-                !coeff_rad=1.                 
+                !coeff_rad=1.
                 !do l=1, nlayer
                 !  do ig=1,ngrid
@@ -1013 +1013 @@
                 !enddo
 
-                ! Option to call scheme once more for clear regions  
+                ! Option to call scheme once more for clear regions
                if(CLFvarying)then
 
@@ -1027 +1027 @@
                                  tau_col1,cloudfrac,totcloudfrac,                    &
                                  clearsky,firstcall,lastcall)
-                  clearsky = .false.  ! just in case.     
+                  clearsky = .false.  ! just in case.
 
                   ! Sum the fluxes and heating rates from cloudy/clear cases
                   do ig=1,ngrid
                      tf=totcloudfrac(ig)
-                     ntf=1.-tf                    
+                     ntf=1.-tf
                      fluxsurf_lw(ig)       = ntf*fluxsurf_lw1(ig)       + tf*fluxsurf_lw(ig)
                      fluxsurf_sw(ig)       = ntf*fluxsurf_sw1(ig)       + tf*fluxsurf_sw(ig)
@@ -1040 +1040 @@
                      fluxabs_sw(ig)        = ntf*fluxabs_sw1(ig)        + tf*fluxabs_sw(ig)
                      tau_col(ig)           = ntf*tau_col1(ig)           + tf*tau_col(ig)
-                    
+
                      zdtlw(ig,1:nlayer) = ntf*zdtlw1(ig,1:nlayer) + tf*zdtlw(ig,1:nlayer)
                      zdtsw(ig,1:nlayer) = ntf*zdtsw1(ig,1:nlayer) + tf*zdtsw(ig,1:nlayer)
 
-                     OSR_nu(ig,1:L_NSPECTV) = ntf*OSR_nu1(ig,1:L_NSPECTV) + tf*OSR_nu(ig,1:L_NSPECTV)       
-                     GSR_nu(ig,1:L_NSPECTV) = ntf*GSR_nu1(ig,1:L_NSPECTV) + tf*GSR_nu(ig,1:L_NSPECTV)                  
-                     OLR_nu(ig,1:L_NSPECTI) = ntf*OLR_nu1(ig,1:L_NSPECTI) + tf*OLR_nu(ig,1:L_NSPECTI)                       
+                     OSR_nu(ig,1:L_NSPECTV) = ntf*OSR_nu1(ig,1:L_NSPECTV) + tf*OSR_nu(ig,1:L_NSPECTV)
+                     GSR_nu(ig,1:L_NSPECTV) = ntf*GSR_nu1(ig,1:L_NSPECTV) + tf*GSR_nu(ig,1:L_NSPECTV)
+                     OLR_nu(ig,1:L_NSPECTI) = ntf*OLR_nu1(ig,1:L_NSPECTI) + tf*OLR_nu(ig,1:L_NSPECTI)
                      if (diagdtau) then
-                       int_dtauv(ig,:,1:L_NSPECTV) = ntf*int_dtauv1(ig,:,1:L_NSPECTV) + tf*int_dtauv(ig,:,1:L_NSPECTV)                       
-                       int_dtaui(ig,:,1:L_NSPECTI) = ntf*int_dtaui1(ig,:,1:L_NSPECTI) + tf*int_dtaui(ig,:,1:L_NSPECTI)                       
+                       int_dtauv(ig,:,1:L_NSPECTV) = ntf*int_dtauv1(ig,:,1:L_NSPECTV) + tf*int_dtauv(ig,:,1:L_NSPECTV)
+                       int_dtaui(ig,:,1:L_NSPECTI) = ntf*int_dtaui1(ig,:,1:L_NSPECTI) + tf*int_dtaui(ig,:,1:L_NSPECTI)
                      endif
-                  enddo                               
+                  enddo
 
                endif ! end of CLFvarying.
@@ -1072 +1072 @@
                ! Net atmospheric radiative heating rate (K.s-1)
                dtrad(1:ngrid,1:nlayer)=zdtsw(1:ngrid,1:nlayer)+zdtlw(1:ngrid,1:nlayer)
-               
-               ! Late initialization of the Ice Spectral Albedo. We needed the visible bands to do that ! 
+
+               ! Late initialization of the Ice Spectral Albedo. We needed the visible bands to do that !
                if (firstcall .and. albedo_spectral_mode) then
                   call spectral_albedo_calc(albedo_snow_SPECTV)
@@ -1079 +1079 @@
 
             elseif(newtonian)then
-            
+
 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 ! II.b Call Newtonian cooling scheme
@@ -1089 +1089 @@
 
             else
-            
+
 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-! II.c Atmosphere has no radiative effect 
+! II.c Atmosphere has no radiative effect
 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
                fluxtop_dn(1:ngrid)  = fract(1:ngrid)*mu0(1:ngrid)*Fat1AU/dist_star**2
@@ -1100 +1100 @@
                print*,'------------WARNING---WARNING------------' ! by MT2015.
                print*,'You are in corrk=false mode, '
-               print*,'and the surface albedo is taken equal to the first visible spectral value'               
-               
+               print*,'and the surface albedo is taken equal to the first visible spectral value'
+
                fluxsurfabs_sw(1:ngrid) = fluxtop_dn(1:ngrid)*(1.-albedo(1:ngrid,1))
                fluxrad_sky(1:ngrid)    = fluxsurfabs_sw(1:ngrid)
@@ -1111 +1111 @@
 
          endif ! of if(mod(icount-1,iradia).eq.0)
-       
+
 
          ! Transformation of the radiative tendencies
@@ -1119 +1119 @@
          fluxrad(1:ngrid)=fluxrad_sky(1:ngrid)-zplanck(1:ngrid)
          pdt(1:ngrid,1:nlayer)=pdt(1:ngrid,1:nlayer)+dtrad(1:ngrid,1:nlayer)
-         
+
          ! Test of energy conservation
          !----------------------------
@@ -1152 +1152 @@
 
       if (calldifv) then
-      
+
          zflubid(1:ngrid)=fluxrad(1:ngrid)+fluxgrd(1:ngrid)
 
          ! JL12 the following if test is temporarily there to allow us to compare the old vdifc with turbdiff.
          if (UseTurbDiff) then
-         
+
             call turbdiff(ngrid,nlayer,nq,rnat,                  &
                           ptimestep,capcal,                      &
@@ -1168 +1168 @@
 
          else
-         
+
             zdh(1:ngrid,1:nlayer)=pdt(1:ngrid,1:nlayer)/zpopsk(1:ngrid,1:nlayer)
- 
+
             call vdifc(ngrid,nlayer,nq,rnat,zpopsk,           &
                        ptimestep,capcal,lwrite,               &
@@ -1199 +1199 @@
 !!         call writediagfi(ngrid,"vdifc_post_zdqsdif"," "," ",2,zdqsdif(1:ngrid,igcm_h2o_vap))
 
-         if (tracer) then 
+         if (tracer) then
            pdq(1:ngrid,1:nlayer,1:nq)=pdq(1:ngrid,1:nlayer,1:nq)+ zdqdif(1:ngrid,1:nlayer,1:nq)
            dqsurf(1:ngrid,1:nq)=dqsurf(1:ngrid,1:nq) + zdqsdif(1:ngrid,1:nq)
@@ -1210 +1210 @@
          !-------------------------
          if(enertest)then
-         
+
             dEzdiff(:,:)=cpp*mass(:,:)*zdtdif(:,:)
             do ig = 1, ngrid
@@ -1216 +1216 @@
                dEzdiff(ig,1)= dEzdiff(ig,1)+ sensibFlux(ig)! subtract flux to the ground
             enddo
-            
+
             call planetwide_sumval(dEdiff(:)*cell_area(:)/totarea_planet,dEtot)
             dEdiffs(:)=capcal(:)*zdtsdif(:)-zflubid(:)-sensibFlux(:)
             call planetwide_sumval(dEdiffs(:)*cell_area(:)/totarea_planet,dEtots)
             call planetwide_sumval(sensibFlux(:)*cell_area(:)/totarea_planet,AtmToSurf_TurbFlux)
-            
+
             if (is_master) then
-            
+
                if (UseTurbDiff) then
                          print*,'In TurbDiff sensible flux (atm=>surf) =',AtmToSurf_TurbFlux,' W m-2'
@@ -1234 +1234 @@
                end if
             endif ! end of 'is_master'
-            
+
          ! JL12 : note that the black body radiative flux emitted by the surface has been updated by the implicit scheme but not given back elsewhere.
          endif ! end of 'enertest'
@@ -1241 +1241 @@
          ! Test water conservation.
          if(watertest.and.water)then
-         
+
             call planetwide_sumval(massarea(:,:)*zdqdif(:,:,igcm_h2o_vap)*ptimestep/totarea_planet,dWtot_tmp)
             call planetwide_sumval(zdqsdif(:,igcm_h2o_vap)*cell_area(:)*ptimestep/totarea_planet,dWtots_tmp)
@@ -1253 +1253 @@
             do ig = 1, ngrid
                vdifcncons(ig)=vdifcncons(ig) + SUM(mass(ig,:)*zdqdif(ig,:,igcm_h2o_ice))
-            enddo            
+            enddo
             call planetwide_maxval(vdifcncons(:),nconsMAX)
 
@@ -1281 +1281 @@
 !   IV. Convection :
 !-------------------
-      
+
 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~
 ! IV.a Thermal plume model :
 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~
-      
+
       IF (calltherm) THEN
-         
+
 ! AB: We need to evaporate ice before calling thermcell_main.
          IF (water) THEN
@@ -1295 +1295 @@
             zqtherm(:,:,:) = pq(:,:,:) + pdq(:,:,:) * ptimestep
          ENDIF
-         
+
          CALL thermcell_main(ngrid, nlayer, nq, ptimestep, firstcall,            &
                              pplay, pplev, pphi, zpopsk,                         &
@@ -1301 +1301 @@
                              zdutherm, zdvtherm, zdttherm, zdqtherm,             &
                              fm, entr, detr, zw2, fraca)
-         
+
          pdu(1:ngrid,1:nlayer) = pdu(1:ngrid,1:nlayer) + zdutherm(1:ngrid,1:nlayer)
          pdv(1:ngrid,1:nlayer) = pdv(1:ngrid,1:nlayer) + zdvtherm(1:ngrid,1:nlayer)
          pdt(1:ngrid,1:nlayer) = pdt(1:ngrid,1:nlayer) + zdttherm(1:ngrid,1:nlayer)
-         
+
          IF (tracer) THEN
             pdq(1:ngrid,1:nlayer,1:nq) = pdq(1:ngrid,1:nlayer,1:nq) + zdqtherm(1:ngrid,1:nlayer,1:nq)
          ENDIF
-         
+
       ENDIF ! end of 'calltherm'
-      
+
 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 ! IV.b Dry convective adjustment :
@@ -1336 +1336 @@
          zdtadj(1:ngrid,1:nlayer) = zdhadj(1:ngrid,1:nlayer)*zpopsk(1:ngrid,1:nlayer) ! for diagnostic only
 
-         if(tracer) then 
-            pdq(1:ngrid,1:nlayer,1:nq) = pdq(1:ngrid,1:nlayer,1:nq) + zdqadj(1:ngrid,1:nlayer,1:nq) 
+         if(tracer) then
+            pdq(1:ngrid,1:nlayer,1:nq) = pdq(1:ngrid,1:nlayer,1:nq) + zdqadj(1:ngrid,1:nlayer,1:nq)
          end if
 
@@ -1356 +1356 @@
             do ig = 1, ngrid
                cadjncons(ig)=cadjncons(ig) + SUM(mass(ig,:)*zdqadj(ig,:,igcm_h2o_ice))
-            enddo            
+            enddo
             call planetwide_maxval(cadjncons(:),nconsMAX)
 
@@ -1363 +1363 @@
                print*,'In convadj MAX non-cons factor          =',nconsMAX,' kg m-2 s-1'
             endif
-            
+
          endif ! end of 'watertest'
-         
+
       endif ! end of 'calladj'
 !----------------------------------------------
@@ -1394 +1394 @@
 
 
-      
+
 !-----------------------------------------------
 !   V. Carbon dioxide condensation-sublimation :
@@ -1400 +1400 @@
 
       if (co2cond) then
-      
+
          if (.not.tracer) then
             print*,'We need a CO2 ice tracer to condense CO2'
@@ -1435 +1435 @@
 
 !---------------------------------------------
-!   VI. Specific parameterizations for tracers 
+!   VI. Specific parameterizations for tracers
 !---------------------------------------------
 
       if (tracer) then
-      
+
   ! ---------------------
   !   VI.1. Water and ice
   ! ---------------------
          if (water) then
-            
+
             ! Water ice condensation in the atmosphere
             if(watercond.and.(RLVTT.gt.1.e-8))then
-               
+
                if ((.not.calltherm).and.moistadjustment) then
                   dqmoist(1:ngrid,1:nlayer,1:nq)=0.
                   dtmoist(1:ngrid,1:nlayer)=0.
-                  
+
                   ! Moist Convective Adjustment.
                   ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -1462 +1462 @@
 !!!!
                   call moistadj(ngrid,nlayer,nq,pt,pq,pdq,pplev,pplay,dtmoist,dqmoist,ptimestep,rneb_man)
-               
+
                   pdq(1:ngrid,1:nlayer,igcm_h2o_vap) = pdq(1:ngrid,1:nlayer,igcm_h2o_vap)     &
                                                      + dqmoist(1:ngrid,1:nlayer,igcm_h2o_vap)
@@ -1468 +1468 @@
                                                      + dqmoist(1:ngrid,1:nlayer,igcm_h2o_ice)
                   pdt(1:ngrid,1:nlayer) = pdt(1:ngrid,1:nlayer)+dtmoist(1:ngrid,1:nlayer)
-                   
+
                   ! Test energy conservation.
                   if(enertest)then
@@ -1475 +1475 @@
                      call planetwide_minval(dtmoist(:,:),dtmoist_min)
                      madjdEz(:,:)=cpp*mass(:,:)*dtmoist(:,:)
-                     
+
                      do ig=1,ngrid
                         madjdE(ig) = cpp*SUM(mass(:,:)*dtmoist(:,:))
                      enddo
-                     
+
                      if (is_master) then
                         print*,'In moistadj atmospheric energy change   =',dEtot,' W m-2'
@@ -1485 +1485 @@
                         print*,'In moistadj MIN atmospheric energy change   =',dtmoist_min*ptimestep,'K/step'
                      endif
-                     
+
                      call planetwide_sumval(massarea(:,:)*dqmoist(:,:,igcm_h2o_vap)*ptimestep/totarea_planet+        &
                                             massarea(:,:)*dqmoist(:,:,igcm_h2o_ice)*ptimestep/totarea_planet,dWtot)
                      if (is_master) print*,'In moistadj atmospheric water change    =',dWtot,' kg m-2'
-                     
+
                   endif ! end of 'enertest'
                else
@@ -1497 +1497 @@
                   dqmoist(:,:,:)=0
                endif ! end of '(.not.calltherm).and.moistadjustment'
-               
+
                ! Large scale condensation/evaporation.
                ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -1519 +1519 @@
                   call planetwide_sumval(massarea(:,:)*dqvaplscale(:,:)*ptimestep/totarea_planet+        &
                                            SUM(massarea(:,:)*dqcldlscale(:,:))*ptimestep/totarea_planet,dWtot)
-                        
+
                   if (is_master) print*,'In largescale atmospheric water change  =',dWtot,' kg m-2'
                endif ! end of 'enertest'
@@ -1526 +1526 @@
                do l = 1, nlayer
                   do ig=1,ngrid
-                     cloudfrac(ig,l)=MAX(rneb_lsc(ig,l),rneb_man(ig,l)) 
+                     cloudfrac(ig,l)=MAX(rneb_lsc(ig,l),rneb_man(ig,l))
                   enddo
                enddo
 
             endif ! end of 'watercond'
-           
+
             ! Water ice / liquid precipitation.
             ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-            zdqrain(1:ngrid,1:nlayer,1:nq) = 0.0  !JL18 need to do that everytimestep if mass redis is on. 
+            zdqrain(1:ngrid,1:nlayer,1:nq) = 0.0  !JL18 need to do that everytimestep if mass redis is on.
 
             if(waterrain)then
@@ -1549 +1549 @@
                                                   + zdqrain(1:ngrid,1:nlayer,igcm_h2o_ice)
                pdt(1:ngrid,1:nlayer) = pdt(1:ngrid,1:nlayer)+zdtrain(1:ngrid,1:nlayer)
-               
+
                dqsurf(1:ngrid,igcm_h2o_vap) = dqsurf(1:ngrid,igcm_h2o_vap)+zdqsrain(1:ngrid)
-               dqsurf(1:ngrid,igcm_h2o_ice) = dqsurf(1:ngrid,igcm_h2o_ice)+zdqssnow(1:ngrid) 
-               
+               dqsurf(1:ngrid,igcm_h2o_ice) = dqsurf(1:ngrid,igcm_h2o_ice)+zdqssnow(1:ngrid)
+
 !!               call writediagfi(ngrid,"rain_post_dqsurf"," "," ",2,dqsurf(1:ngrid,igcm_h2o_vap))
 !!               call writediagfi(ngrid,"rain_post_qsurf"," "," ",2,qsurf(1:ngrid,igcm_h2o_vap))
-                
+
                ! Test energy conservation.
                if(enertest)then
-               
+
                   call planetwide_sumval(cpp*massarea(:,:)*zdtrain(:,:)/totarea_planet,dEtot)
                   if (is_master) print*,'In rain atmospheric T energy change       =',dEtot,' W m-2'
@@ -1568 +1568 @@
                   dVtot = dVtot + dVtot_tmp
                   dEtot = dItot + dVtot
-                  
+
                   if (is_master) then
                      print*,'In rain dItot =',dItot,' W m-2'
@@ -1579 +1579 @@
                           massarea(:,:)*zdqrain(:,:,igcm_h2o_ice)*ptimestep/totarea_planet,dWtot)
                   call planetwide_sumval((zdqsrain(:)+zdqssnow(:))*cell_area(:)*ptimestep/totarea_planet,dWtots)
-                  
+
                   if (is_master) then
                           print*,'In rain atmospheric water change        =',dWtot,' kg m-2'
@@ -1585 +1585 @@
                           print*,'In rain non-cons factor                 =',dWtot+dWtots,' kg m-2'
                   endif
-                  
+
                endif ! end of 'enertest'
 
             end if ! enf of 'waterrain'
-            
+
          end if ! end of 'water'
 
@@ -1648 +1648 @@
 
          !Generic Condensation
-         if (generic_condensation) then 
+         if (generic_condensation) then
             call condensation_generic(ngrid,nlayer,nq,ptimestep,pplev,pplay,   &
                                           pt,pq,pdt,pdq,dt_generic_condensation, &
@@ -1666 +1666 @@
             end if
 
-            if (.not. water) then 
+            if (.not. water) then
                ! Compute GCS (Generic Condensable Specie) cloud fraction. For now we can not have both water cloud fraction and GCS cloud fraction
                ! Water is the priority
@@ -1688 +1688 @@
                      enddo
                   endif
-            end do ! do iq=1,nq loop on tracers 
+            end do ! do iq=1,nq loop on tracers
             endif ! .not. water
 
@@ -1712 +1712 @@
 
             pdt(1:ngrid,1:nlayer) = pdt(1:ngrid,1:nlayer)+zdtrain_generic(1:ngrid,1:nlayer)
-            
+
             dqsurf(1:ngrid,1:nq) = dqsurf(1:ngrid,1:nq)+zdqsrain_generic(1:ngrid,1:nq)
 
@@ -1719 +1719 @@
             ! Test energy conservation.
             if(enertest)then
-               
+
                call planetwide_sumval(cpp*massarea(:,:)*zdtrain_generic(:,:)/totarea_planet,dEtot)
                if (is_master) print*,'In rain_generic atmospheric T energy change       =',dEtot,' W m-2'
@@ -1730 +1730 @@
 
                   call generic_tracer_index(nq,iq,igcm_generic_vap,igcm_generic_ice,call_ice_vap_generic)
-                  
+
                   if (call_ice_vap_generic) then ! to call only one time the ice/vap pair of a tracer
 
@@ -1740 +1740 @@
                      dVtot = dVtot + dVtot_tmp
                      dEtot = dItot + dVtot
-                     
+
                      if (is_master) then
                         print*,'In rain_generic dItot =',dItot,' W m-2'
@@ -1750 +1750 @@
                            massarea(:,:)*dq_rain_generic_cld(:,:,igcm_generic_ice)*ptimestep/totarea_planet,dWtot)
                      call planetwide_sumval((zdqsrain_generic(:,igcm_generic_ice)+zdqssnow_generic(:,igcm_generic_ice))*cell_area(:)*ptimestep/totarea_planet,dWtots)
-                     
+
                      if (is_master) then
                            print*,'In rain_generic atmospheric generic tracer change        =',dWtot,' kg m-2'
@@ -1759 +1759 @@
                   endif
 
-               end do ! do iq=1,nq loop on tracers 
-               
+               end do ! do iq=1,nq loop on tracers
+
             endif ! end of 'enertest'
-         
+
          endif !generic_rain
 
@@ -1772 +1772 @@
 
             if(watertest)then
-           
+
                iq=igcm_h2o_ice
                call planetwide_sumval(massarea(:,:)*pq(:,:,iq)*ptimestep/totarea_planet,dWtot)
@@ -1781 +1781 @@
                endif
             endif
-           
+
             call callsedim(ngrid,nlayer,ptimestep,           &
                           pplev,zzlev,pt,pdt,pq,pdq,        &
@@ -1834 +1834 @@
                zdmassmr_col(ig)=SUM(zdmassmr(ig,1:nlayer))
             enddo
-            
+
             call writediagfi(ngrid,"mass_evap","mass gain"," ",3,zdmassmr)
             call writediagfi(ngrid,"mass_evap_col","mass gain col"," ",2,zdmassmr_col)
@@ -1843 +1843 @@
                                      pu,pv,pdt,zdtsurf,pdq,pdu,pdv,zdmassmr,        &
                                      zdtmr,zdtsurfmr,zdpsrfmr,zdumr,zdvmr,zdqmr,zdqsurfmr)
-         
+
             pdq(1:ngrid,1:nlayer,1:nq) = pdq(1:ngrid,1:nlayer,1:nq) + zdqmr(1:ngrid,1:nlayer,1:nq)
             dqsurf(1:ngrid,1:nq)       = dqsurf(1:ngrid,1:nq) + zdqsurfmr(1:ngrid,1:nq)
@@ -1851 +1851 @@
             pdpsrf(1:ngrid)            = pdpsrf(1:ngrid) + zdpsrfmr(1:ngrid)
             zdtsurf(1:ngrid)           = zdtsurf(1:ngrid) + zdtsurfmr(1:ngrid)
-            
+
          endif
 
@@ -1881 +1881 @@
 !!                                     flux_g, dt_hdiff, &
 !!                                     dt_ekman)
-            
+
             call ocean_slab_noice(icount, ptimestep, knon, knindex, &
                       zdqssnow, tsea_ice, &
@@ -1896 +1896 @@
             call ocean_slab_get_vars(ngrid, tslab, sea_ice, flux_g, &
                       dt_hdiff, dt_ekman)
-  
+
             do ig=1,ngrid
                if (nint(rnat(ig)).eq.1)then
@@ -1919 +1919 @@
 
          if(hydrology)then
-            
+
             call hydrol(ngrid,nq,ptimestep,rnat,tsurf,qsurf,dqsurf,dqs_hyd, &
                         capcal,albedo,albedo_bareground,                    &
@@ -1925 +1925 @@
                         mu0,zdtsurf,zdtsurf_hyd,hice,pctsrf_sic,            &
                         sea_ice)
-                
+
 !!            call writediagfi(ngrid,"hydrol_post0_qsurf"," "," ",2,qsurf(1:ngrid,igcm_h2o_vap))
 
@@ -1931 +1931 @@
             dqsurf(1:ngrid,1:nq) = dqsurf(1:ngrid,1:nq) + dqs_hyd(1:ngrid,1:nq)
 
-            
+
             qsurf(1:ngrid,1:nq) = qsurf(1:ngrid,1:nq) + ptimestep*dqsurf(1:ngrid,1:nq)
 
@@ -1960 +1960 @@
          end if ! of if (hydrology)
 
-         ! Add qsurf to qsurf_hist, which is what we save in diagfi.nc. At the same time, we set the water 
+         ! Add qsurf to qsurf_hist, which is what we save in diagfi.nc. At the same time, we set the water
          ! content of ocean gridpoints back to zero, in order to avoid rounding errors in vdifc, rain.
          qsurf_hist(:,:) = qsurf(:,:)
@@ -1968 +1968 @@
 
 !------------------------------------------------
-!   VII. Surface and sub-surface soil temperature 
+!   VII. Surface and sub-surface soil temperature
 !------------------------------------------------
 
@@ -1981 +1981 @@
            tsurf(ig)=tsurf(ig)+ptimestep*zdtsurf(ig)
          enddo
-   
+
       else
-        tsurf(1:ngrid)=tsurf(1:ngrid)+ptimestep*zdtsurf(1:ngrid) 
+        tsurf(1:ngrid)=tsurf(1:ngrid)+ptimestep*zdtsurf(1:ngrid)
       endif ! end of 'ok_slab_ocean'
 
@@ -1990 +1990 @@
       if (callsoil) then
          call soil(ngrid,nsoilmx,.false.,lastcall,inertiedat,   &
-                   ptimestep,tsurf,tsoil,capcal,fluxgrd)     
+                   ptimestep,tsurf,tsoil,capcal,fluxgrd)
       endif
 
 
       if (ok_slab_ocean) then
-      
+
          do ig=1,ngrid
-         
+
             fluxgrdocean(ig)=fluxgrd(ig)
             if (nint(rnat(ig)).eq.0) then
@@ -2013 +2013 @@
                      capcal(ig)=capcalsno
                   endif
-               endif               
+               endif
             endif
-            
+
          enddo
-         
+
       endif !end of 'ok_slab_ocean'
 
@@ -2023 +2023 @@
       ! Test energy conservation
       if(enertest)then
-         call planetwide_sumval(cell_area(:)*capcal(:)*zdtsurf(:)/totarea_planet,dEtots)         
+         call planetwide_sumval(cell_area(:)*capcal(:)*zdtsurf(:)/totarea_planet,dEtots)
          if (is_master) print*,'Surface energy change                 =',dEtots,' W m-2'
       endif
@@ -2035 +2035 @@
 
 
- 
+
       ! Temperature, zonal and meridional winds.
       zt(1:ngrid,1:nlayer) = pt(1:ngrid,1:nlayer) + pdt(1:ngrid,1:nlayer)*ptimestep
       zu(1:ngrid,1:nlayer) = pu(1:ngrid,1:nlayer) + pdu(1:ngrid,1:nlayer)*ptimestep
       zv(1:ngrid,1:nlayer) = pv(1:ngrid,1:nlayer) + pdv(1:ngrid,1:nlayer)*ptimestep
-      
+
       ! Recast thermal plume vertical velocity array for outputs
       IF (calltherm) THEN
@@ -2113 +2113 @@
             endif
          end do
-                     
+
          if(ngrid.eq.1)then
             DYN=0.0
          endif
-         
+
          if (is_master) then
             print*,'  ISR            ASR            OLR            GND            DYN [W m^-2]'
@@ -2199 +2199 @@
 
          h2otot = icesrf + liqsrf + icecol + vapcol
-         
+
          if (is_master) then
             print*,' Total water amount [kg m^-2]: ',h2otot
@@ -2220 +2220 @@
       ! Calculate RH (Relative Humidity) for diagnostic.
       if(water)then
-      
+
          do l = 1, nlayer
             do ig=1,ngrid
@@ -2241 +2241 @@
 
             call generic_tracer_index(nq,iq,igcm_generic_vap,igcm_generic_ice,call_ice_vap_generic)
-            
+
             if (call_ice_vap_generic) then ! to call only one time the ice/vap pair of a tracer
-      
+
                do l = 1, nlayer
                   do ig=1,ngrid
@@ -2252 +2252 @@
 
             end if
-         
+
          end do ! iq=1,nq
 
@@ -2259 +2259 @@
 
       if (is_master) print*,'--> Ls =',zls*180./pi
-      
-      
+
+
 !----------------------------------------------------------------------
 !        Writing NetCDF file  "RESTARTFI" at the end of the run
@@ -2277 +2277 @@
 
 #ifndef MESOSCALE
-         
+
          if (ngrid.ne.1) then
             write(*,*)'PHYSIQ: for physdem ztime_fin =',ztime_fin
-            
+
             call physdem1("restartfi.nc",nsoilmx,ngrid,nlayer,nq, &
                           ptimestep,ztime_fin,                    &
@@ -2301 +2301 @@
 !  which can later be used to make the statistic files of the run:
 !  if flag "callstats" from callphys.def is .true.)
-         
+
 
          call wstats(ngrid,"ps","Surface pressure","Pa",2,ps)
@@ -2311 +2311 @@
                      "Thermal IR radiative flux to space","W.m-2",2,    &
                      fluxtop_lw)
-                     
+
 !            call wstats(ngrid,"fluxsurf_sw",                               &
 !                        "Solar radiative flux to surface","W.m-2",2,       &
-!                         fluxsurf_sw_tot)                     
+!                         fluxsurf_sw_tot)
 !            call wstats(ngrid,"fluxtop_sw",                                &
 !                        "Solar radiative flux to space","W.m-2",2,         &
@@ -2335 +2335 @@
             do iq=1,nq
                call wstats(ngrid,noms(iq),noms(iq),'kg/kg',3,zq(1,1,iq))
-               call wstats(ngrid,trim(noms(iq))//'_surf',trim(noms(iq))//'_surf',  & 
+               call wstats(ngrid,trim(noms(iq))//'_surf',trim(noms(iq))//'_surf',  &
                            'kg m^-2',2,qsurf(1,iq) )
-               call wstats(ngrid,trim(noms(iq))//'_col',trim(noms(iq))//'_col',    & 
+               call wstats(ngrid,trim(noms(iq))//'_col',trim(noms(iq))//'_col',    &
                           'kg m^-2',2,qcol(1,iq) )
-                          
-!              call wstats(ngrid,trim(noms(iq))//'_reff',                          & 
-!                          trim(noms(iq))//'_reff',                                   & 
+
+!              call wstats(ngrid,trim(noms(iq))//'_reff',                          &
+!                          trim(noms(iq))//'_reff',                                   &
 !                          'm',3,reffrad(1,1,iq))
 
             end do
-            
+
             if (water) then
                vmr=zq(1:ngrid,1:nlayer,igcm_h2o_vap)*mugaz/mmol(igcm_h2o_vap)
                call wstats(ngrid,"vmr_h2ovapor",                             &
-                           "H2O vapour volume mixing ratio","mol/mol",       & 
+                           "H2O vapour volume mixing ratio","mol/mol",       &
                            3,vmr)
             endif
 
-         endif ! end of 'tracer' 
+         endif ! end of 'tracer'
 
          if(watercond.and.CLFvarying)then
@@ -2381 +2381 @@
             call mkstats(ierr)
          endif
-         
+
 
 #ifndef MESOSCALE
-       
+
 !-----------------------------------------------------------------------------------------------------
 !           OUTPUT in netcdf file "DIAGFI.NC", containing any variable for diagnostic
@@ -2427 +2427 @@
          call writediagfi(ngrid,"latentFlux","latent heat flux","w.m^-2",2,zdqsdif(:,igcm_h2o_vap)*RLVTT)
       endif
-      
+
       ! Thermal plume model
       if (calltherm) then
@@ -2442 +2442 @@
          call WRITEDIAGFI(ngrid,"dvgwno","GW non-oro dV","m/s2", 3,d_v_hin)
       endif
-      
+
       ! Total energy balance diagnostics
       if(callrad.and.(.not.newtonian))then
-      
+
          !call writediagfi(ngrid,"ALB","Surface albedo"," ",2,albedo_equivalent)
          !call writediagfi(ngrid,"ALB_1st","First Band Surface albedo"," ",2,albedo(:,1))
@@ -2452 +2452 @@
          call writediagfi(ngrid,"OLR","outgoing longwave rad.","W m-2",2,fluxtop_lw)
          call writediagfi(ngrid,"shad","rings"," ", 2, fract)
-         
+
 !           call writediagfi(ngrid,"ASRcs","absorbed stellar rad (cs).","W m-2",2,fluxabs_sw1)
 !           call writediagfi(ngrid,"OLRcs","outgoing longwave rad (cs).","W m-2",2,fluxtop_lw1)
@@ -2465 +2465 @@
             call writediagfi(ngrid,"GND","heat flux from ground","W m-2",2,fluxgrd)
          endif
-         
+
          call writediagfi(ngrid,"DYN","dynamical heat input","W m-2",2,fluxdyn)
-         
+
       endif ! end of 'callrad'
-        
+
       if(enertest) then
-      
+
          if (calldifv) then
-         
+
             call writediagfi(ngrid,"q2","turbulent kinetic energy","J.kg^-1",3,q2)
             call writediagfi(ngrid,"sensibFlux","sensible heat flux","w.m^-2",2,sensibFlux)
-            
+
 !             call writediagfi(ngrid,"dEzdiff","turbulent diffusion heating (-sensible flux)","w.m^-2",3,dEzdiff)
 !             call writediagfi(ngrid,"dEdiff","integrated turbulent diffusion heating (-sensible flux)","w.m^-2",2,dEdiff)
 !             call writediagfi(ngrid,"dEdiffs","In TurbDiff (correc rad+latent heat) surf nrj change","w.m^-2",2,dEdiffs)
-             
-         endif
-          
+
+         endif
+
          if (corrk) then
             call writediagfi(ngrid,"dEzradsw","radiative heating","w.m^-2",3,dEzradsw)
             call writediagfi(ngrid,"dEzradlw","radiative heating","w.m^-2",3,dEzradlw)
          endif
-          
+
          if(watercond) then
 
-            call writediagfi(ngrid,"lscaledE","heat from largescale","W m-2",2,lscaledE) 
+            call writediagfi(ngrid,"lscaledE","heat from largescale","W m-2",2,lscaledE)
             if ((.not.calltherm).and.moistadjustment) then
                call writediagfi(ngrid,"madjdE","heat from moistadj","W m-2",2,madjdE)
             endif
             call writediagfi(ngrid,"qsatatm","atm qsat"," ",3,qsat)
-             
-!             call writediagfi(ngrid,"lscaledEz","heat from largescale","W m-2",3,lscaledEz) 
-!             call writediagfi(ngrid,"madjdEz","heat from moistadj","W m-2",3,madjdEz)             
+
+!             call writediagfi(ngrid,"lscaledEz","heat from largescale","W m-2",3,lscaledEz)
+!             call writediagfi(ngrid,"madjdEz","heat from moistadj","W m-2",3,madjdEz)
 !             call writediagfi(ngrid,"h2o_max_col","maximum H2O column amount","kg.m^-2",2,H2Omaxcol)
 
@@ -2506 +2506 @@
             call writediagfi(ngrid,"genericconddE","heat from generic condensation","W m-2",2,genericconddE)
             call writediagfi(ngrid,"dt_generic_condensation","heating from generic condensation","K s-1",3,dt_generic_condensation)
-         
-         endif
-          
+
+         endif
+
       endif ! end of 'enertest'
 
         ! Diagnostics of optical thickness
         ! Warning this is exp(-tau), I let you postproc with -log to have tau itself - JVO 19
-        if (diagdtau) then                
+        if (diagdtau) then
           do nw=1,L_NSPECTV
             write(str2,'(i2.2)') nw
@@ -2530 +2530 @@
         call writediagfi(ngrid,"dtrad","radiative heating","K s-1",3,dtrad)
         call writediagfi(ngrid,"zdtdyn","Dyn. heating","T s-1",3,zdtdyn)
-        
+
         ! For Debugging.
         !call writediagfi(ngrid,'rnat','Terrain type',' ',2,real(rnat))
         !call writediagfi(ngrid,'pphi','Geopotential',' ',3,pphi)
-        
+
 
       ! Output aerosols.
@@ -2548 +2548 @@
       ! Output tracers.
       if (tracer) then
-      
+
          do iq=1,nq
             call writediagfi(ngrid,noms(iq),noms(iq),'kg/kg',3,zq(1,1,iq))
-            call writediagfi(ngrid,trim(noms(iq))//'_surf',trim(noms(iq))//'_surf',  & 
+            call writediagfi(ngrid,trim(noms(iq))//'_surf',trim(noms(iq))//'_surf',  &
                              'kg m^-2',2,qsurf_hist(1,iq) )
-            call writediagfi(ngrid,trim(noms(iq))//'_col',trim(noms(iq))//'_col',    & 
+            call writediagfi(ngrid,trim(noms(iq))//'_col',trim(noms(iq))//'_col',    &
                             'kg m^-2',2,qcol(1,iq) )
-!          call writediagfi(ngrid,trim(noms(iq))//'_surf',trim(noms(iq))//'_surf',  & 
-!                          'kg m^-2',2,qsurf(1,iq) )                          
+!          call writediagfi(ngrid,trim(noms(iq))//'_surf',trim(noms(iq))//'_surf',  &
+!                          'kg m^-2',2,qsurf(1,iq) )
 
             if(watercond.or.CLFvarying)then
@@ -2591 +2591 @@
 
          enddo ! end of 'nq' loop
-         
+
        endif ! end of 'tracer'
 
 
       ! Output spectrum.
-      if(specOLR.and.corrk)then      
+      if(specOLR.and.corrk)then
          call writediagspecIR(ngrid,"OLR3D","OLR(lon,lat,band)","W/m^2/cm^-1",3,OLR_nu)
          call writediagspecVI(ngrid,"OSR3D","OSR(lon,lat,band)","W/m^2/cm^-1",3,OSR_nu)
@@ -2628 +2628 @@
       do iq=1,nq
          call generic_tracer_index(nq,iq,igcm_generic_vap,igcm_generic_ice,call_ice_vap_generic)
-                  
+
          if (call_ice_vap_generic) then ! to call only one time the ice/vap pair of a tracer
 
@@ -2645 +2645 @@
             comm_RH(1:ngrid,1:nlayer) = RH_generic(1:ngrid,1:nlayer,iq)
 
-         endif 
-      end do ! do iq=1,nq loop on tracers 
+         endif
+      end do ! do iq=1,nq loop on tracers
 
    else
@@ -2689 +2689 @@
 
 ! XIOS outputs
-#ifdef CPP_XIOS      
+#ifdef CPP_XIOS
       ! Send fields to XIOS: (NB these fields must also be defined as
       ! <field id="..." /> in context_lmdz_physics.xml to be correctly used)
       CALL send_xios_field("ls",zls)
-      
+
       CALL send_xios_field("ps",ps)
       CALL send_xios_field("area",cell_area)
@@ -2701 +2701 @@
       CALL send_xios_field("v",zv)
       CALL send_xios_field("omega",omega)
-      
+
       IF (calltherm) THEN
          CALL send_xios_field('w_plm',zw2_bis)
@@ -2709 +2709 @@
 !         CALL send_xios_field('fraca',fraca)
       ENDIF
-      
+
       IF (water) THEN
          CALL send_xios_field('h2o_vap',zq(:,:,igcm_h2o_vap))
@@ -2725 +2725 @@
       CALL send_xios_field("ASR",fluxabs_sw)
 
-      if (specOLR .and. corrk) then 
+      if (specOLR .and. corrk) then
          call send_xios_field("OLR3D",OLR_nu)
          call send_xios_field("IR_Bandwidth",DWNI)
@@ -2738 +2738 @@
       endif
 #endif
-      
+
       if (check_physics_outputs) then
          ! Check the validity of updated fields at the end of the physics step
@@ -2745 +2745 @@
 
       icount=icount+1
-      
+
       end subroutine physiq
-      
+
    end module physiq_mod