Changeset 1647

Timestamp:

Jan 11, 2017, 3:33:51 PM (8 years ago)

Author:

jvatant

Message:

+ Major clean of the new LMDZ.TITAN from too-generic options and routines (water, co2, ocean, surface type ...)
+ From this revision LMDZ.TITAN begins to be really separated from LMDZ.GENERIC
+ Partial desactivation of aerosols, only the dummy case is still enabled to keep the code running ( new aerosol routines to come in followings commits )

JVO

Location:

trunk/LMDZ.TITAN/libf

Files:

• dynphy_lonlat/phytitan/ini_archive.F (modified) (3 diffs)
• dynphy_lonlat/phytitan/iniphysiq_mod.F90 (modified) (3 diffs)
• dynphy_lonlat/phytitan/lect_start_archive.F (modified) (20 diffs)
• dynphy_lonlat/phytitan/newstart.F (modified) (17 diffs)
• dynphy_lonlat/phytitan/start2archive.F (modified) (14 diffs)
• dynphy_lonlat/phytitan/write_archive.F (modified) (2 diffs)
• phytitan/aeropacity.F90 (modified) (6 diffs)
• phytitan/aerosol_mod.F90 (modified) (1 diff)
• phytitan/albedo_calcv.F90 (deleted)
• phytitan/albedo_snow_calc.F90 (deleted)
• phytitan/calc_cpp_mugaz.F90 (modified) (2 diffs)
• phytitan/calc_rayleigh.F90 (modified) (3 diffs)
• phytitan/callcorrk.F90 (modified) (8 diffs)
• phytitan/callkeys_mod.F90 (modified) (5 diffs)
• phytitan/callsedim.F (modified) (6 diffs)
• phytitan/comsoil_h.F90 (modified) (1 diff)
• phytitan/condense_co2.F90 (deleted)
• phytitan/convadj.F (modified) (12 diffs)
• phytitan/evap.F (deleted)
• phytitan/forceWCfn.F (deleted)
• phytitan/gases_h.F90 (modified) (1 diff)
• phytitan/hydrol.F90 (deleted)
• phytitan/iniaerosol.F (modified) (3 diffs)
• phytitan/inifis_mod.F90 (modified) (8 diffs)
• phytitan/initracer.F (modified) (7 diffs)
• phytitan/interpolateH2He.F90 (deleted)
• phytitan/interpolateH2Ocont_CKD.F90 (deleted)
• phytitan/interpolateH2Ocont_PPC.F90 (deleted)
• phytitan/iostart.F90 (modified) (7 diffs)
• phytitan/largescale.F90 (deleted)
• phytitan/mass_redistribution.F90 (modified) (6 diffs)
• phytitan/moistadj.F90 (deleted)
• phytitan/ocean_slab_mod.F90 (deleted)
• phytitan/optci.F90 (modified) (4 diffs)
• phytitan/optcv.F90 (modified) (3 diffs)
• phytitan/phyetat0.F90 (modified) (5 diffs)
• phytitan/phyredem.F90 (modified) (3 diffs)
• phytitan/physiq_mod.F90 (modified) (49 diffs)
• phytitan/radii_mod.F90 (modified) (6 diffs)
• phytitan/rain.F90 (deleted)
• phytitan/slab_ice_h.F90 (deleted)
• phytitan/su_gases.F90 (modified) (1 diff)
• phytitan/su_watercycle.F90 (deleted)
• phytitan/suaer_corrk.F90 (modified) (5 diffs)
• phytitan/sugas_corrk.F90 (modified) (4 diffs)
• phytitan/surf_heat_transp_mod.F90 (deleted)
• phytitan/surface_nature.F (deleted)
• phytitan/surfdat_h.F90 (modified) (1 diff)
• phytitan/surfini.F (modified) (4 diffs)
• phytitan/totalcloudfrac.F90 (deleted)
• phytitan/tracer_h.F90 (modified) (3 diffs)
• phytitan/turbdiff.F90 (modified) (12 diffs)
• phytitan/vdifc.F (modified) (12 diffs)
• phytitan/watercommon_h.F90 (deleted)
• phytitan/watersat.F90 (deleted)
• phytitan/watersat_grad.F90 (deleted)

trunk/LMDZ.TITAN/libf/dynphy_lonlat/phytitan/ini_archive.F

@@ -36 +36 @@
 
       USE comsoil_h
-      USE slab_ice_h, only: noceanmx
 !      use control_mod
       USE comvert_mod, ONLY: ap,bp,aps,bps,presnivs,pseudoalt
@@ -81 +80 @@
       INTEGER idim_tim
       INTEGER idim_nsoilmx ! "subsurface_layers" dimension ID #
-      INTEGER idim_noceanmx ! "ocean_layers" dimension ID #
       INTEGER nid,nvarid
       real sig_s(llm),s(llm)
@@ -164 +162 @@
       ierr = NF_DEF_DIM (nid, "altitude", llm, idim_llm)
       ierr = NF_DEF_DIM (nid,"subsurface_layers",nsoilmx,idim_nsoilmx)
-      ierr = NF_DEF_DIM (nid,"ocean_layers",noceanmx,idim_noceanmx)
 
       ierr = NF_DEF_DIM (nid,"index", length, idim_index)

trunk/LMDZ.TITAN/libf/dynphy_lonlat/phytitan/iniphysiq_mod.F90

@@ -11 +11 @@
 
 use control_mod, only: nday 
-use surf_heat_transp_mod, only: ini_surf_heat_transp
 use infotrac, only : nqtot ! number of advected tracers
 use planete_mod, only: ini_planete_mod
@@ -53 +52 @@
 integer,intent(in) :: iflag_phys ! type of physics to be called
 
-logical :: ok_slab_ocean
-
   ! the common part for all planetary physics
   !------------------------------------------
@@ -79 +76 @@
 call ini_planete_mod(nlayer,preff,ap,bp)
 
-! for slab ocean, copy over some arrays
-ok_slab_ocean=.false. ! default value
-call getin_p("ok_slab_ocean",ok_slab_ocean)
-if (ok_slab_ocean) then
-  call ini_surf_heat_transp(ip1jm,ip1jmp1,unsairez,fext,unsaire, &
-                            cu,cuvsurcv,cv,cvusurcu,aire,apoln,apols, &
-                            aireu,airev)
-endif
-
 call inifis(klon_omp,nlayer,nqtot,pdayref,punjours,nday,ptimestep, &
             latitude,longitude,cell_area,prad,pg,pr,pcpp)

trunk/LMDZ.TITAN/libf/dynphy_lonlat/phytitan/lect_start_archive.F

@@ -2 +2 @@
      &     date,tsurf,tsoil,emis,q2,
      &     t,ucov,vcov,ps,h,phisold_newgrid,
-     &     q,qsurf,surfith,nid,
-     &     rnat,pctsrf_sic,tslab,tsea_ice,sea_ice)
+     &     q,qsurf,surfith,nid)
 
 !      USE surfdat_h
       USE comsoil_h, ONLY: nsoilmx, layer, mlayer, volcapa, inertiedat
-      USE tracer_h, ONLY: igcm_co2_ice
       USE infotrac, ONLY: tname, nqtot
-      USE slab_ice_h, ONLY: noceanmx
 !      USE control_mod
 ! to use  'getin'
-      USE callkeys_mod, ONLY: ok_slab_ocean
       USE comvert_mod, ONLY: ap,bp,aps,bps,preff
       USE comconst_mod, ONLY: kappa,g,pi
@@ -101 +97 @@
       REAL tsurf(ngrid) ! surface temperature
       REAL tsoil(ngrid,nsoilmx) ! soil temperature
-      REAL co2ice(ngrid) ! CO2 ice layer
       REAL emis(ngrid)
       REAL q2(ngrid,llm+1),qsurf(ngrid,nqtot)
-      REAL tslab(ngrid,noceanmx)
-      REAL rnat(ngrid),pctsrf_sic(ngrid)
-      REAL tsea_ice(ngrid),sea_ice(ngrid)
 c     REAL phisfi(ngrid)
 
@@ -126 +118 @@
       real tsurfS(iip1,jjp1),tsoilS(iip1,jjp1,nsoilmx)
       real inertiedatS(iip1,jjp1,nsoilmx)
-      real co2iceS(iip1,jjp1)
       real emisS(iip1,jjp1)
       REAL q2S(iip1,jjp1,llm+1),qsurfS(iip1,jjp1,nqtot)
-      real tslabS(iip1,jjp1,noceanmx)
-      real pctsrf_sicS(iip1,jjp1),tsea_iceS(iip1,jjp1)
-      real rnatS(iip1,jjp1), sea_iceS(iip1,jjp1)
-
-      real ptotal, co2icetotal
+      
+      real ptotal
 
 c Var intermediaires : vent naturel, mais pas coord scalaire
@@ -156 +144 @@
       real, dimension(:,:,:), allocatable :: inertiedatoldnew
       real, dimension(:,:), allocatable :: psold,phisold
-      real, dimension(:,:), allocatable :: co2iceold
       real, dimension(:,:), allocatable :: tsurfold
       real, dimension(:,:), allocatable :: emisold
       real, dimension(:,:,:,:), allocatable :: qold
-      real, dimension(:,:,:), allocatable :: tslabold
-      real, dimension(:,:), allocatable :: rnatold,pctsrf_sicold
-      real, dimension(:,:), allocatable :: tsea_iceold,sea_iceold
-
 
       real tab_cntrl(100)
 
-      real ptotalold, co2icetotalold
+      real ptotalold
 
       logical :: olddepthdef=.false. ! flag
@@ -300 +283 @@
       allocate(phisold(imold+1,jmold+1))
       allocate(qold(imold+1,jmold+1,lmold,nqtot))
-      allocate(co2iceold(imold+1,jmold+1))
       allocate(tsurfold(imold+1,jmold+1))
       allocate(emisold(imold+1,jmold+1))
@@ -313 +295 @@
       allocate(mlayerold(nsoilold))
       allocate(qsurfold(imold+1,jmold+1,nqtot))
-      allocate(tslabold(imold+1,jmold+1,noceanmx))
-      allocate(rnatold(imold+1,jmold+1))
-      allocate(pctsrf_sicold(imold+1,jmold+1))
-      allocate(tsea_iceold(imold+1,jmold+1))
-      allocate(sea_iceold(imold+1,jmold+1))
 
       allocate(var (imold+1,jmold+1,llm))
@@ -332 +309 @@
 C-----------------------------------------------------------------------
 c 3.1. Lecture du tableau des parametres du run 
-c     (pour  la lecture ulterieure de "ptotalold" et "co2icetotalold")
+c     (pour  la lecture ulterieure de "ptotalold")
 c-----------------------------------------------------------------------
 c
@@ -556 +533 @@
 
 C-----------------------------------------------------------------------
-c   lecture de "ptotalold" et "co2icetotalold"
+c   lecture de "ptotalold"
 c-----------------------------------------------------------------------
       ptotalold = tab_cntrl(tab0+49)
-      co2icetotalold = tab_cntrl(tab0+50)
  
 c-----------------------------------------------------------------------
@@ -636 +612 @@
 
 c-----------------------------------------------------------------------
-c 5.1 Lecture des champs 2D (co2ice, emis,ps,tsurf,Tg[10], qsurf)
+c 5.1 Lecture des champs 2D (emis,ps,tsurf,Tg[10], qsurf)
 c-----------------------------------------------------------------------
  
@@ -642 +618 @@
       count=(/imold+1,jmold+1,1,0/)
        
-! CO2ice is now in qsurf(igcm_co2_ice) ...
-!      ierr = NF_INQ_VARID (nid, "co2ice", nvarid)
-!      IF (ierr .NE. NF_NOERR) THEN
-!         PRINT*, "lect_start_archive: Le champ <co2ice> est absent"
-!         CALL abort
-!      ENDIF
-!#ifdef NC_DOUBLE
-!      ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,count,co2iceold)
-!#else
-!      ierr = NF_GET_VARA_REAL(nid, nvarid,start,count,co2iceold)
-!#endif
-!      IF (ierr .NE. NF_NOERR) THEN
-!         PRINT*, "lect_start_archive: Lecture echouee pour <co2ice>"
-!         PRINT*, NF_STRERROR(ierr)
-!         CALL abort
-!      ENDIF
-c
+
       ierr = NF_INQ_VARID (nid, "emis", nvarid)
       IF (ierr .NE. NF_NOERR) THEN
@@ -718 +678 @@
          CALL abort
       ENDIF
-c
-cc Slab ocean
-      if(ok_slab_ocean) then
-      start=(/1,1,1,memo/)
-      count=(/imold+1,jmold+1,noceanmx,1/)
-
-       ierr=NF_INQ_VARID(nid,"tslab",nvarid)
-       IF (ierr.ne.NF_NOERR) then
-          PRINT*,"lect_start_archive: Cannot find <tslab>"
-       ENDIF
-#ifdef NC_DOUBLE
-      ierr=NF_GET_VARA_DOUBLE(nid,nvarid,start,count,tslabold)
-#else
-      ierr=NF_GET_VARA_REAL(nid,nvarid,start,count,tslabold)
-#endif
-      IF (ierr .NE. NF_NOERR) THEN
-         PRINT*, "lect_start_archive: Lecture echouee pour <tslab>"
-      ENDIF
-
-
-c
-      start=(/1,1,memo,0/)
-      count=(/imold+1,jmold+1,1,0/)
-
-      ierr = NF_INQ_VARID (nid, "rnat", nvarid)
-      IF (ierr .NE. NF_NOERR) THEN
-         PRINT*, "lect_start_archive: Le champ <rnat> est absent"
-      ENDIF
-#ifdef NC_DOUBLE
-      ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,count,rnatold)
-#else
-      ierr = NF_GET_VARA_REAL(nid, nvarid,start,count,rnatold)
-#endif
-      IF (ierr .NE. NF_NOERR) THEN
-         PRINT*, "lect_start_archive: Lecture echouee pour <rnat>"
-      ENDIF
-c
-      ierr = NF_INQ_VARID (nid, "pctsrf_sic", nvarid)
-      IF (ierr .NE. NF_NOERR) THEN
-         PRINT*, "lect_start_archive: Le champ <pctsrf_sic> est absent"
-      ENDIF
-#ifdef NC_DOUBLE
-      ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,count,pctsrf_sicold)
-#else
-      ierr = NF_GET_VARA_REAL(nid, nvarid,start,count,pctsrf_sicold)
-#endif
-      IF (ierr .NE. NF_NOERR) THEN
-         PRINT*, "lect_start_archive: Lecture echouee pour <pctsrf_sic>"
-      ENDIF
-c
-      ierr = NF_INQ_VARID (nid, "tsea_ice", nvarid)
-      IF (ierr .NE. NF_NOERR) THEN
-         PRINT*, "lect_start_archive: Le champ <tsea_ice> est absent"
-      ENDIF
-#ifdef NC_DOUBLE
-      ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,count,tsea_iceold)
-#else
-      ierr = NF_GET_VARA_REAL(nid, nvarid,start,count,tsea_iceold)
-#endif
-      IF (ierr .NE. NF_NOERR) THEN
-         PRINT*, "lect_start_archive: Lecture echouee pour <tsea_ice>"
-      ENDIF
-c
-      ierr = NF_INQ_VARID (nid, "sea_ice", nvarid)
-      IF (ierr .NE. NF_NOERR) THEN
-         PRINT*, "lect_start_archive: Le champ <sea_ice> est absent"
-      ENDIF
-#ifdef NC_DOUBLE
-      ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,count,sea_iceold)
-#else
-      ierr = NF_GET_VARA_REAL(nid, nvarid,start,count,sea_iceold)
-#endif
-      IF (ierr .NE. NF_NOERR) THEN
-         PRINT*, "lect_start_archive: Lecture echouee pour <sea_ice>"
-      ENDIF
-  
-      ENDIF! ok_slab_ocean
+
 c
       write(*,*)"lect_start_archive: rlonuold:"
@@ -813 +697 @@
 
       DO iq=1,nqtot
-          txt=trim(tname(iq))//"_surf"
-          if (txt.eq."h2o_vap") then
-            ! There is no surface tracer for h2o_vap;
-            ! "h2o_ice" should be loaded instead
-            txt="h2o_ice_surf"
-            write(*,*) 'lect_start_archive: loading surface tracer',
-     &                     ' h2o_ice instead of h2o_vap'
-          endif
-
         
         write(*,*) "lect_start_archive: loading tracer ",trim(txt)
@@ -1045 +920 @@
      &                   rlonuold,rlatvold,rlonu,rlatv)
 
-! CO2 ice is now in qsurf(igcm_co2_ice)
-!      call interp_horiz (co2iceold,co2ices,imold,jmold,iim,jjm,1,
-!     &                   rlonuold,rlatvold,rlonu,rlatv)
-
 c Temperature de surface
       call interp_horiz (tsurfold,tsurfs,imold,jmold,iim,jjm,1,
@@ -1088 +959 @@
          END DO
       END DO
-      co2icetotal = 0.
-      if (igcm_co2_ice.ne.0) then
-        ! recast surface CO2 ice on new grid
-        call interp_horiz(qsurfold(1,1,igcm_co2_ice),
-     &                  qsurfs(1,1,igcm_co2_ice),
-     &                  imold,jmold,iim,jjm,1,
-     &                  rlonuold,rlatvold,rlonu,rlatv)
-       DO j=1,jjp1
-         DO i=1,iim
-           !co2icetotal = co2icetotal + co2iceS(i,j)*aire(i,j)
-           co2icetotal=co2icetotal+qsurfS(i,j,igcm_co2_ice)*aire(i,j)
-         END DO
-       END DO
-      else
-        write(*,*) "Warning: No co2_ice surface tracer"
-      endif
+      
 
       write(*,*)
@@ -1110 +966 @@
       write (*,*) 'Ratio new atm./ old atm =', ptotal/ptotalold 
       write(*,*)
-      write(*,*)'Ancienne grille: masse de la glace CO2:',co2icetotalold
-      write(*,*)'Nouvelle grille: masse de la glace CO2:',co2icetotal
-      if (co2icetotalold.ne.0.) then
-      write(*,*)'Ratio new ice./old ice =',co2icetotal/co2icetotalold
-      endif
-      write(*,*)
-
 
       DO j=1,jjp1
@@ -1124 +973 @@
       END DO
 
-      if ( co2icetotalold.gt.0.) then 
-!         DO j=1,jjp1
-!            DO i=1,iip1
-!               co2iceS(i,j)=co2iceS(i,j) * co2icetotalold/co2icetotal
-!            END DO
-!         END DO
-        write(*,*) "Not enforcing conservation of surface CO2 ice"
-        write(*,*) " (should be OK when CO2 ice is a tracer)"
-      end if
 
 c-----------------------------------------------------------------------
@@ -1290 +1130 @@
 
 c-----------------------------------------------------------------------
-c 6.3   Slab Ocean :
-c-----------------------------------------------------------------------
-      call interp_horiz (tslabold,tslabs,imold,jmold,iim,jjm,noceanmx,
-     &                   rlonuold,rlatvold,rlonu,rlatv)
-      call gr_dyn_fi (noceanmx,iim+1,jjm+1,ngrid,tslabs,tslab)
-
-      call interp_horiz (rnatold,rnats,imold,jmold,iim,jjm,1,
-     &                   rlonuold,rlatvold,rlonu,rlatv)
-      call gr_dyn_fi (1,iim+1,jjm+1,ngrid,rnats,rnat)
-
-      call interp_horiz (pctsrf_sicold,pctsrf_sics,imold,jmold,iim,
-     &                   jjm,1,rlonuold,rlatvold,rlonu,rlatv)
-      call gr_dyn_fi (1,iim+1,jjm+1,ngrid,pctsrf_sics,pctsrf_sic)
-
-      call interp_horiz (tsea_iceold,tsea_ices,imold,jmold,iim,jjm,1,
-     &                   rlonuold,rlatvold,rlonu,rlatv)
-      call gr_dyn_fi (1,iim+1,jjm+1,ngrid,tsea_ices,tsea_ice)
-
-      call interp_horiz (sea_iceold,sea_ices,imold,jmold,iim,jjm,1,
-     &                   rlonuold,rlatvold,rlonu,rlatv)
-      call gr_dyn_fi (1,iim+1,jjm+1,ngrid,sea_ices,sea_ice)
-
-c-----------------------------------------------------------------------
-c 6.4 Variable 3d :
+c 6.3 Variable 3d :
 c-----------------------------------------------------------------------
       
@@ -1434 +1251 @@
          end do
       enddo
-      
-!      call gr_dyn_fi (1,iim+1,jjm+1,ngrid,co2ices,co2ice)
-! no need to transfer "co2ice" any more; it is in qsurf(igcm_co2_ice)
 
       endif !! if nqtot .ne. 0
@@ -1470 +1284 @@
       deallocate(phisold)
       deallocate(qold)
-      deallocate(co2iceold)
       deallocate(tsurfold)
       deallocate(emisold)
@@ -1482 +1295 @@
       deallocate(qsurfold)
       deallocate(var,varp1)
-      deallocate(tslabold)
-      deallocate(rnatold)
-      deallocate(pctsrf_sicold)
-      deallocate(tsea_iceold)
-      deallocate(sea_iceold)
 
 !      write(*,*)'lect_start_archive: END'

trunk/LMDZ.TITAN/libf/dynphy_lonlat/phytitan/newstart.F

@@ -19 +19 @@
      &                              is_master
       use infotrac, only: infotrac_init, nqtot, tname
-      USE tracer_h, ONLY: igcm_co2_ice, igcm_h2o_vap, igcm_h2o_ice
       USE comsoil_h, ONLY: nsoilmx, layer, mlayer, inertiedat
       USE surfdat_h, ONLY: phisfi, albedodat,
@@ -28 +27 @@
       use phyredem, only: physdem0, physdem1
       use iostart, only: open_startphy
-      use slab_ice_h, only:noceanmx
       use filtreg_mod, only: inifilr
       USE mod_const_mpi, ONLY: COMM_LMDZ
@@ -95 +93 @@
       REAL tsurf(ngridmx)        ! surface temperature
       REAL,ALLOCATABLE :: tsoil(:,:) ! soil temperature
-!      REAL co2ice(ngridmx)        ! CO2 ice layer
       REAL emis(ngridmx)        ! surface emissivity
       real emisread             ! added by RW
@@ -109 +106 @@
       real mugaz ! molar mass of the atmosphere
 
-      integer ierr 
-
-      REAL rnat(ngridmx)
-      REAL,ALLOCATABLE :: tslab(:,:) ! slab ocean temperature
-      REAL pctsrf_sic(ngridmx) ! sea ice cover
-      REAL tsea_ice(ngridmx) ! temperature sea_ice
-      REAL sea_ice(ngridmx) ! mass of sea ice (kg/m2)
+      integer ierr
 
 c Variables on the new grid along scalar points 
@@ -153 +144 @@
       logical :: flagtset=.false. ,  flagps0=.false.
       real val, val2, val3, val4 ! to store temporary variables
-      real :: iceith=2000 ! thermal inertia of subterranean ice
 
       INTEGER :: itau
@@ -165 +155 @@
 !     added by BC for equilibrium temperature startup
       real teque
-
-!     added by BC for cloud fraction setup
-      REAL hice(ngridmx),cloudfrac(ngridmx,llm)
-      REAL totalfrac(ngridmx)
-
 
 !     added by RW for nuketharsis
@@ -209 +194 @@
       allocate(tsoil(ngridmx,nsoilmx))
       allocate(ith(iip1,jjp1,nsoilmx),ithfi(ngridmx,nsoilmx))
-      allocate(tslab(ngridmx,nsoilmx))
       
 c=======================================================================
@@ -355 +339 @@
         CALL phyetat0 (ngridmx,llm,fichnom,tab0,Lmodif,nsoilmx,
      .        nqtot,day_ini,time,
-     .        tsurf,tsoil,emis,q2,qsurf,   !) ! temporary modif by RDW
-     .        cloudfrac,totalfrac,hice,rnat,pctsrf_sic,tslab,tsea_ice,
-     .        sea_ice)
+     .        tsurf,tsoil,emis,q2,qsurf)
 
         ! copy albedo and soil thermal inertia on (local) physics grid
@@ -549 +531 @@
      &   date,tsurf,tsoil,emis,q2,
      &   t,ucov,vcov,ps,teta,phisold_newgrid,q,qsurf,
-     &   surfith,nid,
-     &   rnat,pctsrf_sic,tslab,tsea_ice,sea_ice)
+     &   surfith,nid)
         write(*,*) "OK, read start_archive file"
         ! copy soil thermal inertia
@@ -784 +765 @@
         else if (trim(modif).eq.'ptot') then
 
-          ! check if we have a co2_ice surface tracer:
-          if (igcm_co2_ice.eq.0) then
-            write(*,*) " No surface CO2 ice !"
-            write(*,*) " only atmospheric pressure will be considered!"
-          endif
 c         calcul de la pression totale glace + atm actuelle
           patm=0.
@@ -800 +776 @@
                 patm = patm + ps(i,j)*aire(i,j)
                 airetot= airetot + aire(i,j)
-                if (igcm_co2_ice.ne.0) then
-                  !pcap = pcap + aire(i,j)*co2ice(ig)*g
-                  pcap = pcap + aire(i,j)*qsurf(ig,igcm_co2_ice)*g
-                endif
              ENDDO
           ENDDO
           ptoto = pcap + patm
 
-          print*,'Current total pressure at surface (co2 ice + atm) ',
+          print*,'Current total pressure at surface (atm) ',
      &       ptoto/airetot
 
@@ -1048 +1020 @@
        else if (trim(modif) .eq. 'wetstart') then
         ! check that there is indeed a water vapor tracer
-        if (igcm_h2o_vap.eq.0) then
+
           write(*,*) "No water vapour tracer! Can't use this option"
           stop
-        endif
-          DO l=1,llm
-            DO j=1,jjp1
-              DO i=1,iip1
-                q(i,j,l,igcm_h2o_vap)=150.e-6 * (rlatu(j)+pi/2.) / pi
-              ENDDO
-            ENDDO
-          ENDDO
-
-         write(*,*) 'Water mass mixing ratio at north pole='
-     *               ,q(1,1,1,igcm_h2o_vap)
-         write(*,*) '---------------------------south pole='
-     *               ,q(1,jjp1,1,igcm_h2o_vap)
 
 c      noglacier : remove tropical water ice (to initialize high res sim)
 c      --------------------------------------------------
         else if (trim(modif) .eq. 'noglacier') then
-           if (igcm_h2o_ice.eq.0) then
+        
              write(*,*) "No water ice tracer! Can't use this option"
              stop
-           endif
-           do ig=1,ngridmx
-             j=(ig-2)/iim +2
-              if(ig.eq.1) j=1
-              write(*,*) 'OK: remove surface ice for |lat|<45'
-              if (abs(rlatu(j)*180./pi).lt.45.) then
-                   qsurf(ig,igcm_h2o_ice)=0.
-              end if
-           end do
+
 
 
@@ -1085 +1036 @@
 c      --------------------------------------------------
         else if (trim(modif) .eq. 'watercapn') then
-           if (igcm_h2o_ice.eq.0) then
+        
              write(*,*) "No water ice tracer! Can't use this option"
              stop
-           endif
-
-           DO j=1,jjp1        
-              DO i=1,iim
-                 ig=1+(j-2)*iim +i
-                 if(j.eq.1) ig=1
-                 if(j.eq.jjp1) ig=ngridmx
-
-                 if (rlatu(j)*180./pi.gt.80.) then
-                    qsurf(ig,igcm_h2o_ice)=3.4e3
-                    !do isoil=1,nsoilmx
-                    !   ith(i,j,isoil) = 18000. ! thermal inertia
-                    !enddo
-                   write(*,*)'     ==> Ice mesh North boundary (deg)= ',
-     &                   rlatv(j-1)*180./pi
-                 end if
-              ENDDO
-           ENDDO
-           CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
-
-c$$$           do ig=1,ngridmx
-c$$$             j=(ig-2)/iim +2
-c$$$              if(ig.eq.1) j=1
-c$$$              if (rlatu(j)*180./pi.gt.80.) then
-c$$$
-c$$$                   qsurf(ig,igcm_h2o_ice)=1.e5
-c$$$                   qsurf(ig,igcm_h2o_vap)=0.0!1.e5
-c$$$
-c$$$                   write(*,*) 'ig=',ig,'    H2O ice mass (kg/m2)= ',
-c$$$     &              qsurf(ig,igcm_h2o_ice)
-c$$$
-c$$$                   write(*,*)'     ==> Ice mesh South boundary (deg)= ',
-c$$$     &              rlatv(j)*180./pi
-c$$$                end if
-c$$$           enddo
 
 c      watercaps : H20 ice on permanent southern cap
 c      -------------------------------------------------
         else if (trim(modif) .eq. 'watercaps') then
-           if (igcm_h2o_ice.eq.0) then
+
               write(*,*) "No water ice tracer! Can't use this option"
               stop
-           endif
-
-           DO j=1,jjp1        
-              DO i=1,iim
-                 ig=1+(j-2)*iim +i
-                 if(j.eq.1) ig=1
-                 if(j.eq.jjp1) ig=ngridmx
-
-                 if (rlatu(j)*180./pi.lt.-80.) then
-                    qsurf(ig,igcm_h2o_ice)=3.4e3
-                    !do isoil=1,nsoilmx
-                    !   ith(i,j,isoil) = 18000. ! thermal inertia
-                    !enddo
-                   write(*,*)'     ==> Ice mesh South boundary (deg)= ',
-     &                   rlatv(j-1)*180./pi
-                 end if
-              ENDDO
-           ENDDO
-           CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
-
-c$$$           do ig=1,ngridmx
-c$$$               j=(ig-2)/iim +2
-c$$$               if(ig.eq.1) j=1
-c$$$               if (rlatu(j)*180./pi.lt.-80.) then
-c$$$                  qsurf(ig,igcm_h2o_ice)=1.e5
-c$$$                  qsurf(ig,igcm_h2o_vap)=0.0 !1.e5
-c$$$
-c$$$                  write(*,*) 'ig=',ig,'   H2O ice mass (kg/m2)= ',
-c$$$     &                 qsurf(ig,igcm_h2o_ice)
-c$$$                  write(*,*)'     ==> Ice mesh North boundary (deg)= ',
-c$$$     &                 rlatv(j-1)*180./pi
-c$$$               end if
-c$$$           enddo
-
 
 c       noacglac : H2O ice across highest terrain
 c       --------------------------------------------
         else if (trim(modif) .eq. 'noacglac') then
-           if (igcm_h2o_ice.eq.0) then
+
              write(*,*) "No water ice tracer! Can't use this option"
              stop
-           endif
-          DO j=1,jjp1        
-             DO i=1,iim
-                ig=1+(j-2)*iim +i
-                if(j.eq.1) ig=1
-                if(j.eq.jjp1) ig=ngridmx
-
-                if(phis(i,j).gt.1000.*g)then
-                    alb(i,j) = 0.5 ! snow value
-                    do isoil=1,nsoilmx
-                       ith(i,j,isoil) = 18000. ! thermal inertia
-                       ! this leads to rnat set to 'ocean' in physiq.F90
-                       ! actually no, because it is soil not surface
-                    enddo
-                endif
-             ENDDO
-          ENDDO
-          CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
-          CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi)
-          CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
-
-
 
 c       oborealis : H2O oceans across Vastitas Borealis
 c       -----------------------------------------------
         else if (trim(modif) .eq. 'oborealis') then
-           if (igcm_h2o_ice.eq.0) then
+
              write(*,*) "No water ice tracer! Can't use this option"
              stop
-           endif
-          DO j=1,jjp1        
-             DO i=1,iim
-                ig=1+(j-2)*iim +i
-                if(j.eq.1) ig=1
-                if(j.eq.jjp1) ig=ngridmx
-
-                if(phis(i,j).lt.-4000.*g)then
-!                if( (phis(i,j).lt.-4000.*g)
-!     &               .and. (rlatu(j)*180./pi.lt.0.) )then ! south hemisphere only
-!                    phis(i,j)=-8000.0*g ! proper ocean
-                    phis(i,j)=-4000.0*g ! scanty ocean
-
-                    alb(i,j) = 0.07 ! oceanic value
-                    do isoil=1,nsoilmx
-                       ith(i,j,isoil) = 18000. ! thermal inertia
-                       ! this leads to rnat set to 'ocean' in physiq.F90
-                       ! actually no, because it is soil not surface
-                    enddo
-                endif
-             ENDDO
-          ENDDO
-          CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
-          CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi)
-          CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
-
+             
 c       iborealis : H2O ice in Northern plains
 c       --------------------------------------
         else if (trim(modif) .eq. 'iborealis') then
-           if (igcm_h2o_ice.eq.0) then
+
              write(*,*) "No water ice tracer! Can't use this option"
              stop
-           endif
-          DO j=1,jjp1        
-             DO i=1,iim
-                ig=1+(j-2)*iim +i
-                if(j.eq.1) ig=1
-                if(j.eq.jjp1) ig=ngridmx
-
-                if(phis(i,j).lt.-4000.*g)then
-                   !qsurf(ig,igcm_h2o_ice)=1.e3
-                   qsurf(ig,igcm_h2o_ice)=241.4 ! to make total 33 kg m^-2
-                endif
-             ENDDO
-          ENDDO
-          CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
-          CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi)
-          CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
-
 
 c       oceanball : H2O liquid everywhere
 c       ----------------------------
         else if (trim(modif) .eq. 'oceanball') then
-           if (igcm_h2o_ice.eq.0) then
+
              write(*,*) "No water ice tracer! Can't use this option"
              stop
-           endif
-          DO j=1,jjp1        
-             DO i=1,iim
-                ig=1+(j-2)*iim +i
-                if(j.eq.1) ig=1
-                if(j.eq.jjp1) ig=ngridmx
-
-                qsurf(ig,igcm_h2o_vap)=0.0    ! for ocean, this is infinite source
-                qsurf(ig,igcm_h2o_ice)=0.0
-                alb(i,j) = 0.07 ! ocean value
-
-                do isoil=1,nsoilmx
-                   ith(i,j,isoil) = 18000. ! thermal inertia
-                   !ith(i,j,isoil) = 50. ! extremely low for test
-                   ! this leads to rnat set to 'ocean' in physiq.F90
-                enddo
-
-             ENDDO
-          ENDDO
-          CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
-          CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi)
-          CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
 
 c       iceball : H2O ice everywhere
 c       ----------------------------
         else if (trim(modif) .eq. 'iceball') then
-           if (igcm_h2o_ice.eq.0) then
+
              write(*,*) "No water ice tracer! Can't use this option"
              stop
-           endif
-          DO j=1,jjp1        
-             DO i=1,iim
-                ig=1+(j-2)*iim +i
-                if(j.eq.1) ig=1
-                if(j.eq.jjp1) ig=ngridmx
-
-                qsurf(ig,igcm_h2o_vap)=-50.    ! for ocean, this is infinite source
-                qsurf(ig,igcm_h2o_ice)=50.     ! == to 0.05 m of oceanic ice
-                alb(i,j) = 0.6 ! ice albedo value
-
-                do isoil=1,nsoilmx
-                   !ith(i,j,isoil) = 18000. ! thermal inertia
-                   ! this leads to rnat set to 'ocean' in physiq.F90
-                enddo
-
-             ENDDO
-          ENDDO
-          CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
-          CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi)
 
 c       supercontinent : H2O ice everywhere
 c       ----------------------------
         else if (trim(modif) .eq. 'supercontinent') then
-          write(*,*) 'Minimum longitude (-180,180)'
-          read(*,*) val
-          write(*,*) 'Maximum longitude (-180,180)'
-          read(*,*) val2
-          write(*,*) 'Minimum latitude (-90,90)'
-          read(*,*) val3
-          write(*,*) 'Maximum latitude (-90,90)'
-          read(*,*) val4
-
-          do j=1,jjp1
-            do i=1,iip1
-              ig=1+(j-2)*iim +i
-              if(j.eq.1) ig=1
-              if(j.eq.jjp1) ig=ngridmx
-
-c             Supercontinent:
-              if (((rlatu(j)*180./pi.le.val4).and.
-     &            (rlatu(j)*180./pi.ge.val3).and.
-     &            (rlonv(i)*180./pi.le.val2).and.
-     &            (rlonv(i)*180./pi.ge.val))) then
-
-                rnat(ig)=1.
-                alb(i,j) = 0.3
-                do isoil=1,nsoilmx
-                  ith(i,j,isoil) = 2000.
-                enddo
-c             Ocean:
-              else
-                rnat(ig)=0.
-                alb(i,j) = 0.07
-                do isoil=1,nsoilmx
-                  ith(i,j,isoil) = 18000.
-                enddo
-              end if
-
-            enddo
-          enddo
-          CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
-          CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi)
+ 
+             write(*,*) "No water ice tracer! Can't use this option"
+             stop
 
 c       isotherm : Isothermal temperatures and no winds
@@ -1445 +1185 @@
 c       ------------------------------------------------
         else if (trim(modif) .eq. 'co2ice=0') then
-         ! check that there is indeed a co2_ice tracer ...
-          if (igcm_co2_ice.ne.0) then
-           do ig=1,ngridmx
-              !co2ice(ig)=0
-              qsurf(ig,igcm_co2_ice)=0
-              emis(ig)=emis(ngridmx/2)
-           end do
-          else
             write(*,*) "Can't remove CO2 ice!! (no co2_ice tracer)"
-          endif
+        
         
 !       therm_ini_s: (re)-set soil thermal inertia to reference surface values
@@ -1484 +1216 @@
 c       ------------------------------------------------
         else if (trim(modif) .eq. 'slab_ocean_0') then
-        write(*,*)'OK: initialisation of slab ocean' 
-
-      DO ig=1, ngridmx
-         rnat(ig)=1.
-         tslab(ig,1)=0.
-         tslab(ig,2)=0.
-         tsea_ice(ig)=0.
-         sea_ice(ig)=0.
-         pctsrf_sic(ig)=0.
-         
-         if(ithfi(ig,1).GT.10000.)then
-           rnat(ig)=0.
-           phisfi(ig)=0.
-           tsea_ice(ig)=273.15-1.8
-           tslab(ig,1)=tsurf(ig)
-           tslab(ig,2)=tsurf(ig)!*2./3.+(273.15-1.8)/3.
-          endif
-
-      ENDDO
-          CALL gr_fi_dyn(1,ngridmx,iip1,jjp1,phisfi,phis)
-
-
+        
+           write(*,*) "No ocean yet on Titan ! Can't use this option"
+           stop
 
         else
@@ -1517 +1230 @@
  999  continue
 
- 
-c=======================================================================
-c   Initialisation for cloud fraction and oceanic ice (added by BC 2010)
-c=======================================================================
-      DO ig=1, ngridmx
-         totalfrac(ig)=0.5
-         DO l=1,llm
-            cloudfrac(ig,l)=0.5
-         ENDDO
-! Ehouarn, march 2012: also add some initialisation for hice
-         hice(ig)=0.0
-      ENDDO
       
-c========================================================
-
-!      DO ig=1,ngridmx
-!         IF(tsurf(ig) .LT. 273.16-1.8) THEN
-!            hice(ig)=(273.16-1.8-tsurf(ig))/(273.16-1.8-240)*1.
-!            hice(ig)=min(hice(ig),1.0)
-!         ENDIF
-!      ENDDO
-
-
-
 
 c=======================================================================
@@ -1647 +1337 @@
       call physdem1("restartfi.nc",nsoilmx,ngridmx,llm,nqtot,
      &                dtphys,real(day_ini),
-     &                tsurf,tsoil,emis,q2,qsurf,
-     &                cloudfrac,totalfrac,hice,
-     &                rnat,pctsrf_sic,tslab,tsea_ice,sea_ice)
+     &                tsurf,tsoil,emis,q2,qsurf)
 
 c=======================================================================

trunk/LMDZ.TITAN/libf/dynphy_lonlat/phytitan/start2archive.F

@@ -25 +25 @@
 !      use control_mod
 !      use comgeomphy, only: initcomgeomphy
-      use slab_ice_h, only: noceanmx
 ! to use  'getin'
       USE ioipsl_getincom
@@ -31 +30 @@
       USE mod_const_mpi, ONLY: COMM_LMDZ
       USE control_mod, only: planet_type
-      USE callkeys_mod, ONLY: ok_slab_ocean
       use filtreg_mod, only: inifilr
       USE comvert_mod, ONLY: ap,bp
@@ -72 +70 @@
       REAL tsurf(ngridmx)        ! Surface temperature
       REAL,ALLOCATABLE :: tsoil(:,:) ! Soil temperature
-      REAL co2ice(ngridmx)        ! CO2 ice layer
       REAL q2(ngridmx,llm+1)
       REAL,ALLOCATABLE :: qsurf(:,:)
@@ -82 +79 @@
       INTEGER*4 day_ini_fi
 
-!     added by FF for cloud fraction setup
-      REAL hice(ngridmx)
-      REAL cloudfrac(ngridmx,llm),totalcloudfrac(ngridmx)
-
-!     added by BC for slab ocean
-      REAL rnat(ngridmx),pctsrf_sic(ngridmx),sea_ice(ngridmx)
-      REAL tslab(ngridmx,noceanmx),tsea_ice(ngridmx)
-
 
 c Variable naturelle / grille scalaire
@@ -97 +86 @@
       REAL,ALLOCATABLE :: tsoilS(:,:)
       REAL,ALLOCATABLE :: ithS(:,:) ! Soil Thermal Inertia
-      REAL co2iceS(ip1jmp1)
       REAL q2S(ip1jmp1,llm+1)
       REAL,ALLOCATABLE :: qsurfS(:,:)
       REAL emisS(ip1jmp1)
-
-!     added by FF for cloud fraction setup
-      REAL hiceS(ip1jmp1)
-      REAL cloudfracS(ip1jmp1,llm),totalcloudfracS(ip1jmp1)
-
-!     added by BC for slab ocean
-      REAL rnatS(ip1jmp1),pctsrf_sicS(ip1jmp1),sea_iceS(ip1jmp1)
-      REAL tslabS(ip1jmp1,noceanmx),tsea_iceS(ip1jmp1)
 
 
@@ -120 +100 @@
       INTEGER Lmodif
 
-      REAL ptotal, co2icetotal
+      REAL ptotal
       REAL timedyn,timefi !fraction du jour dans start, startfi
       REAL date
@@ -236 +216 @@
       CALL phyetat0 (ngridmx,llm,fichnom,0,Lmodif,nsoilmx,nqtot,
      .      day_ini_fi,timefi,
-     .      tsurf,tsoil,emis,q2,qsurf,
-!       change FF 05/2011
-     .       cloudfrac,totalcloudfrac,hice,
-!       change BC 05/2014
-     .       rnat,pctsrf_sic,tslab,tsea_ice,sea_ice)
-
+     .      tsurf,tsoil,emis,q2,qsurf)
 
 
@@ -328 +303 @@
 c
 c tsurf --> tsurfS
-c co2ice --> co2iceS
 c tsoil --> tsoilS
 c emis --> emisS
@@ -337 +311 @@
 
       call gr_fi_dyn(1,ngridmx,iip1,jjp1,tsurf,tsurfS)
-!      call gr_fi_dyn(1,ngridmx,iip1,jjp1,co2ice,co2iceS)
       call gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,tsoil,tsoilS)
       ! Note: thermal inertia "inertiedat" is in comsoil.h
@@ -344 +317 @@
       call gr_fi_dyn(llm+1,ngridmx,iip1,jjp1,q2,q2S)
       call gr_fi_dyn(nqtot,ngridmx,iip1,jjp1,qsurf,qsurfS)
-      call gr_fi_dyn(llm,ngridmx,iip1,jjp1,cloudfrac,cloudfracS)
-      call gr_fi_dyn(1,ngridmx,iip1,jjp1,hice,hiceS)
-      call gr_fi_dyn(1,ngridmx,iip1,jjp1,totalcloudfrac,totalcloudfracS)
-
-      call gr_fi_dyn(1,ngridmx,iip1,jjp1,rnat,rnatS)
-      call gr_fi_dyn(1,ngridmx,iip1,jjp1,pctsrf_sic,pctsrf_sicS)
-      call gr_fi_dyn(1,ngridmx,iip1,jjp1,tsea_ice,tsea_iceS)
-      call gr_fi_dyn(1,ngridmx,iip1,jjp1,sea_ice,sea_iceS)
-      call gr_fi_dyn(noceanmx,ngridmx,iip1,jjp1,tslab,tslabS)
 
 c=======================================================================
@@ -359 +323 @@
 
       ptotal =  0.
-      co2icetotal = 0.
       DO j=1,jjp1
          DO i=1,iim
            ptotal=ptotal+aire(i+(iim+1)*(j-1))*ps(i+(iim+1)*(j-1))/g
-!           co2icetotal = co2icetotal + 
-!     &            co2iceS(i+(iim+1)*(j-1))*aire(i+(iim+1)*(j-1))
          ENDDO
       ENDDO
       write(*,*)'Ancienne grille : masse de l''atm :',ptotal
-!      write(*,*)'Ancienne grille : masse de la glace CO2 :',co2icetotal
-
-c-----------------------------------------------------------------------
-c Passage de "ptotal" et "co2icetotal" par tab_cntrl_fi
+
+c-----------------------------------------------------------------------
+c Passage de "ptotal"  par tab_cntrl_fi
 c-----------------------------------------------------------------------
 
       tab_cntrl_fi(49) = ptotal
-      tab_cntrl_fi(50) = co2icetotal
+      tab_cntrl_fi(50) = 0.
 
 c=======================================================================
@@ -430 +390 @@
 
 c-----------------------------------------------------------------------
-c Ecriture des champs  (co2ice,emis,ps,Tsurf,T,u,v,q2,q,qsurf)
+c Ecriture des champs  (emis,ps,Tsurf,T,u,v,q2,q,qsurf)
 c-----------------------------------------------------------------------
 c ATTENTION: q2 a une couche de plus!!!!
@@ -439 +399 @@
 c-----------------------------------------------------------------------
 
-!      call write_archive(nid,ntime,'co2ice','couche de glace co2',
-!     &  'kg/m2',2,co2iceS)
       call write_archive(nid,ntime,'emis','grd emis',' ',2,emisS)
       call write_archive(nid,ntime,'ps','Psurf','Pa',2,ps)
@@ -497 +455 @@
 ! Note: no need to write volcapa, it is stored in "controle" table
 
-c-----------------------------------------------------------------------
-c Ecriture du champs  cloudfrac,hice,totalcloudfrac
-c-----------------------------------------------------------------------
-      call write_archive(nid,ntime,'hice',
-     &         'Height of oceanic ice','m',2,hiceS)
-      call write_archive(nid,ntime,'totalcloudfrac',
-     &        'Total cloud Fraction','',2,totalcloudfracS)
-      call write_archive(nid,ntime,'cloudfrac'
-     &        ,'Cloud fraction','',3,cloudfracS)
-
-c-----------------------------------------------------------------------
-c Slab ocean
-c-----------------------------------------------------------------------
-      OPEN(99,file='callphys.def',status='old',form='formatted'
-     &     ,iostat=ierr)
-      CLOSE(99)
-
-      IF(ierr.EQ.0) THEN
-
-
-         write(*,*) "Use slab-ocean ?"
-         ok_slab_ocean=.false.         ! default value
-         call getin("ok_slab_ocean",ok_slab_ocean)
-         write(*,*) "ok_slab_ocean = ",ok_slab_ocean
-
-      if(ok_slab_ocean) then
-      call write_archive(nid,ntime,'rnat'
-     &        ,'rnat','',2,rnatS)
-      call write_archive(nid,ntime,'pctsrf_sic'
-     &        ,'pctsrf_sic','',2,pctsrf_sicS)
-      call write_archive(nid,ntime,'sea_ice'
-     &        ,'sea_ice','',2,sea_iceS)
-      call write_archive(nid,ntime,'tslab'
-     &        ,'tslab','',-2,tslabS)
-      call write_archive(nid,ntime,'tsea_ice'
-     &        ,'tsea_ice','',2,tsea_iceS)
-      endif !ok_slab_ocean
-      ENDIF
-c-----------------------------------------------------------------------
+
 c Fin 
 c-----------------------------------------------------------------------

trunk/LMDZ.TITAN/libf/dynphy_lonlat/phytitan/write_archive.F

@@ -33 +33 @@
 
       use comsoil_h, only: nsoilmx
-      use slab_ice_h, only: noceanmx
 
       implicit none
@@ -195 +194 @@
         edges(1)=iip1
         edges(2)=jjp1
-        edges(3)=noceanmx
+        edges(3)=1 ! JVO2017 : was noceanmx before -> set to 1
         edges(4)=1
 #ifdef NC_DOUBLE

trunk/LMDZ.TITAN/libf/phytitan/aeropacity.F90

@@ -1 +1 @@
       Subroutine aeropacity(ngrid,nlayer,nq,pplay,pplev,pq, &
-         aerosol,reffrad,QREFvis3d,QREFir3d,tau_col,cloudfrac,totcloudfrac,clearsky)
+         aerosol,reffrad,QREFvis3d,QREFir3d,tau_col)
 
        use radinc_h, only : L_TAUMAX,naerkind
        use aerosol_mod
-       USE tracer_h, only: noms,rho_co2,rho_ice
+       USE tracer_h, only: noms
        use comcstfi_mod, only: g
-       use callkeys_mod, only: aerofixco2,aerofixh2o,kastprof,cloudlvl, &
-                CLFvarying,CLFfixval,dusttau,                           &
-                pres_bottom_tropo,pres_top_tropo,obs_tau_col_tropo,     &
+       use callkeys_mod, only: pres_bottom_tropo,pres_top_tropo,obs_tau_col_tropo,  &
                 pres_bottom_strato,pres_top_strato,obs_tau_col_strato
                   
@@ -55 +53 @@
       REAL,INTENT(IN) :: QREFir3d(ngrid,nlayer,naerkind)
       REAL,INTENT(OUT):: tau_col(ngrid) !column integrated visible optical depth
-      ! BENJAMIN MODIFS
-      real,intent(in) :: cloudfrac(ngrid,nlayer) ! cloud fraction
-      real,intent(out) :: totcloudfrac(ngrid) ! total cloud fraction
-      logical,intent(in) :: clearsky
 
       real aerosol0
@@ -69 +63 @@
       EXTERNAL CBRT
 
-      INTEGER,SAVE :: i_co2ice=0      ! co2 ice
-      INTEGER,SAVE :: i_h2oice=0      ! water ice
-!$OMP THREADPRIVATE(i_co2ice,i_h2oice)
-      CHARACTER(LEN=20) :: tracername ! to temporarily store text
-
       ! for fixed dust profiles
-      real topdust, expfactor, zp
-      REAL taudusttmp(ngrid) ! Temporary dust opacity used before scaling
-      REAL tauh2so4tmp(ngrid) ! Temporary h2so4 opacity used before scaling
-
-      real CLFtot
+      real expfactor, zp
 
       ! identify tracers
@@ -85 +70 @@
 
         write(*,*) "Tracers found in aeropacity:"
-        do iq=1,nq
-          tracername=noms(iq)
-          if (tracername.eq."co2_ice") then
-            i_co2ice=iq
-          write(*,*) "i_co2ice=",i_co2ice
-
-          endif
-          if (tracername.eq."h2o_ice") then
-            i_h2oice=iq
-            write(*,*) "i_h2oice=",i_h2oice
-          endif
-        enddo
 
         if (noaero) then
@@ -107 +80 @@
         endif
 
-        if ((iaero_co2.ne.0).and.(.not.noaero)) then
-          print*, 'iaero_co2=  ',iaero_co2
-        endif
-        if (iaero_h2o.ne.0) then
-          print*,'iaero_h2o=  ',iaero_h2o    
-        endif
-        if (iaero_dust.ne.0) then
-          print*,'iaero_dust= ',iaero_dust
-        endif
-        if (iaero_h2so4.ne.0) then
-          print*,'iaero_h2so4= ',iaero_h2so4
-        endif
         if (iaero_back2lay.ne.0) then
           print*,'iaero_back2lay= ',iaero_back2lay
@@ -127 +88 @@
 
 
-!     ---------------------------------------------------------
-!==================================================================
-!    CO2 ice aerosols
-!==================================================================
+!     ---------------------------------------------------------   
 
-      if (iaero_co2.ne.0) then
-           iaer=iaero_co2
-!       1. Initialization
-            aerosol(1:ngrid,1:nlayer,iaer)=0.0
-!       2. Opacity calculation
-            if (noaero) then ! aerosol set to zero
-             aerosol(1:ngrid,1:nlayer,iaer)=0.0
-            elseif (aerofixco2.or.(i_co2ice.eq.0)) then !  CO2 ice cloud prescribed
-               aerosol(1:ngrid,1:nlayer,iaer)=1.e-9
-               !aerosol(1:ngrid,12,iaer)=4.0 ! single cloud layer option
-            else
-               DO ig=1, ngrid
-                  DO l=1,nlayer-1 ! to stop the rad tran bug
-
-                     aerosol0 =                         &
-                          (  0.75 * QREFvis3d(ig,l,iaer) /        &
-                          ( rho_co2 * reffrad(ig,l,iaer) )  ) *   &
-                          ( pq(ig,l,i_co2ice) + 1.E-9 ) *         &
-                          ( pplev(ig,l) - pplev(ig,l+1) ) / g
-                     aerosol0           = max(aerosol0,1.e-9)
-                     aerosol0           = min(aerosol0,L_TAUMAX)
-                     aerosol(ig,l,iaer) = aerosol0
-!                     aerosol(ig,l,iaer) = 0.0
-!                     print*, aerosol(ig,l,iaer)
-!        using cloud fraction
-!                     aerosol(ig,l,iaer) = -log(1 - CLF + CLF*exp(-aerosol0/CLF))
-!                     aerosol(ig,l,iaer) = min(aerosol(ig,l,iaer),L_TAUMAX)
-
-
-                  ENDDO
-               ENDDO
-            end if ! if fixed or varying
-      end if ! if CO2 aerosols   
-!==================================================================
-!     Water ice / liquid 
-!==================================================================
-
-      if (iaero_h2o.ne.0) then 
-           iaer=iaero_h2o
-!       1. Initialization
-            aerosol(1:ngrid,1:nlayer,iaer)=0.0
-!       2. Opacity calculation
-            if (aerofixh2o.or.(i_h2oice.eq.0).or.clearsky) then
-               aerosol(1:ngrid,1:nlayer,iaer) =1.e-9
-
-               ! put cloud at cloudlvl
-               if(kastprof.and.(cloudlvl.ne.0.0))then
-                  ig=1
-                  do l=1,nlayer
-                     if(int(cloudlvl).eq.l)then
-                     !if(cloudlvl.gt.(pplay(ig,l)/pplev(ig,1)))then
-                        print*,'Inserting cloud at level ',l
-                        !aerosol(ig,l,iaer)=10.0
-
-                        rho_ice=920.0
-
-                        ! the Kasting approximation
-                        aerosol(ig,l,iaer) =                      &
-                          (  0.75 * QREFvis3d(ig,l,iaer) /        &
-                          ( rho_ice * reffrad(ig,l,iaer) )  ) *   &
-                          !( pq(ig,l,i_h2oice) + 1.E-9 ) *         &
-                          ( 4.0e-4 + 1.E-9 ) *         &
-                          ( pplev(ig,l) - pplev(ig,l+1) ) / g
-
-
-                        open(115,file='clouds.out',form='formatted')
-                        write(115,*) l,aerosol(ig,l,iaer)
-                        close(115)
-
-                        return
-                     endif
-                  end do
-
-                  call abort
-               endif
-
-            else
-
-               do ig=1, ngrid
-                  do l=1,nlayer-1 ! to stop the rad tran bug
-
-                     aerosol(ig,l,iaer) =                                    & !modification by BC
-                          (  0.75 * QREFvis3d(ig,l,iaer) /        &
-                          ( rho_ice * reffrad(ig,l,iaer) )  ) *   &
-                          !  pq(ig,l,i_h2oice) *                   & !JL I dropped the +1e-9 here to have the same
-                          !( pplev(ig,l) - pplev(ig,l+1) ) / g       !   opacity in the clearsky=true and the 
-                                                                     !   clear=false/pq=0 case
-                          ( pq(ig,l,i_h2oice) + 1.E-9 ) *         & ! Doing this makes the code unstable, so I have restored it (RW)
-                          ( pplev(ig,l) - pplev(ig,l+1) ) / g 
-
-                  enddo
-               enddo
-
-               if(CLFvarying)then
-                  call totalcloudfrac(ngrid,nlayer,nq,cloudfrac,totcloudfrac,pplev,pq,aerosol(1,1,iaer))
-                  do ig=1, ngrid
-                     do l=1,nlayer-1 ! to stop the rad tran bug
-                        CLFtot  = max(totcloudfrac(ig),0.01)
-                        aerosol(ig,l,iaer)=aerosol(ig,l,iaer)/CLFtot
-                        aerosol(ig,l,iaer) = max(aerosol(ig,l,iaer),1.e-9)
-                     enddo
-                  enddo
-               else
-                  do ig=1, ngrid
-                     do l=1,nlayer-1 ! to stop the rad tran bug
-                        CLFtot  = CLFfixval
-                        aerosol(ig,l,iaer)=aerosol(ig,l,iaer)/CLFtot
-                        aerosol(ig,l,iaer) = max(aerosol(ig,l,iaer),1.e-9)
-                     enddo
-                  enddo
-              end if!(CLFvarying)
-            endif !(aerofixed.or.(i_h2oice.eq.0).or.clearsky)
-              
-      end if ! End if h2o aerosol
-
-!==================================================================
-!             Dust 
-!==================================================================
-      if (iaero_dust.ne.0) then
-          iaer=iaero_dust
-!         1. Initialization 
-          aerosol(1:ngrid,1:nlayer,iaer)=0.0
-          
-          topdust=30.0 ! km  (used to be 10.0 km) LK
-
-!       2. Opacity calculation
-
-!           expfactor=0.
-           DO l=1,nlayer-1
-             DO ig=1,ngrid
-!             Typical mixing ratio profile
-
-                 zp=(pplev(ig,1)/pplay(ig,l))**(70./topdust)
-                 expfactor=max(exp(0.007*(1.-max(zp,1.))),1.e-3)
-
-!             Vertical scaling function
-              aerosol(ig,l,iaer)= (pplev(ig,l)-pplev(ig,l+1)) &
-               *expfactor
-
-
-             ENDDO
-           ENDDO
-
-!          Rescaling each layer to reproduce the choosen (or assimilated)
-!          dust extinction opacity at visible reference wavelength, which
-!          is scaled to the surface pressure pplev(ig,1)
-
-            taudusttmp(1:ngrid)=0.
-              DO l=1,nlayer
-                DO ig=1,ngrid
-                   taudusttmp(ig) = taudusttmp(ig) &
-                          +  aerosol(ig,l,iaer)
-                ENDDO
-              ENDDO
-            DO l=1,nlayer-1
-               DO ig=1,ngrid
-                  aerosol(ig,l,iaer) = max(1E-20, &
-                          dusttau &
-                       *  pplev(ig,1) / pplev(ig,1) &
-                       *  aerosol(ig,l,iaer) &
-                       /  taudusttmp(ig))
-
-              ENDDO
-            ENDDO
-      end if ! If dust aerosol   
-
-!==================================================================
-!           H2SO4 
-!==================================================================
-! added by LK
-      if (iaero_h2so4.ne.0) then
-         iaer=iaero_h2so4
-
-!       1. Initialization
-         aerosol(1:ngrid,1:nlayer,iaer)=0.0
-
-
-!       2. Opacity calculation
-
-!           expfactor=0.
-         DO l=1,nlayer-1
-            DO ig=1,ngrid
-!              Typical mixing ratio profile
-
-               zp=(pplev(ig,1)/pplay(ig,l))**(70./30) !emulating topdust
-               expfactor=max(exp(0.007*(1.-max(zp,1.))),1.e-3)
-
-!             Vertical scaling function
-               aerosol(ig,l,iaer)= (pplev(ig,l)-pplev(ig,l+1))*expfactor
-
-            ENDDO
-         ENDDO
-         tauh2so4tmp(1:ngrid)=0.
-         DO l=1,nlayer
-            DO ig=1,ngrid
-               tauh2so4tmp(ig) = tauh2so4tmp(ig) + aerosol(ig,l,iaer)
-            ENDDO
-         ENDDO
-         DO l=1,nlayer-1
-            DO ig=1,ngrid
-               aerosol(ig,l,iaer) = max(1E-20, &
-                          1 &
-                       *  pplev(ig,1) / pplev(ig,1) &
-                       *  aerosol(ig,l,iaer) &
-                       /  tauh2so4tmp(ig))
-
-            ENDDO
-         ENDDO
-
-! 1/700. is assuming a "sulfurtau" of 1
-! Sulfur aerosol routine to be improved.
-!                     aerosol0 =                         &
-!                          (  0.75 * QREFvis3d(ig,l,iaer) /        &
-!                          ( rho_h2so4 * reffrad(ig,l,iaer) )  ) *   &
-!                          ( pq(ig,l,i_h2so4) + 1.E-9 ) *         &
-!                          ( pplev(ig,l) - pplev(ig,l+1) ) / g
-!                     aerosol0           = max(aerosol0,1.e-9)
-!                     aerosol0           = min(aerosol0,L_TAUMAX)
-!                     aerosol(ig,l,iaer) = aerosol0
-
-!                  ENDDO
-!               ENDDO
-      end if
- 
-           
-!     ---------------------------------------------------------
 !==================================================================
 !    Two-layer aerosols (unknown composition)

trunk/LMDZ.TITAN/libf/phytitan/aerosol_mod.F90

@@ -8 +8 @@
 !                 corresponding aerosol was not activated in callphys.def
 !                 -- otherwise a value is given in iniaerosol
-      integer :: iaero_co2 = 0 
-      integer :: iaero_h2o = 0
-      integer :: iaero_dust = 0
-      integer :: iaero_h2so4 = 0
       logical :: noaero = .false.
 
 ! two-layer simple aerosol model
       integer :: iaero_back2lay = 0
-!$OMP THREADPRIVATE(iaero_co2,iaero_h2o,iaero_dust,iaero_h2so4,noaero,iaero_back2lay)
+!$OMP THREADPRIVATE(noaero,iaero_back2lay)
       
 !==================================================================

trunk/LMDZ.TITAN/libf/phytitan/calc_cpp_mugaz.F90

@@ -38 +38 @@
          else
             ! all values at 300 K from Engineering Toolbox
-            if(igas.eq.igas_CO2)then
-               mugaz_c = mugaz_c + 44.01*gfrac(igas)
-            elseif(igas.eq.igas_N2)then
+            if(igas.eq.igas_N2)then
                mugaz_c = mugaz_c + 28.01*gfrac(igas)
             elseif(igas.eq.igas_H2)then
                mugaz_c = mugaz_c + 2.01*gfrac(igas)
-            elseif(igas.eq.igas_He)then
-               mugaz_c = mugaz_c + 4.003*gfrac(igas)
-            elseif(igas.eq.igas_H2O)then
-               mugaz_c = mugaz_c + 18.02*gfrac(igas)
-            elseif(igas.eq.igas_SO2)then
-               mugaz_c = mugaz_c + 64.066*gfrac(igas)
-            elseif(igas.eq.igas_H2S)then
-               mugaz_c = mugaz_c + 34.08*gfrac(igas)
             elseif(igas.eq.igas_CH4)then
                mugaz_c = mugaz_c + 16.04*gfrac(igas)
-            elseif(igas.eq.igas_NH3)then
-               mugaz_c = mugaz_c + 17.03*gfrac(igas)
             elseif(igas.eq.igas_C2H6)then 
                ! C2H6 http://encyclopedia.airliquide.com/Encyclopedia.asp?GasID=28
@@ -77 +65 @@
          else
             ! all values at 300 K from Engineering Toolbox
-            if(igas.eq.igas_CO2)then
-               !cpp_c   = cpp_c   + 0.744*gfrac(igas) ! @ ~210 K (better for
-               !Mars conditions) 
-               cpp_c   = cpp_c   + 0.846*gfrac(igas)*44.01/mugaz_c
-            elseif(igas.eq.igas_N2)then
+            if(igas.eq.igas_N2)then
                cpp_c   = cpp_c   + 1.040*gfrac(igas)*28.01/mugaz_c
             elseif(igas.eq.igas_H2)then
                cpp_c   = cpp_c   + 14.31*gfrac(igas)*2.01/mugaz_c
-            elseif(igas.eq.igas_He)then
-               cpp_c   = cpp_c   + 5.19*gfrac(igas)*4.003/mugaz_c
-            elseif(igas.eq.igas_H2O)then
-               cpp_c   = cpp_c   + 1.864*gfrac(igas)*18.02/mugaz_c
-            elseif(igas.eq.igas_SO2)then
-               cpp_c   = cpp_c   + 0.64*gfrac(igas)*64.066/mugaz_c
-            elseif(igas.eq.igas_H2S)then
-               cpp_c   = cpp_c   + 1.003*gfrac(igas)*34.08/mugaz_c ! from wikipedia...
             elseif(igas.eq.igas_CH4)then
                cpp_c   = cpp_c   + 2.226*gfrac(igas)*16.04/mugaz_c
-            elseif(igas.eq.igas_NH3)then
-               cpp_c   = cpp_c   + 2.175*gfrac(igas)*17.03/mugaz_c
-               print*,'WARNING, cpp for NH3 may be for liquid'
             elseif(igas.eq.igas_C2H6)then
                ! C2H6

trunk/LMDZ.TITAN/libf/phytitan/calc_rayleigh.F90

@@ -62 +62 @@
             tauconsti(igas) = 0.0
          else
-            if(igas.eq.igas_CO2) then
-               tauconsti(igas) = (8.7/g)*1.527*scalep/P0
-            elseif(igas.eq.igas_N2)then
+            if(igas.eq.igas_N2)then
                tauconsti(igas) = (9.81/g)*8.569E-3*scalep/(P0/93.0)
-            elseif(igas.eq.igas_H2O)then
-!               tauconsti(igas) = (10.0/g)*9.22E-3*scalep/101325.0 
-               tauconsti(igas) = 1.98017E-10/(g*mugaz*100.) 
             elseif(igas.eq.igas_H2)then
                tauconsti(igas) = (10.0/g)*0.0241*scalep/101325.0
@@ -74 +69 @@
                ! uses n=1.000132 from Optics, Fourth Edition.
                ! was out by a factor 2!
-            elseif(igas.eq.igas_He)then
-               tauconsti(igas) = (10.0/g)*0.00086*scalep/101325.0
             else
                print*,'Error in calc_rayleigh: Gas species not recognised!'
@@ -114 +107 @@
                   tauvari(igas)   = 0.0
                else
-                  if(igas.eq.igas_CO2)then
-                     tauvari(igas) = (1.0+0.013/wl**2)/wl**4
-                  elseif(igas.eq.igas_N2)then
+                  if(igas.eq.igas_N2)then
                      tauvari(igas) = (1.0+0.0113/wl**2+0.00013/wl**4)/wl**4
-                  elseif(igas.eq.igas_H2O)then
-!                     tauvari(igas) = 1.0/wl**4 ! to be improved...
-                     tauvari(igas) = (5.7918E6/(2.38E2-1/wl**2)+1.679E5/(57.36E0-1/wl**2))**2/wl**4 
                   elseif(igas.eq.igas_H2)then
-                     tauvari(igas) = 1.0/wl**4 
-                  elseif(igas.eq.igas_He)then
                      tauvari(igas) = 1.0/wl**4 
                   else

trunk/LMDZ.TITAN/libf/phytitan/callcorrk.F90

@@ -6 +6 @@
           fluxabs_sw,fluxtop_dn,                               &
           OLR_nu,OSR_nu,                                       &
-          tau_col,cloudfrac,totcloudfrac,                      &
-          clearsky,firstcall,lastcall)
+          tau_col,firstcall,lastcall)
 
       use radinc_h
       use radcommon_h
-      use watercommon_h
       use datafile_mod, only: datadir
       use ioipsl_getin_p_mod, only: getin_p
       use gases_h
-      use radii_mod, only : su_aer_radii,co2_reffrad,h2o_reffrad,dust_reffrad,h2so4_reffrad,back2lay_reffrad
-      use aerosol_mod, only : iaero_co2,iaero_h2o,iaero_dust,iaero_h2so4, iaero_back2lay
+      use radii_mod, only : su_aer_radii,back2lay_reffrad
+      use aerosol_mod, only : iaero_back2lay
       USE tracer_h
       use comcstfi_mod, only: pi, mugaz, cpp
-      use callkeys_mod, only: varactive,diurnal,tracer,water,nosurf,varfixed,satval,        &
-                              kastprof,strictboundcorrk,specOLR,CLFvarying
+      use callkeys_mod, only: diurnal,tracer,nosurf,satval,        &
+                              strictboundcorrk,specOLR
 
       implicit none
@@ -62 +60 @@
       REAL,INTENT(IN) :: dist_star                 ! Distance star-planet (AU).
       REAL,INTENT(IN) :: muvar(ngrid,nlayer+1)
-      REAL,INTENT(IN) :: cloudfrac(ngrid,nlayer)   ! Fraction of clouds (%).
-      logical,intent(in) :: clearsky
       logical,intent(in) :: firstcall              ! Signals first call to physics.
       logical,intent(in) :: lastcall               ! Signals last call to physics.
@@ -80 +76 @@
       REAL,INTENT(OUT) :: OSR_nu(ngrid,L_NSPECTV)        ! Outgoing SW radition in each band (Normalized to the band width (W/m2/cm-1).
       REAL,INTENT(OUT) :: tau_col(ngrid)                 ! Diagnostic from aeropacity.
-      REAL,INTENT(OUT) :: albedo_equivalent(ngrid)       ! Spectrally Integrated Albedo. For Diagnostic. By MT2015
-      REAL,INTENT(OUT) :: totcloudfrac(ngrid)            ! Column Fraction of clouds (%).
+      REAL,INTENT(OUT) :: albedo_equivalent(ngrid)       ! Spectrally Integrated Albedo. For Diagnostic. By MT2015  
       
       
-      
-      
-
       ! Globally varying aerosol optical properties on GCM grid ; not needed everywhere so not in radcommon_h.   
       REAL :: QVISsQREF3d(ngrid,nlayer,L_NSPECTV,naerkind)
@@ -230 +222 @@
         
 
-         if((igcm_h2o_vap.eq.0) .and. varactive)then
-            print*,'varactive in callcorrk but no h2o_vap tracer.'
-            stop
-         endif
-
          OLR_nu(:,:) = 0.
          OSR_nu(:,:) = 0.
@@ -278 +265 @@
 
       do iaer=1,naerkind
-
-         if ((iaer.eq.iaero_co2).and.tracer.and.(igcm_co2_ice.gt.0)) then ! Treat condensed co2 particles.
-            call co2_reffrad(ngrid,nlayer,nq,pq,reffrad(1,1,iaero_co2))
-            print*,'Max. CO2 ice particle size = ',maxval(reffrad(1:ngrid,1:nlayer,iaer))/1.e-6,' um'
-            print*,'Min. CO2 ice particle size = ',minval(reffrad(1:ngrid,1:nlayer,iaer))/1.e-6,' um'
-         end if
-         
-         if ((iaer.eq.iaero_h2o).and.water) then ! Treat condensed water particles. To be generalized for other aerosols ...
-            call h2o_reffrad(ngrid,nlayer,pq(1,1,igcm_h2o_ice),pt, &
-                             reffrad(1,1,iaero_h2o),nueffrad(1,1,iaero_h2o))
-            print*,'Max. H2O cloud particle size = ',maxval(reffrad(1:ngrid,1:nlayer,iaer))/1.e-6,' um'
-            print*,'Min. H2O cloud particle size = ',minval(reffrad(1:ngrid,1:nlayer,iaer))/1.e-6,' um'
-         endif
-         
-         if(iaer.eq.iaero_dust)then
-            call dust_reffrad(ngrid,nlayer,reffrad(1,1,iaero_dust))
-            print*,'Dust particle size = ',reffrad(1,1,iaer)/1.e-6,' um'
-         endif
-         
-         if(iaer.eq.iaero_h2so4)then
-            call h2so4_reffrad(ngrid,nlayer,reffrad(1,1,iaero_h2so4))
-            print*,'H2SO4 particle size =',reffrad(1,1,iaer)/1.e-6,' um'
-         endif
-         
+     
           if(iaer.eq.iaero_back2lay)then
             call back2lay_reffrad(ngrid,reffrad(1,1,iaero_back2lay),nlayer,pplev)
@@ -323 +287 @@
       call aeropacity(ngrid,nlayer,nq,pplay,pplev,pq,aerosol,      &
            reffrad,QREFvis3d,QREFir3d,                             & 
-           tau_col,cloudfrac,totcloudfrac,clearsky)                
+           tau_col)                
           
 
@@ -506 +470 @@
 
 
-!-----------------------------------------------------------------------
-!     Water vapour (to be generalised for other gases eventually ...)
-!-----------------------------------------------------------------------
+     do k=1,L_LEVELS
+         qvar(k) = 1.0D-7
+     end do
+
+
+       ! IMPORTANT: Now convert from kg/kg to mol/mol.
+!       do k=1,L_LEVELS
+!          qvar(k) = qvar(k)/(epsi+qvar(k)*(1.-epsi))
+!       end do
+
+       DO l=1,nlayer
+          muvarrad(2*l)   = muvar(ig,nlayer+2-l)
+          muvarrad(2*l+1) = (muvar(ig,nlayer+2-l)+muvar(ig,max(nlayer+1-l,1)))/2
+       END DO
       
-      if(varactive)then
-
-         i_var=igcm_h2o_vap
-         do l=1,nlayer
-            qvar(2*l)   = pq(ig,nlayer+1-l,i_var)
-            qvar(2*l+1) = pq(ig,nlayer+1-l,i_var)    
-!JL13index            qvar(2*l+1) = (pq(ig,nlayer+1-l,i_var)+pq(ig,max(nlayer-l,1),i_var))/2    
-!JL13index            ! Average approximation as for temperature...
-         end do
-         qvar(1)=qvar(2)
-
-      elseif(varfixed)then
-
-         do l=1,nlayer ! Here we will assign fixed water vapour profiles globally.
-            RH = satval * ((pplay(ig,l)/pplev(ig,1) - 0.02) / 0.98)
-            if(RH.lt.0.0) RH=0.0
-            
-            ptemp=pplay(ig,l)
-            Ttemp=pt(ig,l)
-            call watersat(Ttemp,ptemp,qsat)
-
-            !pq_temp(l) = qsat      ! fully saturated everywhere
-            pq_temp(l) = RH * qsat ! ~realistic profile (e.g. 80% saturation at ground)
-         end do
-         
-         do l=1,nlayer
-            qvar(2*l)   = pq_temp(nlayer+1-l)
-            qvar(2*l+1) = (pq_temp(nlayer+1-l)+pq_temp(max(nlayer-l,1)))/2
-         end do
-         
-         qvar(1)=qvar(2)
-
-         ! Lowest layer of atmosphere
-         RH = satval * (1 - 0.02) / 0.98
-         if(RH.lt.0.0) RH=0.0
-
-!         ptemp = pplev(ig,1)
-!         Ttemp = tsurf(ig)
-!         call watersat(Ttemp,ptemp,qsat)
-
-         qvar(2*nlayer+1)= RH * qsat ! ~realistic profile (e.g. 80% saturation at ground)
- 
-      else
-         do k=1,L_LEVELS
-            qvar(k) = 1.0D-7
-         end do
-      end if ! varactive/varfixed
-
-      if(.not.kastprof)then
-         ! IMPORTANT: Now convert from kg/kg to mol/mol.
-         do k=1,L_LEVELS
-            qvar(k) = qvar(k)/(epsi+qvar(k)*(1.-epsi))
-         end do
-      end if
-
-!-----------------------------------------------------------------------
-!     kcm mode only !
-!-----------------------------------------------------------------------
-
-      if(kastprof)then
-
-         ! Initial values equivalent to mugaz.
-         DO l=1,nlayer
-            muvarrad(2*l)   = mugaz
-            muvarrad(2*l+1) = mugaz
-         END DO
-
-         if(ngasmx.gt.1)then
-
-            DO l=1,nlayer
-               muvarrad(2*l)   =  muvar(ig,nlayer+2-l)
-               muvarrad(2*l+1) = (muvar(ig,nlayer+2-l) + &
-                                  muvar(ig,max(nlayer+1-l,1)))/2
-            END DO
-      
-            muvarrad(1) = muvarrad(2)
-            muvarrad(2*nlayer+1) = muvar(ig,1)
-
-            print*,'Recalculating qvar with VARIABLE epsi for kastprof'
-            print*,'Assumes that the variable gas is H2O!!!'
-            print*,'Assumes that there is only one tracer'
-            
-            !i_var=igcm_h2o_vap
-            i_var=1
-            
-            if(nq.gt.1)then
-               print*,'Need 1 tracer only to run kcm1d.e' 
-               stop
-            endif
-            
-            do l=1,nlayer
-               vtmp(l)=pq(ig,l,i_var)/(epsi+pq(ig,l,i_var)*(1.-epsi)) 
-               !vtmp(l)=pq(ig,l,i_var)*muvar(ig,l+1)/mH2O !JL to be changed
-            end do
-
-            do l=1,nlayer
-               qvar(2*l)   = vtmp(nlayer+1-l)
-               qvar(2*l+1) = vtmp(nlayer+1-l)
-!               qvar(2*l+1) = ( vtmp(nlayer+1-l) + vtmp(max(nlayer-l,1)) )/2
-            end do
-            qvar(1)=qvar(2)
-
-            print*,'Warning: reducing qvar in callcorrk.'
-            print*,'Temperature profile no longer consistent ', &
-                   'with saturated H2O. qsat=',satval
-                   
-            do k=1,L_LEVELS
-               qvar(k) = qvar(k)*satval
-            end do
-
-         endif
-      else ! if kastprof
-         DO l=1,nlayer
-            muvarrad(2*l)   = muvar(ig,nlayer+2-l)
-            muvarrad(2*l+1) = (muvar(ig,nlayer+2-l)+muvar(ig,max(nlayer+1-l,1)))/2
-         END DO
-      
-         muvarrad(1) = muvarrad(2)
-         muvarrad(2*nlayer+1)=muvar(ig,1)         
-      endif ! if kastprof
+       muvarrad(1) = muvarrad(2)
+       muvarrad(2*nlayer+1)=muvar(ig,1)         
       
       ! Keep values inside limits for which we have radiative transfer coefficients !!!
@@ -949 +806 @@
 
       ! See physiq.F for explanations about CLFvarying. This is temporary.
-      if (lastcall .and. .not.CLFvarying) then
+      if (lastcall) then
         IF( ALLOCATED( gasi ) ) DEALLOCATE( gasi )
         IF( ALLOCATED( gasv ) ) DEALLOCATE( gasv )

trunk/LMDZ.TITAN/libf/phytitan/callkeys_mod.F90

@@ -4 +4 @@
       logical,save :: callrad,corrk,calldifv,UseTurbDiff
 !$OMP THREADPRIVATE(callrad,corrk,calldifv,UseTurbDiff)
-      logical,save :: calladj,co2cond,callsoil
-!$OMP THREADPRIVATE(calladj,co2cond,callsoil)
+      logical,save :: calladj,callsoil
+!$OMP THREADPRIVATE(calladj,callsoil)
       logical,save :: season,diurnal,tlocked,rings_shadow,lwrite
 !$OMP THREADPRIVATE(season,diurnal,tlocked,rings_shadow,lwrite)
       logical,save :: callstats,calleofdump
 !$OMP THREADPRIVATE(callstats,calleofdump)
-      logical,save :: callgasvis,continuum,H2Ocont_simple,graybody
-!$OMP THREADPRIVATE(callgasvis,continuum,H2Ocont_simple,graybody)
+      logical,save :: callgasvis,continuum,graybody
+!$OMP THREADPRIVATE(callgasvis,continuum,graybody)
       logical,save :: strictboundcorrk                                     
 !$OMP THREADPRIVATE(strictboundcorrk)
@@ -19 +19 @@
       logical,save :: meanOLR
       logical,save :: specOLR
-      logical,save :: kastprof
-!$OMP THREADPRIVATE(enertest,nonideal,meanOLR,kastprof)
+!$OMP THREADPRIVATE(enertest,nonideal,meanOLR)
       logical,save :: newtonian
       logical,save :: check_cpp_match
@@ -29 +28 @@
       logical,save :: stelbbody
       logical,save :: ozone
-      logical,save :: nearco2cond
       logical,save :: tracer
       logical,save :: mass_redistrib
-!$OMP THREADPRIVATE(stelbbody,ozone,nearco2cond,tracer,mass_redistrib)
-      logical,save :: varactive
-      logical,save :: varfixed
-      logical,save :: radfixed
+!$OMP THREADPRIVATE(stelbbody,ozone,tracer,mass_redistrib)
       logical,save :: sedimentation
-!$OMP THREADPRIVATE(varactive,varfixed,radfixed,sedimentation)
-      logical,save :: water,watercond,waterrain
-!$OMP THREADPRIVATE(water,watercond,waterrain)
-      logical,save :: aeroco2,aeroh2o,aeroh2so4,aeroback2lay
-!$OMP THREADPRIVATE(aeroco2,aeroh2o,aeroh2so4,aeroback2lay)
-      logical,save :: aerofixco2,aerofixh2o
-!$OMP THREADPRIVATE(aerofixco2,aerofixh2o)
-      logical,save :: hydrology
-      logical,save :: sourceevol
-      logical,save :: CLFvarying
+!$OMP THREADPRIVATE(sedimentation)
+      logical,save :: aeroback2lay
+!$OMP THREADPRIVATE(aeroback2lay)
       logical,save :: nosurf
       logical,save :: oblate
-!$OMP THREADPRIVATE(hydrology,sourceevol,CLFvarying,nosurf,oblate)
-      logical,save :: ok_slab_ocean
-      logical,save :: ok_slab_sic
-      logical,save :: ok_slab_heat_transp
-      logical,save :: albedo_spectral_mode
-!$OMP THREADPRIVATE(ok_slab_ocean,ok_slab_sic,ok_slab_heat_transp,albedo_spectral_mode)
+!$OMP THREADPRIVATE(nosurf,oblate)
 
       integer,save :: iddist
@@ -63 +46 @@
 
       real,save :: topdustref
-      real,save :: Nmix_co2
-      real,save :: dusttau
       real,save :: Fat1AU
       real,save :: stelTbb
-!$OMP THREADPRIVATE(topdustref,Nmix_co2,dusttau,Fat1AU,stelTbb)
-      real,save :: Tstrat
+!$OMP THREADPRIVATE(topdustref,Fat1AU,stelTbb)
       real,save :: tplanet
       real,save :: obs_tau_col_tropo
       real,save :: obs_tau_col_strato
-!$OMP THREADPRIVATE(Tstrat,tplanet,obs_tau_col_tropo,obs_tau_col_strato)
+!$OMP THREADPRIVATE(tplanet,obs_tau_col_tropo,obs_tau_col_strato)
       real,save :: pres_bottom_tropo
       real,save :: pres_top_tropo
@@ -81 +61 @@
       real,save :: size_strato
       real,save :: satval
-      real,save :: CLFfixval
-      real,save :: n2mixratio
-!$OMP THREADPRIVATE(size_tropo,size_strato,satval,CLFfixval,n2mixratio)
-      real,save :: co2supsat
+!$OMP THREADPRIVATE(size_tropo,size_strato,satval)
       real,save :: pceil
-      real,save :: albedosnow
-      real,save :: albedoco2ice
-      real,save :: maxicethick
-!$OMP THREADPRIVATE(co2supsat,pceil,albedosnow,albedoco2ice,maxicethick)
-      real,save :: Tsaldiff
+!$OMP THREADPRIVATE(pceil)
       real,save :: tau_relax
-      real,save :: cloudlvl
-      real,save :: icetstep
       real,save :: intheat
-!$OMP THREADPRIVATE(Tsaldiff,tau_relax,cloudlvl,icetstep,intheat)
+!$OMP THREADPRIVATE(tau_relax,intheat)
       real,save :: flatten
       real,save :: Rmean

trunk/LMDZ.TITAN/libf/phytitan/callsedim.F

@@ -4 +4 @@
 
       use radinc_h, only : naerkind
-      use radii_mod, only: h2o_reffrad
-      use aerosol_mod, only : iaero_h2o
-      USE tracer_h, only : igcm_co2_ice,igcm_h2o_ice,radius,rho_q
+      USE tracer_h, only : radius, rho_q
       use comcstfi_mod, only: g
-      use callkeys_mod, only : water
 
       IMPLICIT NONE
@@ -62 +59 @@
       real epaisseur (ngrid,nlay) ! Layer thickness (m)
       real wq(ngrid,nlay+1) ! displaced tracer mass (kg.m-2)
-c      real dens(ngrid,nlay) ! Mean density of the ice part. accounting for dust core
 
 
@@ -73 +69 @@
       IF (firstcall) THEN
         firstcall=.false.
-        ! add some tests on presence of required tracers/aerosols:
-        if (water) then
-          if (igcm_h2o_ice.eq.0) then
-            write(*,*) "callsedim error: water=.true.",
-     &                 " but igcm_h2o_ice=0"
-          stop
-          endif
-          if (iaero_h2o.eq.0) then
-            write(*,*) "callsedim error: water=.true.",
-     &                 " but iaero_ho2=0"
-          stop
-          endif
-        endif
       ENDIF ! of IF (firstcall)
       
@@ -106 +89 @@
  
       do iq=1,nq
-       if((radius(iq).gt.1.e-9).and.(iq.ne.igcm_co2_ice)) then   
-!         (no sedim for gases, and co2_ice sedim is done in condense_co2)      
+       if(radius(iq).gt.1.e-9) then   
+!         (no sedim for gases)      
 
 ! store locally updated tracers
@@ -120 +103 @@
 ! Sedimentation 
 !======================================================================
-! Water          
-          if (water.and.(iq.eq.igcm_h2o_ice)) then
-            ! compute radii for h2o_ice 
-             call h2o_reffrad(ngrid,nlay,zqi(1,1,igcm_h2o_ice),zt,
-     &                reffrad(1,1,iaero_h2o),nueffrad(1,1,iaero_h2o))
-            ! call sedimentation for h2o_ice
-             call newsedim(ngrid,nlay,ngrid*nlay,ptimestep,
-     &            pplev,masse,epaisseur,zt,reffrad(1,1,iaero_h2o),
-     &            rho_q(iq),zqi(1,1,igcm_h2o_ice),wq)
-
 ! General Case
-          else 
              call newsedim(ngrid,nlay,1,ptimestep,
      &            pplev,masse,epaisseur,zt,radius(iq),rho_q(iq),
      &            zqi(1,1,iq),wq)
-          endif
 
 !=======================================================================
@@ -152 +123 @@
             ENDDO
           ENDDO
-       endif ! of no gases no co2_ice
+       endif ! of no gases
       enddo ! of do iq=1,nq
       return

trunk/LMDZ.TITAN/libf/phytitan/comsoil_h.F90

@@ -8 +8 @@
 ! Full soil layer depths are set as: layer(k)=lay1_soil*alpha_soil**(k-1) , k=1,nsoil
 ! Mid soil layer depths are set as: mlayer(k)=lay1_soil*alpha_soil**(k-1/2) , k=0,nsoil-1 
-  real,save :: lay1_soil=2.e-4 ! depth (m) of first full soil layer (may be set in callphys.def)
+  real,save :: lay1_soil=2.e-3 ! depth (m) of first full soil layer (may be set in callphys.def)
   real,save :: alpha_soil=2 ! coefficient for soil layer thickness (may be set in callphys.def)
 !$OMP THREADPRIVATE(nsoilmx,lay1_soil,alpha_soil)

trunk/LMDZ.TITAN/libf/phytitan/convadj.F

@@ -8 +8 @@
       USE tracer_h
       use comcstfi_mod, only: g
-      use callkeys_mod, only: tracer,water
+      use callkeys_mod, only: tracer
 
       implicit none
@@ -16 +16 @@
 !     Purpose
 !     -------
-!     Calculates dry convective adjustment. If one tracer is CO2,
-!     we take into account the molecular mass variation (e.g. when
-!     CO2 condenses) to trigger convection (F. Forget 01/2005)
+!     Calculates dry convective adjustment.
 !
 !     Authors
@@ -62 +60 @@
 
 !     Tracers
-      INTEGER iq,ico2
-      save ico2
-!$OMP THREADPRIVATE(ico2)
+      INTEGER iq
       REAL zq(ngrid,nlay,nq), zq2(ngrid,nlay,nq)
-      REAL zqm(nq),zqco2m
-      real m_co2, m_noco2, A , B
+      REAL zqm(nq)
+      real A , B
       save A, B
 !$OMP THREADPRIVATE(A,B)
@@ -73 +69 @@
       real mtot1, mtot2 , mm1, mm2
        integer l1ref, l2ref
-      LOGICAL vtest(ngrid),down,firstcall
+      LOGICAL vtest(ngrid),down, firstcall
       save firstcall
       data firstcall/.true./
@@ -88 +84 @@
 
       IF (firstcall) THEN 
-        ico2=0
-        if (tracer) then
-!     Prepare Special treatment if one of the tracers is CO2 gas
-           do iq=1,nq
-             if (noms(iq).eq."co2") then
-                print*,'dont go there'
-                stop
-                ico2=iq
-                m_co2 = 44.01E-3  ! CO2 molecular mass (kg/mol)   
-                m_noco2 = 33.37E-3  ! Non condensible mol mass (kg/mol)   
-!               Compute A and B coefficient use to compute
-!               mean molecular mass Mair defined by
-!               1/Mair = q(ico2)/m_co2 + (1-q(ico2))/m_noco2
-!               1/Mair = A*q(ico2) + B
-                A =(1/m_co2 - 1/m_noco2)
-                B=1/m_noco2
-             end if
-           enddo
-        endif
         firstcall=.false.
       ENDIF ! of IF (firstcall)
-
+      
       DO l=1,nlay
          DO ig=1,ngrid
@@ -142 +119 @@
       ENDDO
 
-      if (ico2.ne.0) then
-!     Special case if one of the tracers is CO2 gas
-         DO l=1,nlay
-           DO ig=1,ngrid
-             zhc(ig,l) = zh2(ig,l)*(A*zq2(ig,l,ico2)+B)
-           ENDDO
-         ENDDO
-       else
-          CALL scopy(ngrid*nlay,zh2,1,zhc,1)
-       end if
+      CALL scopy(ngrid*nlay,zh2,1,zhc,1)
 
 !     Find out which grid points are convectively unstable
@@ -203 +171 @@
               zdsm = dsig(l2)
               zhm = zh2(i, l2)
-              if(ico2.ne.0) zqco2m = zq2(i,l2,ico2)
 
 !     Test loop downwards
@@ -211 +178 @@
                 zdsm = zdsm + dsig(l)
                 zhm = zhm + sdsig(l) * (zh2(i, l) - zhm) / zsm
-                if(ico2.ne.0) then
-                  zqco2m = 
-     &            zqco2m + dsig(l) * (zq2(i,l,ico2) - zqco2m) / zdsm
-                  zhmc = zhm*(A*zqco2m+B)
-                else 
-                  zhmc = zhm
-                end if
+                zhmc = zhm
  
 !     do we have to extend the column downwards?
@@ -260 +221 @@
               end do
               DO l = l1, l2
-                if(ico2.ne.0) then
-                  zalpha=zalpha+
-     &            ABS(zhc(i,l)/(A+B*zqco2m) -zhm)*dsig(l)
-                else
-                  zalpha=zalpha+ABS(zh2(i,l)-zhm)*dsig(l)
-                endif
+                zalpha=zalpha+ABS(zh2(i,l)-zhm)*dsig(l)
                 zh2(i, l) = zhm
 !     modifs by RDW !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -306 +262 @@
                  end do
               ENDDO
-              if (ico2.ne.0) then
-                DO l=l1, l2
-                  zhc(i,l) = zh2(i,l)*(A*zq2(i,l,ico2)+B)
-                ENDDO
-              end if
 
 
@@ -324 +275 @@
 !     check conservation
          cadjncons=0.0
-         if(water)then
-         do l = 1, nlay
-            masse = (pplev(i,l) - pplev(i,l+1))/g
-            iq    = igcm_h2o_vap
-            cadjncons = cadjncons + 
-     &           masse*(zq2(i,l,iq)-zq(i,l,iq))/ptimestep 
-         end do
-         endif
 
          if(cadjncons.lt.-1.e-6)then
@@ -369 +312 @@
             print*,'zh2(ig,:)        = ',zh2(i,l)
          end do
-         do l = 1, nlay
-            print*,'zq(ig,:,vap)     = ',zq(i,l,igcm_h2o_vap)
-         end do
-         do l = 1, nlay
-            print*,'zq2(ig,:,vap)    = ',zq2(i,l,igcm_h2o_vap)
-         end do
-            print*,'zqm(vap)         = ',zqm(igcm_h2o_vap)
+!         do l = 1, nlay
+!            print*,'zq(ig,:,vap)     = ',zq(i,l,igcm_h2o_vap)
+!         end do
+!         do l = 1, nlay
+!            print*,'zq2(ig,:,vap)    = ',zq2(i,l,igcm_h2o_vap)
+!         end do
+!            print*,'zqm(vap)         = ',zqm(igcm_h2o_vap)
             print*,'jadrs=',jadrs
 

trunk/LMDZ.TITAN/libf/phytitan/gases_h.F90

@@ -19 +19 @@
       ! in analogy with tracer.h ...
       integer :: igas_H2
-      integer :: igas_He
-      integer :: igas_H2O
-      integer :: igas_CO2
-      integer :: igas_CO
       integer :: igas_N2
-      integer :: igas_O2
-      integer :: igas_SO2
-      integer :: igas_H2S
       integer :: igas_CH4
-      integer :: igas_NH3
       integer :: igas_C2H2
       integer :: igas_C2H6
 !!$OMP THREADPRIVATE(ngasmx,vgas,gnom,gfrac,&
-!       !$OMP igas_H2,igas_He,igas_H2O,igas_CO2,igas_CO,igas_N2,&
-!       !$OMP igas_O2,igas_SO2,igas_H2S,igas_CH4,igas_NH3,igas_C2H2,igas_C2H6)
+!       !$OMP igas_H2,igas_N2,igas_CH4,igas_C2H2,igas_C2H6)
 
       end module gases_h

trunk/LMDZ.TITAN/libf/phytitan/iniaerosol.F

@@ -4 +4 @@
       use radinc_h, only: naerkind
       use aerosol_mod
-      use callkeys_mod, only: aeroco2,aeroh2o,dusttau,aeroh2so4,
-     &          aeroback2lay
+      use callkeys_mod, only: aeroback2lay
 
       IMPLICIT NONE
@@ -22 +21 @@
 
       ia=0
-      if (aeroco2) then
-         ia=ia+1
-         iaero_co2=ia
-      endif
-      write(*,*) '--- CO2 aerosol = ', iaero_co2
- 
-      if (aeroh2o) then
-         ia=ia+1
-         iaero_h2o=ia
-      endif
-      write(*,*) '--- H2O aerosol = ', iaero_h2o
-
-      if (dusttau.gt.0) then
-         ia=ia+1
-         iaero_dust=ia
-      endif
-      write(*,*) '--- Dust aerosol = ', iaero_dust
-
-      if (aeroh2so4) then
-         ia=ia+1
-         iaero_h2so4=ia
-      endif
-      write(*,*) '--- H2SO4 aerosol = ', iaero_h2so4
-      
+       
       if (aeroback2lay) then
          ia=ia+1
@@ -59 +35 @@
 
       if (ia.eq.0) then  !For the zero aerosol case. 
+         noaero = .true.
          ia = 1
-         noaero = .true.
-         iaero_co2=ia
       endif
 

trunk/LMDZ.TITAN/libf/phytitan/inifis_mod.F90

@@ -235 +235 @@
      call getin_p("continuum",continuum)
      write(*,*) " continuum = ",continuum
-
-     write(*,*) "use analytic function for H2O continuum ?"
-     H2Ocont_simple=.false. ! default value
-     call getin_p("H2Ocont_simple",H2Ocont_simple)
-     write(*,*) " H2Ocont_simple = ",H2Ocont_simple
  
      write(*,*) "call turbulent vertical diffusion ?"
@@ -255 +250 @@
      call getin_p("calladj",calladj)
      write(*,*) " calladj = ",calladj
-
-     write(*,*) "call CO2 condensation ?"
-     co2cond=.false. ! default value
-     call getin_p("co2cond",co2cond)
-     write(*,*) " co2cond = ",co2cond
-! Test of incompatibility
-     if (co2cond.and.(.not.tracer)) then
-        print*,'We need a CO2 ice tracer to condense CO2'
-        call abort
-     endif 
- 
-     write(*,*) "CO2 supersaturation level ?"
-     co2supsat=1.0 ! default value
-     call getin_p("co2supsat",co2supsat)
-     write(*,*) " co2supsat = ",co2supsat
 
      write(*,*) "Radiative timescale for Newtonian cooling ?"
@@ -313 +293 @@
      write(*,*)" specOLR = ",specOLR
 
-     write(*,*)"Operate in kastprof mode?"
-     kastprof=.false.
-     call getin_p("kastprof",kastprof)
-     write(*,*)" kastprof = ",kastprof
-
      write(*,*)"Uniform absorption in radiative transfer?"
      graybody=.false.
@@ -339 +314 @@
      write(*,*)" alpha_soil = ",alpha_soil
 
-! Slab Ocean 
-     write(*,*) "Use slab-ocean ?"
-     ok_slab_ocean=.false.         ! default value
-     call getin_p("ok_slab_ocean",ok_slab_ocean)
-     write(*,*) "ok_slab_ocean = ",ok_slab_ocean
-     ! Sanity check: for now slab oncean only works in serial mode
-     if (ok_slab_ocean.and.is_parallel) then
-       write(*,*) " Error: slab ocean should only be used in serial mode!"
-       call abort
-     endif
-
-     write(*,*) "Use slab-sea-ice ?"
-     ok_slab_sic=.true.         ! default value
-     call getin_p("ok_slab_sic",ok_slab_sic)
-     write(*,*) "ok_slab_sic = ",ok_slab_sic
-
-     write(*,*) "Use heat transport for the ocean ?"
-     ok_slab_heat_transp=.true.   ! default value
-     call getin_p("ok_slab_heat_transp",ok_slab_heat_transp)
-     write(*,*) "ok_slab_heat_transp = ",ok_slab_heat_transp
-
-
-
-! Test of incompatibility:
-! if kastprof used, we must be in 1D
-     if (kastprof.and.(ngrid.gt.1)) then
-       print*,'kastprof can only be used in 1D!'
-       call abort
-     endif
-
-     write(*,*)"Stratospheric temperature for kastprof mode?"
-     Tstrat=167.0
-     call getin_p("Tstrat",Tstrat)
-     write(*,*)" Tstrat = ",Tstrat
-
      write(*,*)"Remove lower boundary?"
      nosurf=.false.
@@ -441 +381 @@
 ! TRACERS:
 
-     write(*,*)"Varying H2O cloud fraction?"
-     CLFvarying=.false.     ! default value
-     call getin_p("CLFvarying",CLFvarying)
-     write(*,*)" CLFvarying = ",CLFvarying
-
-     write(*,*)"Value of fixed H2O cloud fraction?"
-     CLFfixval=1.0                ! default value
-     call getin_p("CLFfixval",CLFfixval)
-     write(*,*)" CLFfixval = ",CLFfixval
-
-     write(*,*)"fixed radii for Cloud particles?"
-     radfixed=.false. ! default value
-     call getin_p("radfixed",radfixed)
-     write(*,*)" radfixed = ",radfixed
-
-     if(kastprof)then
-        radfixed=.true.
-     endif  
-
-     write(*,*)"Number mixing ratio of CO2 ice particles:"
-     Nmix_co2=1.e6 ! default value
-     call getin_p("Nmix_co2",Nmix_co2)
-     write(*,*)" Nmix_co2 = ",Nmix_co2
-
 !         write(*,*)"Number of radiatively active aerosols:"
 !         naerkind=0. ! default value
@@ -470 +386 @@
 !         write(*,*)" naerkind = ",naerkind
 
-     write(*,*)"Opacity of dust (if used):"
-     dusttau=0. ! default value
-     call getin_p("dusttau",dusttau)
-     write(*,*)" dusttau = ",dusttau
-
-     write(*,*)"Radiatively active CO2 aerosols?"
-     aeroco2=.false.     ! default value
-     call getin_p("aeroco2",aeroco2)
-     write(*,*)" aeroco2 = ",aeroco2
-
-     write(*,*)"Fixed CO2 aerosol distribution?"
-     aerofixco2=.false.     ! default value
-     call getin_p("aerofixco2",aerofixco2)
-     write(*,*)" aerofixco2 = ",aerofixco2
-
-     write(*,*)"Radiatively active water ice?"
-     aeroh2o=.false.     ! default value
-     call getin_p("aeroh2o",aeroh2o)
-     write(*,*)" aeroh2o = ",aeroh2o
-
-     write(*,*)"Fixed H2O aerosol distribution?"
-     aerofixh2o=.false.     ! default value
-     call getin_p("aerofixh2o",aerofixh2o)
-     write(*,*)" aerofixh2o = ",aerofixh2o
-
-     write(*,*)"Radiatively active sulfuric acid aersols?"
-     aeroh2so4=.false.     ! default value
-     call getin_p("aeroh2so4",aeroh2so4)
-     write(*,*)" aeroh2so4 = ",aeroh2so4
-         
+ 
 !=================================
 
@@ -553 +440 @@
 !=================================
 
-     write(*,*)"Cloud pressure level (with kastprof only):"
-     cloudlvl=0. ! default value
-     call getin_p("cloudlvl",cloudlvl)
-     write(*,*)" cloudlvl = ",cloudlvl
-
-     write(*,*)"Is the variable gas species radiatively active?"
-     Tstrat=167.0
-     varactive=.false.
-     call getin_p("varactive",varactive)
-     write(*,*)" varactive = ",varactive
-
-     write(*,*)"Is the variable gas species distribution set?"
-     varfixed=.false.
-     call getin_p("varfixed",varfixed)
-     write(*,*)" varfixed = ",varfixed
-
      write(*,*)"What is the saturation % of the variable species?"
      satval=0.8
@@ -574 +445 @@
      write(*,*)" satval = ",satval
 
-
-! Test of incompatibility:
-! if varactive, then varfixed should be false
-     if (varactive.and.varfixed) then
-       print*,'if varactive, varfixed must be OFF!'
-       stop
-     endif
-
      write(*,*) "Gravitationnal sedimentation ?"
      sedimentation=.false. ! default value
      call getin_p("sedimentation",sedimentation)
-     write(*,*) " sedimentation = ",sedimentation
-
-     write(*,*) "Compute water cycle ?"
-     water=.false. ! default value
-     call getin_p("water",water)
-     write(*,*) " water = ",water
-         
-! Test of incompatibility:
-! if water is true, there should be at least a tracer
-     if (water.and.(.not.tracer)) then
-        print*,'if water is ON, tracer must be ON too!'
-        stop
-     endif
-
-     write(*,*) "Include water condensation ?"
-     watercond=.false. ! default value
-     call getin_p("watercond",watercond)
-     write(*,*) " watercond = ",watercond
-
-! Test of incompatibility:
-! if watercond is used, then water should be used too
-     if (watercond.and.(.not.water)) then
-        print*,'if watercond is used, water should be used too'
-        stop
-     endif
-
-     write(*,*) "Include water precipitation ?"
-     waterrain=.false. ! default value
-     call getin_p("waterrain",waterrain)
-     write(*,*) " waterrain = ",waterrain
-
-     write(*,*) "Include surface hydrology ?"
-     hydrology=.false. ! default value
-     call getin_p("hydrology",hydrology)
-     write(*,*) " hydrology = ",hydrology
-
-     write(*,*) "Evolve surface water sources ?"
-     sourceevol=.false. ! default value
-     call getin_p("sourceevol",sourceevol)
-     write(*,*) " sourceevol = ",sourceevol
-
-     write(*,*) "Ice evolution timestep ?"
-     icetstep=100.0 ! default value
-     call getin_p("icetstep",icetstep)
-     write(*,*) " icetstep = ",icetstep
-         
-     write(*,*) "Spectral Dependant albedo ?"
-     albedo_spectral_mode=.false. ! default value
-     call getin_p("albedo_spectral_mode",albedo_spectral_mode)
-     write(*,*) " albedo_spectral_mode = ",albedo_spectral_mode
-
-     write(*,*) "Snow albedo ?"
-     write(*,*) "If albedo_spectral_mode=.true., then this "
-     write(*,*) "corresponds to the 0.5 microns snow albedo."
-     albedosnow=0.5         ! default value
-     call getin_p("albedosnow",albedosnow)
-     write(*,*) " albedosnow = ",albedosnow
-         
-     write(*,*) "CO2 ice albedo ?"
-     albedoco2ice=0.5       ! default value
-     call getin_p("albedoco2ice",albedoco2ice)
-     write(*,*) " albedoco2ice = ",albedoco2ice
-
-     write(*,*) "Maximum ice thickness ?"
-     maxicethick=2.0         ! default value
-     call getin_p("maxicethick",maxicethick)
-     write(*,*) " maxicethick = ",maxicethick
-
-     write(*,*) "Freezing point of seawater ?"
-     Tsaldiff=-1.8          ! default value
-     call getin_p("Tsaldiff",Tsaldiff)
-     write(*,*) " Tsaldiff = ",Tsaldiff
+     write(*,*) " sedimentation = ",sedimentation       
 
      write(*,*) "Does user want to force cpp and mugaz?"

trunk/LMDZ.TITAN/libf/phytitan/initracer.F

@@ -3 +3 @@
       use surfdat_h
       USE tracer_h
-      USE callkeys_mod, only: water
       IMPLICIT NONE
 c=======================================================================
@@ -25 +24 @@
 
 !      real qsurf(ngrid,nq)       ! tracer on surface (e.g.  kg.m-2)
-!      real co2ice(ngrid)           ! co2 ice mass on surface (e.g.  kg.m-2)
       character(len=20) :: txt ! to store some text
       integer iq,ig,count
@@ -40 +38 @@
 c  alpha_lift(nq)  ! saltation vertical flux/horiz flux ratio (m-1)
 c  alpha_devil(nq) ! lifting coeeficient by dust devil
-c  rho_dust          ! Mars dust density
-c  rho_ice           ! Water ice density
 c  doubleq           ! if method with mass (iq=1) and number(iq=2) mixing ratio
 c  varian            ! Characteristic variance of log-normal distribution
@@ -76 +72 @@
       ! 0. initialize tracer indexes to zero:
       ! NB: igcm_* indexes are commons in 'tracer.h'
-      do iq=1,nq
-        igcm_dustbin(iq)=0
-      enddo
-      igcm_dust_mass=0
-      igcm_dust_number=0
-      igcm_h2o_vap=0
-      igcm_h2o_ice=0
-      igcm_co2=0
       igcm_co=0
       igcm_o=0
@@ -97 +85 @@
       igcm_ar=0
       igcm_ar_n2=0
-      igcm_co2_ice=0
 
       write(*,*) 'initracer: noms() ', noms
@@ -133 +120 @@
 !      endif ! of if (doubleq)
       ! 2. find chemistry and water tracers
-      do iq=1,nq
-        if (noms(iq).eq."co2") then
-          igcm_co2=iq
-          mmol(igcm_co2)=44.
-          count=count+1
-!          write(*,*) 'co2: count=',count
-        endif
-        if (noms(iq).eq."co2_ice") then
-          igcm_co2_ice=iq
-          mmol(igcm_co2_ice)=44.
-          count=count+1
-!          write(*,*) 'co2_ice: count=',count
-        endif
-        if (noms(iq).eq."h2o_vap") then
-          igcm_h2o_vap=iq
-          mmol(igcm_h2o_vap)=18.
-          count=count+1
-!          write(*,*) 'h2o_vap: count=',count
-        endif
-        if (noms(iq).eq."h2o_ice") then
-          igcm_h2o_ice=iq
-          mmol(igcm_h2o_ice)=18.
-          count=count+1
-!          write(*,*) 'h2o_ice: count=',count
-        endif
-      enddo ! of do iq=1,nq
-      
+     
       ! check that we identified all tracers:
       if (count.ne.nq) then
@@ -181 +142 @@
       call zerophys(nq,rho_q)
 
-      rho_dust=2500.  ! Mars dust density (kg.m-3)
-      rho_ice=920.    ! Water ice density (kg.m-3)
-      rho_co2=1620.   ! CO2 ice density (kg.m-3)
-
-
-
-c$$$      if (doubleq) then
-c$$$c       "doubleq" technique 
-c$$$c       -------------------
-c$$$c      (transport of mass and number mixing ratio)
-c$$$c       iq=1: Q mass mixing ratio, iq=2: N number mixing ratio
-c$$$
-c$$$        if( (nq.lt.2).or.(water.and.(nq.lt.3)) ) then
-c$$$          write(*,*)'initracer: nq is too low : nq=', nq
-c$$$          write(*,*)'water= ',water,' doubleq= ',doubleq   
-c$$$        end if
-c$$$
-c$$$        varian=0.637    ! Characteristic variance    
-c$$$        qext(igcm_dust_mass)=3.04   ! reference extinction at 0.67 um for ref dust
-c$$$        qext(igcm_dust_number)=3.04 ! reference extinction at 0.67 um for ref dust
-c$$$        rho_q(igcm_dust_mass)=rho_dust
-c$$$        rho_q(igcm_dust_number)=rho_dust
-c$$$
-c$$$c       Intermediate calcul for computing geometric mean radius r0
-c$$$c       as a function of mass and number mixing ratio Q and N
-c$$$c       (r0 = (r3n_q * Q/ N)
-c$$$        r3n_q = exp(-4.5*varian**2)*(3./4.)/(pi*rho_dust)
-c$$$
-c$$$c       Intermediate calcul for computing effective radius reff
-c$$$c       from geometric mean radius r0
-c$$$c       (reff = ref_r0 * r0)
-c$$$        ref_r0 = exp(2.5*varian**2)
-c$$$        
-c$$$c       lifted dust :
-c$$$c       '''''''''''
-c$$$        reff_lift = 3.e-6      !  Effective radius of lifted dust (m)
-c$$$        alpha_devil(igcm_dust_mass)=9.e-9   !  dust devil lift mass coeff
-c$$$        alpha_lift(igcm_dust_mass)=3.0e-15  !  Lifted mass coeff
-c$$$
-c$$$        r0_lift = reff_lift/ref_r0
-c$$$        alpha_devil(igcm_dust_number)=r3n_q*
-c$$$     &                        alpha_devil(igcm_dust_mass)/r0_lift**3
-c$$$        alpha_lift(igcm_dust_number)=r3n_q*
-c$$$     &                        alpha_lift(igcm_dust_mass)/r0_lift**3
-c$$$
-c$$$c       Not used:
-c$$$        radius(igcm_dust_mass) = 0.
-c$$$        radius(igcm_dust_number) = 0.
-c$$$
-c$$$      else
-
-
-c$$$       if (dustbin.gt.1) then
-c$$$        print*,'ATTENTION:',
-c$$$     $   ' properties of dust need input in initracer !!!'
-c$$$        stop
-c$$$
-c$$$       else if (dustbin.eq.1) then
-c$$$
-c$$$c       This will be used for 1 dust particle size:
-c$$$c       ------------------------------------------
-c$$$        radius(igcm_dustbin(1))=3.e-6
-c$$$        Qext(igcm_dustbin(1))=3.04
-c$$$        alpha_lift(igcm_dustbin(1))=0.0e-6
-c$$$        alpha_devil(igcm_dustbin(1))=7.65e-9
-c$$$        qextrhor(igcm_dustbin(1))=(3./4.)*Qext(igcm_dustbin(1))
-c$$$     &                         /(rho_dust*radius(igcm_dustbin(1)))
-c$$$        rho_q(igcm_dustbin(1))=rho_dust
-c$$$
-c$$$       endif
-c$$$!      end if    ! (doubleq)
-
-c     Initialization for water vapor
-c     ------------------------------
-      if(water) then
-         radius(igcm_h2o_vap)=0.
-         Qext(igcm_h2o_vap)=0.
-         alpha_lift(igcm_h2o_vap) =0.
-         alpha_devil(igcm_h2o_vap)=0.
-         qextrhor(igcm_h2o_vap)= 0.
-
-c       "Dryness coefficient" controlling the evaporation and
-c        sublimation from the ground water ice (close to 1)
-c        HERE, the goal is to correct for the fact
-c        that the simulated permanent water ice polar caps
-c        is larger than the actual cap and the atmospheric
-c        opacity not always realistic.
-
-
-!         if(ngrid.eq.1)
-
-
-!     to be modified for BC+ version?
-
-         !! this is defined in surfdat_h.F90
-         IF (.not.ALLOCATED(dryness)) ALLOCATE(dryness(ngrid))
-         IF (.not.ALLOCATED(watercaptag)) ALLOCATE(watercaptag(ngrid))
-
-         do ig=1,ngrid
-           if (ngrid.ne.1) watercaptag(ig)=.false.
-           dryness(ig) = 1.
-         enddo
-
-
-
-
-!         IF (caps) THEN
-c Perennial H20 north cap defined by watercaptag=true (allows surface to be
-c hollowed by sublimation in vdifc).
-!         do ig=1,ngrid
-!           if (lati(ig)*180./pi.gt.84) then
-!             if (ngrid.ne.1) watercaptag(ig)=.true.
-!             dryness(ig) = 1.
-c Use the following cap definition for high spatial resolution (latitudinal bin <= 5 deg)
-c             if (lati(ig)*180./pi.lt.85.and.long(ig).ge.0) then
-c               if (ngrid.ne.1) watercaptag(ig)=.true.
-c               dryness(ig) = 1.
-c             endif
-c             if (lati(ig)*180./pi.ge.85) then
-c               if (ngrid.ne.1) watercaptag(ig)=.true.
-c               dryness(ig) = 1.
-c             endif
-!           endif  ! (lati>80 deg)
-!         end do ! (ngrid)
-!        ENDIF ! (caps)
-
-!         if(iceparty.and.(nq.ge.2)) then
-
-           radius(igcm_h2o_ice)=3.e-6
-           rho_q(igcm_h2o_ice)=rho_ice
-           Qext(igcm_h2o_ice)=0.
-!           alpha_lift(igcm_h2o_ice) =0.
-!           alpha_devil(igcm_h2o_ice)=0.
-           qextrhor(igcm_h2o_ice)= (3./4.)*Qext(igcm_h2o_ice) 
-     $       / (rho_ice*radius(igcm_h2o_ice))
-
-
-
-!         elseif(iceparty.and.(nq.lt.2)) then
-!            write(*,*) 'nq is too low : nq=', nq
-!            write(*,*) 'water= ',water,' iceparty= ',iceparty
-!         endif
-
-      end if  ! (water)
 
 c     Output for records:

trunk/LMDZ.TITAN/libf/phytitan/iostart.F90

@@ -15 +15 @@
     INTEGER,SAVE :: idim6 ! "nlayer" dimension
     INTEGER,SAVE :: idim7 ! "Time" dimension
-    INTEGER,SAVE :: idim8 ! "ocean_layers" dimension
     INTEGER,SAVE :: timeindex ! current time index (for time-dependent fields)
 !$OMP THREADPRIVATE(idim1,idim2,idim3,idim4,idim5,idim6,idim7,timeindex)
@@ -469 +468 @@
   USE tracer_h, only: nqtot
   USE comsoil_h, only: nsoilmx
-  USE slab_ice_h, only: noceanmx
 
   IMPLICIT NONE
@@ -558 +556 @@
       ENDIF
 
-      ierr=NF90_DEF_DIM(nid_restart,"ocean_layers",noceanmx,idim8)
-      IF (ierr/=NF90_NOERR) THEN
-        write(*,*)'open_restartphy: problem defining oceanic layer dimension '
-        write(*,*)trim(nf90_strerror(ierr))
-        CALL ABORT
-      ENDIF
-
-
       ierr=NF90_ENDDEF(nid_restart)
       IF (ierr/=NF90_NOERR) THEN
@@ -646 +636 @@
   USE mod_grid_phy_lmdz
   USE mod_phys_lmdz_para
-  USE slab_ice_h, only: noceanmx
   IMPLICIT NONE
   CHARACTER(LEN=*),INTENT(IN)    :: field_name
@@ -832 +821 @@
         endif ! of if (.not.present(time))
 
-      ELSE IF (field_size==noceanmx) THEN
-        ! input is a 2D "oceanic field" array
-        if (.not.present(time)) then ! for a time-independent field
-          ierr = NF90_REDEF(nid_restart)
-#ifdef NC_DOUBLE
-          ierr=NF90_DEF_VAR(nid_restart,field_name,NF90_DOUBLE,&
-                            (/idim2,idim8/),nvarid)
-#else
-          ierr=NF90_DEF_VAR(nid_restart,field_name,NF90_FLOAT,&
-                            (/idim2,idim8/),nvarid)
-#endif
-          if (ierr.ne.NF90_NOERR) then
-            write(*,*)"put_field_rgen error: failed to define "//trim(field_name)
-            write(*,*)trim(nf90_strerror(ierr))
-          endif
-          IF (LEN_TRIM(title) > 0) ierr=NF90_PUT_ATT(nid_restart,nvarid,"title",title)
-          ierr = NF90_ENDDEF(nid_restart)
-          ierr = NF90_PUT_VAR(nid_restart,nvarid,field_glo)
-        else
-          ! check if the variable has already been defined:
-          ierr=NF90_INQ_VARID(nid_restart,field_name,nvarid)
-          if (ierr/=NF90_NOERR) then ! variable not found, define it
-            ierr=NF90_REDEF(nid_restart)
-#ifdef NC_DOUBLE
-            ierr=NF90_DEF_VAR(nid_restart,field_name,NF90_DOUBLE,&
-                              (/idim2,idim8,idim7/),nvarid)
-#else
-            ierr=NF90_DEF_VAR(nid_restart,field_name,NF90_FLOAT,&
-                              (/idim2,idim8,idim7/),nvarid)
-#endif
-           if (ierr.ne.NF90_NOERR) then
-              write(*,*)"put_field_rgen error: failed to define "//trim(field_name)
-              write(*,*)trim(nf90_strerror(ierr))
-            endif
-            IF (LEN_TRIM(title) > 0) ierr=NF90_PUT_ATT(nid_restart,nvarid,"title",title)
-            ierr=NF90_ENDDEF(nid_restart)
-          endif
-          ! Write the variable
-          ierr=NF90_PUT_VAR(nid_restart,nvarid,field_glo,&
-                            start=(/1,1,timeindex/))
-
-        endif ! of if (.not.present(time))
-
-
       ELSE
         PRINT *, "Error phyredem(put_field_rgen) : wrong dimension for ",trim(field_name)
@@ -948 +893 @@
   USE comsoil_h, only: nsoilmx
   USE mod_phys_lmdz_para, only: is_master
-  USE slab_ice_h, only: noceanmx
   IMPLICIT NONE
      CHARACTER(LEN=*),INTENT(IN) :: var_name
@@ -1000 +944 @@
         ! We know it is an  "mlayer" kind of 1D array
         idim1d=idim3
-      ELSEIF (var_size==noceanmx) THEN
-        ! We know it is an  "mlayer" kind of 1D array
-        idim1d=idim8
       ELSE 
         PRINT *, "put_var_rgen error : wrong dimension"

trunk/LMDZ.TITAN/libf/phytitan/mass_redistribution.F90

@@ -1 +1 @@
 SUBROUTINE mass_redistribution(ngrid,nlayer,nq,ptimestep,   &
-                       rnat,pcapcal,pplay,pplev,pt,ptsrf,pq,pqs,     &
+                       pcapcal,pplay,pplev,pt,ptsrf,pq,pqs,     &
                        pu,pv,pdt,pdtsrf,pdq,pdu,pdv,pdmassmr,  &
                        pdtmr,pdtsrfmr,pdpsrfmr,pdumr,pdvmr,pdqmr,pdqsmr)
                                                    
-       USE watercommon_h, Only: Tsat_water,RLVTT
        use surfdat_h
        use radcommon_h, only: glat
@@ -10 +9 @@
        USE planete_mod, only: bp
        use comcstfi_mod, only: g
-       USE callkeys_mod, ONLY: water
        
        IMPLICIT NONE
@@ -57 +55 @@
       INTEGER ngrid, nlayer, nq   
       REAL ptimestep
-      REAL pcapcal(ngrid)
-      INTEGER rnat(ngrid)      
+      REAL pcapcal(ngrid)   
       REAL pplay(ngrid,nlayer),pplev(ngrid,nlayer+1)
       REAL pt(ngrid,nlayer),pdt(ngrid,nlayer)
@@ -76 +73 @@
 
 !    Boiling/sublimation
-      REAL Tsat(ngrid),zmassboil(ngrid)
+      REAL zmassboil(ngrid)
 
 !    vertical reorganization of sigma levels
@@ -132 +129 @@
 
 
-      if (water) then
-         do ig=1,ngrid
-            call Tsat_water(pplev(ig,1)+zdmass_sum(ig,1)*g*ptimestep,Tsat(ig))
-         enddo
-               call writediagfi(ngrid,'Tsat','saturation temperature at surface','',2,Tsat)
-         
-         do ig=1,ngrid
-            if (ztsrf(ig).gt.Tsat(ig)) then
-               zmassboil(ig)=(ptsrf(ig)-Tsat(ig))*pcapcal(ig)/RLVTT/ptimestep
-               if ((zmassboil(ig)*ptimestep.gt.pqs(ig,igcm_h2o_vap)).and.(rnat(ig).eq.1)) then
-                  zmassboil(ig)=pqs(ig,igcm_h2o_vap)/ptimestep
-               endif
-               zmassboil(ig)=zmassboil(ig)*0.0 !momentary, should be 1. JL12
-               pdqsmr(ig,igcm_h2o_vap)=-zmassboil(ig)
-               pdtsrfmr(ig)=-zmassboil(ig)*RLVTT/pcapcal(ig)
-               ztsrf(ig)=ptsrf(ig)+pdtsrfmr(ig)*ptimestep
-            else
-               zmassboil(ig)=0.
-               pdtsrfmr(ig)=0.
-            endif
-         enddo
-      endif
-
 !     *************************
 !           UPDATE SURFACE
@@ -219 +193 @@
          zvm(1) = 0.  
          zqm(1,1:nq)=0. ! most tracer do not condense !
-         if (water) zqm(1,igcm_h2o_vap)=1. ! flux is 100% h2o at surface
          
 !        Van Leer scheme:

trunk/LMDZ.TITAN/libf/phytitan/optci.F90

@@ -7 +7 @@
   use gases_h
   use comcstfi_mod, only: g, r, mugaz
-  use callkeys_mod, only: kastprof,continuum,graybody,H2Ocont_simple
+  use callkeys_mod, only: continuum,graybody
   implicit none
 
@@ -77 +77 @@
   integer interm
 
-  !--- Kasting's CIA ----------------------------------------
-  !real*8, parameter :: Ci(L_NSPECTI)=[                         &
-  !     3.8E-5, 1.2E-5, 2.8E-6, 7.6E-7, 4.5E-7, 2.3E-7,    &
-  !     5.4E-7, 1.6E-6, 0.0,                               &
-  !     0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,            & 
-  !     0.0, 4.0E-7, 4.0E-6, 1.4E-5,    &
-  !     1.0E-5, 1.2E-6, 2.0E-7, 5.0E-8, 3.0E-8, 0.0 ] 
-  !real*8, parameter :: Ti(L_NSPECTI)=[ -2.2, -1.9,             &
-  !     -1.7, -1.7, -1.7, -1.7, -1.7, -1.7,                &
-  !     0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, &
-  !     -1.7,-1.7,-1.7,-1.7,-1.7,-1.7,-1.7, -1.7,0.0 ]
-  !----------------------------------------------------------
-
   !! AS: to save time in computing continuum (see bilinearbig)
   IF (.not.ALLOCATED(indi)) THEN
@@ -107 +94 @@
      DPR(k) = PLEV(K)-PLEV(K-1)
 
-     !--- Kasting's CIA ----------------------------------------
-     !dz(k)=dpr(k)*189.02*TMID(K)/(0.03720*PMID(K))
-     ! this is CO2 path length (in cm) as written by Francois
-     ! delta_z = delta_p * R_specific * T / (g * P)
-     ! But Kasting states that W is in units of _atmosphere_ cm
-     ! So we do
-     !dz(k)=dz(k)*(PMID(K)/1013.25)
-     !dz(k)=dz(k)/100.0 ! in m for SI calc
-     !----------------------------------------------------------
-
      ! if we have continuum opacities, we need dz
-     if(kastprof)then
-        dz(k) = dpr(k)*(1000.0d0*8.3145d0/muvar(k))*TMID(K)/(g*PMID(K))
-        U(k)  = Cmk*DPR(k)*mugaz/muvar(k) 
-     else
-        dz(k) = dpr(k)*R*TMID(K)/(glat_ig*PMID(K))*mugaz/muvar(k)
-        U(k)  = Cmk*DPR(k)*mugaz/muvar(k)     ! only Cmk line in optci.F  
+      dz(k) = dpr(k)*R*TMID(K)/(glat_ig*PMID(K))*mugaz/muvar(k)
+      U(k)  = Cmk*DPR(k)*mugaz/muvar(k)     ! only Cmk line in optci.F  
             !JL13 the mugaz/muvar factor takes into account water meanmolecular weight if water is present
-     endif
 
      call tpindex(PMID(K),TMID(K),QVAR(K),pfgasref,tgasref,WREFVAR, &
@@ -185 +157 @@
                        call interpolateN2H2(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.,interm)
                        indi(nw,igas,jgas) = interm
-                    elseif(jgas.eq.igas_He)then 
-                       interm = indi(nw,igas,jgas)
-                       call interpolateH2He(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.,interm)
-                       indi(nw,igas,jgas) = interm
                     endif
                     dtemp = dtemp + dtempc
                  enddo
-
-              elseif(igas.eq.igas_H2O.and.T_cont.gt.200.0)then
-
-                 p_air = dble(PMID(k)*scalep) - p_cont ! note assumes background is air!
-                 if(H2Ocont_simple)then
-                    call interpolateH2Ocont_PPC(wn_cont,T_cont,p_cont,p_air,dtemp,.false.)
-                 else
-                    interm = indi(nw,igas,igas)
-                    call interpolateH2Ocont_CKD(wn_cont,T_cont,p_cont,p_air,dtemp,.false.,interm)
-                    indi(nw,igas,igas) = interm
-                 endif
 
               endif

trunk/LMDZ.TITAN/libf/phytitan/optcv.F90

@@ -7 +7 @@
   use gases_h
   use comcstfi_mod, only: g, r, mugaz
-  use callkeys_mod, only: kastprof,continuum,graybody,H2Ocont_simple,callgasvis
+  use callkeys_mod, only: continuum,graybody,callgasvis
 
   implicit none
@@ -106 +106 @@
 
      ! if we have continuum opacities, we need dz
-     if(kastprof)then
-        dz(k) = dpr(k)*(1000.0d0*8.3145d0/muvar(k))*TMID(K)/(g*PMID(K))
-        U(k)  = Cmk*DPR(k)*mugaz/muvar(k) 
-     else
-        dz(k) = dpr(k)*R*TMID(K)/(glat_ig*PMID(K))*mugaz/muvar(k)
-        U(k)  = Cmk*DPR(k)*mugaz/muvar(k)     ! only Cmk line in optci.F  
-            !JL13 the mugaz/muvar factor takes into account water meanmolecular weight if water is present
-     endif
+
+      dz(k) = dpr(k)*R*TMID(K)/(glat_ig*PMID(K))*mugaz/muvar(k)
+      U(k)  = Cmk*DPR(k)*mugaz/muvar(k)     ! only Cmk line in optci.F  
+          !JL13 the mugaz/muvar factor takes into account water meanmolecular weight if water is present
+     
 
      call tpindex(PMID(K),TMID(K),QVAR(K),pfgasref,tgasref,WREFVAR, &
@@ -186 +183 @@
                        indv(nw,igas,jgas) = interm
                        ! should be irrelevant in the visible
-                    elseif(jgas.eq.igas_He)then 
-                       interm = indv(nw,igas,jgas)
-                       call interpolateH2He(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.,interm)
-                       indv(nw,igas,jgas) = interm
                     endif
                     dtemp = dtemp + dtempc
                  enddo
-
-              elseif(igas.eq.igas_H2O.and.T_cont.gt.200.0)then
-
-                 p_air = dble(PMID(k)*scalep) - p_cont ! note assumes background is air!
-                 if(H2Ocont_simple)then
-                    call interpolateH2Ocont_PPC(wn_cont,T_cont,p_cont,p_air,dtemp,.false.)
-                 else
-                    interm = indv(nw,igas,igas)
-                    call interpolateH2Ocont_CKD(wn_cont,T_cont,p_cont,p_air,dtemp,.false.,interm)
-                    indv(nw,igas,igas) = interm
-                 endif
 
               endif

trunk/LMDZ.TITAN/libf/phytitan/phyetat0.F90

@@ -1 +1 @@
 subroutine phyetat0 (ngrid,nlayer,fichnom,tab0,Lmodif,nsoil,nq, &
                      day_ini,time,tsurf,tsoil, &
-                     emis,q2,qsurf,cloudfrac,totcloudfrac,hice, &
-                     rnat,pctsrf_sic,tslab,tsea_ice,sea_ice)
+                     emis,q2,qsurf)
 
 
@@ -10 +9 @@
                      get_field, get_var, inquire_field, &
                      inquire_dimension, inquire_dimension_length
-  use slab_ice_h, only: noceanmx
-  use callkeys_mod, only: CLFvarying
 
   implicit none
@@ -45 +42 @@
   real,intent(out) :: q2(ngrid,nlayer+1) ! 
   real,intent(out) :: qsurf(ngrid,nq) ! tracers on surface
-! real co2ice(ngrid) ! co2 ice cover
-  real,intent(out) :: cloudfrac(ngrid,nlayer)
-  real,intent(out) :: hice(ngrid), totcloudfrac(ngrid)
-  real,intent(out) :: pctsrf_sic(ngrid),tslab(ngrid,noceanmx)  
-  real,intent(out) :: tsea_ice(ngrid),sea_ice(ngrid)
-  real,intent(out) :: rnat(ngrid)
 
 !======================================================================
@@ -56 +47 @@
 
 !      INTEGER radpas
-!      REAL co2_ppm
 !      REAL solaire
 
@@ -250 +240 @@
              minval(emis), maxval(emis)
 endif
-
-! Cloud fraction (added by BC 2010)
-if (CLFvarying) then
-call get_field("cloudfrac",cloudfrac,found,indextime)
-if (.not.found) then
-  write(*,*) "phyetat0: Failed loading <cloudfrac>"
-  call abort
-else
-  write(*,*) "phyetat0: Cloud fraction <cloudfrac> range:", &
-             minval(cloudfrac), maxval(cloudfrac)
-endif
-else
-cloudfrac(:,:)=0.0
-endif
-
-! Total cloud fraction (added by BC 2010)
-if (CLFvarying) then
-call get_field("totcloudfrac",totcloudfrac,found,indextime)
-if (.not.found) then
-  write(*,*) "phyetat0: Failed loading <totcloudfrac>"
-  call abort
-else
-  write(*,*) "phyetat0: Total cloud fraction <totcloudfrac> range:", &
-             minval(totcloudfrac), maxval(totcloudfrac)
-endif
-else
-totcloudfrac(:)=0.0
-endif
-
-! Height of oceanic ice (added by BC 2010)
-call get_field("hice",hice,found,indextime)
-if (.not.found) then
-  write(*,*) "phyetat0: Failed loading <hice>"
-!  call abort
-      do ig=1,ngrid
-      hice(ig)=0.
-      enddo
-else
-  write(*,*) "phyetat0: Height of oceanic ice <hice> range:", &
-             minval(hice), maxval(hice)
-endif
-
-! SLAB OCEAN (added by BC 2014)
-! nature of the surface
-call get_field("rnat",rnat,found,indextime)
-if (.not.found) then
-  write(*,*) "phyetat0: Failed loading <rnat>"
-      do ig=1,ngrid
-        rnat(ig)=1.
-      enddo
-else
-      do ig=1,ngrid
-        if((nint(rnat(ig)).eq.2).or.(nint(rnat(ig)).eq.0))then
-          rnat(ig)=0.
-        else
-          rnat(ig)=1.
-        endif      
-      enddo
-
-  write(*,*) "phyetat0: Nature of surface <rnat> range:", &
-             minval(rnat), maxval(rnat)
-endif
-! Pourcentage of sea ice cover
-call get_field("pctsrf_sic",pctsrf_sic,found,indextime)
-if (.not.found) then
-  write(*,*) "phyetat0: Failed loading <pctsrf_sic>"
-      do ig=1,ngrid
-      pctsrf_sic(ig)=0.
-      enddo
-else
-  write(*,*) "phyetat0: Pourcentage of sea ice cover <pctsrf_sic> range:", &
-             minval(pctsrf_sic), maxval(pctsrf_sic)
-endif
-! Slab ocean temperature (2 layers)
-call get_field("tslab",tslab,found,indextime)
-if (.not.found) then
-  write(*,*) "phyetat0: Failed loading <tslab>"
-      do ig=1,ngrid
-      do iq=1,noceanmx
-      tslab(ig,iq)=tsurf(ig)
-      enddo
-      enddo
-else
-  write(*,*) "phyetat0: Slab ocean temperature <tslab> range:", &
-             minval(tslab), maxval(tslab)
-endif
-! Oceanic ice temperature
-call get_field("tsea_ice",tsea_ice,found,indextime)
-if (.not.found) then
-  write(*,*) "phyetat0: Failed loading <tsea_ice>"
-      do ig=1,ngrid
-      tsea_ice(ig)=273.15-1.8
-      enddo
-else
-  write(*,*) "phyetat0: Oceanic ice temperature <tsea_ice> range:", &
-             minval(tsea_ice), maxval(tsea_ice)
-endif
-!  Oceanic ice quantity (kg/m^2)
-call get_field("sea_ice",sea_ice,found,indextime)
-if (.not.found) then
-  write(*,*) "phyetat0: Failed loading <sea_ice>"
-      do ig=1,ngrid
-      tsea_ice(ig)=0.
-      enddo
-else
-  write(*,*) "phyetat0: Oceanic ice quantity <sea_ice> range:", &
-             minval(sea_ice), maxval(sea_ice)
-endif
-
-
-
 
 ! pbl wind variance

trunk/LMDZ.TITAN/libf/phytitan/phyredem.F90

@@ -135 +135 @@
 
 subroutine physdem1(filename,nsoil,ngrid,nlay,nq, &
-                    phystep,time,tsurf,tsoil,emis,q2,qsurf, &
-                    cloudfrac,totcloudfrac,hice, &
-                    rnat,pctsrf_sic,tslab,tsea_ice,sea_ice)
+                    phystep,time,tsurf,tsoil,emis,q2,qsurf)
   ! write time-dependent variable to restart file
   use iostart, only : open_restartphy, close_restartphy, & 
                       put_var, put_field
   use tracer_h, only: noms
-  use slab_ice_h, only: noceanmx
-  use callkeys_mod, only: ok_slab_ocean
 
   implicit none
@@ -159 +155 @@
   real,intent(in) :: q2(ngrid,nlay+1)
   real,intent(in) :: qsurf(ngrid,nq)
-  real,intent(in) :: cloudfrac(ngrid,nlay)
-  real,intent(in) :: totcloudfrac(ngrid)
-  real,intent(in) :: hice(ngrid)
-  real,intent(in) :: rnat(ngrid)
-  real,intent(in) :: pctsrf_sic(ngrid)
-  real,intent(in) :: tslab(ngrid,noceanmx)
-  real,intent(in) :: tsea_ice(ngrid)
-  real,intent(in) :: sea_ice(ngrid)
 
   integer :: iq
@@ -189 +177 @@
   call put_field("q2","pbl wind variance",q2)
   
-! cloud fraction and sea ice (NB: these should be optional... to be improved)
-  call put_field("cloudfrac","Cloud fraction",cloudfrac)
-  call put_field("totcloudfrac","Total fraction",totcloudfrac)
-  call put_field("hice","Height of oceanic ice",hice)
-
- !Slab ocean 
- if(ok_slab_ocean) then
-   call put_field("rnat","Nature of surface",rnat)
-   call put_field("pctsrf_sic","Pourcentage sea ice",pctsrf_sic)
-   call put_field("tslab","Temperature slab ocean",tslab)
-   call put_field("tsea_ice","Temperature sea ice",tsea_ice)
-   call put_field("sea_ice","Sea ice mass",sea_ice)
- endif!(ok_slab_ocean)
-
-
 ! tracers
   if (nq>0) then

trunk/LMDZ.TITAN/libf/phytitan/physiq_mod.F90

@@ -15 +15 @@
  
       use radinc_h, only : L_NSPECTI,L_NSPECTV,naerkind
-      use watercommon_h, only : RLVTT, Psat_water,epsi
       use gases_h, only: gnom, gfrac
       use radcommon_h, only: sigma, glat, grav, BWNV
-      use radii_mod, only: h2o_reffrad, co2_reffrad
-      use aerosol_mod, only: iaero_co2, iaero_h2o
-      use surfdat_h, only: phisfi, zmea, zstd, zsig, zgam, zthe, &
-                           dryness, watercaptag
+      use surfdat_h, only: phisfi, zmea, zstd, zsig, zgam, zthe
       use comdiurn_h, only: coslat, sinlat, coslon, sinlon
       use comsaison_h, only: mu0, fract, dist_star, declin, right_ascen
@@ -29 +25 @@
       USE tracer_h, only: noms, mmol, radius, rho_q, qext, &
                           alpha_lift, alpha_devil, qextrhor, &
-                          igcm_h2o_ice, igcm_h2o_vap, igcm_dustbin, &
-                          igcm_co2_ice
+                          igcm_dustbin
       use time_phylmdz_mod, only: ecritphy, iphysiq, nday
       use phyredem, only: physdem0, physdem1
-      use slab_ice_h, only: capcalocean, capcalseaice,capcalsno, &
-                            noceanmx
-      use ocean_slab_mod, only :ocean_slab_init, ocean_slab_ice, &
-                                ini_surf_heat_transp_mod, &
-                                ocean_slab_get_vars,ocean_slab_final
-      use surf_heat_transp_mod,only :init_masquv
       use planetwide_mod, only: planetwide_minval,planetwide_maxval,planetwide_sumval
       use mod_phys_lmdz_para, only : is_master
@@ -81 +70 @@
 !      IV. Dry Convective adjustment :
 !
-!      V. Condensation and sublimation of gases (currently just CO2).
-!
-!      VI. Tracers
-!         VI.1. Water and water ice.
-!         VI.2. Aerosols and particles.
-!         VI.3. Updates (pressure variations, surface budget).
-!         VI.4. Slab Ocean.
-!         VI.5. Surface Tracer Update.
-!
-!      VII. Surface and sub-surface soil temperature.
-!
-!      VIII. Perform diagnostics and write output files.
+!      V. Tracers
+!         V.1. Aerosols and particles.
+!         V.2. Updates (pressure variations, surface budget).
+!         V.3. Surface Tracer Update.
+!
+!      VI. Surface and sub-surface soil temperature.
+!
+!      VII. Perform diagnostics and write output files.
 !
 !
@@ -156 +141 @@
 !           No more ngridmx/nqmx, F90 commons and adaptation to parallel: A. Spiga (2012)
 !           Purge of the code : M. Turbet (2015)
-!==================================================================
+!           Fork for Titan and clean of all too-generic (ocean, water, co2 ...) routines : J. Vatant d'Ollone (2017)
+!============================================================================================
 
 
@@ -212 +198 @@
       real, dimension(:,:),allocatable,save :: albedo              ! Surface Spectral albedo. By MT2015.
       real, dimension(:),allocatable,save :: albedo_equivalent     ! Spectral Mean albedo.
-      real, dimension(:),allocatable,save :: albedo_snow_SPECTV    ! Snow Spectral albedo.
-      real, dimension(:),allocatable,save :: albedo_co2_ice_SPECTV ! CO2 ice Spectral albedo.
-      
-!$OMP THREADPRIVATE(tsurf,tsoil,albedo,albedo_equivalent,albedo_snow_SPECTV,albedo_co2_ice_SPECTV)
+      
+!$OMP THREADPRIVATE(tsurf,tsoil,albedo,albedo_equivalent)
 
       real,dimension(:),allocatable,save :: albedo_bareground ! Bare Ground Albedo. By MT 2015.
-      real,dimension(:),allocatable,save :: rnat              ! Defines the type of the grid (ocean,continent,...). By BC.
-      
-!$OMP THREADPRIVATE(albedo_bareground,rnat)
+      
+!$OMP THREADPRIVATE(albedo_bareground)
 
       real,dimension(:),allocatable,save :: emis        ! Thermal IR surface emissivity.
@@ -276 +259 @@
       real zdtsurf(ngrid)     ! Cumulated tendencies.
       real zdtsurfmr(ngrid)   ! Mass_redistribution routine.
-      real zdtsurfc(ngrid)    ! Condense_co2 routine.
       real zdtsdif(ngrid)     ! Turbdiff/vdifc routines.
-      real zdtsurf_hyd(ngrid) ! Hydrol routine.
             
       ! For Atmospheric Temperatures : (K/s)    
-      real dtlscale(ngrid,nlayer)                             ! Largescale routine.
-      real zdtc(ngrid,nlayer)                                 ! Condense_co2 routine.
-      real zdtdif(ngrid,nlayer)                                      ! Turbdiff/vdifc routines.
+      real zdtdif(ngrid,nlayer)                               ! Turbdiff/vdifc routines.
       real zdtmr(ngrid,nlayer)                                ! Mass_redistribution routine.
-      real zdtrain(ngrid,nlayer)                              ! Rain routine.
-      real dtmoist(ngrid,nlayer)                              ! Moistadj routine.
-      real dt_ekman(ngrid,noceanmx), dt_hdiff(ngrid,noceanmx) ! Slab_ocean routine.
       real zdtsw1(ngrid,nlayer), zdtlw1(ngrid,nlayer)         ! Callcorrk routine.
                               
       ! For Surface Tracers : (kg/m2/s)
       real dqsurf(ngrid,nq)                 ! Cumulated tendencies.
-      real zdqsurfc(ngrid)                  ! Condense_co2 routine.
       real zdqsdif(ngrid,nq)                ! Turbdiff/vdifc routines.
       real zdqssed(ngrid,nq)                ! Callsedim routine.
       real zdqsurfmr(ngrid,nq)              ! Mass_redistribution routine.
-      real zdqsrain(ngrid), zdqssnow(ngrid) ! Rain routine.
-      real dqs_hyd(ngrid,nq)                ! Hydrol routine.
                   
       ! For Tracers : (kg/kg_of_air/s)
-      real zdqc(ngrid,nlayer,nq)      ! Condense_co2 routine.
       real zdqadj(ngrid,nlayer,nq)    ! Convadj routine.
       real zdqdif(ngrid,nlayer,nq)    ! Turbdiff/vdifc routines.
@@ -306 +278 @@
       real zdqsed(ngrid,nlayer,nq)    ! Callsedim routine.
       real zdqmr(ngrid,nlayer,nq)     ! Mass_redistribution routine.
-      real zdqrain(ngrid,nlayer,nq)   ! Rain routine.
-      real dqmoist(ngrid,nlayer,nq)   ! Moistadj routine.
-      real dqvaplscale(ngrid,nlayer)  ! Largescale routine.
-      real dqcldlscale(ngrid,nlayer)  ! Largescale routine.
                         
       ! For Winds : (m/s/s)
@@ -359 +327 @@
       real qcol(ngrid,nq) ! Tracer Column Mass (kg/m2).
 
-!     included by RW for H2O Manabe scheme
-      real rneb_man(ngrid,nlayer) ! H2O cloud fraction (moistadj).
-      real rneb_lsc(ngrid,nlayer) ! H2O cloud fraction (large scale).
-
-
 !     to test energy conservation (RW+JL)
       real mass(ngrid,nlayer),massarea(ngrid,nlayer)
@@ -371 +334 @@
       real dEzRadsw(ngrid,nlayer),dEzRadlw(ngrid,nlayer),dEzdiff(ngrid,nlayer)
       real dEdiffs(ngrid),dEdiff(ngrid)
-      real madjdE(ngrid), lscaledE(ngrid),madjdEz(ngrid,nlayer), lscaledEz(ngrid,nlayer)
       
 !JL12 conservation test for mean flow kinetic energy has been disabled temporarily
 
-      real dtmoist_max,dtmoist_min     
       real dItot, dItot_tmp, dVtot, dVtot_tmp
-
-      real,allocatable,save :: hice(:) ! Oceanic Ice height. by BC
- !$OMP THREADPRIVATE(hice)     
-
-      real h2otot                      ! Total Amount of water. For diagnostic.
-      real icesrf,liqsrf,icecol,vapcol ! Total Amounts of water (ice,liq,vap). For diagnostic.
+      
       real dWtot, dWtot_tmp, dWtots, dWtots_tmp
-      logical,save :: watertest
-!$OMP THREADPRIVATE(watertest)
-
-      real qsat(ngrid,nlayer) ! Water Vapor Volume Mixing Ratio at saturation (kg/kg_of_air).
-      real RH(ngrid,nlayer)   ! Relative humidity.
-      real H2Omaxcol(ngrid)   ! Maximum possible H2O column amount (at 100% saturation) (kg/m2).
-      real psat_tmp
-      
-      logical clearsky ! For double radiative transfer call. By BC
+      
       
       ! For Clear Sky Case.
@@ -400 +348 @@
       real tau_col1(ngrid)                                                   ! For aerosol optical depth diagnostic.
       real OLR_nu1(ngrid,L_NSPECTI), OSR_nu1(ngrid,L_NSPECTV)                ! For Outgoing Radiation diagnostics.
-      real tf, ntf
-
-      real,allocatable,dimension(:,:),save :: cloudfrac  ! Fraction of clouds (%).
-      real,allocatable,dimension(:),save :: totcloudfrac ! Column fraction of clouds (%).
-!$OMP THREADPRIVATE(cloudfrac,totcloudfrac)
-
-      real nconsMAX, vdifcncons(ngrid), cadjncons(ngrid) ! Vdfic water conservation test. By RW
+      real tf, ntf    
 
       real,allocatable,dimension(:,:),save :: qsurf_hist
@@ -418 +360 @@
 
       real reffcol(ngrid,naerkind)
-
-!     Sourceevol for 'accelerated ice evolution'. By RW
-      real,allocatable,dimension(:),save :: ice_initial
-      real delta_ice,ice_tot
-      real,allocatable,dimension(:),save :: ice_min      
-      integer num_run
-      logical,save :: ice_update
-!$OMP THREADPRIVATE(ice_initial,ice_min,ice_update)
-
-!     For slab ocean. By BC
-      real, dimension(:),allocatable,save ::  pctsrf_sic
-      real, dimension(:,:),allocatable,save :: tslab
-      real, dimension(:),allocatable,save ::  tsea_ice
-      real, dimension(:),allocatable,save :: sea_ice
-      real, dimension(:),allocatable,save :: zmasq
-      integer, dimension(:),allocatable,save ::knindex 
-!$OMP THREADPRIVATE(pctsrf_sic,tslab,tsea_ice,sea_ice,zmasq,knindex)
-
-      real :: tsurf2(ngrid)
-      real :: flux_o(ngrid),flux_g(ngrid),fluxgrdocean(ngrid)
-      real :: flux_sens_lat(ngrid)
-      real :: qsurfint(ngrid,nq)
 
 
@@ -457 +377 @@
         ALLOCATE(tsoil(ngrid,nsoilmx))    
         ALLOCATE(albedo(ngrid,L_NSPECTV))
-         ALLOCATE(albedo_equivalent(ngrid))        
-         ALLOCATE(albedo_snow_SPECTV(L_NSPECTV))
-         ALLOCATE(albedo_co2_ice_SPECTV(L_NSPECTV))
-         ALLOCATE(albedo_bareground(ngrid))           
-        ALLOCATE(rnat(ngrid))         
+        ALLOCATE(albedo_equivalent(ngrid))        
+        ALLOCATE(albedo_bareground(ngrid))                
         ALLOCATE(emis(ngrid))   
         ALLOCATE(dtrad(ngrid,nlayer))
@@ -472 +389 @@
         ALLOCATE(ztprevious(ngrid,nlayer))
         ALLOCATE(zuprevious(ngrid,nlayer))
-        ALLOCATE(cloudfrac(ngrid,nlayer))
-        ALLOCATE(totcloudfrac(ngrid))
-        ALLOCATE(hice(ngrid))
         ALLOCATE(qsurf_hist(ngrid,nq))
         ALLOCATE(reffrad(ngrid,nlayer,naerkind))
         ALLOCATE(nueffrad(ngrid,nlayer,naerkind))
-        ALLOCATE(ice_initial(ngrid))
-        ALLOCATE(ice_min(ngrid)) 
         ALLOCATE(fluxsurf_lw(ngrid))
         ALLOCATE(fluxsurf_sw(ngrid))
@@ -493 +405 @@
         ALLOCATE(zdtsw(ngrid,nlayer))
         ALLOCATE(tau_col(ngrid))
-        ALLOCATE(pctsrf_sic(ngrid))
-        ALLOCATE(tslab(ngrid,noceanmx))
-        ALLOCATE(tsea_ice(ngrid))
-        ALLOCATE(sea_ice(ngrid))
-        ALLOCATE(zmasq(ngrid))
-        ALLOCATE(knindex(ngrid))
 
         ! This is defined in comsaison_h
@@ -532 +438 @@
 !        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
          call phyetat0(ngrid,nlayer,"startfi.nc",0,0,nsoilmx,nq,      &
-                       day_ini,time_phys,tsurf,tsoil,emis,q2,qsurf,   &
-                       cloudfrac,totcloudfrac,hice,                   &
-                       rnat,pctsrf_sic,tslab, tsea_ice,sea_ice)
+                       day_ini,time_phys,tsurf,tsoil,emis,q2,qsurf)
 
          if (pday.ne.day_ini) then
@@ -550 +454 @@
          albedo(:,:)=0.0
           albedo_bareground(:)=0.0
-          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) 
+         call surfini(ngrid,nq,qsurf,albedo,albedo_bareground) 
          
 !        Initialize orbital calculation. 
@@ -564 +466 @@
          endif
 
-!        Initialize oceanic variables.
-!        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-         if (ok_slab_ocean)then
-
-           call ocean_slab_init(ngrid,ptimestep, tslab, &
-                                sea_ice, pctsrf_sic)
-
-           call ini_surf_heat_transp_mod()
-           
-           knindex(:) = 0
-
-           do ig=1,ngrid
-              zmasq(ig)=1
-              knindex(ig) = ig
-              if (nint(rnat(ig)).eq.0) then
-                 zmasq(ig)=0
-              endif
-           enddo
-
-           CALL init_masquv(ngrid,zmasq)
-
-         endif ! end of 'ok_slab_ocean'.
-
 
 !        Initialize soil.
@@ -595 +473 @@
             call soil(ngrid,nsoilmx,firstcall,lastcall,inertiedat, &
                       ptimestep,tsurf,tsoil,capcal,fluxgrd)
-
-            if (ok_slab_ocean) then
-               do ig=1,ngrid
-                  if (nint(rnat(ig)).eq.2) then
-                     capcal(ig)=capcalocean
-                     if (pctsrf_sic(ig).gt.0.5) then
-                        capcal(ig)=capcalseaice
-                        if (qsurf(ig,igcm_h2o_ice).gt.0.) then
-                           capcal(ig)=capcalsno
-                        endif
-                     endif
-                  endif
-               enddo
-            endif ! end of 'ok_slab_ocean'.
 
          else ! else of 'callsoil'.
@@ -620 +484 @@
          
          icount=1
-
-!        Decide whether to update ice at end of run.
-!        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-         ice_update=.false.
-         if(sourceevol)then
-!$OMP MASTER
-            open(128,file='num_run',form='formatted', &
-                     status="old",iostat=ierr)
-            if (ierr.ne.0) then
-               write(*,*) "physiq: Error! No num_run file!"
-               write(*,*) " (which is needed for sourceevol option)"
-               stop
-            endif
-            read(128,*) num_run 
-            close(128)
-!$OMP END MASTER
-!$OMP BARRIER
-        
-            if(num_run.ne.0.and.mod(num_run,2).eq.0)then
-               print*,'Updating ice at end of this year!'
-               ice_update=.true.
-               ice_min(:)=1.e4
-            endif
-            
-         endif ! end of 'sourceevol'.
-
-
-         ! Here is defined the type of the surface : Continent or Ocean.
-         ! BC2014 : This is now already done in newstart.
-         ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-         if (.not.ok_slab_ocean) then
-         
-           rnat(:)=1.
-           do ig=1,ngrid
-              if(inertiedat(ig,1).gt.1.E4)then
-                 rnat(ig)=0
-              endif
-           enddo
-
-           print*,'WARNING! Surface type currently decided by surface inertia'
-           print*,'This should be improved e.g. in newstart.F'
            
-         endif ! end of 'ok_slab_ocean'.
-
 
 !        Initialize surface history variable.
@@ -674 +495 @@
          zuprevious(:,:)=pu(:,:)
 
-!        Set temperature just above condensation temperature (for Early Mars)
-!        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-         if(nearco2cond) then
-            write(*,*)' WARNING! Starting at Tcond+1K'
-            do l=1, nlayer
-               do ig=1,ngrid
-                  pdt(ig,l)= ((-3167.8)/(log(.01*pplay(ig,l))-23.23)+4     &
-                      -pt(ig,l)) / ptimestep
-               enddo
-            enddo
-         endif
 
          if(meanOLR)then          
@@ -691 +501 @@
             call system('rm -f h2o_bal.out') ! to record global hydrological balance.
          endif
-
-
-         watertest=.false.        
-         if(water)then ! Initialize variables for water cycle
-            
-            if(enertest)then
-               watertest = .true.
-            endif
-
-            if(ice_update)then
-               ice_initial(:)=qsurf(:,igcm_h2o_ice)
-            endif
-
-         endif
-         
-         call su_watercycle ! even if we don't have a water cycle, we might
-                            ! need epsi for the wvp definitions in callcorrk.F
 
          if (ngrid.ne.1) then ! Note : no need to create a restart file in 1d.
@@ -735 +528 @@
       ! Initialize various variables
       ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-      
-      if ( .not.nearco2cond ) then
-         pdt(1:ngrid,1:nlayer) = 0.0
-      endif      
+
+      pdt(1:ngrid,1:nlayer) = 0.0     
       zdtsurf(1:ngrid)      = 0.0
       pdq(1:ngrid,1:nlayer,1:nq) = 0.0
@@ -746 +537 @@
       pdpsrf(1:ngrid)       = 0.0
       zflubid(1:ngrid)      = 0.0      
-      flux_sens_lat(1:ngrid) = 0.0
       taux(1:ngrid) = 0.0
       tauy(1:ngrid) = 0.0
@@ -870 +660 @@
 ! II.a Call correlated-k radiative transfer scheme
 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-               if(kastprof)then
-                  print*,'kastprof should not = true here'
-                  call abort
-               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)) 
-                  ! take into account water effect on mean molecular weight
-               else
-                  muvar(1:ngrid,1:nlayer+1)=mugaz
-               endif 
-      
-
-               if(ok_slab_ocean) then
-                  tsurf2(:)=tsurf(:) 
-                  do ig=1,ngrid
-                     if (nint(rnat(ig))==0) then
-                        tsurf(ig)=((1.-pctsrf_sic(ig))*tslab(ig,1)**4+pctsrf_sic(ig)*tsea_ice(ig)**4)**0.25
-                     endif
-                  enddo
-               endif !(ok_slab_ocean)
-                
+ 
+               muvar(1:ngrid,1:nlayer+1)=mugaz
+
                ! standard callcorrk
-               clearsky=.false.
                call callcorrk(ngrid,nlayer,pq,nq,qsurf,                           &
                               albedo,albedo_equivalent,emis,mu0,pplev,pplay,pt,   &
@@ -900 +670 @@
                               fluxsurfabs_sw,fluxtop_lw,                          &
                               fluxabs_sw,fluxtop_dn,OLR_nu,OSR_nu,                &
-                              tau_col,cloudfrac,totcloudfrac,                     &
-                              clearsky,firstcall,lastcall)     
-
-                ! Option to call scheme once more for clear regions  
-               if(CLFvarying)then
-
-                  ! ---> PROBLEMS WITH ALLOCATED ARRAYS : temporary solution in callcorrk: do not deallocate if CLFvarying ...
-                  clearsky=.true.
-                  call callcorrk(ngrid,nlayer,pq,nq,qsurf,                           &
-                                 albedo,albedo_equivalent1,emis,mu0,pplev,pplay,pt,   &
-                                 tsurf,fract,dist_star,aerosol,muvar,                &
-                                 zdtlw1,zdtsw1,fluxsurf_lw1,fluxsurf_sw1,            &
-                                 fluxsurfabs_sw1,fluxtop_lw1,                        &
-                                 fluxabs_sw1,fluxtop_dn,OLR_nu1,OSR_nu1,             &
-                                 tau_col1,cloudfrac,totcloudfrac,                    &
-                                 clearsky,firstcall,lastcall)
-                  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                    
-                     fluxsurf_lw(ig)       = ntf*fluxsurf_lw1(ig)       + tf*fluxsurf_lw(ig)
-                     fluxsurf_sw(ig)       = ntf*fluxsurf_sw1(ig)       + tf*fluxsurf_sw(ig)
-                     albedo_equivalent(ig) = ntf*albedo_equivalent1(ig) + tf*albedo_equivalent(ig)
-                     fluxsurfabs_sw(ig)    = ntf*fluxsurfabs_sw1(ig)    + tf*fluxsurfabs_sw(ig)
-                     fluxtop_lw(ig)        = ntf*fluxtop_lw1(ig)        + tf*fluxtop_lw(ig)
-                     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)                       
-                     OLR_nu(ig,1:L_NSPECTI) = ntf*OLR_nu1(ig,1:L_NSPECTI) + tf*OLR_nu(ig,1:L_NSPECTI)                       
-                  enddo                               
-
-               endif ! end of CLFvarying.
-
-               if(ok_slab_ocean) then
-                  tsurf(:)=tsurf2(:)
-               endif
-
+                              tau_col,firstcall,lastcall)
 
                ! Radiative flux from the sky absorbed by the surface (W.m-2).
@@ -957 +685 @@
                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 ! 
-               if (firstcall .and. albedo_spectral_mode) then
-                  call spectral_albedo_calc(albedo_snow_SPECTV)
-               endif
-
             elseif(newtonian)then
             
@@ -1040 +763 @@
          if (UseTurbDiff) then
          
-            call turbdiff(ngrid,nlayer,nq,rnat,                  &
+            call turbdiff(ngrid,nlayer,nq,                       &
                           ptimestep,capcal,lwrite,               &
                           pplay,pplev,zzlay,zzlev,z0,            &
@@ -1046 +769 @@
                           pdt,pdq,zflubid,                       &
                           zdudif,zdvdif,zdtdif,zdtsdif,          &
-                          sensibFlux,q2,zdqdif,zdqevap,zdqsdif,  &
+                          sensibFlux,q2,zdqdif,zdqsdif,          &
                           taux,tauy,lastcall)
 
@@ -1053 +776 @@
             zdh(1:ngrid,1:nlayer)=pdt(1:ngrid,1:nlayer)/zpopsk(1:ngrid,1:nlayer)
  
-            call vdifc(ngrid,nlayer,nq,rnat,zpopsk,           &
+            call vdifc(ngrid,nlayer,nq,zpopsk,                &
                        ptimestep,capcal,lwrite,               &
                        pplay,pplev,zzlay,zzlev,z0,            &
@@ -1072 +795 @@
          pdt(1:ngrid,1:nlayer)=pdt(1:ngrid,1:nlayer)+zdtdif(1:ngrid,1:nlayer)
          zdtsurf(1:ngrid)=zdtsurf(1:ngrid)+zdtsdif(1:ngrid)
-
-
-         if(ok_slab_ocean)then
-            flux_sens_lat(1:ngrid)=(zdtsdif(1:ngrid)*capcal(1:ngrid)-fluxrad(1:ngrid))
-         endif
-
 
          if (tracer) then 
@@ -1116 +833 @@
          endif ! end of 'enertest'
 
-
-         ! 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)
-            do ig = 1, ngrid
-               vdifcncons(ig)=SUM(mass(ig,:)*zdqdif(ig,:,igcm_h2o_vap))
-            enddo
-            call planetwide_sumval(massarea(:,:)*zdqdif(:,:,igcm_h2o_ice)*ptimestep/totarea_planet,dWtot)
-            call planetwide_sumval(zdqsdif(:,igcm_h2o_ice)*cell_area(:)*ptimestep/totarea_planet,dWtots)
-            dWtot = dWtot + dWtot_tmp
-            dWtots = dWtots + dWtots_tmp
-            do ig = 1, ngrid
-               vdifcncons(ig)=vdifcncons(ig) + SUM(mass(ig,:)*zdqdif(ig,:,igcm_h2o_ice))
-            enddo            
-            call planetwide_maxval(vdifcncons(:),nconsMAX)
-
-            if (is_master) then
-               print*,'---------------------------------------------------------------'
-               print*,'In difv atmospheric water change        =',dWtot,' kg m-2'
-               print*,'In difv surface water change            =',dWtots,' kg m-2'
-               print*,'In difv non-cons factor                 =',dWtot+dWtots,' kg m-2'
-               print*,'In difv MAX non-cons factor             =',nconsMAX,' kg m-2 s-1'
-            endif
-
-         endif ! end of 'watertest'
-         !-------------------------
-
       else ! calldifv
 
@@ -1190 +878 @@
          endif
 
-         ! Test water conservation
-         if(watertest)then
-            call planetwide_sumval(massarea(:,:)*zdqadj(:,:,igcm_h2o_vap)*ptimestep/totarea_planet,dWtot_tmp)
-            do ig = 1, ngrid
-               cadjncons(ig)=SUM(mass(ig,:)*zdqadj(ig,:,igcm_h2o_vap))
-            enddo
-            call planetwide_sumval(massarea(:,:)*zdqadj(:,:,igcm_h2o_ice)*ptimestep/totarea_planet,dWtot)
-            dWtot = dWtot + dWtot_tmp
-            do ig = 1, ngrid
-               cadjncons(ig)=cadjncons(ig) + SUM(mass(ig,:)*zdqadj(ig,:,igcm_h2o_ice))
-            enddo            
-            call planetwide_maxval(cadjncons(:),nconsMAX)
-
-            if (is_master) then
-               print*,'In convadj atmospheric water change     =',dWtot,' kg m-2'
-               print*,'In convadj MAX non-cons factor          =',nconsMAX,' kg m-2 s-1'
-            endif
-            
-         endif ! end of 'watertest'
          
       endif ! end of 'calladj'
-
-!-----------------------------------------------
-!   V. Carbon dioxide condensation-sublimation :
-!-----------------------------------------------
-
-      if (co2cond) then
-      
-         if (.not.tracer) then
-            print*,'We need a CO2 ice tracer to condense CO2'
-            call abort
-         endif
-         call condense_co2(ngrid,nlayer,nq,ptimestep,                    &
-                           capcal,pplay,pplev,tsurf,pt,                  &
-                           pdt,zdtsurf,pq,pdq,                           &
-                           qsurf,zdqsurfc,albedo,emis,                   &
-                           albedo_bareground,albedo_co2_ice_SPECTV,      &
-                           zdtc,zdtsurfc,pdpsrf,zdqc)
-
-         pdt(1:ngrid,1:nlayer) = pdt(1:ngrid,1:nlayer)+zdtc(1:ngrid,1:nlayer)
-         zdtsurf(1:ngrid)      = zdtsurf(1:ngrid) + zdtsurfc(1:ngrid)
-
-         pdq(1:ngrid,1:nlayer,1:nq)   = pdq(1:ngrid,1:nlayer,1:nq)+ zdqc(1:ngrid,1:nlayer,1:nq)
-         dqsurf(1:ngrid,igcm_co2_ice) = dqsurf(1:ngrid,igcm_co2_ice) + zdqsurfc(1:ngrid)
-
-         ! test energy conservation
-         if(enertest)then
-            call planetwide_sumval(cpp*massarea(:,:)*zdtc(:,:)/totarea_planet,dEtot)
-            call planetwide_sumval(capcal(:)*zdtsurfc(:)*cell_area(:)/totarea_planet,dEtots)
-            if (is_master) then
-               print*,'In co2cloud atmospheric energy change   =',dEtot,' W m-2'
-               print*,'In co2cloud surface energy change       =',dEtots,' W m-2'
-            endif
-         endif
-
-      endif  ! end of 'co2cond'
-
+      
 
 !---------------------------------------------
-!   VI. Specific parameterizations for tracers 
+!   V. 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
-
-               dqmoist(1:ngrid,1:nlayer,1:nq)=0.
-               dtmoist(1:ngrid,1:nlayer)=0.
-               
-               ! Moist Convective Adjustment.
-               ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-               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)
-               pdq(1:ngrid,1:nlayer,igcm_h2o_ice) = pdq(1:ngrid,1:nlayer,igcm_h2o_ice)     &
-                                                  + 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
-               
-                  call planetwide_sumval(cpp*massarea(:,:)*dtmoist(:,:)/totarea_planet,dEtot)
-                  call planetwide_maxval(dtmoist(:,:),dtmoist_max)
-                  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'
-                     print*,'In moistadj MAX atmospheric energy change   =',dtmoist_max*ptimestep,'K/step'
-                     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'
-               
-
-               ! Large scale condensation/evaporation.
-               ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-               call largescale(ngrid,nlayer,nq,ptimestep,pplev,pplay,pt,pq,pdt,pdq,dtlscale,dqvaplscale,dqcldlscale,rneb_lsc)
-
-               pdt(1:ngrid,1:nlayer) = pdt(1:ngrid,1:nlayer)+dtlscale(1:ngrid,1:nlayer)
-               pdq(1:ngrid,1:nlayer,igcm_h2o_vap) = pdq(1:ngrid,1:nlayer,igcm_h2o_vap)+dqvaplscale(1:ngrid,1:nlayer)
-               pdq(1:ngrid,1:nlayer,igcm_h2o_ice) = pdq(1:ngrid,1:nlayer,igcm_h2o_ice)+dqcldlscale(1:ngrid,1:nlayer)
-
-               ! Test energy conservation.
-               if(enertest)then
-                  lscaledEz(:,:) = cpp*mass(:,:)*dtlscale(:,:)
-                  do ig=1,ngrid
-                     lscaledE(ig) = cpp*SUM(mass(:,:)*dtlscale(:,:))
-                  enddo
-                  call planetwide_sumval(cpp*massarea(:,:)*dtlscale(:,:)/totarea_planet,dEtot)
-
-                  if (is_master) print*,'In largescale atmospheric energy change =',dEtot,' W m-2'
-
-                  ! Test water conservation.
-                  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'
-
-               ! Compute cloud fraction.
-               do l = 1, nlayer
-                  do ig=1,ngrid
-                     cloudfrac(ig,l)=MAX(rneb_lsc(ig,l),rneb_man(ig,l)) 
-                  enddo
-               enddo
-
-            endif ! end of 'watercond'
-           
-            ! Water ice / liquid precipitation.
-            ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-            if(waterrain)then
-
-               zdqrain(1:ngrid,1:nlayer,1:nq) = 0.0
-               zdqsrain(1:ngrid)    = 0.0
-               zdqssnow(1:ngrid)    = 0.0
-
-               call rain(ngrid,nlayer,nq,ptimestep,pplev,pplay,pt,pdt,pq,pdq,            &
-                         zdtrain,zdqrain,zdqsrain,zdqssnow,cloudfrac)
-
-               pdq(1:ngrid,1:nlayer,igcm_h2o_vap) = pdq(1:ngrid,1:nlayer,igcm_h2o_vap) &
-                                                  + zdqrain(1:ngrid,1:nlayer,igcm_h2o_vap)
-               pdq(1:ngrid,1:nlayer,igcm_h2o_ice) = pdq(1:ngrid,1:nlayer,igcm_h2o_ice) &
-                                                  + 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) 
-
-               ! 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'
-                  call planetwide_sumval(massarea(:,:)*zdqrain(:,:,igcm_h2o_ice)/totarea_planet*RLVTT/cpp,dItot_tmp)
-                  call planetwide_sumval(cell_area(:)*zdqssnow(:)/totarea_planet*RLVTT/cpp,dItot)
-                  dItot = dItot + dItot_tmp
-                  call planetwide_sumval(massarea(:,:)*zdqrain(:,:,igcm_h2o_vap)*ptimestep/totarea_planet,dVtot_tmp)
-                  call planetwide_sumval(cell_area(:)*zdqsrain(:)/totarea_planet*RLVTT/cpp,dVtot)
-                  dVtot = dVtot + dVtot_tmp
-                  dEtot = dItot + dVtot
-                  
-                  if (is_master) then
-                     print*,'In rain dItot =',dItot,' W m-2'
-                     print*,'In rain dVtot =',dVtot,' W m-2'
-                     print*,'In rain atmospheric L energy change       =',dEtot,' W m-2'
-                  endif
-
-                  ! Test water conservation
-                  call planetwide_sumval(massarea(:,:)*zdqrain(:,:,igcm_h2o_vap)*ptimestep/totarea_planet+        &
-                          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'
-                          print*,'In rain surface water change            =',dWtots,' kg m-2'
-                          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'
-
+      
   ! -------------------------
-  !   VI.2. Aerosol particles
+  !   V.1. Aerosol particles
   ! -------------------------
 
@@ -1398 +898 @@
             zdqsed(1:ngrid,1:nlayer,1:nq) = 0.0
             zdqssed(1:ngrid,1:nq)  = 0.0
-
-            if(watertest)then
-           
-               iq=igcm_h2o_ice
-               call planetwide_sumval(massarea(:,:)*pq(:,:,iq)*ptimestep/totarea_planet,dWtot)
-               call planetwide_sumval(massarea(:,:)*pdq(:,:,iq)*ptimestep/totarea_planet,dWtots)
-               if (is_master) then
-                        print*,'Before sedim pq  =',dWtot,' kg m-2'
-                  print*,'Before sedim pdq =',dWtots,' kg m-2'
-               endif
-            endif
            
             call callsedim(ngrid,nlayer,ptimestep,           &
@@ -1414 +903 @@
                           zdqsed,zdqssed,nq)
 
-            if(watertest)then
-               iq=igcm_h2o_ice
-               call planetwide_sumval(massarea(:,:)*pq(:,:,iq)*ptimestep/totarea_planet,dWtot)
-               call planetwide_sumval(massarea(:,:)*pdq(:,:,iq)*ptimestep/totarea_planet,dWtots)
-               if (is_master) then
-                        print*,'After sedim pq  =',dWtot,' kg m-2'
-                        print*,'After sedim pdq =',dWtots,' kg m-2'
-                 endif
-            endif
-
             ! Whether it falls as rain or snow depends only on the surface temperature
             pdq(1:ngrid,1:nlayer,1:nq) = pdq(1:ngrid,1:nlayer,1:nq) + zdqsed(1:ngrid,1:nlayer,1:nq)
             dqsurf(1:ngrid,1:nq) = dqsurf(1:ngrid,1:nq) + zdqssed(1:ngrid,1:nq)
 
-            ! Test water conservation
-            if(watertest)then
-               call planetwide_sumval(massarea(:,:)*(zdqsed(:,:,igcm_h2o_vap)+zdqsed(:,:,igcm_h2o_ice))*ptimestep/totarea_planet,dWtot)
-               call planetwide_sumval((zdqssed(:,igcm_h2o_vap)+zdqssed(:,igcm_h2o_ice))*cell_area(:)*ptimestep/totarea_planet,dWtots)
-               if (is_master) then
-                        print*,'In sedim atmospheric ice change         =',dWtot,' kg m-2'
-                        print*,'In sedim surface ice change             =',dWtots,' kg m-2'
-                        print*,'In sedim non-cons factor                =',dWtot+dWtots,' kg m-2'
-               endif
-            endif
-
          end if ! end of 'sedimentation'
 
 
   ! ---------------
-  !   VI.3. Updates
+  !   V.2. Updates
   ! ---------------
 
@@ -1449 +917 @@
          if(mass_redistrib) then
 
-            zdmassmr(1:ngrid,1:nlayer) = mass(1:ngrid,1:nlayer) *     &
-               (   zdqevap(1:ngrid,1:nlayer)                          &
-                 + zdqrain(1:ngrid,1:nlayer,igcm_h2o_vap)             &
-                 + dqmoist(1:ngrid,1:nlayer,igcm_h2o_vap)             &
-                 + dqvaplscale(1:ngrid,1:nlayer) )
+            zdmassmr(1:ngrid,1:nlayer) = mass(1:ngrid,1:nlayer) * zdqevap(1:ngrid,1:nlayer)
 
             do ig = 1, ngrid
@@ -1464 +928 @@
 
             call mass_redistribution(ngrid,nlayer,nq,ptimestep,                     &
-                                     rnat,capcal,pplay,pplev,pt,tsurf,pq,qsurf,     &
+                                     capcal,pplay,pplev,pt,tsurf,pq,qsurf,          &
                                      pu,pv,pdt,zdtsurf,pdq,pdu,pdv,zdmassmr,        &
                                      zdtmr,zdtsurfmr,zdpsrfmr,zdumr,zdvmr,zdqmr,zdqsurfmr)
@@ -1478 +942 @@
          endif
 
-  ! ------------------
-  !   VI.4. Slab Ocean
-  ! ------------------
-
-         if (ok_slab_ocean)then
-
-            do ig=1,ngrid
-               qsurfint(:,igcm_h2o_ice)=qsurf(:,igcm_h2o_ice)
-            enddo
-
-            call ocean_slab_ice(ptimestep,                          &
-                                ngrid, knindex, tsea_ice, fluxrad,  &
-                                zdqssnow, qsurf(:,igcm_h2o_ice),    &
-                              - zdqsdif(:,igcm_h2o_vap),            &
-                                flux_sens_lat,tsea_ice, pctsrf_sic, &
-                                taux,tauy,icount)
-
-
-            call ocean_slab_get_vars(ngrid,tslab,      &
-                                     sea_ice, flux_o,  &
-                                     flux_g, dt_hdiff, &
-                                     dt_ekman)
-   
-            do ig=1,ngrid
-               if (nint(rnat(ig)).eq.1)then
-                  tslab(ig,1) = 0.
-                  tslab(ig,2) = 0.
-                  tsea_ice(ig) = 0.
-                  sea_ice(ig) = 0.
-                  pctsrf_sic(ig) = 0.
-                  qsurf(ig,igcm_h2o_ice) = qsurfint(ig,igcm_h2o_ice)
-               endif
-            enddo
-
-         endif ! end of 'ok_slab_ocean'
-
   ! -----------------------------
-  !   VI.5. Surface Tracer Update
+  !   V.3. Surface Tracer Update
   ! -----------------------------
 
-         if(hydrology)then
-            
-            call hydrol(ngrid,nq,ptimestep,rnat,tsurf,qsurf,dqsurf,dqs_hyd, &
-                        capcal,albedo,albedo_bareground,                    &
-                        albedo_snow_SPECTV,albedo_co2_ice_SPECTV,           &
-                        mu0,zdtsurf,zdtsurf_hyd,hice,pctsrf_sic,            &
-                        sea_ice)
-
-            zdtsurf(1:ngrid)     = zdtsurf(1:ngrid) + zdtsurf_hyd(1:ngrid)
-            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)
-
-            ! Test energy conservation
-            if(enertest)then
-               call planetwide_sumval(cell_area(:)*capcal(:)*zdtsurf_hyd(:)/totarea_planet,dEtots)
-               if (is_master) print*,'In hydrol surface energy change     =',dEtots,' W m-2'
-            endif
-
-            ! test water conservation
-            if(watertest)then
-               call planetwide_sumval(dqs_hyd(:,igcm_h2o_ice)*cell_area(:)*ptimestep/totarea_planet,dWtots)
-               if (is_master) print*,'In hydrol surface ice change            =',dWtots,' kg m-2'
-                  call planetwide_sumval(dqs_hyd(:,igcm_h2o_vap)*cell_area(:)*ptimestep/totarea_planet,dWtots)
-               if (is_master) then
-                  print*,'In hydrol surface water change          =',dWtots,' kg m-2'
-                  print*,'---------------------------------------------------------------'
-               endif
-            endif
-
-         else ! of if (hydrology)
-
-            qsurf(1:ngrid,1:nq) = qsurf(1:ngrid,1:nq) + ptimestep*dqsurf(1:ngrid,1:nq)
-
-         end if ! of if (hydrology)
+        qsurf(1:ngrid,1:nq) = qsurf(1:ngrid,1:nq) + ptimestep*dqsurf(1:ngrid,1:nq)
 
          ! Add qsurf to qsurf_hist, which is what we save in diagfi.nc. At the same time, we set the water 
@@ -1558 +952 @@
          qsurf_hist(:,:) = qsurf(:,:)
 
-         if(ice_update)then
-            ice_min(1:ngrid)=min(ice_min(1:ngrid),qsurf(1:ngrid,igcm_h2o_ice))
-         endif
-
       endif! end of if 'tracer'
 
 
 !------------------------------------------------
-!   VII. Surface and sub-surface soil temperature 
+!   VI. Surface and sub-surface soil temperature 
 !------------------------------------------------
 
 
       ! Increment surface temperature
-      if(ok_slab_ocean)then
-         do ig=1,ngrid
-           if (nint(rnat(ig)).eq.0)then
-             zdtsurf(ig)= (tslab(ig,1)              &
-             + pctsrf_sic(ig)*(tsea_ice(ig)-tslab(ig,1))-tsurf(ig))/ptimestep
-           endif
-           tsurf(ig)=tsurf(ig)+ptimestep*zdtsurf(ig)
-         enddo
-   
-      else
-        tsurf(1:ngrid)=tsurf(1:ngrid)+ptimestep*zdtsurf(1:ngrid) 
-      endif ! end of 'ok_slab_ocean'
-
+
+      tsurf(1:ngrid)=tsurf(1:ngrid)+ptimestep*zdtsurf(1:ngrid) 
 
       ! Compute soil temperatures and subsurface heat flux.
@@ -1592 +971 @@
 
 
-      if (ok_slab_ocean) then
-      
-         do ig=1,ngrid
-         
-            fluxgrdocean(ig)=fluxgrd(ig)
-            if (nint(rnat(ig)).eq.0) then
-               capcal(ig)=capcalocean
-               fluxgrd(ig)=0.
-               fluxgrdocean(ig)=pctsrf_sic(ig)*flux_g(ig)+(1-pctsrf_sic(ig))*(dt_hdiff(ig,1)+dt_ekman(ig,1))
-               do iq=1,nsoilmx
-                  tsoil(ig,iq)=tsurf(ig)
-               enddo
-               if (pctsrf_sic(ig).gt.0.5) then
-                  capcal(ig)=capcalseaice
-                  if (qsurf(ig,igcm_h2o_ice).gt.0.) then
-                     capcal(ig)=capcalsno
-                  endif
-               endif               
-            endif
-            
-         enddo
-         
-      endif !end of 'ok_slab_ocean'
-
-
       ! Test energy conservation
       if(enertest)then
@@ -1625 +979 @@
 
 !---------------------------------------------------
-!   VIII. Perform diagnostics and write output files
+!   VII. Perform diagnostics and write output files
 !---------------------------------------------------
 
@@ -1757 +1111 @@
          ! Generalised for arbitrary aerosols now. By LK
          reffcol(1:ngrid,1:naerkind)=0.0
-         if(co2cond.and.(iaero_co2.ne.0))then
-            call co2_reffrad(ngrid,nlayer,nq,zq,reffrad(1,1,iaero_co2))
-            do ig=1,ngrid
-               reffcol(ig,iaero_co2) = SUM(zq(ig,1:nlayer,igcm_co2_ice)*reffrad(ig,1:nlayer,iaero_co2)*mass(ig,1:nlayer))
-            enddo
-         endif
-         if(water.and.(iaero_h2o.ne.0))then
-            call h2o_reffrad(ngrid,nlayer,zq(1,1,igcm_h2o_ice),zt, &
-                             reffrad(1,1,iaero_h2o),nueffrad(1,1,iaero_h2o))
-            do ig=1,ngrid
-               reffcol(ig,iaero_h2o) = SUM(zq(ig,1:nlayer,igcm_h2o_ice)*reffrad(ig,1:nlayer,iaero_h2o)*mass(ig,1:nlayer))
-            enddo
-         endif
 
       endif ! end of 'tracer'
-
-
-      ! Test for water conservation.
-      if(water)then
-
-         call planetwide_sumval(cell_area(:)*qsurf_hist(:,igcm_h2o_ice)/totarea_planet,icesrf)
-         call planetwide_sumval(cell_area(:)*qsurf_hist(:,igcm_h2o_vap)/totarea_planet,liqsrf)
-         call planetwide_sumval(cell_area(:)*qcol(:,igcm_h2o_ice)/totarea_planet,icecol)
-         call planetwide_sumval(cell_area(:)*qcol(:,igcm_h2o_vap)/totarea_planet,vapcol)
-
-         h2otot = icesrf + liqsrf + icecol + vapcol
-         
-         if (is_master) then
-            print*,' Total water amount [kg m^-2]: ',h2otot
-            print*,' Surface ice    Surface liq.   Atmos. con.     Atmos. vap. [kg m^-2] '
-            print*, icesrf,liqsrf,icecol,vapcol
-         endif
-
-         if(meanOLR .and. is_master)then
-            if((ngrid.gt.1) .or. (mod(icount-1,ecritphy).eq.0))then
-               ! to record global water balance
-               open(98,file="h2o_bal.out",form='formatted',position='append')
-               write(98,*) zday,icesrf,liqsrf,icecol,vapcol
-               close(98)
-            endif
-         endif
-
-      endif
-
-
-      ! Calculate RH (Relative Humidity) for diagnostic.
-      if(water)then
-      
-         do l = 1, nlayer
-            do ig=1,ngrid
-               call Psat_water(zt(ig,l),pplay(ig,l),psat_tmp,qsat(ig,l))
-               RH(ig,l) = zq(ig,l,igcm_h2o_vap) / qsat(ig,l)
-            enddo
-         enddo
-
-         ! Compute maximum possible H2O column amount (100% saturation).
-         do ig=1,ngrid
-            H2Omaxcol(ig) = SUM( qsat(ig,:) * mass(ig,:))
-         enddo
-
-      endif ! end of 'water'
 
 
@@ -1838 +1133 @@
          ztime_fin = ptime + ptimestep/(float(iphysiq)*daysec)
 
-         ! Update surface ice distribution to iterate to steady state if requested
-         if(ice_update)then
-
-            do ig=1,ngrid
-
-               delta_ice = (qsurf(ig,igcm_h2o_ice)-ice_initial(ig))
-               
-               ! add multiple years of evolution
-               qsurf_hist(ig,igcm_h2o_ice) = qsurf_hist(ig,igcm_h2o_ice) + delta_ice*icetstep 
-
-               ! if ice has gone -ve, set to zero
-               if(qsurf_hist(ig,igcm_h2o_ice).lt.0.0)then
-                  qsurf_hist(ig,igcm_h2o_ice) = 0.0 
-               endif
-
-               ! if ice is seasonal, set to zero (NEW)
-               if(ice_min(ig).lt.0.01)then
-                  qsurf_hist(ig,igcm_h2o_ice) = 0.0 
-               endif
-
-            enddo
-
-            ! enforce ice conservation
-            ice_tot= SUM(qsurf_hist(:,igcm_h2o_ice)*cell_area(:) )/SUM(cell_area(:))
-            qsurf_hist(:,igcm_h2o_ice) = qsurf_hist(:,igcm_h2o_ice)*(icesrf/ice_tot)
-
-         endif
-
          if (ngrid.ne.1) then
             write(*,*)'PHYSIQ: for physdem ztime_fin =',ztime_fin
@@ -1871 +1138 @@
             call physdem1("restartfi.nc",nsoilmx,ngrid,nlayer,nq, &
                           ptimestep,ztime_fin,                    &
-                          tsurf,tsoil,emis,q2,qsurf_hist,         &
-                          cloudfrac,totcloudfrac,hice,            &
-                          rnat,pctsrf_sic,tslab,tsea_ice,sea_ice)
-         endif
-
-         if(ok_slab_ocean) then
-            call ocean_slab_final!(tslab, seaice)
-         end if
-
+                          tsurf,tsoil,emis,q2,qsurf_hist)
+         endif
          
       endif ! end of 'lastcall'
@@ -1937 +1197 @@
 
             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",       & 
-                           3,vmr)
-            endif
 
          endif ! end of 'tracer' 
-
-         if(watercond.and.CLFvarying)then
-            call wstats(ngrid,"rneb_man","H2O cloud fraction (conv)"," ",3,rneb_man)
-            call wstats(ngrid,"rneb_lsc","H2O cloud fraction (large scale)"," ",3,rneb_lsc)
-            call wstats(ngrid,"CLF","H2O cloud fraction"," ",3,cloudfrac)
-            call wstats(ngrid,"CLFt","H2O column cloud fraction"," ",2,totcloudfrac)
-            call wstats(ngrid,"RH","relative humidity"," ",3,RH)
-         endif
-
-         if (ok_slab_ocean) then
-            call wstats(ngrid,"dt_hdiff1","dt_hdiff1","K/s",2,dt_hdiff(:,1))
-            call wstats(ngrid,"dt_hdiff2","dt_hdiff2","K/s",2,dt_hdiff(:,2))
-            call wstats(ngrid,"dt_ekman1","dt_ekman1","K/s",2,dt_ekman(:,1))
-            call wstats(ngrid,"dt_ekman2","dt_ekman2","K/s",2,dt_ekman(:,2))
-            call wstats(ngrid,"tslab1","tslab1","K",2,tslab(:,1))
-            call wstats(ngrid,"tslab2","tslab2","K",2,tslab(:,2))
-            call wstats(ngrid,"pctsrf_sic","pct ice/sea","",2,pctsrf_sic)
-            call wstats(ngrid,"tsea_ice","tsea_ice","K",2,tsea_ice)
-            call wstats(ngrid,"sea_ice","sea ice","kg/m2",2,sea_ice)
-            call wstats(ngrid,"rnat","nature of the surface","",2,rnat)
-         endif
 
          if(lastcall) then
@@ -2003 +1235 @@
 !        call writediagsoil(ngrid,"temp","temperature","K",3,tsoil)
 
-      ! Oceanic layers
-      if(ok_slab_ocean) then
-         call writediagfi(ngrid,"pctsrf_sic","pct ice/sea","",2,pctsrf_sic)
-         call writediagfi(ngrid,"tsea_ice","tsea_ice","K",2,tsea_ice)
-         call writediagfi(ngrid,"sea_ice","sea ice","kg/m2",2,sea_ice)
-         call writediagfi(ngrid,"tslab1","tslab1","K",2,tslab(:,1))
-         call writediagfi(ngrid,"tslab2","tslab2","K",2,tslab(:,2))
-         call writediagfi(ngrid,"dt_hdiff1","dt_hdiff1","K/s",2,dt_hdiff(:,1))
-         call writediagfi(ngrid,"dt_hdiff2","dt_hdiff2","K/s",2,dt_hdiff(:,2))
-         call writediagfi(ngrid,"dt_ekman1","dt_ekman1","K/s",2,dt_ekman(:,1))
-         call writediagfi(ngrid,"dt_ekman2","dt_ekman2","K/s",2,dt_ekman(:,2))
-         call writediagfi(ngrid,"rnat","nature of the surface","",2,rnat)
-         call writediagfi(ngrid,"sensibFlux","sensible heat flux","w.m^-2",2,sensibFlux)
-         call writediagfi(ngrid,"latentFlux","latent heat flux","w.m^-2",2,zdqsdif(:,igcm_h2o_vap)*RLVTT)
-      endif
       
 
@@ -2037 +1254 @@
 !           call writediagfi(ngrid,"fluxsurflwcs","lw back radiation (cs).","W m-2",2,fluxsurf_lw1)
 
-         if(ok_slab_ocean) then
-            call writediagfi(ngrid,"GND","heat flux from ground","W m-2",2,fluxgrdocean)
-         else
-            call writediagfi(ngrid,"GND","heat flux from ground","W m-2",2,fluxgrd)
-         endif
+
+         call writediagfi(ngrid,"GND","heat flux from ground","W m-2",2,fluxgrd)
          
          call writediagfi(ngrid,"DYN","dynamical heat input","W m-2",2,fluxdyn)
@@ -2065 +1279 @@
          endif
           
-         if(watercond) then
-
-            call writediagfi(ngrid,"lscaledE","heat from largescale","W m-2",2,lscaledE) 
-            call writediagfi(ngrid,"madjdE","heat from moistadj","W m-2",2,madjdE)
-            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,"h2o_max_col","maximum H2O column amount","kg.m^-2",2,H2Omaxcol)
-
-         endif
-          
       endif ! end of 'enertest'
 
@@ -2087 +1289 @@
         
         ! For Debugging.
-        !call writediagfi(ngrid,'rnat','Terrain type',' ',2,real(rnat))
         !call writediagfi(ngrid,'pphi','Geopotential',' ',3,pphi)
         
-
-      ! Output aerosols.
-      if (igcm_co2_ice.ne.0.and.iaero_co2.ne.0) &
-         call writediagfi(ngrid,'CO2ice_reff','CO2ice_reff','m',3,reffrad(1,1,iaero_co2))
-      if (igcm_h2o_ice.ne.0.and.iaero_h2o.ne.0) &
-         call writediagfi(ngrid,'H2Oice_reff','H2Oice_reff','m',3,reffrad(:,:,iaero_h2o))
-      if (igcm_co2_ice.ne.0.and.iaero_co2.ne.0) &
-         call writediagfi(ngrid,'CO2ice_reffcol','CO2ice_reffcol','um kg m^-2',2,reffcol(1,iaero_co2))
-      if (igcm_h2o_ice.ne.0.and.iaero_h2o.ne.0) &
-         call writediagfi(ngrid,'H2Oice_reffcol','H2Oice_reffcol','um kg m^-2',2,reffcol(1,iaero_h2o))
 
       ! Output tracers.
@@ -2114 +1305 @@
 !                          'kg m^-2',2,qsurf(1,iq) )                          
 
-            if(watercond.or.CLFvarying)then
-               call writediagfi(ngrid,"rneb_man","H2O cloud fraction (conv)"," ",3,rneb_man)
-               call writediagfi(ngrid,"rneb_lsc","H2O cloud fraction (large scale)"," ",3,rneb_lsc)
-               call writediagfi(ngrid,"CLF","H2O cloud fraction"," ",3,cloudfrac)
-               call writediagfi(ngrid,"CLFt","H2O column cloud fraction"," ",2,totcloudfrac)
-               call writediagfi(ngrid,"RH","relative humidity"," ",3,RH)
-            endif
-
-            if(waterrain)then
-               call writediagfi(ngrid,"rain","rainfall","kg m-2 s-1",2,zdqsrain)
-               call writediagfi(ngrid,"snow","snowfall","kg m-2 s-1",2,zdqssnow)
-            endif
-
-            if((hydrology).and.(.not.ok_slab_ocean))then
-               call writediagfi(ngrid,"hice","oceanic ice height","m",2,hice)
-            endif
-
-            if(ice_update)then
-               call writediagfi(ngrid,"ice_min","min annual ice","m",2,ice_min)
-               call writediagfi(ngrid,"ice_ini","initial annual ice","m",2,ice_initial)
-            endif
-
             do ig=1,ngrid
                if(tau_col(ig).gt.1.e3)then

trunk/LMDZ.TITAN/libf/phytitan/radii_mod.F90

@@ -3 +3 @@
 !==================================================================
 !  module to centralize the radii calculations for aerosols
-! OK for water but should be extended to other aerosols (CO2,...)
 !==================================================================
       
-!     water cloud optical properties
-
-      use callkeys_mod, only: radfixed,Nmix_co2,                    &
-                pres_bottom_tropo,pres_top_tropo,size_tropo,        &
-                pres_bottom_strato,size_strato
-      
-      real, save ::  rad_h2o
-      real, save ::  rad_h2o_ice
-      real, save ::  Nmix_h2o
-      real, save ::  Nmix_h2o_ice
-!$OMP THREADPRIVATE(rad_h2o,rad_h2o_ice,Nmix_h2o,Nmix_h2o_ice)
-      real, parameter ::  coef_chaud=0.13
-      real, parameter ::  coef_froid=0.09
-
+      use callkeys_mod, only: pres_bottom_tropo,pres_top_tropo, &
+                size_tropo,pres_bottom_strato,size_strato
 
 contains
@@ -38 +25 @@
       use ioipsl_getin_p_mod, only: getin_p
       use radinc_h, only: naerkind
-      use aerosol_mod, only: iaero_back2lay, iaero_co2, iaero_dust, &
-                             iaero_h2o, iaero_h2so4
+      use aerosol_mod, only: iaero_back2lay
       Implicit none
 
@@ -59 +45 @@
 !     .def file. To be improved!
 
-            if(iaer.eq.iaero_co2)then ! CO2 ice
+
+!     WARNING : Titan adapt. (J. Vatant d'Ollone, 2017)
+!            - ONLY THE NO AEROSOL CASE FOR NOW SINCE WE COMPUTE THEM ANOTHER WAY !
+!            - This routine is just here to keep the code running without unplugging all (yet)
+!            - There's only the dummy aerosol case on iaer = 1      
+            if(iaer.eq.1)then
                reffrad(1:ngrid,1:nlayer,iaer) = 1.e-4
-               nueffrad(1:ngrid,1:nlayer,iaer) = 0.1 
-            endif
-
-            if(iaer.eq.iaero_h2o)then ! H2O ice
-               reffrad(1:ngrid,1:nlayer,iaer) = 1.e-5
-               nueffrad(1:ngrid,1:nlayer,iaer) = 0.1 
-            endif
-
-            if(iaer.eq.iaero_dust)then ! dust
-               reffrad(1:ngrid,1:nlayer,iaer) = 1.e-5
-               nueffrad(1:ngrid,1:nlayer,iaer) = 0.1 
-            endif
- 
-            if(iaer.eq.iaero_h2so4)then ! H2O ice
-               reffrad(1:ngrid,1:nlayer,iaer) = 1.e-6
                nueffrad(1:ngrid,1:nlayer,iaer) = 0.1 
             endif
@@ -83 +59 @@
                nueffrad(1:ngrid,1:nlayer,iaer) = 0.1 
             endif
-
-
 
             if(iaer.gt.5)then
@@ -95 +69 @@
          enddo
 
-
-         if (radfixed) then
-
-            write(*,*)"radius of H2O water particles:"
-            rad_h2o=13. ! default value
-            call getin_p("rad_h2o",rad_h2o)
-            write(*,*)" rad_h2o = ",rad_h2o
-
-            write(*,*)"radius of H2O ice particles:"
-            rad_h2o_ice=35. ! default value
-            call getin_p("rad_h2o_ice",rad_h2o_ice)
-            write(*,*)" rad_h2o_ice = ",rad_h2o_ice
-
-         else
-
-            write(*,*)"Number mixing ratio of H2O water particles:"
-            Nmix_h2o=1.e6 ! default value
-            call getin_p("Nmix_h2o",Nmix_h2o)
-            write(*,*)" Nmix_h2o = ",Nmix_h2o
-
-            write(*,*)"Number mixing ratio of H2O ice particles:"
-            Nmix_h2o_ice=Nmix_h2o ! default value
-            call getin_p("Nmix_h2o_ice",Nmix_h2o_ice)
-            write(*,*)" Nmix_h2o_ice = ",Nmix_h2o_ice
-         endif
-
       print*,'exit su_aer_radii'
 
@@ -126 +74 @@
 !==================================================================
 
-
-!==================================================================
-   subroutine h2o_reffrad(ngrid,nlayer,pq,pt,reffrad,nueffrad)
-!==================================================================
-!     Purpose
-!     -------
-!     Compute the effective radii of liquid and icy water particles
-!
-!     Authors
-!     -------
-!     Jeremy Leconte (2012)
-!
-!==================================================================
-      use watercommon_h, Only: T_h2O_ice_liq,T_h2O_ice_clouds,rhowater,rhowaterice
-      use comcstfi_mod, only: pi
-      Implicit none
-
-      integer,intent(in) :: ngrid
-      integer,intent(in) :: nlayer
-
-      real, intent(in) :: pq(ngrid,nlayer) !water ice mixing ratios (kg/kg)
-      real, intent(in) :: pt(ngrid,nlayer) !temperature (K)
-      real, intent(out) :: reffrad(ngrid,nlayer)      !aerosol radii
-      real, intent(out) :: nueffrad(ngrid,nlayer) ! dispersion      
-
-      integer :: ig,l
-      real zfice ,zrad,zrad_liq,zrad_ice
-      real,external :: CBRT            
-      
-
-      if (radfixed) then
-         do l=1,nlayer
-            do ig=1,ngrid
-               zfice = 1.0 - (pt(ig,l)-T_h2O_ice_clouds) / (T_h2O_ice_liq-T_h2O_ice_clouds)
-               zfice = MIN(MAX(zfice,0.0),1.0)
-               reffrad(ig,l)= rad_h2o * (1.-zfice) + rad_h2o_ice * zfice
-               nueffrad(ig,l) = coef_chaud * (1.-zfice) + coef_froid * zfice
-            enddo
-         enddo
-      else
-         do l=1,nlayer
-            do ig=1,ngrid
-               zfice = 1.0 - (pt(ig,l)-T_h2O_ice_clouds) / (T_h2O_ice_liq-T_h2O_ice_clouds)
-               zfice = MIN(MAX(zfice,0.0),1.0)
-               zrad_liq  = CBRT( 3*pq(ig,l)/(4*Nmix_h2o*pi*rhowater) )
-               zrad_ice  = CBRT( 3*pq(ig,l)/(4*Nmix_h2o_ice*pi*rhowaterice) )
-               nueffrad(ig,l) = coef_chaud * (1.-zfice) + coef_froid * zfice
-               zrad = zrad_liq * (1.-zfice) + zrad_ice * zfice
-
-               reffrad(ig,l) = min(max(zrad,1.e-6),1000.e-6)
-               enddo
-            enddo      
-      end if
-
-   end subroutine h2o_reffrad
-!==================================================================
-
-
-!==================================================================
-   subroutine h2o_cloudrad(ngrid,nlayer,pql,reffliq,reffice)
-!==================================================================
-!     Purpose
-!     -------
-!     Compute the effective radii of liquid and icy water particles
-!
-!     Authors
-!     -------
-!     Jeremy Leconte (2012)
-!
-!==================================================================
-      use watercommon_h, Only: rhowater,rhowaterice
-      use comcstfi_mod, only: pi
-      Implicit none
-
-      integer,intent(in) :: ngrid
-      integer,intent(in) :: nlayer
-
-      real, intent(in) :: pql(ngrid,nlayer) !condensed water mixing ratios (kg/kg)
-      real, intent(out) :: reffliq(ngrid,nlayer),reffice(ngrid,nlayer)     !liquid and ice water particle radii (m)
-
-      real,external :: CBRT            
-      integer :: i,k
-
-      if (radfixed) then
-         reffliq(1:ngrid,1:nlayer)= rad_h2o
-         reffice(1:ngrid,1:nlayer)= rad_h2o_ice
-      else
-         do k=1,nlayer
-           do i=1,ngrid
-             reffliq(i,k) = CBRT(3*pql(i,k)/(4*Nmix_h2o*pi*rhowater))
-             reffliq(i,k) = min(max(reffliq(i,k),1.e-6),1000.e-6)
-           
-             reffice(i,k) = CBRT(3*pql(i,k)/(4*Nmix_h2o_ice*pi*rhowaterice))
-             reffice(i,k) = min(max(reffice(i,k),1.e-6),1000.e-6)
-           enddo
-         enddo
-      endif
-
-   end subroutine h2o_cloudrad
-!==================================================================
-
-
-
-!==================================================================
-   subroutine co2_reffrad(ngrid,nlayer,nq,pq,reffrad)
-!==================================================================
-!     Purpose
-!     -------
-!     Compute the effective radii of co2 ice particles
-!
-!     Authors
-!     -------
-!     Jeremy Leconte (2012)
-!
-!==================================================================
-      USE tracer_h, only:igcm_co2_ice,rho_co2
-      use comcstfi_mod, only: pi
-      Implicit none
-
-      integer,intent(in) :: ngrid,nlayer,nq
-
-      real, intent(in) :: pq(ngrid,nlayer,nq) !tracer mixing ratios (kg/kg)
-      real, intent(out) :: reffrad(ngrid,nlayer)      !co2 ice particles radii (m)
-
-      integer :: ig,l
-      real :: zrad   
-      real,external :: CBRT            
-            
-      
-
-      if (radfixed) then
-         reffrad(1:ngrid,1:nlayer) = 5.e-5 ! CO2 ice
-      else
-         do l=1,nlayer
-            do ig=1,ngrid
-               zrad = CBRT( 3*pq(ig,l,igcm_co2_ice)/(4*Nmix_co2*pi*rho_co2) )
-               reffrad(ig,l) = min(max(zrad,1.e-6),100.e-6)
-            enddo
-         enddo      
-      end if
-
-   end subroutine co2_reffrad
-!==================================================================
-
-
-
-!==================================================================
-   subroutine dust_reffrad(ngrid,nlayer,reffrad)
-!==================================================================
-!     Purpose
-!     -------
-!     Compute the effective radii of dust particles
-!
-!     Authors
-!     -------
-!     Jeremy Leconte (2012)
-!
-!==================================================================
-      Implicit none
-
-      integer,intent(in) :: ngrid
-      integer,intent(in) :: nlayer
-
-      real, intent(out) :: reffrad(ngrid,nlayer)      !dust particles radii (m)
-            
-      reffrad(1:ngrid,1:nlayer) = 2.e-6 ! dust
-
-   end subroutine dust_reffrad
-!==================================================================
-
-
-!==================================================================
-   subroutine h2so4_reffrad(ngrid,nlayer,reffrad)
-!==================================================================
-!     Purpose
-!     -------
-!     Compute the effective radii of h2so4 particles
-!
-!     Authors
-!     -------
-!     Jeremy Leconte (2012)
-!
-!==================================================================
-      Implicit none
-
-      integer,intent(in) :: ngrid
-      integer,intent(in) :: nlayer
-
-      real, intent(out) :: reffrad(ngrid,nlayer)      !h2so4 particle radii (m)
-                
-      reffrad(1:ngrid,1:nlayer) = 1.e-6 ! h2so4
-
-   end subroutine h2so4_reffrad
-!==================================================================
 
 !==================================================================

trunk/LMDZ.TITAN/libf/phytitan/su_gases.F90

@@ -73 +73 @@
            igas_H2=igas
            count=count+1
-        elseif (trim(gnom(igas)).eq."He_" .or. trim(gnom(igas)).eq."He") then
-           igas_He=igas
-           count=count+1
-        elseif (trim(gnom(igas)).eq."H2O") then
-           igas_H2O=igas
-           count=count+1
-        elseif (trim(gnom(igas)).eq."CO2") then
-           igas_CO2=igas
-           count=count+1
-        elseif (trim(gnom(igas)).eq."CO_" .or. trim(gnom(igas)).eq."CO") then
-           igas_CO=igas
-           count=count+1
         elseif (trim(gnom(igas)).eq."N2_" .or. trim(gnom(igas)).eq."N2") then
            igas_N2=igas
            count=count+1
-        elseif (trim(gnom(igas)).eq."O2_" .or. trim(gnom(igas)).eq."O2") then
-           igas_O2=igas
-           count=count+1
-        elseif (trim(gnom(igas)).eq."SO2") then
-           igas_SO2=igas
-           count=count+1
-        elseif (trim(gnom(igas)).eq."H2S") then
-           igas_H2S=igas
-           count=count+1
         elseif (trim(gnom(igas)).eq."CH4") then
            igas_CH4=igas
-           count=count+1
-        elseif (trim(gnom(igas)).eq."NH3") then
-           igas_NH3=igas
            count=count+1
         elseif (trim(gnom(igas)).eq."C2H6") then

trunk/LMDZ.TITAN/libf/phytitan/suaer_corrk.F90

@@ -42 +42 @@
 !     MODIF R. Wordsworth (2009)
 !     - generalisation to arbitrary spectral bands 
-!
-!==================================================================
+!     
+!     WARNING : Titan adapt. (J. Vatant d'Ollone, 2017)
+!            - ONLY THE NO AEROSOL CASE FOR NOW SINCE WE COMPUTE THEM ANOTHER WAY !
+!            - This routine is just here to keep the code running without unplugging all (yet)                           
+!================================================================================================
 
 !     Optical properties (read in external ASCII files)
@@ -106 +109 @@
       endif
       do iaer=1,naerkind
-       if (iaer.eq.iaero_co2) then
+       if (iaer.eq.1) then ! JVO 2017 : Now iaer = 1 is always dummy co2 for noaero case, since we don't use aerosols
         forgetit=.true.
-          if (.not.noaero) then
-              print*, 'naerkind= co2', iaer
-          end if
 !     visible
         if(forgetit)then
@@ -130 +130 @@
 ! JB thinks this shouldn't matter too much, but it needs to be 
 ! fixed / decided for the final version.
-
-       if (iaer.eq.iaero_h2o) then
-        print*, 'naerkind= h2o', iaer
-
-!     visible
-         file_id(iaer,1) = 'optprop_icevis_n50.dat'
-         lamrefvis(iaer)=1.5E-6   ! 1.5um OMEGA/MEx
-!     infrared
-         file_id(iaer,2) = 'optprop_iceir_n50.dat'
-         lamrefir(iaer)=12.1E-6   ! 825cm-1 TES/MGS
-       endif
-
-       if (iaer.eq.iaero_dust) then
-        print*, 'naerkind= dust', iaer
-
-!     visible
-         file_id(iaer,1) = 'optprop_dustvis_n50.dat'
-         !lamrefvis(iaer)=1.5E-6   ! 1.5um OMEGA/MEx
-         lamrefvis(iaer)=0.67e-6
-!     infrared
-         file_id(iaer,2) = 'optprop_dustir_n50.dat'
-         !lamrefir(iaer)=12.1E-6   ! 825cm-1 TES/MGS
-         lamrefir(iaer)=9.3E-6
-       endif 
-
-       if (iaer.eq.iaero_h2so4) then
-         print*, 'naerkind= h2so4', iaer
-
-!     visible
-         file_id(iaer,1) = 'optprop_h2so4vis_n20.dat'
-         !lamrefir(iaer)=  doesn't exist?
-         lamrefvis(iaer)=1.5E-6   ! no idea, must find
-!     infrared
-         file_id(iaer,2) = 'optprop_h2so4ir_n20.dat'
-         !lamrefir(iaer)=12.1E-6   ! 825cm-1 TES/MGS
-         lamrefir(iaer)=9.3E-6 ! no idea, must find
-! added by LK
-       endif
 
        if (iaer.eq.iaero_back2lay) then
@@ -328 +290 @@
 
       jfile = jfile+1
-      IF ((idomain.NE.iaero_co2).OR.(iaer.NE.iaero_co2)) THEN
+      IF ((idomain.NE.1).OR.(iaer.NE.1)) THEN
          endwhile = .true.
       ENDIF
@@ -346 +308 @@
 
 !     1.5 If Francois' data, convert wvl to metres
-       if(iaer.eq.iaero_co2)then
+       if(iaer.eq.1)then
          if(forgetit)then
          !   print*,'please sort out forgetit for naerkind>1'

trunk/LMDZ.TITAN/libf/phytitan/sugas_corrk.F90

@@ -29 +29 @@
       use gases_h
       use ioipsl_getin_p_mod, only: getin_p
-      use callkeys_mod, only: varactive,varfixed,graybody,callgasvis,&
-                continuum,H2Ocont_simple
+      use callkeys_mod, only: graybody,callgasvis, continuum
       implicit none
 
@@ -86 +85 @@
       endif 
 
-      if(ngas.eq.1 .and. (varactive.or.varfixed))then
-         print*,'You have varactive/fixed=.true. but the database [',  &
-                     corrkdir(1:LEN_TRIM(corrkdir)),           &
-                '] has no variable species, exiting.'
-         call abort
-      elseif(ngas.gt.5 .or. ngas.lt.1)then
+      if(ngas.gt.5 .or. ngas.lt.1)then
          print*,ngas,' species in database [',               &
                      corrkdir(1:LEN_TRIM(corrkdir)),           &
@@ -112 +106 @@
       read(111,*) wrefvar
       close(111)
-
-      if(L_REFVAR.gt.1 .and. (.not.varactive) .and. (.not.varfixed))then
-         print*,'You have varactive and varfixed=.false. and the database [', &
-                     corrkdir(1:LEN_TRIM(corrkdir)),           &
-                '] has a variable species.'
-         call abort
-      endif
 
       ! Check that gastype and gnom match
@@ -638 +625 @@
                   dummy = -9999
                   call interpolateN2H2(592.D+0,278.15D+0,200000.D+0,10000.D+0,testcont,.true.,dummy)
-               elseif (jgas .eq. igas_He) then
-                  dummy = -9999
-                  call interpolateH2He(500.D+0,250.D+0,200000.D+0,10000.D+0,testcont,.true.,dummy)
                endif
             enddo
-
-         elseif (igas .eq. igas_H2O) then
-
-            ! H2O is special 
-            if(H2Ocont_simple)then
-               call interpolateH2Ocont_PPC(990.D+0,296.D+0,683.2D+0*2,0.D+0,testcont,.true.)
-            else
-               dummy = -9999
-               call interpolateH2Ocont_CKD(990.D+0,296.D+0,683.2D+0*2,0.D+0,testcont,.true.,dummy)
-            endif
 
          endif

trunk/LMDZ.TITAN/libf/phytitan/surfdat_h.F90

@@ -17 +17 @@
 !$OMP THREADPRIVATE(zmea,zstd,zsig,zgam,zthe)
 
-       real,allocatable,dimension(:) :: dryness  !"Dryness coefficient" for grnd water ice sublimation
-                                                 ! AS: previously in tracer.h. it is more logical here.
-
-       logical,allocatable,dimension(:) :: watercaptag !! was in watercap.h
-!$OMP THREADPRIVATE(dryness,watercaptag)
-
        end module surfdat_h
 

trunk/LMDZ.TITAN/libf/phytitan/surfini.F

@@ -1 @@
-      SUBROUTINE surfini(ngrid,nq,qsurf,albedo,albedo_bareground,
-     &                   albedo_snow_SPECTV,albedo_co2_ice_SPECTV)
+      SUBROUTINE surfini(ngrid,nq,qsurf,albedo,albedo_bareground)
 
       USE surfdat_h, only: albedodat
-      USE tracer_h, only: igcm_co2_ice
       use planetwide_mod, only: planetwide_maxval, planetwide_minval
       use radinc_h, only : L_NSPECTV
-      use callkeys_mod, only : albedosnow, albedoco2ice
 
       IMPLICIT NONE
@@ -26 +23 @@
       REAL,INTENT(OUT) :: albedo(ngrid,L_NSPECTV)
       REAL,INTENT(OUT) :: albedo_bareground(ngrid)
-      REAL,INTENT(OUT) :: albedo_snow_SPECTV(L_NSPECTV)
-      REAL,INTENT(OUT) :: albedo_co2_ice_SPECTV(L_NSPECTV)
       REAL,INTENT(IN) :: qsurf(ngrid,nq) ! tracer on surface (kg/m2)
 
@@ -35 +30 @@
 c=======================================================================
 
-      ! Step 1 : Initialisation of the Spectral Albedos.
-      DO nw=1,L_NSPECTV
-         albedo_snow_SPECTV(nw)=albedosnow
-         albedo_co2_ice_SPECTV(nw)=albedoco2ice
-      ENDDO
-
-
-      ! Step 2 : We get the bare ground albedo from the start files.
+      ! We get the bare ground albedo from the start files.
       DO ig=1,ngrid
          albedo_bareground(ig)=albedodat(ig)
@@ -54 +42 @@
       write(*,*) 'surfini: maximum bare ground albedo',max_albedo
 
-
-      ! Step 3 : We modify the albedo considering some CO2 at the surface. We dont take into account water ice (this is made in hydrol after the first timestep) ...
-      if (igcm_co2_ice.ne.0) then
-         DO ig=1,ngrid
-            IF (qsurf(ig,igcm_co2_ice) .GT. 1.) THEN ! This was changed by MT2015. Condition for ~1mm of CO2 ice deposit.
-               DO nw=1,L_NSPECTV
-                  albedo(ig,nw)=albedo_co2_ice_SPECTV(nw)
-               ENDDO
-            END IF   
-         ENDDO   
-      else
-         write(*,*) "surfini: No CO2 ice tracer on surface  ..."
-         write(*,*) "         and therefore no albedo change."
-      endif     
+    
       call planetwide_minval(albedo,min_albedo)
       call planetwide_maxval(albedo,max_albedo)

trunk/LMDZ.TITAN/libf/phytitan/tracer_h.F90

@@ -20 +20 @@
       real r3n_q     ! used to compute r0 from number and mass mixing ratio
       real rho_dust     ! Mars dust density (kg.m-3)
-      real rho_ice     ! Water ice density (kg.m-3)
-      real rho_co2     ! CO2 ice density (kg.m-3)
       real ref_r0        ! for computing reff=ref_r0*r0 (in log.n. distribution)
 !$OMP THREADPRIVATE(noms,mmol,radius,rho_q,qext,alpha_lift,alpha_devil,qextrhor, &
-        !$OMP varian,r3n_q,rho_dust,rho_ice,rho_co2,ref_r0)
+        !$OMP varian,r3n_q,rho_dust,rho_ice,ref_r0)
 
 ! tracer indexes: these are initialized in initracer and should be 0 if the
@@ -33 +31 @@
       integer :: igcm_dust_mass   ! dust mass mixing ratio (for transported dust)
       integer :: igcm_dust_number ! dust number mixing ratio (transported dust)
-      ! water
-      integer :: igcm_h2o_vap ! water vapour
-      integer :: igcm_h2o_ice ! water ice
       ! chemistry:
-      integer :: igcm_co2
       integer :: igcm_co
       integer :: igcm_o
@@ -52 +46 @@
       ! other tracers
       integer :: igcm_ar_n2 ! for simulations using co2 +neutral gaz
-      integer :: igcm_co2_ice ! CO2 ice 
-!$OMP THREADPRIVATE(igcm_dustbin,igcm_dust_mass,igcm_dust_number,igcm_h2o_vap,igcm_h2o_ice, &
-        !$OMP igcm_co2,igcm_co,igcm_o,igcm_o1d,igcm_o2,igcm_o3,igcm_h,igcm_h2,igcm_oh,      &
-        !$OMP igcm_ho2,igcm_h2o2,igcm_n2,igcm_ar,igcm_ar_n2,igcm_co2_ice)
+!$OMP THREADPRIVATE(igcm_dustbin,igcm_dust_mass,igcm_dust_number, &
+        !$OMP igcm_co,igcm_o,igcm_o1d,igcm_o2,igcm_o3,igcm_h,igcm_h2,igcm_oh,       &
+        !$OMP igcm_ho2,igcm_h2o2,igcm_n2,igcm_ar,igcm_ar_n2)
 
        end module tracer_h

trunk/LMDZ.TITAN/libf/phytitan/turbdiff.F90

@@ -1 @@
-      subroutine turbdiff(ngrid,nlay,nq,rnat,          &
+      subroutine turbdiff(ngrid,nlay,nq,               &
           ptimestep,pcapcal,lecrit,                    &   
           pplay,pplev,pzlay,pzlev,pz0,                 &
@@ -5 +5 @@
           pdtfi,pdqfi,pfluxsrf,            &
           Pdudif,pdvdif,pdtdif,pdtsrf,sensibFlux,pq2,  &
-          pdqdif,pdqevap,pdqsdif,flux_u,flux_v,lastcall)
-
-      use watercommon_h, only : RLVTT, T_h2O_ice_liq, RCPD, mx_eau_sol,Psat_water
+          pdqdif,pdqsdif,flux_u,flux_v,lastcall)
+
       use radcommon_h, only : sigma, glat
-      use surfdat_h, only: dryness
-      use tracer_h, only: igcm_h2o_vap, igcm_h2o_ice
       use comcstfi_mod, only: rcp, g, r, cpp
-      use callkeys_mod, only: water,tracer,nosurf
+      use callkeys_mod, only: tracer,nosurf
 
       implicit none
@@ -62 +59 @@
       REAL,INTENT(IN) :: pcapcal(ngrid)
       REAL,INTENT(INOUT) :: pq2(ngrid,nlay+1)
-      REAL,INTENT(OUT) :: flux_u(ngrid),flux_v(ngrid)
-      REAL,INTENT(IN) :: rnat(ngrid)      
+      REAL,INTENT(OUT) :: flux_u(ngrid),flux_v(ngrid)  
       LOGICAL,INTENT(IN) :: lastcall ! not used
 
@@ -109 +105 @@
       INTEGER iq
       REAL zq(ngrid,nlay,nq)
-      REAL zqnoevap(ngrid,nlay) !special for water case to compute where evaporated water goes.
-      REAL pdqevap(ngrid,nlay) !special for water case to compute where evaporated water goes.
       REAL zdmassevap(ngrid)
       REAL rho(ngrid)         ! near-surface air density
-      REAL qsat(ngrid),psat(ngrid)
       REAL kmixmin
 
-!     Variables added for implicit latent heat inclusion
-!     --------------------------------------------------
-      real dqsat(ngrid), qsat_temp1, qsat_temp2,psat_temp
-
-      integer, save :: ivap, iliq, iliq_surf,iice_surf ! also make liq for clarity on surface...
-!$OMP THREADPRIVATE(ivap,iliq,iliq_surf,iice_surf)
 
       real, parameter :: karman=0.4
@@ -133 +120 @@
 !     --------------
 
-      IF (firstcall) THEN
-
-         if(water)then
-             ivap=igcm_h2o_vap
-             iliq=igcm_h2o_ice
-             iliq_surf=igcm_h2o_vap
-             iice_surf=igcm_h2o_ice ! simply to make the code legible               
-                                  ! to be generalised
-         else
-             ivap=0
-             iliq=0
-             iliq_surf=0
-             iice_surf=0 ! simply to make the code legible                        
-         endif
+      IF (firstcall) THEN         
+
          sensibFlux(:)=0.
 
@@ -202 +177 @@
             ENDDO
          ENDDO
-         if (water) then
-            DO ilev=1,nlay
-               DO ig=1,ngrid
-                  zqnoevap(ig,ilev)=pq(ig,ilev,ivap) + pdqfi(ig,ilev,ivap)*ptimestep
-               ENDDO
-            ENDDO
-         Endif
       end if
 
@@ -405 +373 @@
 !       flux (if any) from subsurface
 
-      if(.not.water) then
 
          DO ig=1,ngrid
@@ -427 +394 @@
          ENDDO
 
-      endif                     ! not water
 
 !-----------------------------------------------------------------------
@@ -446 +412 @@
 !     -----------------------
          do iq=1,nq 
-
-            if (iq.ne.ivap) then
 
                DO ig=1,ngrid
@@ -462 +426 @@
                   ENDDO
                ENDDO
-
-               if ((water).and.(iq.eq.iliq)) then
-                  ! special case for condensed water tracer: do not include
-                  ! h2o ice tracer from surface (which is set when handling
-                  ! h2o vapour case (see further down).
-                  ! zb(ig,1)=0 if iq ne ivap
-                  DO ig=1,ngrid
-                     z1(ig)=1./(zmass(ig,1)+zfluxq(ig,2)*(1.-zdq(ig,2)))
-                     zcq(ig,1)=(zmass(ig,1)*zq(ig,1,iq)+zfluxq(ig,2)*zcq(ig,2))*z1(ig)
-                  ENDDO
-               else             ! general case
-                  do ig=1,ngrid
-                     z1(ig)=1./(zmass(ig,1)+zfluxq(ig,2)*(1.-zdq(ig,2)))
-                     zcq(ig,1)=(zmass(ig,1)*zq(ig,1,iq)+zfluxq(ig,2)*zcq(ig,2)+(-pdqsdif(ig,iq))*ptimestep)*z1(ig)
-                          ! tracer flux from surface
-                          ! currently pdqsdif always zero here,
-                          ! so last line is superfluous
-                  enddo
-               endif            ! of if (water.and.(iq.eq.igcm_h2o_ice))
-
+               
+               do ig=1,ngrid
+                  z1(ig)=1./(zmass(ig,1)+zfluxq(ig,2)*(1.-zdq(ig,2)))
+                  zcq(ig,1)=(zmass(ig,1)*zq(ig,1,iq)+zfluxq(ig,2)*zcq(ig,2)+(-pdqsdif(ig,iq))*ptimestep)*z1(ig)
+                  ! tracer flux from surface
+                  ! currently pdqsdif always zero here,
+                  ! so last line is superfluous
+               enddo
 
 !     Starting upward calculations for simple tracer mixing (e.g., dust)
@@ -493 +445 @@
                   end do
                end do
-
-            endif               ! if (iq.ne.ivap)
-
-!     Calculate temperature tendency including latent heat term
-!     and assuming an infinite source of water on the ground
-!     ------------------------------------------------------------------
-
-            if (water.and.(iq.eq.ivap)) then 
-            
-               ! compute evaporation efficiency
-               do ig=1,ngrid
-                  if(nint(rnat(ig)).eq.1)then
-                     dryness(ig)=pqsurf(ig,iliq_surf)+pqsurf(ig,iice_surf) 
-                     dryness(ig)=MIN(1.,2*dryness(ig)/mx_eau_sol)
-                     dryness(ig)=MAX(0.,dryness(ig))
-                  endif
-               enddo
-
-               do ig=1,ngrid
-                ! Calculate the value of qsat at the surface (water)
-                call Psat_water(ptsrf(ig),pplev(ig,1),psat(ig),qsat(ig))
-                call Psat_water(ptsrf(ig)-0.0001,pplev(ig,1),psat_temp,qsat_temp1)
-                call Psat_water(ptsrf(ig)+0.0001,pplev(ig,1),psat_temp,qsat_temp2)
-                dqsat(ig)=(qsat_temp2-qsat_temp1)/0.0002
-                ! calculate dQsat / dT by finite differences
-                ! we cannot use the updated temperature value yet...
-               enddo
-
-! coefficients for q
-
-               do ig=1,ngrid
-                  z1(ig)=1./(zmass(ig,nlay)+zfluxq(ig,nlay))
-                  zcq(ig,nlay)=zmass(ig,nlay)*zq(ig,nlay,iq)*z1(ig)
-                  zdq(ig,nlay)=zfluxq(ig,nlay)*z1(ig)
-               enddo 
-          
-               do ilay=nlay-1,2,-1
-                  do ig=1,ngrid
-                     z1(ig)=1./(zmass(ig,ilay)+zfluxq(ig,ilay)+zfluxq(ig,ilay+1)*(1.-zdq(ig,ilay+1)))
-                     zcq(ig,ilay)=(zmass(ig,ilay)*zq(ig,ilay,iq)+zfluxq(ig,ilay+1)*zcq(ig,ilay+1))*z1(ig)
-                     zdq(ig,ilay)=zfluxq(ig,ilay)*z1(ig)
-                  enddo
-               enddo
-
-               do ig=1,ngrid
-                  z1(ig)=1./(zmass(ig,1)+zfluxq(ig,1)*dryness(ig)+zfluxq(ig,2)*(1.-zdq(ig,2)))
-                  zcq(ig,1)=(zmass(ig,1)*zq(ig,1,iq)+zfluxq(ig,2)*zcq(ig,2))*z1(ig)
-                  zdq(ig,1)=dryness(ig)*zfluxq(ig,1)*z1(ig)
-               enddo
-
-              do ig=1,ngrid
-!calculation of surface temperature
-                  zdq0(ig) = dqsat(ig)
-                  zcq0(ig) = qsat(ig)-dqsat(ig)*ptsrf(ig)
-
-                  z1(ig) = pcapcal(ig)*ptsrf(ig) +cpp*zfluxt(ig,1)*zct(ig,1)*zovExner(ig,1)   &
-                      + zdplanck(ig)*ptsrf(ig) + pfluxsrf(ig)*ptimestep                       &
-                      + zfluxq(ig,1)*dryness(ig)*RLVTT*((zdq(ig,1)-1.0)*zcq0(ig)+zcq(ig,1))
-
-                  z2(ig) = pcapcal(ig) + cpp*zfluxt(ig,1)*(1.-zovExner(ig,1)*zdt(ig,1))       &
-                      + zdplanck(ig)+zfluxq(ig,1)*dryness(ig)*RLVTT*zdq0(ig)*(1.0-zdq(ig,1))
-
-                  ztsrf(ig) = z1(ig) / z2(ig)
-
-! calculation of qs and q1
-                  zq0(ig)     = zcq0(ig)+zdq0(ig)*ztsrf(ig)
-                  zq(ig,1,iq) = zcq(ig,1)+zdq(ig,1)*zq0(ig)
-
-! calculation of evaporation              
-                  dqsdif_total(ig)=zfluxq(ig,1)*dryness(ig)*(zq(ig,1,ivap)-zq0(ig))
-
-!     --------------------------------------------------------
-!     Now check if we've taken too much water from the surface
-!     This can only occur on the continent 
-!     If we do, we recompute Tsurf, T1 and q1 accordingly
-                  if((-dqsdif_total(ig).gt.(pqsurf(ig,iice_surf)+pqsurf(ig,iliq_surf))).and.rnat(ig).eq.1)then
-                      !water flux * ptimestep
-                      dqsdif_total(ig)=-(pqsurf(ig,iice_surf)+pqsurf(ig,iliq_surf))
-
-                      !recompute surface temperature  
-                      z1(ig) = pcapcal(ig)*ptsrf(ig) +cpp*zfluxq(ig,1)*zct(ig,1)*zovExner(ig,1)   &
-                        + zdplanck(ig)*ptsrf(ig) + pfluxsrf(ig)*ptimestep                       &
-                        + RLVTT*dqsdif_total(ig)
-                      z2(ig) = pcapcal(ig) + cpp*zfluxq(ig,1)*(1.-zovExner(ig,1)*zdt(ig,1))       &
-                        + zdplanck(ig)
-                      ztsrf(ig) = z1(ig) / z2(ig)
-
-                      !recompute q1 with new water flux from surface  
-                      zq(ig,1,iq) = (zmass(ig,1)*(pq(ig,1,iq)+ptimestep*pdqfi(ig,1,iq))  &
-                                            +zfluxq(ig,2)*zcq(ig,2)-dqsdif_total(ig))     &
-                                 / (zmass(ig,1)+(1.-zdq(ig,2))*zfluxq(ig,2))                 
-                  end if
-                  
-! calculation surface T tendency  and T(1)           
-                  pdtsrf(ig) = (ztsrf(ig) - ptsrf(ig))/ptimestep
-                  zt(ig,1)   = zct(ig,1) + zdt(ig,1)*ztsrf(ig)                
-               enddo
-
-
-! recalculate temperature and q(vap) to top of atmosphere, starting from ground
-               do ilay=2,nlay
-                  do ig=1,ngrid
-                     zq(ig,ilay,iq)=zcq(ig,ilay)+zdq(ig,ilay)*zq(ig,ilay-1,iq)
-                     zt(ig,ilay)=zct(ig,ilay)+zdt(ig,ilay)*zt(ig,ilay-1)
-                  end do
-               end do
-
-
-               do ig=1,ngrid
-!     --------------------------------------------------------------------------
-!     On the ocean, if T > 0 C then the vapour tendency must replace the ice one
-!     The surface vapour tracer is actually liquid. To make things difficult.
-
-                  if (nint(rnat(ig)).eq.0) then ! unfrozen ocean
-                     
-                     pdqsdif(ig,iliq_surf)=dqsdif_total(ig)/ptimestep
-                     pdqsdif(ig,iice_surf)=0.0
-
-                  elseif (nint(rnat(ig)).eq.1) then ! (continent)
-!     If water is evaporating / subliming, we take it from ice before liquid
-!     -- is this valid??
-                     if(dqsdif_total(ig).lt.0)then
-                        if (-dqsdif_total(ig).gt.pqsurf(ig,iice_surf))then
-                           pdqsdif(ig,iice_surf) = -pqsurf(ig,iice_surf)/ptimestep ! removes all the ice!
-                           pdqsdif(ig,iliq_surf) = dqsdif_total(ig)/ptimestep- pdqsdif(ig,iice_surf) ! take the remainder from the liquid instead
-                        else               
-                           pdqsdif(ig,iice_surf)=dqsdif_total(ig)/ptimestep
-                           pdqsdif(ig,iliq_surf)=0.
-                        end if
-                     else !dqsdif_total(ig).ge.0
-                        !If water vapour is condensing, we must decide whether it forms ice or liquid.
-                        if(ztsrf(ig).gt.T_h2O_ice_liq)then
-                           pdqsdif(ig,iice_surf)=0.0
-                           pdqsdif(ig,iliq_surf)=dqsdif_total(ig)/ptimestep
-                        else
-                           pdqsdif(ig,iice_surf)=dqsdif_total(ig)/ptimestep
-                           pdqsdif(ig,iliq_surf)=0.0
-                        endif               
-                     endif
-
-                  elseif (nint(rnat(ig)).eq.2) then ! (continental glaciers)
-                     pdqsdif(ig,iliq_surf)=0.0
-                     pdqsdif(ig,iice_surf)=dqsdif_total(ig)/ptimestep
-
-                  endif !rnat
-               end do            ! of DO ig=1,ngrid
-
-           endif                ! if (water et iq=ivap)
         end do                  ! of do iq=1,nq
-
-        if (water) then  ! special case where we recompute water mixing without any evaporation.
-                         !    The difference with the first calculation then tells us where evaporated water has gone
-
-            DO ig=1,ngrid
-               z1(ig)=1./(zmass(ig,nlay)+zfluxq(ig,nlay))
-               zcq(ig,nlay)=zmass(ig,nlay)*zqnoevap(ig,nlay)*z1(ig)
-               zdq(ig,nlay)=zfluxq(ig,nlay)*z1(ig)
-            ENDDO 
-            
-            DO ilay=nlay-1,2,-1
-               DO ig=1,ngrid
-                  z1(ig)=1./(zmass(ig,ilay)+zfluxq(ig,ilay)+zfluxq(ig,ilay+1)*(1.-zdq(ig,ilay+1)))
-                  zcq(ig,ilay)=(zmass(ig,ilay)*zqnoevap(ig,ilay)+zfluxq(ig,ilay+1)*zcq(ig,ilay+1))*z1(ig)
-                  zdq(ig,ilay)=zfluxq(ig,ilay)*z1(ig)
-               ENDDO
-            ENDDO
-
-            do ig=1,ngrid
-               z1(ig)=1./(zmass(ig,1)+zfluxq(ig,2)*(1.-zdq(ig,2)))
-               zcq(ig,1)=(zmass(ig,1)*zqnoevap(ig,1)+zfluxq(ig,2)*zcq(ig,2))*z1(ig)
-            enddo
-
-!     Starting upward calculations for simple tracer mixing (e.g., dust)
-            do ig=1,ngrid
-               zqnoevap(ig,1)=zcq(ig,1)
-            end do
-
-            do ilay=2,nlay
-               do ig=1,ngrid
-                  zqnoevap(ig,ilay)=zcq(ig,ilay)+zdq(ig,ilay)*zqnoevap(ig,ilay-1)
-               end do
-            end do
-
-         endif               ! if water
-        
-        
+                       
       endif                     ! tracer
 
@@ -706 +474 @@
             enddo
          enddo
-         if (water) then
-            do ilev = 1, nlay
-               do ig=1,ngrid
-                  pdqevap(ig,ilev)=(zq(ig,ilev,ivap)-zqnoevap(ig,ilev))/ptimestep
-               enddo
-            enddo
-            do ig=1,ngrid
-               zdmassevap(ig)=SUM(pdqevap(ig,:)*zmass(ig,:))*ptimestep
-            end do          
-         endif
       endif 
-
-      if(water)then
-         call writediagfi(ngrid,'beta','Dryness coefficient',' ',2,dryness)
-         if (tracer) then
-            call writediagfi(ngrid,'dqevap','evaporated water vapor specific concentration','s-1',3,pdqevap)
-         endif
-      endif
 
 !      if(lastcall)then

trunk/LMDZ.TITAN/libf/phytitan/vdifc.F

@@ -1 @@
-      subroutine vdifc(ngrid,nlay,nq,rnat,ppopsk,         
+      subroutine vdifc(ngrid,nlay,nq,ppopsk,         
      &     ptimestep,pcapcal,lecrit,                        
      &     pplay,pplev,pzlay,pzlev,pz0,
@@ -7 +7 @@
      &     pdqdif,pdqsdif,lastcall)
 
-      use watercommon_h, only : RLVTT, T_h2O_ice_liq, RCPD, mx_eau_sol
       use radcommon_h, only : sigma
       USE surfdat_h
       USE tracer_h
       use comcstfi_mod, only: g, r, cpp, rcp
-      use callkeys_mod, only: water,tracer,nosurf
+      use callkeys_mod, only: tracer,nosurf
 
       implicit none
@@ -53 +52 @@
       REAL pdtsrf(ngrid),sensibFlux(ngrid),pcapcal(ngrid)
       REAL pq2(ngrid,nlay+1)
-      
-      real rnat(ngrid)      
+           
 
 !     Arguments added for condensation
@@ -116 +114 @@
       REAL kmixmin
 
-!     Variables added for implicit latent heat inclusion
-!     --------------------------------------------------
-      real latconst, dqsat(ngrid), qsat_temp1, qsat_temp2
-      real z1_Tdry(ngrid), z2_Tdry(ngrid)
-      real z1_T(ngrid), z2_T(ngrid)
-      real zb_T(ngrid)
-      real zc_T(ngrid,nlay)
-      real zd_T(ngrid,nlay)
-      real lat1(ngrid), lat2(ngrid)
-      real surfh2otot
-      logical surffluxdiag
-      integer isub ! sub-loop for precision
-
-      integer ivap, iice ! also make liq for clarity on surface...
-      save ivap, iice
-!$OMP THREADPRIVATE(ivap,iice)
-
       real, parameter :: karman=0.4
       real cd0, roughratio
 
-      logical forceWC
       real masse, Wtot, Wdiff
 
       real dqsdif_total(ngrid) 
       real zq0(ngrid) 
-
-      forceWC=.true.
-!      forceWC=.false.
 
 
@@ -155 +132 @@
 !     PRINT*,'In vdifc: Tcond(P=1mb)=',tcond1mb,' Lcond=',latcond
 !     PRINT*,'          acond,bcond',acond,bcond
-
-         if(water)then
-!                iliq=igcm_h2o_vap
-                ivap=igcm_h2o_vap
-                iice=igcm_h2o_ice ! simply to make the code legible               
-                                  ! to be generalised later
-         endif
 
          firstcall=.false.
@@ -404 +374 @@
 !       flux (if any) from subsurface
 
-      if(.not.water) then
 
          DO ig=1,ngrid
@@ -427 +396 @@
          ENDDO
 
-      endif                     ! not water
 
 !-----------------------------------------------------------------------
@@ -446 +414 @@
 !     -----------------------
          do iq=1,nq 
-
-            if (iq.ne.igcm_h2o_vap) then
 
                DO ig=1,ngrid
@@ -465 +431 @@
                ENDDO
 
-               if ((water).and.(iq.eq.iice)) then
-                  ! special case for water ice tracer: do not include
-                  ! h2o ice tracer from surface (which is set when handling
-                  ! h2o vapour case (see further down).
-                  ! zb(ig,1)=0 if iq ne ivap
-                  DO ig=1,ngrid
-                     z1(ig)=1./(za(ig,1)+
-     &                    zb(ig,2)*(1.-zdq(ig,2)))
-                     zcq(ig,1)=(za(ig,1)*zq(ig,1,iq)+
-     &                    zb(ig,2)*zcq(ig,2))*z1(ig)
-                  ENDDO
-               else             ! general case
+
+
                   DO ig=1,ngrid
                      z1(ig)=1./(za(ig,1)+
@@ -487 +443 @@
                           ! so last line is superfluous
                   enddo
-               endif            ! of if (water.and.(iq.eq.igcm_h2o_ice))
 
 
@@ -502 +457 @@
                end do
 
-            endif               ! if (iq.ne.igcm_h2o_vap)
-
-!     Calculate temperature tendency including latent heat term
-!     and assuming an infinite source of water on the ground
-!     ------------------------------------------------------------------
-
-            if (water.and.(iq.eq.igcm_h2o_vap)) then 
-            
-               ! compute evaporation efficiency
-               do ig = 1, ngrid
-                  if(nint(rnat(ig)).eq.1)then
-                     dryness(ig)=pqsurf(ig,ivap)+pqsurf(ig,iice) 
-                     dryness(ig)=MIN(1.,2*dryness(ig)/mx_eau_sol)
-                     dryness(ig)=MAX(0.,dryness(ig))
-                  endif
-               enddo
-
-               do ig=1,ngrid
-
-                ! Calculate the value of qsat at the surface (water)
-                call watersat(ptsrf(ig),pplev(ig,1),qsat(ig))
-                call watersat(ptsrf(ig)-0.0001,pplev(ig,1),qsat_temp1)
-                call watersat(ptsrf(ig)+0.0001,pplev(ig,1),qsat_temp2)
-                dqsat(ig)=(qsat_temp2-qsat_temp1)/0.0002
-                ! calculate dQsat / dT by finite differences
-                ! we cannot use the updated temperature value yet...
- 
-               enddo
-
-! coefficients for q
-
-               do ig=1,ngrid
-                  z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
-                  zcq(ig,nlay)=za(ig,nlay)*zq(ig,nlay,iq)*z1(ig)
-                  zdq(ig,nlay)=zb(ig,nlay)*z1(ig)
-               enddo 
-            
-               do ilay=nlay-1,2,-1
-                  do ig=1,ngrid
-                     z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+
-     $                    zb(ig,ilay+1)*(1.-zdq(ig,ilay+1)))
-                     zcq(ig,ilay)=(za(ig,ilay)*zq(ig,ilay,iq)+
-     $                    zb(ig,ilay+1)*zcq(ig,ilay+1))*z1(ig)
-                     zdq(ig,ilay)=zb(ig,ilay)*z1(ig)
-                  enddo
-               enddo
-
-               do ig=1,ngrid
-                  z1(ig)=1./(za(ig,1)+zb(ig,1)*dryness(ig)+
-     $                 zb(ig,2)*(1.-zdq(ig,2)))
-                  zcq(ig,1)=(za(ig,1)*zq(ig,1,iq)+
-     $                 zb(ig,2)*zcq(ig,2))*z1(ig)
-                  zdq(ig,1)=dryness(ig)*zb(ig,1)*z1(ig)
-               enddo
-
-! calculation of h0 and h1
-               do ig=1,ngrid
-                  zdq0(ig) = dqsat(ig)
-                  zcq0(ig) = qsat(ig)-dqsat(ig)*ptsrf(ig)
-
-                  z1(ig) = pcapcal(ig)*ptsrf(ig) +cpp*zb(ig,1)*zc(ig,1)
-     &                 + zdplanck(ig)*ptsrf(ig) + pfluxsrf(ig)*ptimestep 
-     &                 +  zb(ig,1)*dryness(ig)*RLVTT*
-     &                 ((zdq(ig,1)-1.0)*zcq0(ig)+zcq(ig,1))
-
-                  z2(ig) = pcapcal(ig) + cpp*zb(ig,1)*(1.-zd(ig,1))
-     &                 +zdplanck(ig)
-     &                 +zb(ig,1)*dryness(ig)*RLVTT*zdq0(ig)*
-     &                 (1.0-zdq(ig,1))
-
-                  ztsrf2(ig) = z1(ig) / z2(ig)
-                  pdtsrf(ig) = (ztsrf2(ig) - ptsrf(ig))/ptimestep
-                  zh(ig,1)   = zc(ig,1) + zd(ig,1)*ztsrf2(ig)
-               enddo
-
-! calculation of qs and q1
-               do ig=1,ngrid
-                  zq0(ig)     = zcq0(ig)+zdq0(ig)*ztsrf2(ig)
-                  zq(ig,1,iq) = zcq(ig,1)+zdq(ig,1)*zq0(ig)
-               enddo
-
-! calculation of evaporation              
-               do ig=1,ngrid 
-                  evap(ig)= zb(ig,1)*dryness(ig)*(zq(ig,1,ivap)-zq0(ig))
-                  dqsdif_total(ig)=evap(ig)
-               enddo
-
-! recalculate temperature and q(vap) to top of atmosphere, starting from ground
-               do ilay=2,nlay
-                  do ig=1,ngrid
-                     zq(ig,ilay,iq)=zcq(ig,ilay)
-     &                    +zdq(ig,ilay)*zq(ig,ilay-1,iq)
-                     zh(ig,ilay)=zc(ig,ilay)+zd(ig,ilay)*zh(ig,ilay-1)
-                  end do
-               end do
-
-               do ig=1,ngrid
-
-!     --------------------------------------------------------------------------
-!     On the ocean, if T > 0 C then the vapour tendency must replace the ice one
-!     The surface vapour tracer is actually liquid. To make things difficult.
-
-                  if (nint(rnat(ig)).eq.0) then ! unfrozen ocean
-                     
-                     pdqsdif(ig,ivap)=dqsdif_total(ig)/ptimestep
-                     pdqsdif(ig,iice)=0.0
-
-
-                  elseif (nint(rnat(ig)).eq.1) then ! (continent)
-
-!     --------------------------------------------------------
-!     Now check if we've taken too much water from the surface
-!     This can only occur on the continent 
-
-!     If water is evaporating / subliming, we take it from ice before liquid
-!     -- is this valid??
-                     if(dqsdif_total(ig).lt.0)then
-                        pdqsdif(ig,iice)=dqsdif_total(ig)/ptimestep
-                        pdqsdif(ig,iice)=max(-pqsurf(ig,iice)/ptimestep
-     &                       ,pdqsdif(ig,iice))
-                     endif
-                     ! sublimation only greater than qsurf(ice)
-                     ! ----------------------------------------
-                     ! we just convert some liquid to vapour too
-                     ! if latent heats are the same, no big deal
-                     if (-dqsdif_total(ig).gt.pqsurf(ig,iice))then            
-                       pdqsdif(ig,iice) = -pqsurf(ig,iice)/ptimestep ! removes all the ice!
-                       pdqsdif(ig,ivap) = dqsdif_total(ig)/ptimestep
-     &                       - pdqsdif(ig,iice) ! take the remainder from the liquid instead
-                       pdqsdif(ig,ivap) = max(-pqsurf(ig,ivap)/ptimestep
-     &                       ,pdqsdif(ig,ivap))
-                    endif
-
-                 endif          ! if (rnat.ne.1)
-
-!     If water vapour is condensing, we must decide whether it forms ice or liquid.
-                 if(dqsdif_total(ig).gt.0)then ! a bug was here!
-                    if(ztsrf2(ig).gt.T_h2O_ice_liq)then
-                       pdqsdif(ig,iice)=0.0
-                       pdqsdif(ig,ivap)=dqsdif_total(ig)/ptimestep
-                    else
-                       pdqsdif(ig,iice)=dqsdif_total(ig)/ptimestep
-                       pdqsdif(ig,ivap)=0.0
-                    endif
-                 endif
-
-              end do            ! of DO ig=1,ngrid
-           endif                ! if (water et iq=ivap)
+
         end do                  ! of do iq=1,nq
       endif                     ! traceur
@@ -685 +493 @@
          enddo
 
-         if(water.and.forceWC)then ! force water conservation in model
-                                   ! we calculate the difference and add it to the ground
-                                   ! this is ugly and should be improved in the future
-            do ig=1,ngrid
-               Wtot=0.0
-               do ilay=1,nlay
-                  masse = (pplev(ig,ilay) - pplev(ig,ilay+1))/g
-!                  Wtot=Wtot+masse*(zq(ig,ilay,iice)-
-!     &                 (pq(ig,ilay,iice)+pdqfi(ig,ilay,iice)*ptimestep))
-                  Wtot=Wtot+masse*(zq(ig,ilay,ivap)-
-     &                 (pq(ig,ilay,ivap)+pdqfi(ig,ilay,ivap)*ptimestep))
-               enddo
-               Wdiff=Wtot/ptimestep+pdqsdif(ig,ivap)+pdqsdif(ig,iice)
-
-               if(ztsrf2(ig).gt.T_h2O_ice_liq)then
-                  pdqsdif(ig,ivap)=pdqsdif(ig,ivap)-Wdiff
-               else
-                  pdqsdif(ig,iice)=pdqsdif(ig,iice)-Wdiff
-               endif
-            enddo
-
-         endif
-
       endif 
-
-      if(water)then
-      call writediagfi(ngrid,'beta','Dryness coefficient',' ',2,dryness)
-      endif
 
 !      if(lastcall)then