Changeset 1297

Timestamp:

Jun 18, 2014, 4:33:56 AM (11 years ago)

Author:

bclmd

Message:

LMDZ.GENERIC: slab ocean added

Location:

trunk/LMDZ.GENERIC/libf/phystd

Files:

• callkeys.h (modified) (2 diffs)
• comsoil_h.F90 (modified) (1 diff)
• hydrol.F90 (modified) (15 diffs)
• inifis.F (modified) (1 diff)
• initracer.F (modified) (1 diff)
• iostart.F90 (modified) (7 diffs)
• lect_start_archive.F (modified) (11 diffs)
• newstart.F (modified) (13 diffs)
• phyetat0.F90 (modified) (6 diffs)
• phyredem.F90 (modified) (4 diffs)
• physiq.F90 (modified) (28 diffs)
• soil_settings.F (modified) (10 diffs)
• start2archive.F (modified) (7 diffs)
• sugas_corrk.F90 (modified) (2 diffs)
• turbdiff.F90 (modified) (6 diffs)
• vdifc.F (modified) (4 diffs)

trunk/LMDZ.GENERIC/libf/phystd/callkeys.h

@@ -26 +26 @@
      &   , hydrology, sourceevol                                        &
      &   , CLFvarying                                                   &
-     &   , strictboundcorrk                                                                                     
+     &   , strictboundcorrk                                             &                                        
+     &   , ok_slab_ocean                                                &
+     &   , ok_slab_sic                                                  &
+     &   , ok_slab_heat_transp                                          
+
 
       COMMON/callkeys_i/iaervar,iddist,iradia,startype
@@ -73 +77 @@
       logical nosurf
       logical oblate
+      logical ok_slab_ocean
+      logical ok_slab_sic
+      logical ok_slab_heat_transp
 
       integer iddist

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

@@ -5 +5 @@
 !integer, parameter :: nsoilmx = 18 ! for z1=0.0002 m, depth = 18 m => mars case 
 !integer, parameter :: nsoilmx = 13 ! for z1=0.03 m, depth = 104.8 m => earth case
-  integer, parameter :: nsoilmx = 18 
+  integer, parameter :: nsoilmx = 18
 
   real,save,allocatable,dimension(:) :: layer      ! soil layer depths

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

@@ -1 +1 @@
 subroutine hydrol(ngrid,nq,ptimestep,rnat,tsurf,          &
      qsurf,dqsurf,dqs_hyd,pcapcal,                &
-     albedo0,albedo,mu0,pdtsurf,pdtsurf_hyd,hice)
+     albedo0,albedo,mu0,pdtsurf,pdtsurf_hyd,hice, &
+     pctsrf_sic,sea_ice)
 
   use ioipsl_getincom 
@@ -9 +10 @@
   USE comgeomfi_h
   USE tracer_h
+  use slab_ice_h
 
   implicit none
@@ -46 +48 @@
 !     Inputs
 !     ------
-      real albedoocean
-      parameter (albedoocean=0.07)
       real albedoice
       save albedoice
@@ -62 +62 @@
 !     Arguments
 !     ---------
-      integer rnat(ngrid) ! I changed this to integer (RW)
+      real rnat(ngrid) ! I changed this to integer (RW)
       real,dimension(:),allocatable,save :: runoff
       real totalrunoff, tsea, oceanarea
@@ -73 +73 @@
       real hice(ngrid)
       real albedo0(ngrid), albedo(ngrid)
+      real pctsrf_sic(ngrid), sea_ice(ngrid)
 
       real oceanarea2
@@ -91 +92 @@
       real zqsurf(ngrid,nq)
       real ztsurf(ngrid)
+      real albedo_sic, alb_ice
+      real zfra
 
       integer, save :: ivap, iliq, iice
@@ -134 +137 @@
          oceanarea=0.
          do ig=1,ngrid
-            if(rnat(ig).eq.0)then
+            if(nint(rnat(ig)).eq.0)then
                oceanarea=oceanarea+area(ig)
             endif
@@ -175 +178 @@
 !     Ocean
 !     -----
-         if(rnat(ig).eq.0)then
+         if(nint(rnat(ig)).eq.0)then
 
 !     re-calculate oceanic albedo
-            if(diurnal.and.oceanalbvary)then
-               fauxo      = ( 1.47 - ACOS( mu0(ig) ) )/0.15 ! where does this come from (Benjamin)?
-               albedo(ig) = 1.1*( .03 + .630/( 1. + fauxo*fauxo))
-               albedo(ig) = MAX(MIN(albedo(ig),0.60),0.04)
-            else
-               albedo(ig) = albedoocean ! modif Benjamin
-            end if
+!            if(diurnal.and.oceanalbvary)then
+!               fauxo      = ( 1.47 - ACOS( mu0(ig) ) )/0.15 ! where does this come from (Benjamin)?
+!               albedo(ig) = 1.1*( .03 + .630/( 1. + fauxo*fauxo))
+!               albedo(ig) = MAX(MIN(albedo(ig),0.60),0.04)
+!            else
+               albedo(ig) = alb_ocean ! modif Benjamin
+!            end if
+
+
+          if(ok_slab_ocean) then
+
+! Fraction neige (hauteur critique 45kg/m2~15cm)
+            zfra = MAX(0.0,MIN(1.0,zqsurf(ig,iice)/45.0))
+! Albedo glace
+            alb_ice=alb_ice_max-(alb_ice_max-alb_ice_min) &
+            *exp(-sea_ice(ig)/h_alb_ice)
+! Albedo glace+neige
+            albedo_sic= albedosnow*zfra + alb_ice*(1.0-zfra)
+
+! Albedo final
+            albedo(ig) = pctsrf_sic(ig)*albedo_sic + (1.-pctsrf_sic(ig))*alb_ocean
+!     oceanic ice height, just for diagnostics
+            hice(ig)    = MIN(10.,sea_ice(ig)/rhowater)
+          else !ok_slab_ocean
+
 
 !     calculate oceanic ice height including the latent heat of ice formation
@@ -214 +235 @@
 
                pdtsurf_hyd(ig) = -hicebis(ig)*RLFTT*rhowater/pcapcal(ig)/ptimestep                
-               albedo(ig)      = albedoocean
+               albedo(ig)      = alb_ocean
 
             endif
@@ -221 +242 @@
             zqsurf(ig,iice) = hice(ig)*rhowater
 
+          endif!(ok_slab_ocean)
+
 
 !     Continent
 !     ---------
-         elseif (rnat(ig).eq.1) then
+         elseif (nint(rnat(ig)).eq.1) then
 
 !     melt the snow
@@ -298 +321 @@
          oceanarea2=0.       
          DO ig=1,ngrid
-            if((rnat(ig).eq.0).and.(hice(ig).eq.0.))then
+            if((nint(rnat(ig)).eq.0).and.(hice(ig).eq.0.))then
                oceanarea2=oceanarea2+area(ig)*pcapcal(ig)
             end if
@@ -305 +328 @@
          tsea=0.
          DO ig=1,ngrid
-            if((rnat(ig).eq.0).and.(hice(ig).eq.0.))then       
+            if((nint(rnat(ig)).eq.0).and.(hice(ig).eq.0.))then       
                tsea=tsea+ztsurf(ig)*area(ig)*pcapcal(ig)/oceanarea2
             end if
@@ -311 +334 @@
 
          DO ig=1,ngrid
-            if((rnat(ig).eq.0).and.(hice(ig).eq.0))then
+            if((nint(rnat(ig)).eq.0).and.(hice(ig).eq.0))then
                pdtsurf_hyd(ig) = pdtsurf_hyd(ig) + (tsea-ztsurf(ig))/oceantime
             end if
@@ -326 +349 @@
          totalrunoff=0.
          do ig=1,ngrid
-            if (rnat(ig).eq.1) then
+            if (nint(rnat(ig)).eq.1) then
                totalrunoff = totalrunoff + area(ig)*runoff(ig)
             endif
@@ -332 +355 @@
          
          do ig=1,ngrid
-            if (rnat(ig).eq.0) then
+            if (nint(rnat(ig)).eq.0) then
                zqsurf(ig,iliq) = zqsurf(ig,iliq) + &
                     totalrunoff/oceanarea

trunk/LMDZ.GENERIC/libf/phystd/inifis.F

@@ -329 +329 @@
          call getin("graybody",graybody)
          write(*,*)" graybody = ",graybody
+
+! Slab Ocean 
+         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
+
+         write(*,*) "Use slab-sea-ice ?"
+         ok_slab_sic=.true.         ! default value
+         call getin("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("ok_slab_heat_transp",ok_slab_heat_transp)
+         write(*,*) "ok_slab_heat_transp = ",ok_slab_heat_transp
+
+
 
 ! Test of incompatibility:

trunk/LMDZ.GENERIC/libf/phystd/initracer.F

@@ -172 +172 @@
           write(*,*)'      ',iq,' ',trim(noms(iq))
         enddo
-        stop
+!        stop
       else
         write(*,*) "initracer: found all expected tracers, namely:"

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

@@ -14 +14 @@
     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)
     INTEGER,PARAMETER :: length=100 ! size of tab_cntrl array
@@ -466 +467 @@
   USE infotrac, only: nqtot
   USE comsoil_h, only: nsoilmx
+  USE slab_ice_h, only: noceanmx
+
   IMPLICIT NONE
     CHARACTER(LEN=*),INTENT(IN) :: filename
@@ -552 +555 @@
       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
@@ -632 +643 @@
   USE mod_grid_phy_lmdz
   USE mod_phys_lmdz_para
+  USE slab_ice_h, only: noceanmx
   IMPLICIT NONE
   CHARACTER(LEN=*),INTENT(IN)    :: field_name
@@ -817 +829 @@
         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)
@@ -889 +945 @@
   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
@@ -939 +996 @@
         ! 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.GENERIC/libf/phystd/lect_start_archive.F

@@ -1 +1 @@
       SUBROUTINE lect_start_archive(date,tsurf,tsoil,emis,q2,
      &     t,ucov,vcov,ps,h,phisold_newgrid,
-     &     q,qsurf,surfith,nid)
+     &     q,qsurf,surfith,nid,
+     &     rnat,pctsrf_sic,tslab,tsea_ice,sea_ice)
 
 !      USE surfdat_h
@@ -7 +8 @@
       USE tracer_h, ONLY: igcm_co2_ice
       USE infotrac, ONLY: tname, nqtot
+      USE slab_ice_h, ONLY: noceanmx
 !      USE control_mod
+! to use  'getin'
 
 c=======================================================================
@@ -39 +42 @@
 #include "netcdf.inc"
 !#include"advtrac.h"
+#include "callkeys.h"
 c=======================================================================
 c   Declarations
@@ -100 +104 @@
       REAL emis(ngridmx)
       REAL q2(ngridmx,nlayermx+1),qsurf(ngridmx,nqtot)
+      REAL tslab(ngridmx,noceanmx)
+      REAL rnat(ngridmx),pctsrf_sic(ngridmx)
+      REAL tsea_ice(ngridmx),sea_ice(ngridmx)
 c     REAL phisfi(ngridmx)
 
@@ -122 +129 @@
       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
@@ -150 +160 @@
       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)
@@ -299 +313 @@
       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))
@@ -700 +719 @@
       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:"
      &           ,rlonuold," rlatvold:",rlatvold
@@ -971 +1067 @@
       call gr_dyn_fi (1,iim+1,jjm+1,ngridmx,emiss,emis)
 c     write(46,*) 'emis',emis
+
+
+
 c-----------------------------------------------------------------------
 c 6.1.2 Traitement special de la pression au sol :
@@ -1190 +1289 @@
        call gr_dyn_fi (nsoilmx,iim+1,jjm+1,ngridmx,tsoilS,tsoil)
 
-
-
-c-----------------------------------------------------------------------
-c 6.3 Variable 3d :
+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 (1,iim+1,jjm+1,ngridmx,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,ngridmx,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,ngridmx,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,ngridmx,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,ngridmx,sea_ices,sea_ice)
+
+c-----------------------------------------------------------------------
+c 6.4 Variable 3d :
 c-----------------------------------------------------------------------
       
@@ -1362 +1482 @@
       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.GENERIC/libf/phystd/newstart.F

@@ -28 +28 @@
       use iostart, only: open_startphy
       use comgeomphy, only: initcomgeomphy
+      use slab_ice_h, only:noceanmx
       implicit none
 
@@ -108 +109 @@
 
       integer ierr 
+
+      REAL rnat(ngridmx)
+      REAL tslab(ngridmx,nsoilmx) ! 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)
 
 c Variables on the new grid along scalar points 
@@ -167 +174 @@
       REAL totalfrac(ngridmx)
 
+
 !     added by RW for nuketharsis
       real fact1
@@ -186 +194 @@
       CALL init_phys_lmdz(iim,jjp1,llm,1,(/(jjm-1)*iim+2/))
       call initcomgeomphy
-
 ! Load tracer number and names:
       call iniadvtrac(nqtot,numvanle)
@@ -196 +203 @@
 c   Choice of the start file(s) to use
 c=======================================================================
-
       write(*,*) 'From which kind of files do you want to create new',
      .  'start and startfi files'
@@ -288 +294 @@
 c Lecture du tableau des parametres du run (pour la dynamique)
 c-----------------------------------------------------------------------
-
       if (choix_1.eq.0) then
 
@@ -339 +344 @@
      .        day_ini,time,
      .        tsurf,tsoil,emis,q2,qsurf,   !) ! temporary modif by RDW
-     .        cloudfrac,totalfrac,hice)
+     .        cloudfrac,totalfrac,hice,rnat,pctsrf_sic,tslab,tsea_ice,
+     .        sea_ice)
 
         ! copy albedo and soil thermal inertia on (local) physics grid
@@ -529 +535 @@
         write(*,*) 'Reading file START_ARCHIVE'
         CALL lect_start_archive(date,tsurf,tsoil,emis,q2,
-     .   t,ucov,vcov,ps,teta,phisold_newgrid,q,qsurf,
-     &   surfith,nid)
+     &   t,ucov,vcov,ps,teta,phisold_newgrid,q,qsurf,
+     &   surfith,nid,
+     &   rnat,pctsrf_sic,tslab,tsea_ice,sea_ice)
         write(*,*) "OK, read start_archive file"
         ! copy soil thermal inertia
@@ -589 +596 @@
       write(*,*) 'therm_ini_s : Set soil thermal inertia to reference sur
      &face values'
+      write(*,*) 'slab_ocean_0 : initialisation of slab 
+     $ocean variables'
 
         write(*,*)
@@ -1373 +1382 @@
 !       ----------------------------------------------------------------------
 
-        else if (trim(modif).eq.'therm_ini_s') then
+        else if (trim(modif) .eq. 'therm_ini_s') then
 !          write(*,*)"surfithfi(1):",surfithfi(1)
           do isoil=1,nsoilmx
@@ -1394 +1403 @@
 
 
+
+c       slab_ocean_initialisation
+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)
+
+
+
         else
           write(*,*) '       Unknown (misspelled?) option!!!'
         end if ! of if (trim(modif) .eq. '...') elseif ...
-        
-
-
-
-
-
-
-
 
 
@@ -1433 +1461 @@
 !         ENDIF
 !      ENDDO
+
+
+
 
 c=======================================================================
@@ -1540 +1571 @@
      &                dtphys,real(day_ini),
      &                tsurf,tsoil,emis,q2,qsurf,
-     &                cloudfrac,totalfrac,hice)
+     &                cloudfrac,totalfrac,hice,
+     &                rnat,pctsrf_sic,tslab,tsea_ice,sea_ice)
 
 c=======================================================================

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

@@ -1 +1 @@
 subroutine phyetat0 (ngrid,fichnom,tab0,Lmodif,nsoil,nq, &
                      day_ini,time,tsurf,tsoil, &
-                     emis,q2,qsurf,cloudfrac,totcloudfrac,hice)
+                     emis,q2,qsurf,cloudfrac,totcloudfrac,hice, &
+                     rnat,pctsrf_sic,tslab,tsea_ice,sea_ice)
+
 
   USE infotrac, ONLY: tname
@@ -8 +10 @@
                      get_field, get_var, inquire_field, &
                      inquire_dimension, inquire_dimension_length
+  use slab_ice_h, only: noceanmx
 
   implicit none
@@ -47 +50 @@
   real,intent(out) :: cloudfrac(ngrid,nlayermx)
   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)
 
 !======================================================================
@@ -89 +95 @@
 IF (.not. ALLOCATED(zthe)) ALLOCATE(zthe(ngrid))
 
+
 ! open physics initial state file:
 call open_startphy(fichnom)
@@ -270 +277 @@
 if (.not.found) then
   write(*,*) "phyetat0: Failed loading <hice>"
-  call abort
+!  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
@@ -308 +388 @@
 endif ! of if (nq.ge.1)
 
+
 ! Call to soil_settings, in order to read soil temperatures,
 ! as well as thermal inertia and volumetric heat capacity
-
 call soil_settings(nid_start,ngrid,nsoil,tsurf,tsoil,indextime)
-
 !
 ! close file:

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

@@ -135 +135 @@
 subroutine physdem1(filename,nsoil,ngrid,nlay,nq, &
                     phystep,time,tsurf,tsoil,emis,q2,qsurf, &
-                    cloudfrac,totcloudfrac,hice)
+                    cloudfrac,totcloudfrac,hice, &
+                    rnat,pctsrf_sic,tslab,tsea_ice,sea_ice)
   ! write time-dependent variable to restart file
   use iostart, only : open_restartphy, close_restartphy, & 
                       put_var, put_field
   use infotrac, only: tname
+  use slab_ice_h, only: noceanmx
+
   implicit none
 !======================================================================
@@ -145 +148 @@
 #include "comcstfi.h"
 #include "planete.h"
+#include "callkeys.h"
 !======================================================================
   character(len=*),intent(in) :: filename
@@ -161 +165 @@
   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
@@ -187 +196 @@
   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

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

@@ -26 +26 @@
       use control_mod, only: ecritphy, iphysiq, nday
       use phyredem, only: physdem0, physdem1
+      use slab_ice_h
+      use ocean_slab_mod, only :ocean_slab_init, ocean_slab_ice, &
+                                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
+
 
       implicit none
@@ -184 +189 @@
 
       real,dimension(:),allocatable,save :: albedo0 ! Surface albedo
-      integer,dimension(:),allocatable,save :: rnat ! added by BC
+      real,dimension(:),allocatable,save :: rnat ! added by BC
 
       real,dimension(:),allocatable,save :: emis ! Thermal IR surface emissivity
@@ -276 +281 @@
       real zdh(ngrid,nlayermx)
       real gmplanet
+      real taux(ngrid),tauy(ngrid)
 
       integer length
@@ -431 +437 @@
       real, dimension(:), allocatable :: tmp_fract ! day fraction of the time interval 
       real, dimension(:), allocatable :: tmp_mu0 ! equivalent solar angle
+
+!     included by BC for slab ocean
+      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 
+
+      real :: tsurf2(ngrid)
+      real :: flux_o(ngrid),flux_g(ngrid),fluxgrdocean(ngrid)
+      real :: dt_ekman(ngrid,noceanmx),dt_hdiff(ngrid,noceanmx)
+      real :: flux_sens_lat(ngrid)
+      real :: qsurfint(ngrid,nq)
+
+
 !=======================================================================
 
@@ -475 +497 @@
         ALLOCATE(zdtsw(ngrid,nlayermx))
         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))
+!        ALLOCATE(qsurfint(ngrid,nqmx))
+
 
         !! this is defined in comsaison_h
@@ -513 +543 @@
          call phyetat0(ngrid,"startfi.nc",0,0,nsoilmx,nq,   &
                day_ini,time_phys,tsurf,tsoil,emis,q2,qsurf,   &
-               cloudfrac,totcloudfrac,hice)
+               cloudfrac,totcloudfrac,hice,  &
+               rnat,pctsrf_sic,tslab, tsea_ice,sea_ice)
 
          if (pday.ne.day_ini) then
@@ -537 +568 @@
          endif
 
+!        Initialize oceanic variables
+!        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+         if (ok_slab_ocean)then
+
+           call ocean_slab_init(ngrid,ptimestep, tslab, &
+                sea_ice, pctsrf_sic)
+
+            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(zmasq)
+
+
+         endif
+
+
 !        initialize soil 
 !        ~~~~~~~~~~~~~~~
@@ -542 +598 @@
             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
+
+
+
+
          else
             print*,'WARNING! Thermal conduction in the soil turned off'
@@ -572 +646 @@
          endif
 
+!  In newstart now, will have to be remove (BC 2014)
 !        define surface as continent or ocean
 !        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-         rnat(:)=1
-         do ig=1,ngrid
-!            if(iceball.or.oceanball.or.(inertiedat(ig,1).gt.1.E4))then
-            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'
-
+         if (.not.ok_slab_ocean) then
+           rnat(:)=1.
+           do ig=1,ngrid
+!             if(iceball.or.oceanball.or.(inertiedat(ig,1).gt.1.E4))then
+             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!(.not.ok_slab_ocean)
 
 !        initialise surface history variable 
@@ -644 +720 @@
 !     Initialize various variables
 !     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
+      
       pdu(1:ngrid,1:nlayermx) = 0.0
       pdv(1:ngrid,1:nlayermx) = 0.0
@@ -655 +731 @@
       zdtsurf(1:ngrid)      = 0.0
       dqsurf(1:ngrid,1:nq)= 0.0
+      flux_sens_lat(1:ngrid) = 0.0
+      taux(1:ngrid) = 0.0
+      tauy(1:ngrid) = 0.0
+
+
+
 
       zday=pday+ptime ! compute time, in sols (and fraction thereof)
@@ -674 +756 @@
         else if (flatten .eq. 0.0 .or. J2 .eq. 0.0 .or. Rmean .eq. 0.0 .or. MassPlanet .eq. 0.0) then
         print*,'I need values for flatten, J2, Rmean and MassPlanet to compute glat (else set oblate=.false.)'
+
         call abort
         else
@@ -817 +900 @@
       
 !             standard callcorrk
-              clearsky=.false.
-              call callcorrk(ngrid,nlayer,pq,nq,qsurf,                   &
+
+
+
+
+              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)
+                
+
+                clearsky=.false.
+                call callcorrk(ngrid,nlayer,pq,nq,qsurf,                   &
                       albedo,emis,mu0,pplev,pplay,pt,                    &
                       tsurf,fract,dist_star,aerosol,muvar,               &
@@ -827 +924 @@
 
 !             Option to call scheme once more for clear regions
-              if(CLFvarying)then
+                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,                  &
+                   clearsky=.true.
+                   call callcorrk(ngrid,nlayer,pq,nq,qsurf,                  &
                       albedo,emis,mu0,pplev,pplay,pt,                      &
                       tsurf,fract,dist_star,aerosol,muvar,                 &
@@ -839 +936 @@
                       tau_col1,cloudfrac,totcloudfrac,                     &
                       clearsky,firstcall,lastcall)
-                 clearsky = .false.  ! just in case.     
+                   clearsky = .false.  ! just in case.     
+
 
                  ! Sum the fluxes and heating rates from cloudy/clear cases
-                 do ig=1,ngrid
-                    tf=totcloudfrac(ig)
-                    ntf=1.-tf
+                   do ig=1,ngrid
+                      tf=totcloudfrac(ig)
+                      ntf=1.-tf
                     
                     fluxsurf_lw(ig) = ntf*fluxsurf_lw1(ig) + tf*fluxsurf_lw(ig)
@@ -861 +959 @@
 
               endif !CLFvarying
+
+              if(ok_slab_ocean) then
+                tsurf(:)=tsurf2(:)
+              endif!(ok_slab_ocean)
+
 
 !             Radiative flux from the sky absorbed by the surface (W.m-2)
@@ -958 +1061 @@
               zdum1,zdum2,pdt,pdq,zflubid,           &
               zdudif,zdvdif,zdtdif,zdtsdif,          &
-              sensibFlux,q2,zdqdif,zdqevap,zdqsdif,lastcall)
+              sensibFlux,q2,zdqdif,zdqevap,zdqsdif,  &
+              taux,tauy,lastcall)
 
          else
@@ -970 +1074 @@
               zdum1,zdum2,zdh,pdq,zflubid,           &
               zdudif,zdvdif,zdhdif,zdtsdif,          &
-              sensibFlux,q2,zdqdif,zdqsdif,lastcall)
+              sensibFlux,q2,zdqdif,zdqsdif,          &
+              taux,tauy,lastcall)
 
            zdtdif(1:ngrid,1:nlayermx)=zdhdif(1:ngrid,1:nlayermx)*zpopsk(1:ngrid,1:nlayermx) ! for diagnostic only
@@ -981 +1086 @@
          pdt(1:ngrid,1:nlayermx)=pdt(1:ngrid,1:nlayermx)+zdtdif(1:ngrid,1:nlayermx)
          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 
            pdq(1:ngrid,1:nlayermx,1:nq)=pdq(1:ngrid,1:nlayermx,1:nq)+ zdqdif(1:ngrid,1:nlayermx,1:nq)
@@ -1403 +1515 @@
 
 
+!   7e. Ocean
+!     ---------
+
+!       ---------------------------------
+!       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! (ok_slab_ocean)
 
 !       ---------------------------------
@@ -1409 +1559 @@
 
          if(hydrology)then
-
+            
             call hydrol(ngrid,nq,ptimestep,rnat,tsurf,qsurf,dqsurf,dqs_hyd,    &
-               capcal,albedo0,albedo,mu0,zdtsurf,zdtsurf_hyd,hice)
+               capcal,albedo0,albedo,mu0,zdtsurf,zdtsurf_hyd,hice,pctsrf_sic,  &
+               sea_ice)
             ! note: for now, also changes albedo in the subroutine
 
@@ -1462 +1613 @@
 
 !     Increment surface temperature
-      tsurf(1:ngrid)=tsurf(1:ngrid)+ptimestep*zdtsurf(1:ngrid) 
+
+      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!(ok_slab_ocean)
 
 !     Compute soil temperatures and subsurface heat flux
       if (callsoil) then
          call soil(ngrid,nsoilmx,.false.,lastcall,inertiedat,   &
-                ptimestep,tsurf,tsoil,capcal,fluxgrd)
+                ptimestep,tsurf,tsoil,capcal,fluxgrd)     
       endif
+
+
+      if (ok_slab_ocean) then
+            do ig=1,ngrid
+               fluxgrdocean(ig)=fluxgrd(ig)
+               if (nint(rnat(ig)).eq.0) then
+               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 !(ok_slab_ocean)
 
 !-------------------------
@@ -1744 +1928 @@
                       ptimestep,ztime_fin, &
                       tsurf,tsoil,emis,q2,qsurf_hist, &
-                      cloudfrac,totcloudfrac,hice)
+                      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
+
          
          endif ! of if (lastcall)
@@ -1778 +1968 @@
             call wstats(ngrid,"ASR","absorbed stellar rad.","W m-2",2,fluxabs_sw)
             call wstats(ngrid,"OLR","outgoing longwave rad.","W m-2",2,fluxtop_lw)
-
+            call wstats(ngrid,"ALB","Surface albedo"," ",2,albedo)
+            call wstats(ngrid,"p","Pressure","Pa",3,pplay)
             call wstats(ngrid,"temp","Atmospheric temperature","K",3,zt)
             call wstats(ngrid,"u","Zonal (East-West) wind","m.s-1",3,zu)
@@ -1806 +1997 @@
            endif !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! (ok_slab_ocean)
+
             if(lastcall) then
               write (*,*) "Writing stats..."
@@ -1833 +2047 @@
 !        call writediagsoil(ngrid,"tsurf","Surface temperature","K",2,tsurf)
 !        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! (ok_slab_ocean)
 
 !     Total energy balance diagnostics
@@ -1847 +2077 @@
 !           call writediagfi(ngrid,"fluxsurfswcs","sw surface flux (cs).","W m-2",2,fluxsurf_sw1)
 !           call writediagfi(ngrid,"fluxsurflwcs","lw back radiation (cs).","W m-2",2,fluxsurf_lw1)
-           call writediagfi(ngrid,"GND","heat flux from ground","W m-2",2,fluxgrd)
+           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!(ok_slab_ocean)
            call writediagfi(ngrid,"DYN","dynamical heat input","W m-2",2,fluxdyn)
          endif
@@ -1917 +2151 @@
           endif
 
-          if(hydrology)then
+          if((hydrology).and.(.not.ok_slab_ocean))then
              call writediagfi(ngrid,"hice","oceanic ice height","m",2,hice)
           endif

trunk/LMDZ.GENERIC/libf/phystd/soil_settings.F

@@ -65 +65 @@
       real,dimension(:,:),allocatable :: oldinertiedat
       real,dimension(:,:),allocatable :: oldtsoil
-      
       ! for interpolation
       real,dimension(:),allocatable :: oldgrid
@@ -77 +76 @@
 
       logical :: found,ok
-      
-!======================================================================
-
+!======================================================================
 ! 1. Depth coordinate
 ! -------------------
-
 ! 1.1 Start by reading how many layers of soil there are
-
         dimlen=inquire_dimension_length("subsurface_layers")
 
@@ -94 +89 @@
           interpol=.true.
         endif
-
 ! 1.2 Find out the # of dimensions <inertiedat> was defined as using
-
       ndims=inquire_field_ndims("inertiedat")
-
 ! 1.3 Read depths values or set olddepthdef flag and values
       if (ndims.eq.1) then ! we know that there is none
@@ -125 +117 @@
         else ! put values in mlayer
           call get_var("soildepth",mlayer,found)
-          if (.not.found) then
+          print*,"mlayer",mlayer
+          if (.not.found) then
             write(*,*)'soil_settings: Problem while reading <soildepth>'
           endif
         endif !of if (interpol)
       endif !of if (ndims.eq.1)
-
 ! 1.4 Build mlayer(), if necessary
       if (interpol) then
@@ -139 +131 @@
         do iloop=0,nsoil-1
           mlayer(iloop)=lay1*(alpha**(iloop-0.5))
-        enddo
+        enddo
       endif
-
 ! 1.5 Build layer(); following the same law as mlayer()
       ! Assuming layer distribution follows mid-layer law:
@@ -150 +141 @@
         layer(iloop)=lay1*(alpha**(iloop-1))
       enddo
-
 ! 2. Volumetric heat capacity (note: it is declared in comsoil_h)
 ! ---------------------------
@@ -168 +158 @@
 ! 3. Thermal inertia (note: it is declared in comsoil_h)
 ! ------------------
-
 ! 3.1 Look (again) for thermal inertia data (to reset nvarid)
 
@@ -219 +208 @@
 ! 4. Read soil temperatures
 ! -------------------------
-
 !      ierr=nf90_inq_varid(nid,"tsoil",nvarid)
       ok=inquire_field("tsoil")
@@ -255 +243 @@
 ! 5. If necessary, interpolate soil temperatures and thermal inertias
 ! -------------------------------------------------------------------
-
       if (olddepthdef) then
       ! tsoil needs to be interpolated, but not therm_i
@@ -339 +326 @@
       xmax = MAXVAL(tsoil)
       write(*,*)'Soil temperature <tsoil>:',xmin,xmax
-
       end

trunk/LMDZ.GENERIC/libf/phystd/start2archive.F

@@ -24 +24 @@
 !      use control_mod
       use comgeomphy, only: initcomgeomphy
+      use slab_ice_h, only: noceanmx
+! to use  'getin'
+      USE ioipsl_getincom
+
       implicit none
 
@@ -41 +45 @@
 !#include"advtrac.h"
 #include "netcdf.inc"
-
+#include "callkeys.h"
 c-----------------------------------------------------------------------
 c   Declarations
@@ -78 +82 @@
       REAL cloudfrac(ngridmx,nlayermx),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
@@ -93 +101 @@
       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)
+
 
 c Variables intermediaires : vent naturel, mais pas coord scalaire
@@ -185 +198 @@
       Lmodif=0
 
+
       CALL phyetat0 (ngridmx,fichnom,0,Lmodif,nsoilmx,nqtot,day_ini_fi,
      .      timefi,
      .      tsurf,tsoil,emis,q2,qsurf,
 !       change FF 05/2011
-     .       cloudfrac,totalcloudfrac,hice)
+     .       cloudfrac,totalcloudfrac,hice,
+!       change BC 05/2014
+     .       rnat,pctsrf_sic,tslab,tsea_ice,sea_ice)
+
+
 
 
@@ -303 +321 @@
       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=======================================================================
 c Info pour controler
@@ -455 +479 @@
       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','',3,tslabS)
+      call write_archive(nid,ntime,'tsea_ice'
+     &        ,'tsea_ice','',2,tsea_iceS)
+      endif !ok_slab_ocean
+      ENDIF
 c-----------------------------------------------------------------------
 c Fin 

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

@@ -607 +607 @@
 !=======================================================================
 !     Initialise the continuum absorption data
-
+      if(continuum)then
       do igas=1,ngasmx
 
@@ -644 +644 @@
 
       enddo
+      endif
 
       print*,'----------------------------------------------------'

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

@@ -5 +5 @@
           pdufi,pdvfi,pdtfi,pdqfi,pfluxsrf,            &
           Pdudif,pdvdif,pdtdif,pdtsrf,sensibFlux,pq2,  &
-          pdqdif,pdqevap,pdqsdif,lastcall)
+          pdqdif,pdqevap,pdqsdif,flux_u,flux_v,lastcall)
 
       use watercommon_h, only : RLVTT, T_h2O_ice_liq, RCPD, mx_eau_sol,Psat_water
@@ -65 +65 @@
       REAL,INTENT(IN) :: pcapcal(ngrid)
       REAL,INTENT(INOUT) :: pq2(ngrid,nlay+1)
-      
-      integer,intent(in) :: rnat(ngrid)      
+      REAL,INTENT(OUT) :: flux_u(ngrid),flux_v(ngrid)
+      REAL,INTENT(IN) :: rnat(ngrid)      
       LOGICAL,INTENT(IN) :: lastcall ! not used
 
@@ -334 +334 @@
       ENDDO
 
+!     Calcul of wind stress
+
+      DO ig=1,ngrid
+         flux_u(ig) = zfluxv(ig,1)/ptimestep*zu(ig,1)
+         flux_v(ig) = zfluxv(ig,1)/ptimestep*zv(ig,1)
+      ENDDO
 
 
@@ -499 +505 @@
                ! compute evaporation efficiency
                do ig=1,ngrid
-                  if(rnat(ig).eq.1)then
+                  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)
@@ -601 +607 @@
 !     The surface vapour tracer is actually liquid. To make things difficult.
 
-                  if (rnat(ig).eq.0) then ! unfrozen ocean
+                  if (nint(rnat(ig)).eq.0) then ! unfrozen ocean
                      
                      pdqsdif(ig,iliq_surf)=dqsdif_total(ig)/ptimestep
                      pdqsdif(ig,iice_surf)=0.0
 
-                  elseif (rnat(ig).eq.1) then ! (continent)
+                  elseif (nint(rnat(ig)).eq.1) then ! (continent)
 !     If water is evaporating / subliming, we take it from ice before liquid
 !     -- is this valid??
@@ -628 +634 @@
                      endif
 
-                  elseif (rnat(ig).eq.2) then ! (continental glaciers)
+                  elseif (nint(rnat(ig)).eq.2) then ! (continental glaciers)
                      pdqsdif(ig,iliq_surf)=0.0
                      pdqsdif(ig,iice_surf)=dqsdif_total(ig)/ptimestep

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

@@ -57 +57 @@
       REAL pq2(ngrid,nlay+1)
       
-      integer rnat(ngrid)      
+      real rnat(ngrid)      
 
 !     Arguments added for condensation
@@ -512 +512 @@
                ! compute evaporation efficiency
                do ig = 1, ngrid
-                  if(rnat(ig).eq.1)then
+                  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)
@@ -604 +604 @@
 !     The surface vapour tracer is actually liquid. To make things difficult.
 
-                  if (rnat(ig).eq.0) then ! unfrozen ocean
+                  if (nint(rnat(ig)).eq.0) then ! unfrozen ocean
                      
                      pdqsdif(ig,ivap)=dqsdif_total(ig)/ptimestep
@@ -610 +610 @@
 
 
-                  elseif (rnat(ig).eq.1) then ! (continent)
+                  elseif (nint(rnat(ig)).eq.1) then ! (continent)
 
 !     --------------------------------------------------------