source: LMDZ5/trunk/libf/dynlonlat_phylonlat/phylmd/ce0l.F90 @ 2331

PROGRAM ce0l
  !
  ! Purpose: Calls etat0, creates initial states and limit_netcdf
  !
  ! interbar=.T. for barycentric interpolation inter_barxy
  ! extrap  =.T. for data extrapolation, like for the SSTs when file does not
  !                  contain ocean points only.
  ! oldice  =.T. for old-style ice, obtained using grille_m (grid_atob).
  ! masque is created in etat0, passed to limit to ensure consistancy.

  USE control_mod, only: DAY_STEP, DAYREF, NSPLIT_PHYS
  USE etat0dyn, only: etat0dyn_netcdf
  USE netcdf, ONLY: NF90_OPEN, NF90_NOWRITE, NF90_CLOSE, NF90_NOERR
  USE ioipsl, ONLY: ioconf_calendar, getin, flininfo, flinopen, flinget, flinclo

  USE etat0phys, only: etat0phys_netcdf
  USE dimphy, only: KLON
  USE infotrac, only: TYPE_TRAC, infotrac_init
  USE test_disvert_m, ONLY: test_disvert

  IMPLICIT NONE

  ! Local variables:
  include "dimensions.h"
  include "paramet.h"
  include "comgeom.h"
  include "comconst.h"
  include "comvert.h"
  include "iniprint.h"
  include "temps.h"
  include "logic.h"
  REAL               :: masque(iip1, jjp1)             !--- CONTINENTAL MASK
  REAL               :: phis  (iip1, jjp1)             !--- GROUND GEOPOTENTIAL
  CHARACTER(LEN=256) :: modname, fmt, calnd           !--- CALENDAR TYPE
  LOGICAL            :: use_filtre_fft
  LOGICAL, PARAMETER :: interbar=.TRUE., extrap=.FALSE., oldice=.FALSE.

  !--- Local variables for ocean mask reading:
  INTEGER            :: nid_o2a, iml_omask, jml_omask, j
  INTEGER            :: fid, iret, llm_tmp, ttm_tmp, itaul(1)
  REAL, ALLOCATABLE  :: lon_omask(:, :), dlon_omask(:), ocemask(:, :)
  REAL, ALLOCATABLE  :: lat_omask(:, :), dlat_omask(:), ocetmp (:, :)
  REAL               :: date, lev(1)

  !----------------------------------------------------------------------
  modname="ce0l"

  !--- Constants
  pi     = 4. * ATAN(1.)
  rad    = 6371229.
  daysec = 86400.
  omeg   = 2.*pi/daysec
  g      = 9.8
  kappa  = 0.2857143
  cpp    = 1004.70885
  jmp1   = jjm + 1
  preff   = 101325.
  pa      = 50000.

  CALL conf_gcm( 99, .TRUE. )

  dtvr = daysec/REAL(day_step)
  WRITE(lunout, *)'dtvr', dtvr

  CALL iniconst()
  CALL inigeom()

#ifdef CPP_IOIPSL
  calnd='gregorian'
  SELECT CASE(calend)
  CASE('earth_360d')
     CALL ioconf_calendar('360d')
     calnd='with 360 days/year'
  CASE('earth_365d')
     CALL ioconf_calendar('noleap')
     calnd='with no leap year'
  CASE('earth_366d')
     CALL ioconf_calendar('366d')
     calnd='with leap years only'
  CASE('gregorian')
     CALL ioconf_calendar('gregorian')
  CASE('standard')
     CALL ioconf_calendar('gregorian')
  CASE('julian')
     CALL ioconf_calendar('julian')
     calnd='julian'
  CASE('proleptic_gregorian')
     CALL ioconf_calendar('gregorian')
     !--- DC Bof...  => IOIPSL a mettre a jour: proleptic_gregorian /= gregorian
  CASE DEFAULT
     CALL abort_gcm('ce0l', 'Bad choice for calendar', 1)
  END SELECT
  WRITE(lunout, *)'CHOSEN CALENDAR: Earth '//TRIM(calnd)
#endif

  use_filtre_fft=.FALSE.
  CALL getin('use_filtre_fft', use_filtre_fft)
  IF(use_filtre_fft) THEN
     WRITE(lunout, *)"FFT filter not available for sequential dynamics."
     WRITE(lunout, *)"Your setting of variable use_filtre_fft is not used."
  ENDIF

  !--- LAND MASK. TWO CASES:
  !   1) read from ocean model    file "o2a.nc"    (coupled runs)
  !   2) computed from topography file "Relief.nc" (masque(:, :)=-99999.)
  ! Coupled simulations (case 1) use the ocean model mask to compute the
  ! weights to ensure ocean fractions are the same for atmosphere and ocean.

  IF(NF90_OPEN("o2a.nc", NF90_NOWRITE, nid_o2a)/=NF90_NOERR) THEN
     WRITE(lunout, *)'BEWARE !! No ocean mask "o2a.nc" file found'
     WRITE(lunout, *)'Forced run.'
     masque(:, :)=-99999.
  ELSE
     iret=NF90_CLOSE(nid_o2a)
     WRITE(lunout, *)'BEWARE !! Ocean mask "o2a.nc" file found'
     WRITE(lunout, *)'Coupled run.'
     CALL flininfo("o2a.nc", iml_omask, jml_omask, llm_tmp, ttm_tmp, nid_o2a)
     IF(iml_omask/=iim .OR.jml_omask/=jjp1) THEN
        WRITE(lunout, *)'Mismatching dimensions for ocean mask'
        WRITE(lunout, *)'iim  = ', iim , ' iml_omask = ', iml_omask
        WRITE(lunout, *)'jjp1 = ', jjp1, ' jml_omask = ', jml_omask
        CALL abort_gcm(modname, '', 1)
     END IF
     ALLOCATE(ocemask(iim, jjp1), lon_omask(iim, jjp1), dlon_omask(iim ))
     ALLOCATE(ocetmp (iim, jjp1), lat_omask(iim, jjp1), dlat_omask(jjp1))
     CALL flinopen("o2a.nc", .FALSE., iim, jjp1, llm_tmp, lon_omask, &
          lat_omask, lev, ttm_tmp, itaul, date, dt, fid)
     CALL flinget(fid, "OceMask",    iim, jjp1, llm_tmp, ttm_tmp, 1, 1, ocetmp)
     CALL flinclo(fid)
     dlon_omask(1:iim ) = lon_omask(1:iim, 1)
     dlat_omask(1:jjp1) = lat_omask(1, 1:jjp1)
     ocemask = ocetmp
     IF(dlat_omask(1)<dlat_omask(jml_omask)) THEN
        DO j=1, jjp1
           ocemask(:, j) = ocetmp(:, jjp1-j+1)
        END DO
     END IF
     DEALLOCATE(ocetmp, lon_omask, lat_omask, dlon_omask, dlat_omask)
     IF(prt_level>=1) THEN
        WRITE(fmt, "(i4, 'i1)')")iim 
        fmt='('//ADJUSTL(fmt)
        WRITE(lunout, *)'OCEAN MASK :'
        WRITE(lunout, fmt) NINT(ocemask)
     END IF
     masque(1:iim, :)=1.-ocemask(:, :)
     masque(iip1 , :)=masque(1, :)
     DEALLOCATE(ocemask)
  END IF
  phis(:, :)=-99999.

  CALL Init_Phys_lmdz(iim, jjp1, llm, 1, (/(jjm-1)*iim+2/))
  WRITE(lunout, *)'---> klon=', klon

  call infotrac_init
  CALL iniphysiq(iim, jjm, llm, daysec, dayref, dtphys / nsplit_phys, rlatu, &
       rlonv, aire, cu, cv, rad, g, r, cpp, iflag_phys)

  IF(pressure_exner) CALL test_disvert

  IF (type_trac == 'inca') THEN
#ifdef INCA
     CALL init_const_lmdz(nbtr, anneeref, dayref, iphysiq, day_step, nday)
     CALL init_inca_para(iim, jjm+1, klon, 1, klon_mpi_para_nb, 0)
     WRITE(lunout, *)'nbtr =' , nbtr 
#endif
  END IF
  IF(ok_etat0) THEN
     WRITE(lunout, '(//)')
     WRITE(lunout, *) '  ************************  '
     WRITE(lunout, *) '  ***  etat0phy_netcdf ***  '
     WRITE(lunout, *) '  ************************  '
     WRITE(lunout, '(//)')
     WRITE(lunout, *) ' interbar = ', interbar
     CALL etat0phys_netcdf(interbar, masque, phis)
  END IF

  IF(ok_etat0) THEN
     WRITE(lunout, '(//)')
     WRITE(lunout, *) '  ************************  '
     WRITE(lunout, *) '  ***  etat0dyn_netcdf ***  '
     WRITE(lunout, *) '  ************************  '
     WRITE(lunout, '(//)')
     WRITE(lunout, *) ' interbar = ', interbar
     CALL etat0dyn_netcdf(interbar, masque, phis)
  END IF

  IF(ok_limit) THEN
     WRITE(lunout, '(//)')
     WRITE(lunout, *) '  *********************  '
     WRITE(lunout, *) '  ***  Limit_netcdf ***  '
     WRITE(lunout, *) '  *********************  '
     WRITE(lunout, '(//)')
     CALL limit_netcdf(interbar, extrap, oldice, masque)
  END IF

  WRITE(lunout, '(//)')
  WRITE(lunout, *) '  ***************************  '
  WRITE(lunout, *) '  ***  grilles_gcm_netcdf ***  '
  WRITE(lunout, *) '  ***************************  '
  WRITE(lunout, '(//)')
  CALL grilles_gcm_netcdf_sub(masque, phis)

END PROGRAM ce0l