source: LMDZ4/branches/LMDZ4V5.0-dev/libf/dyn3dpar/etat0_netcdf.F @ 5305

!
! $Id: etat0_netcdf.F 1379 2010-05-06 12:19:18Z abarral $
!
c
c
      SUBROUTINE etat0_netcdf (interbar, masque)
#ifdef CPP_EARTH        
      USE startvar
      USE ioipsl
      USE dimphy
      USE control_mod
      USE infotrac
      USE fonte_neige_mod
      USE pbl_surface_mod
      USE phys_state_var_mod
      USE filtreg_mod
      use regr_lat_time_climoz_m, only: regr_lat_time_climoz
      use conf_phys_m, only: conf_phys
!     For parameterization of ozone chemistry:
      use regr_lat_time_coefoz_m, only: regr_lat_time_coefoz
      use press_coefoz_m, only: press_coefoz
      use regr_pr_o3_m, only: regr_pr_o3
#endif
!#endif of #ifdef CPP_EARTH
      use netcdf, only: nf90_open, NF90_NOWRITE, nf90_close
      !
      IMPLICIT NONE
      !
#include "dimensions.h"
#include "paramet.h"
      !
      !
!      INTEGER, PARAMETER :: KIDIA=1, KFDIA=iim*(jjm-1)+2, 
!     .KLON=KFDIA-KIDIA+1,KLEV=llm
      !
#ifdef CPP_EARTH    
#include "comgeom2.h"
#include "comvert.h"
#include "comconst.h"
#include "indicesol.h"
#include "dimsoil.h"
#include "temps.h"
#endif
!#endif of #ifdef CPP_EARTH
      ! arguments:
      LOGICAL interbar
      REAL :: masque(iip1,jjp1)

#ifdef CPP_EARTH
      ! local variables:
      REAL :: latfi(klon), lonfi(klon)
      REAL :: orog(iip1,jjp1), rugo(iip1,jjp1)
      REAL :: psol(iip1, jjp1), phis(iip1, jjp1)
      REAL :: p3d(iip1, jjp1, llm+1)
      REAL :: uvent(iip1, jjp1, llm)
      REAL :: vvent(iip1, jjm, llm)
      REAL :: t3d(iip1, jjp1, llm), tpot(iip1, jjp1, llm)
      REAL :: qsat(iip1, jjp1, llm)
      REAL,ALLOCATABLE :: q3d(:, :, :,:)
      REAL :: tsol(klon), qsol(klon), sn(klon)
!!      REAL :: tsolsrf(klon,nbsrf)
      real qsolsrf(klon,nbsrf),snsrf(klon,nbsrf) 
      REAL :: albe(klon,nbsrf), evap(klon,nbsrf)
      REAL :: alblw(klon,nbsrf)
      REAL :: tsoil(klon,nsoilmx,nbsrf) 
      REAL :: frugs(klon,nbsrf), agesno(klon,nbsrf)
      REAL :: rugmer(klon)
      REAL :: qd(iip1, jjp1, llm)
      REAL :: run_off_lic_0(klon)
      ! declarations pour lecture glace de mer
      REAL :: rugv(klon)
      INTEGER :: iml_lic, jml_lic, llm_tmp, ttm_tmp, iret
      INTEGER :: itaul(1), fid
      REAL :: lev(1), date
      REAL, ALLOCATABLE, DIMENSION(:,:) :: lon_lic, lat_lic
      REAL, ALLOCATABLE, DIMENSION(:)  :: dlon_lic, dlat_lic
      REAL, ALLOCATABLE, DIMENSION (:,:) :: fraclic
      REAL :: flic_tmp(iip1, jjp1)
      REAL :: champint(iim, jjp1)
      !

      CHARACTER(len=80) :: varname
      !
      INTEGER :: i,j, ig, l, ji,ii1,ii2
      REAL :: xpi
      !
      REAL :: alpha(iip1,jjp1,llm),beta(iip1,jjp1,llm)
      REAL :: pk(iip1,jjp1,llm), pls(iip1,jjp1,llm), pks(ip1jmp1)
      REAL :: workvar(iip1,jjp1,llm)
      !
      REAL ::  prefkap, unskap
      !
      real :: time_step,t_ops,t_wrt

#include "comdissnew.h"
#include "serre.h"
#include "clesphys.h"

      INTEGER  ::        longcles
      PARAMETER      ( longcles  = 20 )
      REAL :: clesphy0 ( longcles       )
      REAL :: p(iip1,jjp1,llm)
      INTEGER :: itau, iday
      REAL :: masse(iip1,jjp1,llm)
      REAL :: xpn,xps,xppn(iim),xpps(iim)
      real :: time
      REAL :: phi(ip1jmp1,llm)
      REAL :: pbaru(ip1jmp1,llm),pbarv(ip1jm,llm)
      REAL :: w(ip1jmp1,llm)
      REAL ::phystep
CC      REAL :: rugsrel(iip1*jjp1)
      REAL :: fder(klon)
!!      real zrel(iip1*jjp1),chmin,chmax

!!      CHARACTER(len=80) :: visu_file
      INTEGER :: visuid

! pour la lecture du fichier masque ocean
      integer :: nid_o2a
      logical :: couple = .false.
      INTEGER :: iml_omask, jml_omask
      REAL, ALLOCATABLE, DIMENSION(:,:) :: lon_omask, lat_omask
      REAL, ALLOCATABLE, DIMENSION(:)  :: dlon_omask, dlat_omask
      REAL, ALLOCATABLE, DIMENSION (:,:) :: ocemask, ocetmp
      real, dimension(klon) :: ocemask_fi
      integer :: isst(klon-2)
      real zx_tmp_2d(iim,jjp1)

      REAL :: dummy

      logical              :: ok_newmicro
      integer              :: iflag_radia
      logical              :: ok_journe, ok_mensuel, ok_instan, ok_hf
      logical              :: ok_LES
      LOGICAL              :: ok_ade, ok_aie, aerosol_couple, new_aod
      INTEGER              :: flag_aerosol
      REAL                 :: bl95_b0, bl95_b1
      real                 :: fact_cldcon, facttemps,ratqsbas,ratqshaut
      real                 :: tau_ratqs
      integer              :: iflag_cldcon
      integer              :: iflag_ratqs
      integer :: iflag_coupl
      integer :: iflag_clos
      integer :: iflag_wake
      integer :: iflag_thermals,nsplit_thermals
      real    :: tau_thermals
      integer :: iflag_thermals_ed,iflag_thermals_optflux
      REAL      :: solarlong0
      real :: seuil_inversion
      real :: alp_offset
      logical found

      integer  read_climoz ! read ozone climatology
C     Allowed values are 0, 1 and 2
C     0: do not read an ozone climatology
C     1: read a single ozone climatology that will be used day and night
C     2: read two ozone climatologies, the average day and night
C     climatology and the daylight climatology

      !
      !   Constantes 
      !
      pi     = 4. * ATAN(1.)
      rad    = 6371229.
      omeg   = 4.* ASIN(1.)/(24.*3600.)
      g      = 9.8
      daysec = 86400.
      kappa  = 0.2857143
      cpp    = 1004.70885
      !
      preff     = 101325.
      pa        =  50000.
      unskap = 1./kappa
      !
      jmp1    = jjm + 1
      !
      !    Construct a grid
      !

!      CALL defrun_new(99,.TRUE.,clesphy0)
      CALL conf_gcm( 99, .TRUE. , clesphy0 )
      call conf_phys(  ok_journe, ok_mensuel, ok_instan, ok_hf, ok_LES, &
     &                 solarlong0,seuil_inversion,                      &
     &                 fact_cldcon, facttemps,ok_newmicro,iflag_radia,  &
     &                 iflag_cldcon,                                    &
     &                 iflag_ratqs,ratqsbas,ratqshaut,tau_ratqs,        &
     &                 ok_ade, ok_aie, aerosol_couple,                  &
     &                 flag_aerosol, new_aod,                           &
     &                 bl95_b0, bl95_b1,                                &
     &                 iflag_thermals,nsplit_thermals,tau_thermals,     &
     &                 iflag_thermals_ed,iflag_thermals_optflux,        &
     &                 iflag_coupl,iflag_clos,iflag_wake, read_climoz,  &
     &                 alp_offset)

! co2_ppm0 : initial value of atmospheric CO2 from .def file (co2_ppm value)
      co2_ppm0 = co2_ppm

      dtvr   = daysec/REAL(day_step)
      print*,'dtvr',dtvr

      CALL iniconst()
      CALL inigeom()

! Initialisation pour traceurs
      call infotrac_init
      ALLOCATE(q3d(iip1, jjp1, llm, nqtot))

      CALL inifilr()
      CALL phys_state_var_init(read_climoz)
      !
      latfi(1) = ASIN(1.0)
      DO j = 2, jjm
        DO i = 1, iim
          latfi((j-2)*iim+1+i)=  rlatu(j)
        ENDDO
      ENDDO
      latfi(klon) = - ASIN(1.0)
      !
      lonfi(1) = 0.0
      DO j = 2, jjm
        DO i = 1, iim
          lonfi((j-2)*iim+1+i) =  rlonv(i)
        ENDDO
      ENDDO
      lonfi(klon) = 0.0
      !
      xpi = 2.0 * ASIN(1.0)
      DO ig = 1, klon
        latfi(ig) = latfi(ig) * 180.0 / xpi
        lonfi(ig) = lonfi(ig) * 180.0 / xpi
      ENDDO
      !
      rlat(1) = ASIN(1.0)
      DO j = 2, jjm
        DO i = 1, iim
          rlat((j-2)*iim+1+i)=  rlatu(j)
        ENDDO
      ENDDO
      rlat(klon) = - ASIN(1.0)
      !
      rlon(1) = 0.0
      DO j = 2, jjm
        DO i = 1, iim
          rlon((j-2)*iim+1+i) =  rlonv(i)
        ENDDO
      ENDDO
      rlon(klon) = 0.0
      !
      xpi = 2.0 * ASIN(1.0)
      DO ig = 1, klon
        rlat(ig) = rlat(ig) * 180.0 / xpi
        rlon(ig) = rlon(ig) * 180.0 / xpi
      ENDDO
      !
      


C
C En cas de simulation couplee, lecture du masque ocean issu du modele ocean
C utilise pour calculer les poids et pour assurer l'adequation entre les
C fractions d'ocean vu par l'atmosphere et l'ocean. Sinon, on cree le masque 
C a partir du fichier relief
C

      write(*,*)'Essai de lecture masque ocean'
      iret = nf90_open("o2a.nc", NF90_NOWRITE, nid_o2a)
      if (iret .ne. 0) then
        write(*,*)'ATTENTION!! pas de fichier o2a.nc trouve'
        write(*,*)'Run force'
        varname = 'masque'
        masque(:,:) = 0.0
        CALL startget_phys2d(varname, iip1, jjp1, rlonv, rlatu, masque,
     $       0.0, jjm ,rlonu,rlatv , interbar )
        WRITE(*,*) 'MASQUE construit : Masque'
        WRITE(*,'(97I1)') nINT(masque(:,:))
        call gr_dyn_fi(1, iip1, jjp1, klon, masque, zmasq)
        WHERE (zmasq(1 : klon) .LT. EPSFRA)
            zmasq(1 : klon) = 0.
        END WHERE 
        WHERE (1. - zmasq(1 : klon) .LT. EPSFRA)
            zmasq(1 : klon) = 1.
        END WHERE 
      else
        couple = .true.
        iret = nf90_close(nid_o2a)
        call flininfo("o2a.nc", iml_omask, jml_omask, llm_tmp, ttm_tmp
     $    , nid_o2a)
        if (iml_omask /= iim .or. jml_omask /= jjp1) then
          write(*,*)'Dimensions non compatibles pour masque ocean'
          write(*,*)'iim = ',iim,' iml_omask = ',iml_omask
          write(*,*)'jjp1 = ',jjp1,' jml_omask = ',jml_omask
          stop
        endif
        ALLOCATE(lat_omask(iml_omask, jml_omask), stat=iret)
        ALLOCATE(lon_omask(iml_omask, jml_omask), stat=iret)
        ALLOCATE(dlon_omask(iml_omask), stat=iret)
        ALLOCATE(dlat_omask(jml_omask), stat=iret)
        ALLOCATE(ocemask(iml_omask, jml_omask), stat=iret)
        ALLOCATE(ocetmp(iml_omask, jml_omask), stat=iret)
        CALL flinopen("o2a.nc", .FALSE., iml_omask, jml_omask, llm_tmp
     $    , lon_omask, lat_omask, lev, ttm_tmp, itaul, date, dt, fid)
        CALL flinget(fid, 'OceMask', iml_omask, jml_omask, llm_tmp, 
     $      ttm_tmp, 1, 1, ocetmp)
        CALL flinclo(fid)
        dlon_omask(1 : iml_omask) = lon_omask(1 : iml_omask, 1)
        dlat_omask(1 : jml_omask) = lat_omask(1 , 1 : jml_omask)
        ocemask = ocetmp
        if (dlat_omask(1) < dlat_omask(jml_omask)) then
          do j = 1, jml_omask
            ocemask(:,j) = ocetmp(:,jml_omask-j+1)
          enddo
        endif 
C
C passage masque ocean a la grille physique
C
        write(*,*)'ocemask '
        write(*,'(96i1)')int(ocemask)
        ocemask_fi(1) = ocemask(1,1)
        do j = 2, jjm
          do i = 1, iim
            ocemask_fi((j-2)*iim + i + 1) = ocemask(i,j)
          enddo
        enddo
        ocemask_fi(klon) = ocemask(1,jjp1)
        zmasq = 1. - ocemask_fi
      endif

      call gr_fi_dyn(1, klon, iip1, jjp1, zmasq, masque)

      varname = 'relief'
      ! This line needs to be replaced by a call to restget to get the values in the restart file
      orog(:,:) = 0.0
       CALL startget_phys2d(varname, iip1, jjp1, rlonv, rlatu, orog,
     $     0.0 , jjm ,rlonu,rlatv , interbar, masque )
      !
      WRITE(*,*) 'OUT OF GET VARIABLE : Relief'
!      WRITE(*,'(49I1)') INT(orog(:,:))
      !
      varname = 'rugosite'
      ! This line needs to be replaced by a call to restget to get the values in the restart file
      rugo(:,:) = 0.0
       CALL startget_phys2d(varname, iip1, jjp1, rlonv, rlatu, rugo,
     $     0.0 , jjm, rlonu,rlatv , interbar )
      !
      WRITE(*,*) 'OUT OF GET VARIABLE : Rugosite' 
!      WRITE(*,'(49I1)') INT(rugo(:,:)*10)
      !
C
C on initialise les sous surfaces
C
      pctsrf=0.
c
      varname = 'psol'
      psol(:,:) = 0.0
      CALL startget_phys2d(varname, iip1, jjp1, rlonv, rlatu, psol,
     $     0.0 , jjm ,rlonu,rlatv , interbar )
      !
      !  Compute here the pressure on the intermediate levels. One would expect that this is available in the GCM 
      !  anyway.
      !
!      WRITE(*,*) 'PSOL :', psol(10,20)
!      WRITE(*,*) ap(:), bp(:)
      CALL pression(ip1jmp1, ap, bp, psol, p3d)
!      WRITE(*,*) 'P3D :', p3d(10,20,:)
      CALL exner_hyb(ip1jmp1, psol, p3d, alpha, beta, pks, pk, workvar)
!      WRITE(*,*) 'PK:', pk(10,20,:)
      !
      !
      !
      prefkap =  preff  ** kappa
!      WRITE(*,*) 'unskap, cpp,  preff :', unskap, cpp,  preff
      DO l = 1, llm
        DO j=1,jjp1
          DO i =1, iip1
            pls(i,j,l) = preff * ( pk(i,j,l)/cpp) ** unskap
           ENDDO
        ENDDO
      ENDDO
      !
!      WRITE(*,*) 'PLS :', pls(10,20,:)
      !
      varname = 'surfgeo'
      phis(:,:) = 0.0
      CALL startget_phys2d(varname, iip1, jjp1, rlonv, rlatu, phis,
     $     0.0 , jjm ,rlonu,rlatv, interbar )
      !
      varname = 'u'
      uvent(:,:,:) = 0.0
      CALL startget_dyn(varname, rlonu, rlatu, pls, workvar, uvent, 0.,
     $     rlonv, rlatv, interbar )
      !  
      varname = 'v'
      vvent(:,:,:) = 0.0
      CALL startget_dyn(varname, rlonv, rlatv, pls(:, :jjm, :),
     . workvar(:, :jjm, :), vvent, 0., rlonu, rlatu(:jjm), interbar )
      !
      varname = 't'
      t3d(:,:,:) = 0.0
      CALL startget_dyn(varname, rlonv, rlatu, pls, workvar, t3d, 0.,
     $     rlonu, rlatv , interbar )
      !
      WRITE(*,*) 'T3D min,max:',minval(t3d(:,:,:)),
     .                          maxval(t3d(:,:,:))
      varname = 'tpot'
      tpot(:,:,:) = 0.0
      CALL startget_dyn(varname, rlonv, rlatu, pls, pk, tpot, 0., rlonu,
     $     rlatv, interbar)
      !
      WRITE(*,*) 'T3D min,max:',minval(t3d(:,:,:)),
     .                          maxval(t3d(:,:,:))
      WRITE(*,*) 'PLS min,max:',minval(pls(:,:,:)),
     .                          maxval(pls(:,:,:))

c Calcul de l'humidite a saturation
      print*,'avant q_sat'
      call q_sat(llm*jjp1*iip1,t3d,pls,qsat)
      print*,'apres q_sat'

      WRITE(*,*) 'QSAT min,max:',minval(qsat(:,:,:)),
     .                           maxval(qsat(:,:,:))
      !
CC      WRITE(*,*) 'QSAT :', qsat(10,20,:)
      !
      varname = 'q'
      qd(:,:,:) = 0.0
      q3d(:,:,:,:) = 0.0
      WRITE(*,*) 'QSAT min,max:',minval(qsat(:,:,:)),
     .                           maxval(qsat(:,:,:))
      CALL startget_dyn(varname, rlonv, rlatu, pls, qsat, qd, 0., rlonu,
     $     rlatv , interbar )
      q3d(:,:,:,1) = qd(:,:,:)
      !
!     Parameterization of ozone chemistry:
C     Look for ozone tracer:
      i = 1
      DO
         found = tname(i)=="O3" .OR. tname(i)=="o3"
         if (found .or. i == nqtot) exit
         i = i + 1
      end do
      if (found) then
         call regr_lat_time_coefoz
         call press_coefoz
         call regr_pr_o3(p3d, q3d(:, :, :, i))
C     Convert from mole fraction to mass fraction:
         q3d(:, :, :, i) = q3d(:, :, :, i)  * 48. / 29.
      end if

!     Ozone climatology:
      if (read_climoz >= 1) call regr_lat_time_climoz(read_climoz)

      varname = 'tsol'
      ! This line needs to be replaced by a call to restget to get the values in the restart file
      tsol(:) = 0.0
      CALL startget_phys1d(varname, iip1, jjp1, rlonv, rlatu, klon,
     $     tsol, 0.0, jjm, rlonu, rlatv , interbar )
      !
      WRITE(*,*) 'TSOL construit :'
!      WRITE(*,'(48I3)') INT(TSOL(2:klon)-273)
      !
      varname = 'qsol'
      qsol(:) = 0.0
      CALL startget_phys1d(varname, iip1, jjp1, rlonv, rlatu, klon,
     $     qsol, 0.0, jjm, rlonu, rlatv , interbar )
      !
      varname = 'snow'
      sn(:) = 0.0
      CALL startget_phys1d(varname, iip1, jjp1, rlonv, rlatu, klon, sn,
     $     0.0, jjm, rlonu, rlatv , interbar )
      !
      varname = 'rads'
      radsol(:) = 0.0
      CALL startget_phys1d(varname,iip1,jjp1,rlonv,rlatu,klon,radsol,
     $     0.0, jjm, rlonu, rlatv , interbar )
      !
      varname = 'rugmer'
      rugmer(:) = 0.0
      CALL startget_phys1d(varname,iip1,jjp1,rlonv,rlatu,klon,rugmer,
     $     0.0, jjm, rlonu, rlatv , interbar )
      !
!      varname = 'agesno'
!      agesno(:) = 0.0
!      CALL startget_phys1d(varname,iip1,jjp1,rlonv,rlatu,klon,agesno,0.0,
!     .     jjm, rlonu, rlatv , interbar )

      varname = 'zmea'
      zmea(:) = 0.0
      CALL startget_phys1d(varname,iip1,jjp1,rlonv,rlatu,klon,zmea,0.0,
     .     jjm, rlonu, rlatv , interbar )

      varname = 'zstd'
      zstd(:) = 0.0
      CALL startget_phys1d(varname,iip1,jjp1,rlonv,rlatu,klon,zstd,0.0,
     .     jjm, rlonu, rlatv , interbar )
      varname = 'zsig'
      zsig(:) = 0.0
      CALL startget_phys1d(varname,iip1,jjp1,rlonv,rlatu,klon,zsig,0.0,
     .     jjm, rlonu, rlatv , interbar )
      varname = 'zgam'
      zgam(:) = 0.0
      CALL startget_phys1d(varname,iip1,jjp1,rlonv,rlatu,klon,zgam,0.0,
     .     jjm, rlonu, rlatv , interbar )
      varname = 'zthe'
      zthe(:) = 0.0
      CALL startget_phys1d(varname,iip1,jjp1,rlonv,rlatu,klon,zthe,0.0,
     .     jjm, rlonu, rlatv , interbar )
      varname = 'zpic'
      zpic(:) = 0.0
      CALL startget_phys1d(varname,iip1,jjp1,rlonv,rlatu,klon,zpic,0.0,
     .     jjm, rlonu, rlatv , interbar )
      varname = 'zval'
      zval(:) = 0.0
      CALL startget_phys1d(varname,iip1,jjp1,rlonv,rlatu,klon,zval,0.0,
     .     jjm, rlonu, rlatv , interbar )
c
cc      rugsrel(:) = 0.0
cc      IF(ok_orodr)  THEN
cc        DO i = 1, iip1* jjp1
cc         rugsrel(i) = MAX( 1.e-05, zstd(i)* zsig(i) /2. )
cc        ENDDO
cc      ENDIF


C
C lecture du fichier glace de terre pour fixer la fraction de terre 
C et de glace de terre
C
      CALL flininfo("landiceref.nc", iml_lic, jml_lic,llm_tmp, ttm_tmp
     $    , fid)
      ALLOCATE(lat_lic(iml_lic, jml_lic), stat=iret)
      ALLOCATE(lon_lic(iml_lic, jml_lic), stat=iret)
      ALLOCATE(dlon_lic(iml_lic), stat=iret)
      ALLOCATE(dlat_lic(jml_lic), stat=iret)
      ALLOCATE(fraclic(iml_lic, jml_lic), stat=iret)
      CALL flinopen("landiceref.nc", .FALSE., iml_lic, jml_lic, llm_tmp
     $    , lon_lic, lat_lic, lev, ttm_tmp, itaul, date, dt, fid)
      CALL flinget(fid, 'landice', iml_lic, jml_lic, llm_tmp, ttm_tmp
     $    , 1, 1, fraclic)
      CALL flinclo(fid)
C
C interpolation sur la grille T du modele
C
      WRITE(*,*) 'dimensions de landice iml_lic, jml_lic : ', 
     $    iml_lic, jml_lic
c
C sil les coordonnees sont en degres, on les transforme
C
      IF( MAXVAL( lon_lic(:,:) ) .GT. 2.0 * asin(1.0) )  THEN
          lon_lic(:,:) = lon_lic(:,:) * 2.0* ASIN(1.0) / 180.
      ENDIF 
      IF( maxval( lat_lic(:,:) ) .GT. 2.0 * asin(1.0)) THEN 
          lat_lic(:,:) = lat_lic(:,:) * 2.0 * asin(1.0) / 180.
      ENDIF 

      dlon_lic(1 : iml_lic) = lon_lic(1 : iml_lic, 1)
      dlat_lic(1 : jml_lic) = lat_lic(1 , 1 : jml_lic) 
C
      CALL grille_m(iml_lic, jml_lic, dlon_lic, dlat_lic, fraclic
     $    ,iim, jjp1,
     $    rlonv, rlatu, flic_tmp(1 : iim, 1 : jjp1))
cx$$$      flic_tmp(1 : iim, 1 : jjp1) = champint(1: iim, 1 : jjp1)
      flic_tmp(iip1, 1 : jjp1) = flic_tmp(1 , 1 : jjp1)
C
C passage sur la grille physique
C
      CALL gr_dyn_fi(1, iip1, jjp1, klon, flic_tmp,
     $    pctsrf(1:klon, is_lic))
C adequation avec le maque terre/mer
c      zmasq(157) = 0.
      WHERE (pctsrf(1 : klon, is_lic) .LT. EPSFRA ) 
          pctsrf(1 : klon, is_lic) = 0. 
      END WHERE
      WHERE (zmasq( 1 : klon) .LT. EPSFRA) 
          pctsrf(1 : klon, is_lic) = 0.
      END WHERE 
      pctsrf(1 : klon, is_ter) = zmasq(1 : klon)
      DO ji = 1, klon
        IF (zmasq(ji) .GT. EPSFRA) THEN 
            IF ( pctsrf(ji, is_lic) .GE. zmasq(ji)) THEN
                pctsrf(ji, is_lic) = zmasq(ji)
                pctsrf(ji, is_ter) = 0.
            ELSE 
                pctsrf(ji,is_ter) = zmasq(ji) - pctsrf(ji, is_lic)
                IF (pctsrf(ji,is_ter) .LT. EPSFRA) THEN
                    pctsrf(ji,is_ter) = 0.
                    pctsrf(ji, is_lic) = zmasq(ji)
                ENDIF 
            ENDIF 
        ENDIF 
      END DO 
C
C sous surface ocean et glace de mer (pour demarrer on met glace de mer a 0)
C
      pctsrf(1 : klon, is_oce) = (1. - zmasq(1 : klon))


      WHERE (pctsrf(1 : klon, is_oce) .LT. EPSFRA)
          pctsrf(1 : klon, is_oce) = 0.
      END WHERE 

      if (couple) pctsrf(1 : klon, is_oce) = ocemask_fi(1 : klon)

      isst = 0
      where (pctsrf(2:klon-1,is_oce) >0.) isst = 1
C
C verif que somme des sous surface = 1
C
      ji=count( (abs( sum(pctsrf(1 : klon, 1 : nbsrf),dim=2))-1.0) 
     $    .GT. EPSFRA)
      IF (ji .NE. 0) THEN
          WRITE(*,*) 'pb repartition sous maille pour ',ji,' points'
      ENDIF 

!      where (pctsrf(1:klon, is_ter) >= .5) 
!        pctsrf(1:klon, is_ter) = 1.
!        pctsrf(1:klon, is_oce) = 0.
!        pctsrf(1:klon, is_sic) = 0.
!        pctsrf(1:klon, is_lic) = 0.
!        zmasq = 1.
!      endwhere
!      where (pctsrf(1:klon, is_lic) >= .5) 
!        pctsrf(1:klon, is_ter) = 0.
!        pctsrf(1:klon, is_oce) = 0.
!        pctsrf(1:klon, is_sic) = 0.
!        pctsrf(1:klon, is_lic) = 1.
!        zmasq = 1.
!      endwhere
!      where (pctsrf(1:klon, is_oce) >= .5) 
!        pctsrf(1:klon, is_ter) = 0.
!        pctsrf(1:klon, is_oce) = 1.
!        pctsrf(1:klon, is_sic) = 0.
!        pctsrf(1:klon, is_lic) = 0.
!        zmasq = 0.
!      endwhere
!      where (pctsrf(1:klon, is_sic) >= .5) 
!        pctsrf(1:klon, is_ter) = 0.
!        pctsrf(1:klon, is_oce) = 0.
!        pctsrf(1:klon, is_sic) = 1.
!        pctsrf(1:klon, is_lic) = 0.
!        zmasq = 0.
!      endwhere
!      call gr_fi_dyn(1, klon, iip1, jjp1, zmasq, masque)
C
C verif que somme des sous surface = 1
C
!      ji=count( (abs( sum(pctsrf(1 : klon, 1 : nbsrf), dim = 2)) - 1.0 ) 
!     $    .GT. EPSFRA)
!      IF (ji .NE. 0) THEN
!          WRITE(*,*) 'pb repartition sous maille pour ',ji,' points'
!     ENDIF 

      CALL gr_fi_ecrit(1,klon,iim,jjp1,zmasq,zx_tmp_2d)
      write(*,*)'zmasq = '
      write(*,'(96i1)')nint(zx_tmp_2d)
      call gr_fi_dyn(1, klon, iip1, jjp1, zmasq, masque)
      WRITE(*,*) 'MASQUE construit : Masque'
      WRITE(*,'(97I1)') nINT(masque(:,:))



C Calcul intermediaire
c 
      CALL massdair( p3d, masse  )
c

      print *,' ALPHAX ',alphax

      DO  l = 1, llm
        DO  i    = 1, iim
          xppn(i) = aire( i, 1   ) * masse(  i     ,  1   , l )
          xpps(i) = aire( i,jjp1 ) * masse(  i     , jjp1 , l )
        ENDDO
          xpn      = SUM(xppn)/apoln
          xps      = SUM(xpps)/apols
        DO i   = 1, iip1
          masse(   i   ,   1     ,  l )   = xpn
          masse(   i   ,   jjp1  ,  l )   = xps
        ENDDO
      ENDDO
      q3d(iip1,:,:,:) = q3d(1,:,:,:)
      phis(iip1,:) = phis(1,:)

C Ecriture
      CALL inidissip( lstardis, nitergdiv, nitergrot, niterh   ,
     *                tetagdiv, tetagrot , tetatemp              )
      print*,'sortie inidissip'
      itau = 0
      itau_dyn = 0
      itau_phy = 0
      iday = dayref +itau/day_step
      time = real(itau-(iday-dayref)*day_step)/day_step
c     
      IF(time.GT.1)  THEN
       time = time - 1
       iday = iday + 1
      ENDIF
      day_ref = dayref
      annee_ref = anneeref

      CALL geopot  ( ip1jmp1, tpot  , pk , pks,  phis  , phi   )
      print*,'sortie geopot'
      
      CALL caldyn0 ( itau,uvent,vvent,tpot,psol,masse,pk,phis ,
     *                phi,w, pbaru,pbarv,time+iday-dayref   )
       print*,'sortie caldyn0'     
      CALL dynredem0("start.nc",dayref,phis)
      print*,'sortie dynredem0'
      CALL dynredem1("start.nc",0.0,vvent,uvent,tpot,q3d,masse ,
     .                            psol)
      print*,'sortie dynredem1' 
C
C Ecriture etat initial physique
C
      write(*,*)'phystep ',dtvr,iphysiq,nbapp_rad
      phystep   = dtvr * REAL(iphysiq)
      radpas    = NINT (86400./phystep/ REAL(nbapp_rad) )
      write(*,*)'phystep =', phystep, radpas
cIM : lecture de co2_ppm & solaire ds physiq.def
c     co2_ppm   = 348.0
c     solaire   = 1365.0

c
c Initialisation 
c tsol, qsol, sn,albe, evap,tsoil,rain_fall, snow_fall,solsw, sollw,frugs
c
      ftsol(:,is_ter) = tsol
      ftsol(:,is_lic) = tsol
      ftsol(:,is_oce) = tsol
      ftsol(:,is_sic) = tsol
      snsrf(:,is_ter) = sn
      snsrf(:,is_lic) = sn
      snsrf(:,is_oce) = sn
      snsrf(:,is_sic) = sn
      falb1(:,is_ter) = 0.08
      falb1(:,is_lic) = 0.6
      falb1(:,is_oce) = 0.5
      falb1(:,is_sic) = 0.6
      falb2 = falb1
      evap(:,:) = 0.
      qsolsrf(:,is_ter) = 150
      qsolsrf(:,is_lic) = 150
      qsolsrf(:,is_oce) = 150.
      qsolsrf(:,is_sic) = 150.
      do i = 1, nbsrf
        do j = 1, nsoilmx
          tsoil(:,j,i) = tsol
        enddo
      enddo
      rain_fall = 0.; snow_fall = 0.
      solsw = 165.
      sollw = -53.
      t_ancien = 273.15
      q_ancien = 0.
      agesno = 0.
c
      frugs(1:klon,is_oce) = rugmer(1:klon)
      frugs(1:klon,is_ter) = MAX(1.0e-05, zstd(1:klon)*zsig(1:klon)/2.0)
      frugs(1:klon,is_lic) = MAX(1.0e-05, zstd(1:klon)*zsig(1:klon)/2.0)
      frugs(1:klon,is_sic) = 0.001
      fder = 0.0
      clwcon = 0.0
      rnebcon = 0.0
      ratqs = 0.0
      run_off_lic_0 = 0.0 
      rugoro = 0.0

c
c Avant l'appel a phyredem, on initialize les modules de surface
c avec les valeurs qui vont etre ecrit dans startphy.nc
c
      dummy = 1.0
      pbl_tke(:,:,:) = 1.e-8 
      zmax0(:) = 40.
      f0(:) = 1.e-5
      ema_work1(:,:) = 0.
      ema_work2(:,:) = 0.
      wake_deltat(:,:) = 0.
      wake_deltaq(:,:) = 0.
      wake_s(:) = 0.
      wake_cstar(:) = 0.
      wake_fip(:) = 0.
      wake_pe = 0.
      fm_therm = 0.
      entr_therm = 0.
      detr_therm = 0.

      call fonte_neige_init(run_off_lic_0)
      call pbl_surface_init(qsol, fder, snsrf, qsolsrf,
     $     evap, frugs, agesno, tsoil)

      call phyredem("startphy.nc")



C     Sortie Visu pour les champs dynamiques
cc      if (1.eq.0 ) then
cc      print*,'sortie visu'
cc      time_step = 1.
cc      t_ops = 2.
cc      t_wrt = 2.
cc      itau = 2.
cc      visu_file='Etat0_visu.nc'
cc      CALL initdynav(visu_file,dayref,anneeref,time_step,
cc     .              t_ops, t_wrt, visuid)
cc      CALL writedynav(visuid, itau,vvent ,
cc     .                uvent,tpot,pk,phi,q3d,masse,psol,phis)
cc      else
         print*,'CCCCCCCCCCCCCCCCCC REACTIVER SORTIE VISU DANS ETAT0'
cc      endif
      print*,'entree histclo'
      CALL histclo

#endif 
!#endif of #ifdef CPP_EARTH
      RETURN
      !
      END SUBROUTINE etat0_netcdf