Ignore:
Timestamp:
Feb 21, 2024, 5:45:11 PM (10 months ago)
Author:
afalco
Message:

Pluto PCM:
1D functional
AF

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/LMDZ.GENERIC/libf/phystd/physiq_mod.F90

    r3204 r3232  
    11      module physiq_mod
    2      
     2
    33      implicit none
    4      
     4
    55      contains
    6      
     6
    77      subroutine physiq(ngrid,nlayer,nq,   &
    88                  firstcall,lastcall,      &
     
    1414
    1515!!
    16       use write_field_phy, only: Writefield_phy 
     16      use write_field_phy, only: Writefield_phy
    1717!!
    18       use ioipsl_getin_p_mod, only: getin_p 
     18      use ioipsl_getin_p_mod, only: getin_p
    1919      use radinc_h, only : L_NSPECTI,L_NSPECTV,naerkind, corrkdir, banddir
    2020      use watercommon_h, only : RLVTT, Psat_water,epsi,su_watercycle, RV, T_h2o_ice_liq
     
    8080      use turb_mod, only : q2,sensibFlux,turb_resolved
    8181      use mass_redistribution_mod, only: mass_redistribution
    82       use condensation_generic_mod, only: condensation_generic 
     82      use condensation_generic_mod, only: condensation_generic
    8383      use datafile_mod, only: datadir
    8484#ifndef MESOSCALE
     
    9898#endif
    9999
    100 #ifdef CPP_XIOS     
     100#ifdef CPP_XIOS
    101101      use xios_output_mod, only: initialize_xios_output, &
    102102                                 update_xios_timestep, &
     
    109109
    110110!==================================================================
    111 !     
     111!
    112112!     Purpose
    113113!     -------
     
    199199!           Frederic Hourdin        15/10/93
    200200!           Francois Forget        1994
    201 !           Christophe Hourdin        02/1997 
     201!           Christophe Hourdin        02/1997
    202202!           Subroutine completely rewritten by F. Forget (01/2000)
    203203!           Water ice clouds: Franck Montmessin (update 06/2003)
    204204!           Radiatively active tracers: J.-B. Madeleine (10/2008-06/2009)
    205205!           New correlated-k radiative scheme: R. Wordsworth (2009)
    206 !           Many specifically Martian subroutines removed: R. Wordsworth (2009)       
     206!           Many specifically Martian subroutines removed: R. Wordsworth (2009)
    207207!           Improved water cycle: R. Wordsworth / B. Charnay (2010)
    208208!           To F90: R. Wordsworth (2010)
     
    229229      integer,intent(in) :: nlayer            ! Number of atmospheric layers.
    230230      integer,intent(in) :: nq                ! Number of tracers.
    231      
     231
    232232      logical,intent(in) :: firstcall ! Signals first call to physics.
    233233      logical,intent(in) :: lastcall  ! Signals last call to physics.
    234      
     234
    235235      real,intent(in) :: pday                  ! Number of elapsed sols since reference Ls=0.
    236236      real,intent(in) :: ptime                 ! "Universal time", given as fraction of sol (e.g.: 0.5 for noon).
     
    265265! -----------------
    266266
    267 !     Aerosol (dust or ice) extinction optical depth  at reference wavelength 
     267!     Aerosol (dust or ice) extinction optical depth  at reference wavelength
    268268!     for the "naerkind" optically active aerosols:
    269269
     
    276276      integer l,ig,ierr,iq,nw,isoil,iesp, igcm_generic_vap, igcm_generic_ice
    277277      logical call_ice_vap_generic ! to call only one time the ice/vap pair of a tracer
    278      
     278
    279279      real zls                       ! Solar longitude (radians).
    280280      real zlss                      ! Sub solar point longitude (radians).
     
    284284
    285285! VARIABLES for the thermal plume model
    286      
     286
    287287      real f(ngrid)                       ! Mass flux norm
    288288      real fm(ngrid,nlayer+1)             ! Mass flux
     
    297297      real zw2(ngrid,nlayer+1)            ! Vertical speed
    298298      real zw2_bis(ngrid,nlayer)          ! Recasted zw2
    299      
     299
    300300! TENDENCIES due to various processes :
    301301
     
    306306      real zdtsdif(ngrid)     ! Turbdiff/vdifc routines.
    307307      real zdtsurf_hyd(ngrid) ! Hydrol routine.
    308            
    309       ! For Atmospheric Temperatures : (K/s)   
     308
     309      ! For Atmospheric Temperatures : (K/s)
    310310      real dtlscale(ngrid,nlayer)                             ! Largescale routine.
    311       real dt_generic_condensation(ngrid,nlayer)              ! condensation_generic routine. 
     311      real dt_generic_condensation(ngrid,nlayer)              ! condensation_generic routine.
    312312      real zdtc(ngrid,nlayer)                                 ! Condense_co2 routine.
    313313      real zdtdif(ngrid,nlayer)                               ! Turbdiff/vdifc routines.
     
    320320      real zdtsw1(ngrid,nlayer), zdtlw1(ngrid,nlayer)         ! Callcorrk routine.
    321321      real zdtchim(ngrid,nlayer)                              ! Calchim routine.
    322                              
     322
    323323      ! For Surface Tracers : (kg/m2/s)
    324324      real dqsurf(ngrid,nq)                 ! Cumulated tendencies.
     
    332332      real reevap_precip(ngrid)             ! re-evaporation flux of precipitation (integrated over the atmospheric column)
    333333      real reevap_precip_generic(ngrid,nq)
    334                  
     334
    335335      ! For Tracers : (kg/kg_of_air/s)
    336336      real zdqc(ngrid,nlayer,nq)      ! Condense_co2 routine.
     
    355355      REAL array_zero1(ngrid)
    356356      REAL array_zero2(ngrid,nlayer)
    357                        
     357
    358358      ! For Winds : (m/s/s)
    359359      real zdvadj(ngrid,nlayer), zduadj(ngrid,nlayer)       ! Convadj routine.
     
    363363      real zdhdif(ngrid,nlayer)                             ! Turbdiff/vdifc routines.
    364364      real zdhadj(ngrid,nlayer)                             ! Convadj routine.
    365      
     365
    366366      ! For Pressure and Mass :
    367367      real zdmassmr(ngrid,nlayer) ! Atmospheric Mass tendency for mass_redistribution (kg_of_air/m2/s).
     
    369369      real zdpsrfmr(ngrid)        ! Pressure tendency for mass_redistribution routine (Pa/s).
    370370
    371      
    372    
     371
     372
    373373! Local variables for LOCAL CALCULATIONS:
    374374! ---------------------------------------
     
    395395      real vmr(ngrid,nlayer)                        ! volume mixing ratio
    396396      real time_phys
    397      
     397
    398398      real ISR,ASR,OLR,GND,DYN,GSR,Ts1,Ts2,Ts3,TsS ! for Diagnostic.
    399      
     399
    400400      real qcol(ngrid,nq) ! Tracer Column Mass (kg/m2).
    401401
     
    410410      real,save :: dEtotSW, dEtotsSW, dEtotLW, dEtotsLW
    411411!$OMP THREADPRIVATE(dEtotSW, dEtotsSW, dEtotLW, dEtotsLW)
    412      
     412
    413413!JL12 conservation test for mean flow kinetic energy has been disabled temporarily
    414414
    415       real dtmoist_max,dtmoist_min     
     415      real dtmoist_max,dtmoist_min
    416416      real dItot, dItot_tmp, dVtot, dVtot_tmp
    417417
     
    437437
    438438      logical clearsky ! For double radiative transfer call. By BC
    439      
     439
    440440      ! For Clear Sky Case.
    441441      real fluxsurfabs_sw1(ngrid)  ! For SW/LW flux.
     
    466466      real :: flux_sens_lat(ngrid)
    467467      real :: qsurfint(ngrid,nq)
    468 #ifdef MESOSCALE 
     468#ifdef MESOSCALE
    469469      REAL :: lsf_dt(nlayer)
    470470      REAL :: lsf_dq(nlayer)
     
    474474      logical, save ::  check_physics_inputs=.false.
    475475      logical, save ::  check_physics_outputs=.false.
    476 !$OPM THREADPRIVATE(check_physics_inputs,check_physics_outputs)     
     476!$OPM THREADPRIVATE(check_physics_inputs,check_physics_outputs)
    477477
    478478      ! Misc
    479479      character*2 :: str2
    480       character(len=10) :: tmp1 
     480      character(len=10) :: tmp1
    481481      character(len=10) :: tmp2
    482482!==================================================================================================
     
    534534!        Read 'startfi.nc' file.
    535535!        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    536 #ifndef MESOSCALE 
     536#ifndef MESOSCALE
    537537         call phyetat0(startphy_file,                                 &
    538538                       ngrid,nlayer,"startfi.nc",0,0,nsoilmx,nq,      &
     
    541541                       rnat,pctsrf_sic,tslab, tsea_ice,sea_ice)
    542542
    543 #else   
     543#else
    544544
    545545         day_ini = pday
     
    555555             tsoil(1:ngrid,isoil)=tsurf(1:ngrid)
    556556           enddo
    557            if (is_master) write(*,*) "Physiq: initializing day_ini to pdat !"
     557           if (is_master) write(*,*) "Physiq: initializing day_ini to pday !"
    558558           day_ini=pday
    559559         endif
     
    569569         write (*,*) 'In physiq day_ini =', day_ini
    570570
    571 !        Initialize albedo calculation. 
     571!        Initialize albedo calculation.
    572572!        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    573573         albedo(:,:)=0.0
     
    575575         albedo_snow_SPECTV(:)=0.0
    576576         albedo_co2_ice_SPECTV(:)=0.0
    577          call surfini(ngrid,nq,qsurf,albedo,albedo_bareground,albedo_snow_SPECTV,albedo_co2_ice_SPECTV) 
    578 
    579 !        Initialize orbital calculation. 
     577         call surfini(ngrid,nq,qsurf,albedo,albedo_bareground,albedo_snow_SPECTV,albedo_co2_ice_SPECTV)
     578
     579!        Initialize orbital calculation.
    580580!        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    581581         call iniorbit(apoastr,periastr,year_day,peri_day,obliquit)
     
    596596
    597597!!           call ini_surf_heat_transp_mod()
    598            
     598
    599599           knindex(:) = 0
    600600           knon = 0
     
    620620!        ~~~~~~~~~~~~~~~~
    621621         if (callsoil) then
    622          
     622
    623623            call soil(ngrid,nsoilmx,firstcall,lastcall,inertiedat, &
    624624                      ptimestep,tsurf,tsoil,capcal,fluxgrd)
     
    639639
    640640         else ! else of 'callsoil'.
    641          
     641
    642642            print*,'WARNING! Thermal conduction in the soil turned off'
    643643            capcal(:)=1.e6
    644644            fluxgrd(:)=intheat
    645645            print*,'Flux from ground = ',intheat,' W m^-2'
    646            
     646
    647647         endif ! end of 'callsoil'.
    648          
     648
    649649         icount=1
    650650
     
    653653         ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    654654         if (.not.ok_slab_ocean) then
    655          
     655
    656656           rnat(:)=1.
    657657           if (.not.nosurf) then ! inertiedat only defined if there is a surface
     
    691691         endif
    692692
    693          if(meanOLR)then         
    694             call system('rm -f rad_bal.out') ! to record global radiative balance.           
    695             call system('rm -f tem_bal.out') ! to record global mean/max/min temperatures.           
     693         if(meanOLR)then
     694            call system('rm -f rad_bal.out') ! to record global radiative balance.
     695            call system('rm -f tem_bal.out') ! to record global mean/max/min temperatures.
    696696            call system('rm -f h2o_bal.out') ! to record global hydrological balance.
    697697         endif
    698698
    699699
    700          watertest=.false.       
     700         watertest=.false.
    701701         if(water)then ! Initialize variables for water cycle
    702            
     702
    703703            if(enertest)then
    704704               watertest = .true.
     
    711711         metallicity=0.0 ! default value --- is not used here but necessary to call function Psat_generic
    712712         call getin_p("metallicity",metallicity) ! --- is not used here but necessary to call function Psat_generic
    713          
     713
    714714!        Set some parameters for the thermal plume model
    715715!        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     
    731731
    732732!!         call WriteField_phy("post_physdem_qsurf",qsurf(1:ngrid,igcm_h2o_vap),1)
    733 #endif         
    734          if (corrk) then 
    735             ! We initialise the spectral grid here instead of 
    736             ! at firstcall of callcorrk so we can output XspecIR, XspecVI 
     733#endif
     734         if (corrk) then
     735            ! We initialise the spectral grid here instead of
     736            ! at firstcall of callcorrk so we can output XspecIR, XspecVI
    737737            ! when using Dynamico
    738738            print*, "physiq_mod: Correlated-k data base folder:",trim(datadir)
     
    757757                                     year_day,presnivs,pseudoalt,WNOI,WNOV)
    758758#endif
    759          
     759
    760760!!         call WriteField_phy("post_xios_qsurf",qsurf(1:ngrid,igcm_h2o_vap),1)
    761761
     
    764764
    765765!!      call WriteField_phy("post_firstcall_qsurf",qsurf(1:ngrid,igcm_h2o_vap),1)
    766 !!      call writediagfi(ngrid,"firstcall_post_qsurf"," "," ",2,qsurf(1:ngrid,igcm_h2o_vap)) 
     766!!      call writediagfi(ngrid,"firstcall_post_qsurf"," "," ",2,qsurf(1:ngrid,igcm_h2o_vap))
    767767
    768768      if (check_physics_inputs) then
     
    778778! ------------------------------------------------------
    779779
    780 #ifdef CPP_XIOS     
     780#ifdef CPP_XIOS
    781781      ! update XIOS time/calendar
    782782      call update_xios_timestep
    783 #endif     
     783#endif
    784784
    785785      ! Initialize various variables
    786786      ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    787      
     787
    788788      if ( .not.nearco2cond ) then
    789789         pdt(1:ngrid,1:nlayer) = 0.0
    790       endif     
     790      endif
    791791      zdtsurf(1:ngrid)      = 0.0
    792792      pdq(1:ngrid,1:nlayer,1:nq) = 0.0
     
    795795      pdv(1:ngrid,1:nlayer) = 0.0
    796796      pdpsrf(1:ngrid)       = 0.0
    797       zflubid(1:ngrid)      = 0.0     
     797      zflubid(1:ngrid)      = 0.0
    798798      flux_sens_lat(1:ngrid) = 0.0
    799799      taux(1:ngrid) = 0.0
     
    811811
    812812      call orbite(zls,dist_star,declin,right_ascen)
    813      
     813
    814814      if (tlocked) then
    815815              zlss=Mod(-(2.*pi*(zday/year_day)*nres - right_ascen),2.*pi)
     
    818818      else if(diurnal .eqv. .false.) then
    819819              zlss=9999.
    820       endif 
     820      endif
    821821
    822822
    823823      ! Compute variations of g with latitude (oblate case).
    824824      ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    825       if (oblate .eqv. .false.) then     
    826          glat(:) = g         
    827       else if (flatten .eq. 0.0 .or. J2 .eq. 0.0 .or. Rmean .eq. 0.0 .or. MassPlanet .eq. 0.0) then     
     825      if (oblate .eqv. .false.) then
     826         glat(:) = g
     827      else if (flatten .eq. 0.0 .or. J2 .eq. 0.0 .or. Rmean .eq. 0.0 .or. MassPlanet .eq. 0.0) then
    828828         print*,'I need values for flatten, J2, Rmean and MassPlanet to compute glat (else set oblate=.false.)'
    829829         call abort
     
    831831         gmplanet = MassPlanet*grav*1e24
    832832         do ig=1,ngrid
    833             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)))) 
     833            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))))
    834834         end do
    835835      endif
     
    838838      ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    839839      zzlay(1:ngrid,1:nlayer)=pphi(1:ngrid,1:nlayer)
    840       do l=1,nlayer         
     840      do l=1,nlayer
    841841         zzlay(1:ngrid,l)= zzlay(1:ngrid,l)/glat(1:ngrid)
    842842      enddo
     
    850850            zzlev(ig,l)=(z1*zzlay(ig,l-1)+z2*zzlay(ig,l))/(z1+z2)
    851851         enddo
    852       enddo     
     852      enddo
    853853
    854854      !Altitude of top interface (nlayer+1), using the thicknesss of the level below the top one. LT22
    855      
     855
    856856      zzlev(1:ngrid,nlayer+1) = 2*zzlev(1:ngrid,nlayer)-zzlev(1:ngrid,nlayer-1)
    857857
     
    859859      ! Note : Potential temperature calculation may not be the same in physiq and dynamic...
    860860      ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    861       do l=1,nlayer         
     861      do l=1,nlayer
    862862         do ig=1,ngrid
    863863            zpopsk(ig,l)=(pplay(ig,l)/pplev(ig,1))**rcp
     
    894894      if(photochem) then
    895895         call concentrations(ngrid,nlayer,nq,pplay,pt,pdt,pq,pdq,ptimestep)
    896       endif 
     896      endif
    897897#endif
    898898
     
    928928         endif
    929929
    930       endif 
     930      endif
    931931
    932932      if (callrad) then
    933933         if( mod(icount-1,iradia).eq.0.or.lastcall) then
    934            
    935             ! Eclipse incoming sunlight (e.g. Saturn ring shadowing).       
     934
     935            ! Eclipse incoming sunlight (e.g. Saturn ring shadowing).
    936936            if(rings_shadow) then
    937937                call call_rings(ngrid, ptime, pday, diurnal)
    938             endif   
     938            endif
    939939
    940940!!            call writediagfi(ngrid,"corrk_pre_dqsurf"," "," ",2,dqsurf(1:ngrid,igcm_h2o_vap))
    941 !!            call writediagfi(ngrid,"corrk_pre_qsurf"," "," ",2,qsurf(1:ngrid,igcm_h2o_vap)) 
     941!!            call writediagfi(ngrid,"corrk_pre_qsurf"," "," ",2,qsurf(1:ngrid,igcm_h2o_vap))
    942942
    943943
    944944            if (corrk) then
    945            
     945
    946946! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    947947! II.a Call correlated-k radiative transfer scheme
     
    952952               endif
    953953               if(water) then
    954                   muvar(1:ngrid,1:nlayer)=mugaz/(1.e0+(1.e0/epsi-1.e0)*pq(1:ngrid,1:nlayer,igcm_h2o_vap)) 
    955                   muvar(1:ngrid,nlayer+1)=mugaz/(1.e0+(1.e0/epsi-1.e0)*pq(1:ngrid,nlayer,igcm_h2o_vap)) 
     954                  muvar(1:ngrid,1:nlayer)=mugaz/(1.e0+(1.e0/epsi-1.e0)*pq(1:ngrid,1:nlayer,igcm_h2o_vap))
     955                  muvar(1:ngrid,nlayer+1)=mugaz/(1.e0+(1.e0/epsi-1.e0)*pq(1:ngrid,nlayer,igcm_h2o_vap))
    956956                  ! take into account water effect on mean molecular weight
    957957               elseif(generic_condensation) then
     
    959959
    960960                     call generic_tracer_index(nq,iq,igcm_generic_vap,igcm_generic_ice,call_ice_vap_generic)
    961                      
     961
    962962                     if (call_ice_vap_generic) then ! to call only one time the ice/vap pair of a tracer
    963963
    964964                        epsi_generic=constants_epsi_generic(iq)
    965965
    966                         muvar(1:ngrid,1:nlayer)=mugaz/(1.e0+(1.e0/epsi_generic-1.e0)*pq(1:ngrid,1:nlayer,igcm_generic_vap)) 
    967                         muvar(1:ngrid,nlayer+1)=mugaz/(1.e0+(1.e0/epsi_generic-1.e0)*pq(1:ngrid,nlayer,igcm_generic_vap)) 
     966                        muvar(1:ngrid,1:nlayer)=mugaz/(1.e0+(1.e0/epsi_generic-1.e0)*pq(1:ngrid,1:nlayer,igcm_generic_vap))
     967                        muvar(1:ngrid,nlayer+1)=mugaz/(1.e0+(1.e0/epsi_generic-1.e0)*pq(1:ngrid,nlayer,igcm_generic_vap))
    968968
    969969                     endif
    970                   end do ! do iq=1,nq loop on tracers 
     970                  end do ! do iq=1,nq loop on tracers
    971971               ! take into account generic condensable specie (GCS) effect on mean molecular weight
    972                
     972
    973973               else
    974974                  muvar(1:ngrid,1:nlayer+1)=mugaz
    975                endif 
    976      
     975               endif
     976
    977977
    978978               if(ok_slab_ocean) then
    979                   tsurf2(:)=tsurf(:) 
     979                  tsurf2(:)=tsurf(:)
    980980                  do ig=1,ngrid
    981981                     if (nint(rnat(ig))==0) then
     
    984984                  enddo
    985985               endif !(ok_slab_ocean)
    986                
     986
    987987               ! standard callcorrk
    988988               clearsky=.false.
     
    995995                              int_dtaui,int_dtauv,                                &
    996996                              tau_col,cloudfrac,totcloudfrac,                     &
    997                               clearsky,firstcall,lastcall)     
     997                              clearsky,firstcall,lastcall)
    998998
    999999                !GG (feb2021): Option to "artificially" decrease the raditive time scale in
    10001000                !the deep atmosphere  press > 0.1 bar. Suggested by MT.
    1001                 !! COEFF_RAD to be "tuned" to facilitate convergence of tendency 
    1002        
     1001                !! COEFF_RAD to be "tuned" to facilitate convergence of tendency
     1002
    10031003                !coeff_rad=0.   ! 0 values, it doesn't accelerate the convergence
    10041004                !coeff_rad=0.5
    1005                 !coeff_rad=1.                 
     1005                !coeff_rad=1.
    10061006                !do l=1, nlayer
    10071007                !  do ig=1,ngrid
     
    10131013                !enddo
    10141014
    1015                 ! Option to call scheme once more for clear regions 
     1015                ! Option to call scheme once more for clear regions
    10161016               if(CLFvarying)then
    10171017
     
    10271027                                 tau_col1,cloudfrac,totcloudfrac,                    &
    10281028                                 clearsky,firstcall,lastcall)
    1029                   clearsky = .false.  ! just in case.     
     1029                  clearsky = .false.  ! just in case.
    10301030
    10311031                  ! Sum the fluxes and heating rates from cloudy/clear cases
    10321032                  do ig=1,ngrid
    10331033                     tf=totcloudfrac(ig)
    1034                      ntf=1.-tf                   
     1034                     ntf=1.-tf
    10351035                     fluxsurf_lw(ig)       = ntf*fluxsurf_lw1(ig)       + tf*fluxsurf_lw(ig)
    10361036                     fluxsurf_sw(ig)       = ntf*fluxsurf_sw1(ig)       + tf*fluxsurf_sw(ig)
     
    10401040                     fluxabs_sw(ig)        = ntf*fluxabs_sw1(ig)        + tf*fluxabs_sw(ig)
    10411041                     tau_col(ig)           = ntf*tau_col1(ig)           + tf*tau_col(ig)
    1042                    
     1042
    10431043                     zdtlw(ig,1:nlayer) = ntf*zdtlw1(ig,1:nlayer) + tf*zdtlw(ig,1:nlayer)
    10441044                     zdtsw(ig,1:nlayer) = ntf*zdtsw1(ig,1:nlayer) + tf*zdtsw(ig,1:nlayer)
    10451045
    1046                      OSR_nu(ig,1:L_NSPECTV) = ntf*OSR_nu1(ig,1:L_NSPECTV) + tf*OSR_nu(ig,1:L_NSPECTV)       
    1047                      GSR_nu(ig,1:L_NSPECTV) = ntf*GSR_nu1(ig,1:L_NSPECTV) + tf*GSR_nu(ig,1:L_NSPECTV)                 
    1048                      OLR_nu(ig,1:L_NSPECTI) = ntf*OLR_nu1(ig,1:L_NSPECTI) + tf*OLR_nu(ig,1:L_NSPECTI)                       
     1046                     OSR_nu(ig,1:L_NSPECTV) = ntf*OSR_nu1(ig,1:L_NSPECTV) + tf*OSR_nu(ig,1:L_NSPECTV)
     1047                     GSR_nu(ig,1:L_NSPECTV) = ntf*GSR_nu1(ig,1:L_NSPECTV) + tf*GSR_nu(ig,1:L_NSPECTV)
     1048                     OLR_nu(ig,1:L_NSPECTI) = ntf*OLR_nu1(ig,1:L_NSPECTI) + tf*OLR_nu(ig,1:L_NSPECTI)
    10491049                     if (diagdtau) then
    1050                        int_dtauv(ig,:,1:L_NSPECTV) = ntf*int_dtauv1(ig,:,1:L_NSPECTV) + tf*int_dtauv(ig,:,1:L_NSPECTV)                       
    1051                        int_dtaui(ig,:,1:L_NSPECTI) = ntf*int_dtaui1(ig,:,1:L_NSPECTI) + tf*int_dtaui(ig,:,1:L_NSPECTI)                       
     1050                       int_dtauv(ig,:,1:L_NSPECTV) = ntf*int_dtauv1(ig,:,1:L_NSPECTV) + tf*int_dtauv(ig,:,1:L_NSPECTV)
     1051                       int_dtaui(ig,:,1:L_NSPECTI) = ntf*int_dtaui1(ig,:,1:L_NSPECTI) + tf*int_dtaui(ig,:,1:L_NSPECTI)
    10521052                     endif
    1053                   enddo                               
     1053                  enddo
    10541054
    10551055               endif ! end of CLFvarying.
     
    10721072               ! Net atmospheric radiative heating rate (K.s-1)
    10731073               dtrad(1:ngrid,1:nlayer)=zdtsw(1:ngrid,1:nlayer)+zdtlw(1:ngrid,1:nlayer)
    1074                
    1075                ! Late initialization of the Ice Spectral Albedo. We needed the visible bands to do that ! 
     1074
     1075               ! Late initialization of the Ice Spectral Albedo. We needed the visible bands to do that !
    10761076               if (firstcall .and. albedo_spectral_mode) then
    10771077                  call spectral_albedo_calc(albedo_snow_SPECTV)
     
    10791079
    10801080            elseif(newtonian)then
    1081            
     1081
    10821082! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    10831083! II.b Call Newtonian cooling scheme
     
    10891089
    10901090            else
    1091            
     1091
    10921092! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    1093 ! II.c Atmosphere has no radiative effect 
     1093! II.c Atmosphere has no radiative effect
    10941094! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    10951095               fluxtop_dn(1:ngrid)  = fract(1:ngrid)*mu0(1:ngrid)*Fat1AU/dist_star**2
     
    11001100               print*,'------------WARNING---WARNING------------' ! by MT2015.
    11011101               print*,'You are in corrk=false mode, '
    1102                print*,'and the surface albedo is taken equal to the first visible spectral value'               
    1103                
     1102               print*,'and the surface albedo is taken equal to the first visible spectral value'
     1103
    11041104               fluxsurfabs_sw(1:ngrid) = fluxtop_dn(1:ngrid)*(1.-albedo(1:ngrid,1))
    11051105               fluxrad_sky(1:ngrid)    = fluxsurfabs_sw(1:ngrid)
     
    11111111
    11121112         endif ! of if(mod(icount-1,iradia).eq.0)
    1113        
     1113
    11141114
    11151115         ! Transformation of the radiative tendencies
     
    11191119         fluxrad(1:ngrid)=fluxrad_sky(1:ngrid)-zplanck(1:ngrid)
    11201120         pdt(1:ngrid,1:nlayer)=pdt(1:ngrid,1:nlayer)+dtrad(1:ngrid,1:nlayer)
    1121          
     1121
    11221122         ! Test of energy conservation
    11231123         !----------------------------
     
    11521152
    11531153      if (calldifv) then
    1154      
     1154
    11551155         zflubid(1:ngrid)=fluxrad(1:ngrid)+fluxgrd(1:ngrid)
    11561156
    11571157         ! JL12 the following if test is temporarily there to allow us to compare the old vdifc with turbdiff.
    11581158         if (UseTurbDiff) then
    1159          
     1159
    11601160            call turbdiff(ngrid,nlayer,nq,rnat,                  &
    11611161                          ptimestep,capcal,                      &
     
    11681168
    11691169         else
    1170          
     1170
    11711171            zdh(1:ngrid,1:nlayer)=pdt(1:ngrid,1:nlayer)/zpopsk(1:ngrid,1:nlayer)
    1172  
     1172
    11731173            call vdifc(ngrid,nlayer,nq,rnat,zpopsk,           &
    11741174                       ptimestep,capcal,lwrite,               &
     
    11991199!!         call writediagfi(ngrid,"vdifc_post_zdqsdif"," "," ",2,zdqsdif(1:ngrid,igcm_h2o_vap))
    12001200
    1201          if (tracer) then 
     1201         if (tracer) then
    12021202           pdq(1:ngrid,1:nlayer,1:nq)=pdq(1:ngrid,1:nlayer,1:nq)+ zdqdif(1:ngrid,1:nlayer,1:nq)
    12031203           dqsurf(1:ngrid,1:nq)=dqsurf(1:ngrid,1:nq) + zdqsdif(1:ngrid,1:nq)
     
    12101210         !-------------------------
    12111211         if(enertest)then
    1212          
     1212
    12131213            dEzdiff(:,:)=cpp*mass(:,:)*zdtdif(:,:)
    12141214            do ig = 1, ngrid
     
    12161216               dEzdiff(ig,1)= dEzdiff(ig,1)+ sensibFlux(ig)! subtract flux to the ground
    12171217            enddo
    1218            
     1218
    12191219            call planetwide_sumval(dEdiff(:)*cell_area(:)/totarea_planet,dEtot)
    12201220            dEdiffs(:)=capcal(:)*zdtsdif(:)-zflubid(:)-sensibFlux(:)
    12211221            call planetwide_sumval(dEdiffs(:)*cell_area(:)/totarea_planet,dEtots)
    12221222            call planetwide_sumval(sensibFlux(:)*cell_area(:)/totarea_planet,AtmToSurf_TurbFlux)
    1223            
     1223
    12241224            if (is_master) then
    1225            
     1225
    12261226               if (UseTurbDiff) then
    12271227                         print*,'In TurbDiff sensible flux (atm=>surf) =',AtmToSurf_TurbFlux,' W m-2'
     
    12341234               end if
    12351235            endif ! end of 'is_master'
    1236            
     1236
    12371237         ! JL12 : note that the black body radiative flux emitted by the surface has been updated by the implicit scheme but not given back elsewhere.
    12381238         endif ! end of 'enertest'
     
    12411241         ! Test water conservation.
    12421242         if(watertest.and.water)then
    1243          
     1243
    12441244            call planetwide_sumval(massarea(:,:)*zdqdif(:,:,igcm_h2o_vap)*ptimestep/totarea_planet,dWtot_tmp)
    12451245            call planetwide_sumval(zdqsdif(:,igcm_h2o_vap)*cell_area(:)*ptimestep/totarea_planet,dWtots_tmp)
     
    12531253            do ig = 1, ngrid
    12541254               vdifcncons(ig)=vdifcncons(ig) + SUM(mass(ig,:)*zdqdif(ig,:,igcm_h2o_ice))
    1255             enddo           
     1255            enddo
    12561256            call planetwide_maxval(vdifcncons(:),nconsMAX)
    12571257
     
    12811281!   IV. Convection :
    12821282!-------------------
    1283      
     1283
    12841284! ~~~~~~~~~~~~~~~~~~~~~~~~~~
    12851285! IV.a Thermal plume model :
    12861286! ~~~~~~~~~~~~~~~~~~~~~~~~~~
    1287      
     1287
    12881288      IF (calltherm) THEN
    1289          
     1289
    12901290! AB: We need to evaporate ice before calling thermcell_main.
    12911291         IF (water) THEN
     
    12951295            zqtherm(:,:,:) = pq(:,:,:) + pdq(:,:,:) * ptimestep
    12961296         ENDIF
    1297          
     1297
    12981298         CALL thermcell_main(ngrid, nlayer, nq, ptimestep, firstcall,            &
    12991299                             pplay, pplev, pphi, zpopsk,                         &
     
    13011301                             zdutherm, zdvtherm, zdttherm, zdqtherm,             &
    13021302                             fm, entr, detr, zw2, fraca)
    1303          
     1303
    13041304         pdu(1:ngrid,1:nlayer) = pdu(1:ngrid,1:nlayer) + zdutherm(1:ngrid,1:nlayer)
    13051305         pdv(1:ngrid,1:nlayer) = pdv(1:ngrid,1:nlayer) + zdvtherm(1:ngrid,1:nlayer)
    13061306         pdt(1:ngrid,1:nlayer) = pdt(1:ngrid,1:nlayer) + zdttherm(1:ngrid,1:nlayer)
    1307          
     1307
    13081308         IF (tracer) THEN
    13091309            pdq(1:ngrid,1:nlayer,1:nq) = pdq(1:ngrid,1:nlayer,1:nq) + zdqtherm(1:ngrid,1:nlayer,1:nq)
    13101310         ENDIF
    1311          
     1311
    13121312      ENDIF ! end of 'calltherm'
    1313      
     1313
    13141314! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    13151315! IV.b Dry convective adjustment :
     
    13361336         zdtadj(1:ngrid,1:nlayer) = zdhadj(1:ngrid,1:nlayer)*zpopsk(1:ngrid,1:nlayer) ! for diagnostic only
    13371337
    1338          if(tracer) then 
    1339             pdq(1:ngrid,1:nlayer,1:nq) = pdq(1:ngrid,1:nlayer,1:nq) + zdqadj(1:ngrid,1:nlayer,1:nq) 
     1338         if(tracer) then
     1339            pdq(1:ngrid,1:nlayer,1:nq) = pdq(1:ngrid,1:nlayer,1:nq) + zdqadj(1:ngrid,1:nlayer,1:nq)
    13401340         end if
    13411341
     
    13561356            do ig = 1, ngrid
    13571357               cadjncons(ig)=cadjncons(ig) + SUM(mass(ig,:)*zdqadj(ig,:,igcm_h2o_ice))
    1358             enddo           
     1358            enddo
    13591359            call planetwide_maxval(cadjncons(:),nconsMAX)
    13601360
     
    13631363               print*,'In convadj MAX non-cons factor          =',nconsMAX,' kg m-2 s-1'
    13641364            endif
    1365            
     1365
    13661366         endif ! end of 'watertest'
    1367          
     1367
    13681368      endif ! end of 'calladj'
    13691369!----------------------------------------------
     
    13941394
    13951395
    1396      
     1396
    13971397!-----------------------------------------------
    13981398!   V. Carbon dioxide condensation-sublimation :
     
    14001400
    14011401      if (co2cond) then
    1402      
     1402
    14031403         if (.not.tracer) then
    14041404            print*,'We need a CO2 ice tracer to condense CO2'
     
    14351435
    14361436!---------------------------------------------
    1437 !   VI. Specific parameterizations for tracers 
     1437!   VI. Specific parameterizations for tracers
    14381438!---------------------------------------------
    14391439
    14401440      if (tracer) then
    1441      
     1441
    14421442  ! ---------------------
    14431443  !   VI.1. Water and ice
    14441444  ! ---------------------
    14451445         if (water) then
    1446            
     1446
    14471447            ! Water ice condensation in the atmosphere
    14481448            if(watercond.and.(RLVTT.gt.1.e-8))then
    1449                
     1449
    14501450               if ((.not.calltherm).and.moistadjustment) then
    14511451                  dqmoist(1:ngrid,1:nlayer,1:nq)=0.
    14521452                  dtmoist(1:ngrid,1:nlayer)=0.
    1453                  
     1453
    14541454                  ! Moist Convective Adjustment.
    14551455                  ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     
    14621462!!!!
    14631463                  call moistadj(ngrid,nlayer,nq,pt,pq,pdq,pplev,pplay,dtmoist,dqmoist,ptimestep,rneb_man)
    1464                
     1464
    14651465                  pdq(1:ngrid,1:nlayer,igcm_h2o_vap) = pdq(1:ngrid,1:nlayer,igcm_h2o_vap)     &
    14661466                                                     + dqmoist(1:ngrid,1:nlayer,igcm_h2o_vap)
     
    14681468                                                     + dqmoist(1:ngrid,1:nlayer,igcm_h2o_ice)
    14691469                  pdt(1:ngrid,1:nlayer) = pdt(1:ngrid,1:nlayer)+dtmoist(1:ngrid,1:nlayer)
    1470                    
     1470
    14711471                  ! Test energy conservation.
    14721472                  if(enertest)then
     
    14751475                     call planetwide_minval(dtmoist(:,:),dtmoist_min)
    14761476                     madjdEz(:,:)=cpp*mass(:,:)*dtmoist(:,:)
    1477                      
     1477
    14781478                     do ig=1,ngrid
    14791479                        madjdE(ig) = cpp*SUM(mass(:,:)*dtmoist(:,:))
    14801480                     enddo
    1481                      
     1481
    14821482                     if (is_master) then
    14831483                        print*,'In moistadj atmospheric energy change   =',dEtot,' W m-2'
     
    14851485                        print*,'In moistadj MIN atmospheric energy change   =',dtmoist_min*ptimestep,'K/step'
    14861486                     endif
    1487                      
     1487
    14881488                     call planetwide_sumval(massarea(:,:)*dqmoist(:,:,igcm_h2o_vap)*ptimestep/totarea_planet+        &
    14891489                                            massarea(:,:)*dqmoist(:,:,igcm_h2o_ice)*ptimestep/totarea_planet,dWtot)
    14901490                     if (is_master) print*,'In moistadj atmospheric water change    =',dWtot,' kg m-2'
    1491                      
     1491
    14921492                  endif ! end of 'enertest'
    14931493               else
     
    14971497                  dqmoist(:,:,:)=0
    14981498               endif ! end of '(.not.calltherm).and.moistadjustment'
    1499                
     1499
    15001500               ! Large scale condensation/evaporation.
    15011501               ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     
    15191519                  call planetwide_sumval(massarea(:,:)*dqvaplscale(:,:)*ptimestep/totarea_planet+        &
    15201520                                           SUM(massarea(:,:)*dqcldlscale(:,:))*ptimestep/totarea_planet,dWtot)
    1521                        
     1521
    15221522                  if (is_master) print*,'In largescale atmospheric water change  =',dWtot,' kg m-2'
    15231523               endif ! end of 'enertest'
     
    15261526               do l = 1, nlayer
    15271527                  do ig=1,ngrid
    1528                      cloudfrac(ig,l)=MAX(rneb_lsc(ig,l),rneb_man(ig,l)) 
     1528                     cloudfrac(ig,l)=MAX(rneb_lsc(ig,l),rneb_man(ig,l))
    15291529                  enddo
    15301530               enddo
    15311531
    15321532            endif ! end of 'watercond'
    1533            
     1533
    15341534            ! Water ice / liquid precipitation.
    15351535            ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    1536             zdqrain(1:ngrid,1:nlayer,1:nq) = 0.0  !JL18 need to do that everytimestep if mass redis is on. 
     1536            zdqrain(1:ngrid,1:nlayer,1:nq) = 0.0  !JL18 need to do that everytimestep if mass redis is on.
    15371537
    15381538            if(waterrain)then
     
    15491549                                                  + zdqrain(1:ngrid,1:nlayer,igcm_h2o_ice)
    15501550               pdt(1:ngrid,1:nlayer) = pdt(1:ngrid,1:nlayer)+zdtrain(1:ngrid,1:nlayer)
    1551                
     1551
    15521552               dqsurf(1:ngrid,igcm_h2o_vap) = dqsurf(1:ngrid,igcm_h2o_vap)+zdqsrain(1:ngrid)
    1553                dqsurf(1:ngrid,igcm_h2o_ice) = dqsurf(1:ngrid,igcm_h2o_ice)+zdqssnow(1:ngrid) 
    1554                
     1553               dqsurf(1:ngrid,igcm_h2o_ice) = dqsurf(1:ngrid,igcm_h2o_ice)+zdqssnow(1:ngrid)
     1554
    15551555!!               call writediagfi(ngrid,"rain_post_dqsurf"," "," ",2,dqsurf(1:ngrid,igcm_h2o_vap))
    15561556!!               call writediagfi(ngrid,"rain_post_qsurf"," "," ",2,qsurf(1:ngrid,igcm_h2o_vap))
    1557                
     1557
    15581558               ! Test energy conservation.
    15591559               if(enertest)then
    1560                
     1560
    15611561                  call planetwide_sumval(cpp*massarea(:,:)*zdtrain(:,:)/totarea_planet,dEtot)
    15621562                  if (is_master) print*,'In rain atmospheric T energy change       =',dEtot,' W m-2'
     
    15681568                  dVtot = dVtot + dVtot_tmp
    15691569                  dEtot = dItot + dVtot
    1570                  
     1570
    15711571                  if (is_master) then
    15721572                     print*,'In rain dItot =',dItot,' W m-2'
     
    15791579                          massarea(:,:)*zdqrain(:,:,igcm_h2o_ice)*ptimestep/totarea_planet,dWtot)
    15801580                  call planetwide_sumval((zdqsrain(:)+zdqssnow(:))*cell_area(:)*ptimestep/totarea_planet,dWtots)
    1581                  
     1581
    15821582                  if (is_master) then
    15831583                          print*,'In rain atmospheric water change        =',dWtot,' kg m-2'
     
    15851585                          print*,'In rain non-cons factor                 =',dWtot+dWtots,' kg m-2'
    15861586                  endif
    1587                  
     1587
    15881588               endif ! end of 'enertest'
    15891589
    15901590            end if ! enf of 'waterrain'
    1591            
     1591
    15921592         end if ! end of 'water'
    15931593
     
    16481648
    16491649         !Generic Condensation
    1650          if (generic_condensation) then 
     1650         if (generic_condensation) then
    16511651            call condensation_generic(ngrid,nlayer,nq,ptimestep,pplev,pplay,   &
    16521652                                          pt,pq,pdt,pdq,dt_generic_condensation, &
     
    16661666            end if
    16671667
    1668             if (.not. water) then 
     1668            if (.not. water) then
    16691669               ! Compute GCS (Generic Condensable Specie) cloud fraction. For now we can not have both water cloud fraction and GCS cloud fraction
    16701670               ! Water is the priority
     
    16881688                     enddo
    16891689                  endif
    1690             end do ! do iq=1,nq loop on tracers 
     1690            end do ! do iq=1,nq loop on tracers
    16911691            endif ! .not. water
    16921692
     
    17121712
    17131713            pdt(1:ngrid,1:nlayer) = pdt(1:ngrid,1:nlayer)+zdtrain_generic(1:ngrid,1:nlayer)
    1714            
     1714
    17151715            dqsurf(1:ngrid,1:nq) = dqsurf(1:ngrid,1:nq)+zdqsrain_generic(1:ngrid,1:nq)
    17161716
     
    17191719            ! Test energy conservation.
    17201720            if(enertest)then
    1721                
     1721
    17221722               call planetwide_sumval(cpp*massarea(:,:)*zdtrain_generic(:,:)/totarea_planet,dEtot)
    17231723               if (is_master) print*,'In rain_generic atmospheric T energy change       =',dEtot,' W m-2'
     
    17301730
    17311731                  call generic_tracer_index(nq,iq,igcm_generic_vap,igcm_generic_ice,call_ice_vap_generic)
    1732                  
     1732
    17331733                  if (call_ice_vap_generic) then ! to call only one time the ice/vap pair of a tracer
    17341734
     
    17401740                     dVtot = dVtot + dVtot_tmp
    17411741                     dEtot = dItot + dVtot
    1742                      
     1742
    17431743                     if (is_master) then
    17441744                        print*,'In rain_generic dItot =',dItot,' W m-2'
     
    17501750                           massarea(:,:)*dq_rain_generic_cld(:,:,igcm_generic_ice)*ptimestep/totarea_planet,dWtot)
    17511751                     call planetwide_sumval((zdqsrain_generic(:,igcm_generic_ice)+zdqssnow_generic(:,igcm_generic_ice))*cell_area(:)*ptimestep/totarea_planet,dWtots)
    1752                      
     1752
    17531753                     if (is_master) then
    17541754                           print*,'In rain_generic atmospheric generic tracer change        =',dWtot,' kg m-2'
     
    17591759                  endif
    17601760
    1761                end do ! do iq=1,nq loop on tracers 
    1762                
     1761               end do ! do iq=1,nq loop on tracers
     1762
    17631763            endif ! end of 'enertest'
    1764          
     1764
    17651765         endif !generic_rain
    17661766
     
    17721772
    17731773            if(watertest)then
    1774            
     1774
    17751775               iq=igcm_h2o_ice
    17761776               call planetwide_sumval(massarea(:,:)*pq(:,:,iq)*ptimestep/totarea_planet,dWtot)
     
    17811781               endif
    17821782            endif
    1783            
     1783
    17841784            call callsedim(ngrid,nlayer,ptimestep,           &
    17851785                          pplev,zzlev,pt,pdt,pq,pdq,        &
     
    18341834               zdmassmr_col(ig)=SUM(zdmassmr(ig,1:nlayer))
    18351835            enddo
    1836            
     1836
    18371837            call writediagfi(ngrid,"mass_evap","mass gain"," ",3,zdmassmr)
    18381838            call writediagfi(ngrid,"mass_evap_col","mass gain col"," ",2,zdmassmr_col)
     
    18431843                                     pu,pv,pdt,zdtsurf,pdq,pdu,pdv,zdmassmr,        &
    18441844                                     zdtmr,zdtsurfmr,zdpsrfmr,zdumr,zdvmr,zdqmr,zdqsurfmr)
    1845          
     1845
    18461846            pdq(1:ngrid,1:nlayer,1:nq) = pdq(1:ngrid,1:nlayer,1:nq) + zdqmr(1:ngrid,1:nlayer,1:nq)
    18471847            dqsurf(1:ngrid,1:nq)       = dqsurf(1:ngrid,1:nq) + zdqsurfmr(1:ngrid,1:nq)
     
    18511851            pdpsrf(1:ngrid)            = pdpsrf(1:ngrid) + zdpsrfmr(1:ngrid)
    18521852            zdtsurf(1:ngrid)           = zdtsurf(1:ngrid) + zdtsurfmr(1:ngrid)
    1853            
     1853
    18541854         endif
    18551855
     
    18811881!!                                     flux_g, dt_hdiff, &
    18821882!!                                     dt_ekman)
    1883            
     1883
    18841884            call ocean_slab_noice(icount, ptimestep, knon, knindex, &
    18851885                      zdqssnow, tsea_ice, &
     
    18961896            call ocean_slab_get_vars(ngrid, tslab, sea_ice, flux_g, &
    18971897                      dt_hdiff, dt_ekman)
    1898  
     1898
    18991899            do ig=1,ngrid
    19001900               if (nint(rnat(ig)).eq.1)then
     
    19191919
    19201920         if(hydrology)then
    1921            
     1921
    19221922            call hydrol(ngrid,nq,ptimestep,rnat,tsurf,qsurf,dqsurf,dqs_hyd, &
    19231923                        capcal,albedo,albedo_bareground,                    &
     
    19251925                        mu0,zdtsurf,zdtsurf_hyd,hice,pctsrf_sic,            &
    19261926                        sea_ice)
    1927                
     1927
    19281928!!            call writediagfi(ngrid,"hydrol_post0_qsurf"," "," ",2,qsurf(1:ngrid,igcm_h2o_vap))
    19291929
     
    19311931            dqsurf(1:ngrid,1:nq) = dqsurf(1:ngrid,1:nq) + dqs_hyd(1:ngrid,1:nq)
    19321932
    1933            
     1933
    19341934            qsurf(1:ngrid,1:nq) = qsurf(1:ngrid,1:nq) + ptimestep*dqsurf(1:ngrid,1:nq)
    19351935
     
    19601960         end if ! of if (hydrology)
    19611961
    1962          ! Add qsurf to qsurf_hist, which is what we save in diagfi.nc. At the same time, we set the water 
     1962         ! Add qsurf to qsurf_hist, which is what we save in diagfi.nc. At the same time, we set the water
    19631963         ! content of ocean gridpoints back to zero, in order to avoid rounding errors in vdifc, rain.
    19641964         qsurf_hist(:,:) = qsurf(:,:)
     
    19681968
    19691969!------------------------------------------------
    1970 !   VII. Surface and sub-surface soil temperature 
     1970!   VII. Surface and sub-surface soil temperature
    19711971!------------------------------------------------
    19721972
     
    19811981           tsurf(ig)=tsurf(ig)+ptimestep*zdtsurf(ig)
    19821982         enddo
    1983    
     1983
    19841984      else
    1985         tsurf(1:ngrid)=tsurf(1:ngrid)+ptimestep*zdtsurf(1:ngrid) 
     1985        tsurf(1:ngrid)=tsurf(1:ngrid)+ptimestep*zdtsurf(1:ngrid)
    19861986      endif ! end of 'ok_slab_ocean'
    19871987
     
    19901990      if (callsoil) then
    19911991         call soil(ngrid,nsoilmx,.false.,lastcall,inertiedat,   &
    1992                    ptimestep,tsurf,tsoil,capcal,fluxgrd)     
     1992                   ptimestep,tsurf,tsoil,capcal,fluxgrd)
    19931993      endif
    19941994
    19951995
    19961996      if (ok_slab_ocean) then
    1997      
     1997
    19981998         do ig=1,ngrid
    1999          
     1999
    20002000            fluxgrdocean(ig)=fluxgrd(ig)
    20012001            if (nint(rnat(ig)).eq.0) then
     
    20132013                     capcal(ig)=capcalsno
    20142014                  endif
    2015                endif               
     2015               endif
    20162016            endif
    2017            
     2017
    20182018         enddo
    2019          
     2019
    20202020      endif !end of 'ok_slab_ocean'
    20212021
     
    20232023      ! Test energy conservation
    20242024      if(enertest)then
    2025          call planetwide_sumval(cell_area(:)*capcal(:)*zdtsurf(:)/totarea_planet,dEtots)         
     2025         call planetwide_sumval(cell_area(:)*capcal(:)*zdtsurf(:)/totarea_planet,dEtots)
    20262026         if (is_master) print*,'Surface energy change                 =',dEtots,' W m-2'
    20272027      endif
     
    20352035
    20362036
    2037  
     2037
    20382038      ! Temperature, zonal and meridional winds.
    20392039      zt(1:ngrid,1:nlayer) = pt(1:ngrid,1:nlayer) + pdt(1:ngrid,1:nlayer)*ptimestep
    20402040      zu(1:ngrid,1:nlayer) = pu(1:ngrid,1:nlayer) + pdu(1:ngrid,1:nlayer)*ptimestep
    20412041      zv(1:ngrid,1:nlayer) = pv(1:ngrid,1:nlayer) + pdv(1:ngrid,1:nlayer)*ptimestep
    2042      
     2042
    20432043      ! Recast thermal plume vertical velocity array for outputs
    20442044      IF (calltherm) THEN
     
    21132113            endif
    21142114         end do
    2115                      
     2115
    21162116         if(ngrid.eq.1)then
    21172117            DYN=0.0
    21182118         endif
    2119          
     2119
    21202120         if (is_master) then
    21212121            print*,'  ISR            ASR            OLR            GND            DYN [W m^-2]'
     
    21992199
    22002200         h2otot = icesrf + liqsrf + icecol + vapcol
    2201          
     2201
    22022202         if (is_master) then
    22032203            print*,' Total water amount [kg m^-2]: ',h2otot
     
    22202220      ! Calculate RH (Relative Humidity) for diagnostic.
    22212221      if(water)then
    2222      
     2222
    22232223         do l = 1, nlayer
    22242224            do ig=1,ngrid
     
    22412241
    22422242            call generic_tracer_index(nq,iq,igcm_generic_vap,igcm_generic_ice,call_ice_vap_generic)
    2243            
     2243
    22442244            if (call_ice_vap_generic) then ! to call only one time the ice/vap pair of a tracer
    2245      
     2245
    22462246               do l = 1, nlayer
    22472247                  do ig=1,ngrid
     
    22522252
    22532253            end if
    2254          
     2254
    22552255         end do ! iq=1,nq
    22562256
     
    22592259
    22602260      if (is_master) print*,'--> Ls =',zls*180./pi
    2261      
    2262      
     2261
     2262
    22632263!----------------------------------------------------------------------
    22642264!        Writing NetCDF file  "RESTARTFI" at the end of the run
     
    22772277
    22782278#ifndef MESOSCALE
    2279          
     2279
    22802280         if (ngrid.ne.1) then
    22812281            write(*,*)'PHYSIQ: for physdem ztime_fin =',ztime_fin
    2282            
     2282
    22832283            call physdem1("restartfi.nc",nsoilmx,ngrid,nlayer,nq, &
    22842284                          ptimestep,ztime_fin,                    &
     
    23012301!  which can later be used to make the statistic files of the run:
    23022302!  if flag "callstats" from callphys.def is .true.)
    2303          
     2303
    23042304
    23052305         call wstats(ngrid,"ps","Surface pressure","Pa",2,ps)
     
    23112311                     "Thermal IR radiative flux to space","W.m-2",2,    &
    23122312                     fluxtop_lw)
    2313                      
     2313
    23142314!            call wstats(ngrid,"fluxsurf_sw",                               &
    23152315!                        "Solar radiative flux to surface","W.m-2",2,       &
    2316 !                         fluxsurf_sw_tot)                     
     2316!                         fluxsurf_sw_tot)
    23172317!            call wstats(ngrid,"fluxtop_sw",                                &
    23182318!                        "Solar radiative flux to space","W.m-2",2,         &
     
    23352335            do iq=1,nq
    23362336               call wstats(ngrid,noms(iq),noms(iq),'kg/kg',3,zq(1,1,iq))
    2337                call wstats(ngrid,trim(noms(iq))//'_surf',trim(noms(iq))//'_surf',  & 
     2337               call wstats(ngrid,trim(noms(iq))//'_surf',trim(noms(iq))//'_surf',  &
    23382338                           'kg m^-2',2,qsurf(1,iq) )
    2339                call wstats(ngrid,trim(noms(iq))//'_col',trim(noms(iq))//'_col',    & 
     2339               call wstats(ngrid,trim(noms(iq))//'_col',trim(noms(iq))//'_col',    &
    23402340                          'kg m^-2',2,qcol(1,iq) )
    2341                          
    2342 !              call wstats(ngrid,trim(noms(iq))//'_reff',                          & 
    2343 !                          trim(noms(iq))//'_reff',                                   & 
     2341
     2342!              call wstats(ngrid,trim(noms(iq))//'_reff',                          &
     2343!                          trim(noms(iq))//'_reff',                                   &
    23442344!                          'm',3,reffrad(1,1,iq))
    23452345
    23462346            end do
    2347            
     2347
    23482348            if (water) then
    23492349               vmr=zq(1:ngrid,1:nlayer,igcm_h2o_vap)*mugaz/mmol(igcm_h2o_vap)
    23502350               call wstats(ngrid,"vmr_h2ovapor",                             &
    2351                            "H2O vapour volume mixing ratio","mol/mol",       & 
     2351                           "H2O vapour volume mixing ratio","mol/mol",       &
    23522352                           3,vmr)
    23532353            endif
    23542354
    2355          endif ! end of 'tracer' 
     2355         endif ! end of 'tracer'
    23562356
    23572357         if(watercond.and.CLFvarying)then
     
    23812381            call mkstats(ierr)
    23822382         endif
    2383          
     2383
    23842384
    23852385#ifndef MESOSCALE
    2386        
     2386
    23872387!-----------------------------------------------------------------------------------------------------
    23882388!           OUTPUT in netcdf file "DIAGFI.NC", containing any variable for diagnostic
     
    24272427         call writediagfi(ngrid,"latentFlux","latent heat flux","w.m^-2",2,zdqsdif(:,igcm_h2o_vap)*RLVTT)
    24282428      endif
    2429      
     2429
    24302430      ! Thermal plume model
    24312431      if (calltherm) then
     
    24422442         call WRITEDIAGFI(ngrid,"dvgwno","GW non-oro dV","m/s2", 3,d_v_hin)
    24432443      endif
    2444      
     2444
    24452445      ! Total energy balance diagnostics
    24462446      if(callrad.and.(.not.newtonian))then
    2447      
     2447
    24482448         !call writediagfi(ngrid,"ALB","Surface albedo"," ",2,albedo_equivalent)
    24492449         !call writediagfi(ngrid,"ALB_1st","First Band Surface albedo"," ",2,albedo(:,1))
     
    24522452         call writediagfi(ngrid,"OLR","outgoing longwave rad.","W m-2",2,fluxtop_lw)
    24532453         call writediagfi(ngrid,"shad","rings"," ", 2, fract)
    2454          
     2454
    24552455!           call writediagfi(ngrid,"ASRcs","absorbed stellar rad (cs).","W m-2",2,fluxabs_sw1)
    24562456!           call writediagfi(ngrid,"OLRcs","outgoing longwave rad (cs).","W m-2",2,fluxtop_lw1)
     
    24652465            call writediagfi(ngrid,"GND","heat flux from ground","W m-2",2,fluxgrd)
    24662466         endif
    2467          
     2467
    24682468         call writediagfi(ngrid,"DYN","dynamical heat input","W m-2",2,fluxdyn)
    2469          
     2469
    24702470      endif ! end of 'callrad'
    2471        
     2471
    24722472      if(enertest) then
    2473      
     2473
    24742474         if (calldifv) then
    2475          
     2475
    24762476            call writediagfi(ngrid,"q2","turbulent kinetic energy","J.kg^-1",3,q2)
    24772477            call writediagfi(ngrid,"sensibFlux","sensible heat flux","w.m^-2",2,sensibFlux)
    2478            
     2478
    24792479!             call writediagfi(ngrid,"dEzdiff","turbulent diffusion heating (-sensible flux)","w.m^-2",3,dEzdiff)
    24802480!             call writediagfi(ngrid,"dEdiff","integrated turbulent diffusion heating (-sensible flux)","w.m^-2",2,dEdiff)
    24812481!             call writediagfi(ngrid,"dEdiffs","In TurbDiff (correc rad+latent heat) surf nrj change","w.m^-2",2,dEdiffs)
    2482              
    2483          endif
    2484          
     2482
     2483         endif
     2484
    24852485         if (corrk) then
    24862486            call writediagfi(ngrid,"dEzradsw","radiative heating","w.m^-2",3,dEzradsw)
    24872487            call writediagfi(ngrid,"dEzradlw","radiative heating","w.m^-2",3,dEzradlw)
    24882488         endif
    2489          
     2489
    24902490         if(watercond) then
    24912491
    2492             call writediagfi(ngrid,"lscaledE","heat from largescale","W m-2",2,lscaledE) 
     2492            call writediagfi(ngrid,"lscaledE","heat from largescale","W m-2",2,lscaledE)
    24932493            if ((.not.calltherm).and.moistadjustment) then
    24942494               call writediagfi(ngrid,"madjdE","heat from moistadj","W m-2",2,madjdE)
    24952495            endif
    24962496            call writediagfi(ngrid,"qsatatm","atm qsat"," ",3,qsat)
    2497              
    2498 !             call writediagfi(ngrid,"lscaledEz","heat from largescale","W m-2",3,lscaledEz) 
    2499 !             call writediagfi(ngrid,"madjdEz","heat from moistadj","W m-2",3,madjdEz)             
     2497
     2498!             call writediagfi(ngrid,"lscaledEz","heat from largescale","W m-2",3,lscaledEz)
     2499!             call writediagfi(ngrid,"madjdEz","heat from moistadj","W m-2",3,madjdEz)
    25002500!             call writediagfi(ngrid,"h2o_max_col","maximum H2O column amount","kg.m^-2",2,H2Omaxcol)
    25012501
     
    25062506            call writediagfi(ngrid,"genericconddE","heat from generic condensation","W m-2",2,genericconddE)
    25072507            call writediagfi(ngrid,"dt_generic_condensation","heating from generic condensation","K s-1",3,dt_generic_condensation)
    2508          
    2509          endif
    2510          
     2508
     2509         endif
     2510
    25112511      endif ! end of 'enertest'
    25122512
    25132513        ! Diagnostics of optical thickness
    25142514        ! Warning this is exp(-tau), I let you postproc with -log to have tau itself - JVO 19
    2515         if (diagdtau) then               
     2515        if (diagdtau) then
    25162516          do nw=1,L_NSPECTV
    25172517            write(str2,'(i2.2)') nw
     
    25302530        call writediagfi(ngrid,"dtrad","radiative heating","K s-1",3,dtrad)
    25312531        call writediagfi(ngrid,"zdtdyn","Dyn. heating","T s-1",3,zdtdyn)
    2532        
     2532
    25332533        ! For Debugging.
    25342534        !call writediagfi(ngrid,'rnat','Terrain type',' ',2,real(rnat))
    25352535        !call writediagfi(ngrid,'pphi','Geopotential',' ',3,pphi)
    2536        
     2536
    25372537
    25382538      ! Output aerosols.
     
    25482548      ! Output tracers.
    25492549      if (tracer) then
    2550      
     2550
    25512551         do iq=1,nq
    25522552            call writediagfi(ngrid,noms(iq),noms(iq),'kg/kg',3,zq(1,1,iq))
    2553             call writediagfi(ngrid,trim(noms(iq))//'_surf',trim(noms(iq))//'_surf',  & 
     2553            call writediagfi(ngrid,trim(noms(iq))//'_surf',trim(noms(iq))//'_surf',  &
    25542554                             'kg m^-2',2,qsurf_hist(1,iq) )
    2555             call writediagfi(ngrid,trim(noms(iq))//'_col',trim(noms(iq))//'_col',    & 
     2555            call writediagfi(ngrid,trim(noms(iq))//'_col',trim(noms(iq))//'_col',    &
    25562556                            'kg m^-2',2,qcol(1,iq) )
    2557 !          call writediagfi(ngrid,trim(noms(iq))//'_surf',trim(noms(iq))//'_surf',  & 
    2558 !                          'kg m^-2',2,qsurf(1,iq) )                         
     2557!          call writediagfi(ngrid,trim(noms(iq))//'_surf',trim(noms(iq))//'_surf',  &
     2558!                          'kg m^-2',2,qsurf(1,iq) )
    25592559
    25602560            if(watercond.or.CLFvarying)then
     
    25912591
    25922592         enddo ! end of 'nq' loop
    2593          
     2593
    25942594       endif ! end of 'tracer'
    25952595
    25962596
    25972597      ! Output spectrum.
    2598       if(specOLR.and.corrk)then     
     2598      if(specOLR.and.corrk)then
    25992599         call writediagspecIR(ngrid,"OLR3D","OLR(lon,lat,band)","W/m^2/cm^-1",3,OLR_nu)
    26002600         call writediagspecVI(ngrid,"OSR3D","OSR(lon,lat,band)","W/m^2/cm^-1",3,OSR_nu)
     
    26282628      do iq=1,nq
    26292629         call generic_tracer_index(nq,iq,igcm_generic_vap,igcm_generic_ice,call_ice_vap_generic)
    2630                  
     2630
    26312631         if (call_ice_vap_generic) then ! to call only one time the ice/vap pair of a tracer
    26322632
     
    26452645            comm_RH(1:ngrid,1:nlayer) = RH_generic(1:ngrid,1:nlayer,iq)
    26462646
    2647          endif 
    2648       end do ! do iq=1,nq loop on tracers 
     2647         endif
     2648      end do ! do iq=1,nq loop on tracers
    26492649
    26502650   else
     
    26892689
    26902690! XIOS outputs
    2691 #ifdef CPP_XIOS     
     2691#ifdef CPP_XIOS
    26922692      ! Send fields to XIOS: (NB these fields must also be defined as
    26932693      ! <field id="..." /> in context_lmdz_physics.xml to be correctly used)
    26942694      CALL send_xios_field("ls",zls)
    2695      
     2695
    26962696      CALL send_xios_field("ps",ps)
    26972697      CALL send_xios_field("area",cell_area)
     
    27012701      CALL send_xios_field("v",zv)
    27022702      CALL send_xios_field("omega",omega)
    2703      
     2703
    27042704      IF (calltherm) THEN
    27052705         CALL send_xios_field('w_plm',zw2_bis)
     
    27092709!         CALL send_xios_field('fraca',fraca)
    27102710      ENDIF
    2711      
     2711
    27122712      IF (water) THEN
    27132713         CALL send_xios_field('h2o_vap',zq(:,:,igcm_h2o_vap))
     
    27252725      CALL send_xios_field("ASR",fluxabs_sw)
    27262726
    2727       if (specOLR .and. corrk) then 
     2727      if (specOLR .and. corrk) then
    27282728         call send_xios_field("OLR3D",OLR_nu)
    27292729         call send_xios_field("IR_Bandwidth",DWNI)
     
    27382738      endif
    27392739#endif
    2740      
     2740
    27412741      if (check_physics_outputs) then
    27422742         ! Check the validity of updated fields at the end of the physics step
     
    27452745
    27462746      icount=icount+1
    2747      
     2747
    27482748      end subroutine physiq
    2749      
     2749
    27502750   end module physiq_mod
Note: See TracChangeset for help on using the changeset viewer.