trunk/dyn3d/etat0.f

module etat0_mod

  use indicesol, only: nbsrf
  use dimphy, only: klon

  IMPLICIT NONE

  REAL pctsrf(klon, nbsrf)
  ! ("pctsrf(i, :)" is the composition of the surface at horizontal
  ! position "i")

  private nbsrf, klon

contains

  SUBROUTINE etat0

    ! From "etat0_netcdf.F", version 1.3, 2005/05/25 13:10:09

    use caldyn0_m, only: caldyn0
    use comconst, only: cpp, kappa, iniconst
    use comgeom, only: rlatu, rlonv, rlonu, rlatv, aire_2d, apoln, apols, &
         cu_2d, cv_2d, inigeom
    use conf_gcm_m, only: dayref, anneeref
    use dimens_m, only: iim, jjm, llm, nqmx
    use dimphy, only: zmasq
    use dimsoil, only: nsoilmx
    use disvert_m, only: ap, bp, preff, pa, disvert
    use dynredem0_m, only: dynredem0
    use dynredem1_m, only: dynredem1
    use exner_hyb_m, only: exner_hyb
    use geopot_m, only: geopot
    use grid_atob, only: grille_m
    use grid_change, only: init_dyn_phy, dyn_phy
    use histclo_m, only: histclo
    use indicesol, only: is_oce, is_sic, is_ter, is_lic, epsfra
    use iniadvtrac_m, only: iniadvtrac
    use inifilr_m, only: inifilr
    use massdair_m, only: massdair
    use netcdf, only: nf90_nowrite
    use netcdf95, only: nf95_close, nf95_get_var, nf95_gw_var, &
         nf95_inq_varid, nf95_open
    use nr_util, only: pi, assert
    use paramet_m, only: ip1jm, ip1jmp1
    use phyredem_m, only: phyredem
    use q_sat_m, only: q_sat
    use regr_lat_time_coefoz_m, only: regr_lat_time_coefoz
    use regr_pr_o3_m, only: regr_pr_o3
    use serre, only: alphax
    use startdyn, only: start_init_dyn
    USE start_init_orog_m, only: start_init_orog, mask
    use start_init_phys_m, only: start_init_phys
    use start_inter_3d_m, only: start_inter_3d
    use temps, only: itau_phy, annee_ref, day_ref

    ! Variables local to the procedure:

    REAL latfi(klon), lonfi(klon)
    ! (latitude and longitude of a point of the scalar grid identified
    ! by a simple index, in degrees)

    REAL, dimension(iim + 1, jjm + 1, llm):: ucov, t3d, teta
    REAL vcov(iim + 1, jjm, llm)

    REAL q(iim + 1, jjm + 1, llm, nqmx)
    ! (mass fractions of trace species
    ! "q(i, j, l)" is at longitude "rlonv(i)", latitude "rlatu(j)"
    ! and pressure level "pls(i, j, l)".)

    real qsat(iim + 1, jjm + 1, llm) ! mass fraction of saturating water vapor
    REAL tsol(klon), qsol(klon), sn(klon)
    REAL tsolsrf(klon, nbsrf), qsolsrf(klon, nbsrf), snsrf(klon, nbsrf) 
    REAL albe(klon, nbsrf), evap(klon, nbsrf)
    REAL alblw(klon, nbsrf)
    REAL tsoil(klon, nsoilmx, nbsrf) 
    REAL radsol(klon), rain_fall(klon), snow_fall(klon)
    REAL solsw(klon), sollw(klon), fder(klon)
    !IM "slab" ocean
    REAL tslab(klon)
    real seaice(klon) ! kg m-2
    REAL frugs(klon, nbsrf), agesno(klon, nbsrf)
    REAL rugmer(klon)
    REAL phis(iim + 1, jjm + 1) ! surface geopotential, in m2 s-2
    real, dimension(iim + 1, jjm + 1):: zmea_2d, zstd_2d, zsig_2d, zgam_2d
    real, dimension(iim + 1, jjm + 1):: zthe_2d, zpic_2d, zval_2d
    real, dimension(iim + 1, jjm + 1):: tsol_2d, qsol_2d, ps
    REAL zmea(klon), zstd(klon)
    REAL zsig(klon), zgam(klon)
    REAL zthe(klon)
    REAL zpic(klon), zval(klon)
    REAL t_ancien(klon, llm), q_ancien(klon, llm)
    REAL run_off_lic_0(klon)
    real clwcon(klon, llm), rnebcon(klon, llm), ratqs(klon, llm)

    ! D\'eclarations pour lecture glace de mer :
    INTEGER iml_lic, jml_lic
    INTEGER ncid, varid
    REAL, pointer:: dlon_lic(:), dlat_lic(:)
    REAL, ALLOCATABLE:: fraclic(:, :) ! fraction land ice
    REAL flic_tmp(iim + 1, jjm + 1) ! fraction land ice temporary

    INTEGER l, ji

    REAL pk(iim + 1, jjm + 1, llm) ! fonction d'Exner aux milieux des couches 
    real pks(iim + 1, jjm + 1)

    REAL masse(iim + 1, jjm + 1, llm)
    REAL phi(iim + 1, jjm + 1, llm)
    REAL pbaru(ip1jmp1, llm), pbarv(ip1jm, llm)
    REAL w(iim + 1, jjm + 1, llm)

    real sig1(klon, llm) ! section adiabatic updraft
    real w01(klon, llm) ! vertical velocity within adiabatic updraft

    real pls(iim + 1, jjm + 1, llm)
    ! (pressure at mid-layer of LMDZ grid, in Pa)
    ! "pls(i, j, l)" is at longitude "rlonv(i)", latitude "rlatu(j)",
    ! for layer "l")

    REAL p3d(iim + 1, jjm + 1, llm+1) ! pressure at layer interfaces, in Pa
    ! ("p3d(i, j, l)" is at longitude "rlonv(i)", latitude "rlatu(j)",
    ! for interface "l")

    !---------------------------------

    print *, "Call sequence information: etat0"

    CALL iniconst

    ! Construct a grid:

    pa = 5e4
    CALL disvert
    CALL inigeom
    CALL inifilr

    latfi(1) = 90.
    latfi(2:klon-1) = pack(spread(rlatu(2:jjm), 1, iim), .true.) * 180. / pi
    ! (with conversion to degrees)
    latfi(klon) = - 90.

    lonfi(1) = 0.
    lonfi(2:klon-1) = pack(spread(rlonv(:iim), 2, jjm - 1), .true.) * 180. / pi
    ! (with conversion to degrees)
    lonfi(klon) = 0.

    call start_init_orog(phis, zmea_2d, zstd_2d, zsig_2d, zgam_2d, zthe_2d, &
         zpic_2d, zval_2d) ! also compute "mask"
    call init_dyn_phy ! define the mask "dyn_phy" for distinct grid points
    zmasq = pack(mask, dyn_phy)
    PRINT *, 'Masque construit'

    call start_init_phys(tsol_2d, qsol_2d)
    CALL start_init_dyn(tsol_2d, phis, ps)

    ! Compute pressure on intermediate levels:
    forall(l = 1: llm + 1) p3d(:, :, l) = ap(l) + bp(l) * ps
    CALL exner_hyb(ps, p3d, pks, pk)
    call assert(MINVAL(pk) /= MAXVAL(pk), '"pk" should not be constant')

    pls = preff * (pk / cpp)**(1. / kappa)
    PRINT *, "minval(pls) = ", minval(pls)
    print *, "maxval(pls) = ", maxval(pls)

    call start_inter_3d('U', rlonv, rlatv, pls, ucov)
    forall (l = 1: llm) ucov(:iim, :, l) = ucov(:iim, :, l) * cu_2d(:iim, :)
    ucov(iim+1, :, :) = ucov(1, :, :)

    call start_inter_3d('V', rlonu, rlatu(:jjm), pls(:, :jjm, :), vcov)
    forall (l = 1: llm) vcov(:iim, :, l) = vcov(:iim, :, l) * cv_2d(:iim, :)
    vcov(iim + 1, :, :) = vcov(1, :, :)

    call start_inter_3d('TEMP', rlonu, rlatv, pls, t3d)
    PRINT *, 'minval(t3d) = ', minval(t3d)
    print *, "maxval(t3d) = ", maxval(t3d)

    teta(:iim, :, :) = t3d(:iim, :, :) * cpp / pk(:iim, :, :)
    teta(iim + 1, :, :) = teta(1, :, :)
    DO l = 1, llm
       teta(:, 1, l) = SUM(aire_2d(:, 1) * teta(:, 1, l)) / apoln
       teta(:, jjm + 1, l) = SUM(aire_2d(:, jjm + 1) * teta(:, jjm + 1, l)) &
            / apols
    ENDDO

    ! Calcul de l'humidit\'e \`a saturation :
    qsat = q_sat(t3d, pls)
    PRINT *, "minval(qsat) = ", minval(qsat)
    print *, "maxval(qsat) = ", maxval(qsat)
    IF (MINVAL(qsat) == MAXVAL(qsat)) stop '"qsat" should not be constant'

    ! Water vapor:
    call start_inter_3d('R', rlonu, rlatv, pls, q(:, :, :, 1))
    q(:, :, :, 1) = 0.01 * q(:, :, :, 1) * qsat
    WHERE (q(:, :, :, 1) < 0.) q(:, :, :, 1) = 1E-10
    DO l = 1, llm
       q(:, 1, l, 1) = SUM(aire_2d(:, 1) * q(:, 1, l, 1)) / apoln
       q(:, jjm + 1, l, 1) &
            = SUM(aire_2d(:, jjm + 1) * q(:, jjm + 1, l, 1)) / apols
    ENDDO

    q(:, :, :, 2:4) = 0. ! liquid water, radon and lead

    if (nqmx >= 5) then
       ! Ozone:
       call regr_lat_time_coefoz
       call regr_pr_o3(p3d, q(:, :, :, 5))
       ! Convert from mole fraction to mass fraction:
       q(:, :, :, 5) = q(:, :, :, 5) * 48. / 29.
    end if

    tsol = pack(tsol_2d, dyn_phy)
    qsol = pack(qsol_2d, dyn_phy)
    sn = 0. ! snow
    radsol = 0.
    tslab = 0. ! IM "slab" ocean
    seaice = 0.
    rugmer = 0.001
    zmea = pack(zmea_2d, dyn_phy)
    zstd = pack(zstd_2d, dyn_phy)
    zsig = pack(zsig_2d, dyn_phy)
    zgam = pack(zgam_2d, dyn_phy)
    zthe = pack(zthe_2d, dyn_phy)
    zpic = pack(zpic_2d, dyn_phy)
    zval = pack(zval_2d, dyn_phy)

    ! On initialise les sous-surfaces.
    ! Lecture du fichier glace de terre pour fixer la fraction de terre 
    ! et de glace de terre :

    call nf95_open("landiceref.nc", nf90_nowrite, ncid)

    call nf95_inq_varid(ncid, 'longitude', varid)
    call nf95_gw_var(ncid, varid, dlon_lic)
    iml_lic = size(dlon_lic)

    call nf95_inq_varid(ncid, 'latitude', varid)
    call nf95_gw_var(ncid, varid, dlat_lic)
    jml_lic = size(dlat_lic)

    call nf95_inq_varid(ncid, 'landice', varid)
    ALLOCATE(fraclic(iml_lic, jml_lic))
    call nf95_get_var(ncid, varid, fraclic)

    call nf95_close(ncid)

    ! Interpolation sur la grille T du mod\`ele :
    PRINT *, 'Dimensions de "landiceref.nc"'
    print *, "iml_lic = ", iml_lic
    print *, "jml_lic = ", jml_lic

    ! Si les coordonn\'ees sont en degr\'es, on les transforme :
    IF (MAXVAL(dlon_lic) > pi) THEN
       dlon_lic = dlon_lic * pi / 180.
    ENDIF
    IF (maxval(dlat_lic) > pi) THEN 
       dlat_lic = dlat_lic * pi/ 180.
    ENDIF

    flic_tmp(:iim, :) = grille_m(dlon_lic, dlat_lic, fraclic, rlonv(:iim), &
         rlatu)
    flic_tmp(iim + 1, :) = flic_tmp(1, :)

    deallocate(dlon_lic, dlat_lic) ! pointers

    ! Passage sur la grille physique
    pctsrf = 0.
    pctsrf(:, is_lic) = pack(flic_tmp, dyn_phy)
    ! Ad\'equation avec le maque terre/mer
    WHERE (pctsrf(:, is_lic) < EPSFRA) pctsrf(:, is_lic) = 0.
    WHERE (zmasq < EPSFRA) pctsrf(:, is_lic) = 0.
    pctsrf(:, is_ter) = zmasq
    where (zmasq > EPSFRA)
       where (pctsrf(:, is_lic) >= zmasq)
          pctsrf(:, is_lic) = zmasq
          pctsrf(:, is_ter) = 0.
       elsewhere
          pctsrf(:, is_ter) = zmasq - pctsrf(:, is_lic)
          where (pctsrf(:, is_ter) < EPSFRA)
             pctsrf(:, is_ter) = 0.
             pctsrf(:, is_lic) = zmasq
          end where
       end where
    end where

    ! Sous-surface oc\'ean et glace de mer (pour d\'emarrer on met glace
    ! de mer \`a 0) :
    pctsrf(:, is_oce) = 1. - zmasq
    WHERE (pctsrf(:, is_oce) < EPSFRA) pctsrf(:, is_oce) = 0.

    ! V\'erification que somme des sous-surfaces vaut 1 :
    ji = count(abs(sum(pctsrf, dim = 2) - 1.) > EPSFRA)
    IF (ji /= 0) then
       PRINT *, 'Bad surface percentages for ', ji, 'points'
    end IF

    ! Calcul interm\'ediaire :
    CALL massdair(p3d, masse)

    print *, 'ALPHAX = ', alphax

    forall (l = 1:llm)
       masse(:, 1, l) = SUM(aire_2d(:iim, 1) * masse(:iim, 1, l)) / apoln
       masse(:, jjm + 1, l) = &
            SUM(aire_2d(:iim, jjm + 1) * masse(:iim, jjm + 1, l)) / apols
    END forall

    ! Initialisation pour traceurs:
    call iniadvtrac
    itau_phy = 0
    day_ref = dayref
    annee_ref = anneeref

    CALL geopot(teta, pk , pks, phis, phi)
    CALL caldyn0(ucov, vcov, teta, ps, masse, pk, phis, phi, w, pbaru, &
         pbarv)
    CALL dynredem0("start.nc", dayref, phis)
    CALL dynredem1("start.nc", vcov, ucov, teta, q, masse, ps, itau=0)

    ! Initialisations :
    tsolsrf(:, is_ter) = tsol
    tsolsrf(:, is_lic) = tsol
    tsolsrf(:, is_oce) = tsol
    tsolsrf(:, is_sic) = tsol
    snsrf(:, is_ter) = sn
    snsrf(:, is_lic) = sn
    snsrf(:, is_oce) = sn
    snsrf(:, is_sic) = sn
    albe(:, is_ter) = 0.08
    albe(:, is_lic) = 0.6
    albe(:, is_oce) = 0.5
    albe(:, is_sic) = 0.6
    alblw = albe
    evap = 0.
    qsolsrf(:, is_ter) = 150.
    qsolsrf(:, is_lic) = 150.
    qsolsrf(:, is_oce) = 150.
    qsolsrf(:, is_sic) = 150.
    tsoil = spread(spread(tsol, 2, nsoilmx), 3, nbsrf)
    rain_fall = 0.
    snow_fall = 0.
    solsw = 165.
    sollw = -53.
    t_ancien = 273.15
    q_ancien = 0.
    agesno = 0.
    !IM "slab" ocean
    tslab = tsolsrf(:, is_oce)
    seaice = 0.

    frugs(:, is_oce) = rugmer
    frugs(:, is_ter) = MAX(1e-5, zstd * zsig / 2)
    frugs(:, is_lic) = MAX(1e-5, zstd * zsig / 2)
    frugs(:, is_sic) = 0.001
    fder = 0.
    clwcon = 0.
    rnebcon = 0.
    ratqs = 0.
    run_off_lic_0 = 0.
    sig1 = 0.
    w01 = 0.

    call phyredem("startphy.nc", latfi, lonfi, pctsrf, &
         tsolsrf, tsoil, tslab, seaice, qsolsrf, qsol, snsrf, albe, alblw, &
         evap, rain_fall, snow_fall, solsw, sollw, fder, radsol, frugs, &
         agesno, zmea, zstd, zsig, zgam, zthe, zpic, zval, &
         t_ancien, q_ancien, rnebcon, ratqs, clwcon, run_off_lic_0, sig1, w01)
    CALL histclo

  END SUBROUTINE etat0

end module etat0_mod
1	guez	3	module etat0_mod
2
3			use indicesol, only: nbsrf
4			use dimphy, only: klon
5
6			IMPLICIT NONE
7
8			REAL pctsrf(klon, nbsrf)
9	guez	15	! ("pctsrf(i, :)" is the composition of the surface at horizontal
10	guez	78	! position "i")
11	guez	3
12			private nbsrf, klon
13
14			contains
15
16			SUBROUTINE etat0
17
18	guez	90	! From "etat0_netcdf.F", version 1.3, 2005/05/25 13:10:09
19	guez	3
20	guez	27	use caldyn0_m, only: caldyn0
21	guez	79	use comconst, only: cpp, kappa, iniconst
22	guez	27	use comgeom, only: rlatu, rlonv, rlonu, rlatv, aire_2d, apoln, apols, &
23	guez	78	cu_2d, cv_2d, inigeom
24	guez	79	use conf_gcm_m, only: dayref, anneeref
25	guez	27	use dimens_m, only: iim, jjm, llm, nqmx
26	guez	3	use dimphy, only: zmasq
27			use dimsoil, only: nsoilmx
28	guez	79	use disvert_m, only: ap, bp, preff, pa, disvert
29	guez	27	use dynredem0_m, only: dynredem0
30			use dynredem1_m, only: dynredem1
31			use exner_hyb_m, only: exner_hyb
32	guez	43	use geopot_m, only: geopot
33	guez	3	use grid_atob, only: grille_m
34			use grid_change, only: init_dyn_phy, dyn_phy
35	guez	61	use histclo_m, only: histclo
36	guez	27	use indicesol, only: is_oce, is_sic, is_ter, is_lic, epsfra
37	guez	18	use iniadvtrac_m, only: iniadvtrac
38	guez	54	use inifilr_m, only: inifilr
39	guez	67	use massdair_m, only: massdair
40	guez	48	use netcdf, only: nf90_nowrite
41	guez	68	use netcdf95, only: nf95_close, nf95_get_var, nf95_gw_var, &
42			nf95_inq_varid, nf95_open
43	guez	90	use nr_util, only: pi, assert
44	guez	27	use paramet_m, only: ip1jm, ip1jmp1
45			use phyredem_m, only: phyredem
46			use q_sat_m, only: q_sat
47	guez	7	use regr_lat_time_coefoz_m, only: regr_lat_time_coefoz
48			use regr_pr_o3_m, only: regr_pr_o3
49	guez	27	use serre, only: alphax
50	guez	68	use startdyn, only: start_init_dyn
51	guez	78	USE start_init_orog_m, only: start_init_orog, mask
52	guez	43	use start_init_phys_m, only: start_init_phys
53	guez	54	use start_inter_3d_m, only: start_inter_3d
54	guez	28	use temps, only: itau_phy, annee_ref, day_ref
55	guez	3
56			! Variables local to the procedure:
57
58			REAL latfi(klon), lonfi(klon)
59			! (latitude and longitude of a point of the scalar grid identified
60	guez	90	! by a simple index, in degrees)
61	guez	3
62	guez	73	REAL, dimension(iim + 1, jjm + 1, llm):: ucov, t3d, teta
63	guez	65	REAL vcov(iim + 1, jjm, llm)
64	guez	3
65	guez	68	REAL q(iim + 1, jjm + 1, llm, nqmx)
66	guez	3	! (mass fractions of trace species
67	guez	68	! "q(i, j, l)" is at longitude "rlonv(i)", latitude "rlatu(j)"
68	guez	3	! and pressure level "pls(i, j, l)".)
69
70			real qsat(iim + 1, jjm + 1, llm) ! mass fraction of saturating water vapor
71			REAL tsol(klon), qsol(klon), sn(klon)
72			REAL tsolsrf(klon, nbsrf), qsolsrf(klon, nbsrf), snsrf(klon, nbsrf)
73			REAL albe(klon, nbsrf), evap(klon, nbsrf)
74			REAL alblw(klon, nbsrf)
75			REAL tsoil(klon, nsoilmx, nbsrf)
76			REAL radsol(klon), rain_fall(klon), snow_fall(klon)
77			REAL solsw(klon), sollw(klon), fder(klon)
78			!IM "slab" ocean
79			REAL tslab(klon)
80			real seaice(klon) ! kg m-2
81			REAL frugs(klon, nbsrf), agesno(klon, nbsrf)
82			REAL rugmer(klon)
83	guez	78	REAL phis(iim + 1, jjm + 1) ! surface geopotential, in m2 s-2
84			real, dimension(iim + 1, jjm + 1):: zmea_2d, zstd_2d, zsig_2d, zgam_2d
85	guez	3	real, dimension(iim + 1, jjm + 1):: zthe_2d, zpic_2d, zval_2d
86	guez	73	real, dimension(iim + 1, jjm + 1):: tsol_2d, qsol_2d, ps
87	guez	3	REAL zmea(klon), zstd(klon)
88			REAL zsig(klon), zgam(klon)
89			REAL zthe(klon)
90			REAL zpic(klon), zval(klon)
91	guez	78	REAL t_ancien(klon, llm), q_ancien(klon, llm)
92	guez	3	REAL run_off_lic_0(klon)
93			real clwcon(klon, llm), rnebcon(klon, llm), ratqs(klon, llm)
94	guez	68
95	guez	90	! D\'eclarations pour lecture glace de mer :
96	guez	68	INTEGER iml_lic, jml_lic
97			INTEGER ncid, varid
98			REAL, pointer:: dlon_lic(:), dlat_lic(:)
99	guez	3	REAL, ALLOCATABLE:: fraclic(:, :) ! fraction land ice
100	guez	68	REAL flic_tmp(iim + 1, jjm + 1) ! fraction land ice temporary
101	guez	3
102			INTEGER l, ji
103
104			REAL pk(iim + 1, jjm + 1, llm) ! fonction d'Exner aux milieux des couches
105			real pks(iim + 1, jjm + 1)
106
107			REAL masse(iim + 1, jjm + 1, llm)
108	guez	70	REAL phi(iim + 1, jjm + 1, llm)
109	guez	3	REAL pbaru(ip1jmp1, llm), pbarv(ip1jm, llm)
110	guez	91	REAL w(iim + 1, jjm + 1, llm)
111	guez	3
112	guez	72	real sig1(klon, llm) ! section adiabatic updraft
113			real w01(klon, llm) ! vertical velocity within adiabatic updraft
114
115	guez	97	real pls(iim + 1, jjm + 1, llm)
116			! (pressure at mid-layer of LMDZ grid, in Pa)
117			! "pls(i, j, l)" is at longitude "rlonv(i)", latitude "rlatu(j)",
118			! for layer "l")
119
120			REAL p3d(iim + 1, jjm + 1, llm+1) ! pressure at layer interfaces, in Pa
121			! ("p3d(i, j, l)" is at longitude "rlonv(i)", latitude "rlatu(j)",
122			! for interface "l")
123
124	guez	3	!---------------------------------
125
126			print *, "Call sequence information: etat0"
127
128	guez	79	CALL iniconst
129	guez	3
130	guez	36	! Construct a grid:
131
132	guez	3	pa = 5e4
133	guez	79	CALL disvert
134	guez	3	CALL inigeom
135			CALL inifilr
136
137			latfi(1) = 90.
138			latfi(2:klon-1) = pack(spread(rlatu(2:jjm), 1, iim), .true.) * 180. / pi
139			! (with conversion to degrees)
140			latfi(klon) = - 90.
141
142			lonfi(1) = 0.
143			lonfi(2:klon-1) = pack(spread(rlonv(:iim), 2, jjm - 1), .true.) * 180. / pi
144			! (with conversion to degrees)
145			lonfi(klon) = 0.
146
147	guez	78	call start_init_orog(phis, zmea_2d, zstd_2d, zsig_2d, zgam_2d, zthe_2d, &
148			zpic_2d, zval_2d) ! also compute "mask"
149	guez	3	call init_dyn_phy ! define the mask "dyn_phy" for distinct grid points
150	guez	15	zmasq = pack(mask, dyn_phy)
151	guez	3	PRINT *, 'Masque construit'
152
153	guez	43	call start_init_phys(tsol_2d, qsol_2d)
154	guez	78	CALL start_init_dyn(tsol_2d, phis, ps)
155	guez	3
156			! Compute pressure on intermediate levels:
157	guez	73	forall(l = 1: llm + 1) p3d(:, :, l) = ap(l) + bp(l) * ps
158			CALL exner_hyb(ps, p3d, pks, pk)
159	guez	90	call assert(MINVAL(pk) /= MAXVAL(pk), '"pk" should not be constant')
160	guez	3
161	guez	36	pls = preff * (pk / cpp)**(1. / kappa)
162			PRINT *, "minval(pls) = ", minval(pls)
163			print *, "maxval(pls) = ", maxval(pls)
164	guez	3
165	guez	65	call start_inter_3d('U', rlonv, rlatv, pls, ucov)
166			forall (l = 1: llm) ucov(:iim, :, l) = ucov(:iim, :, l) * cu_2d(:iim, :)
167			ucov(iim+1, :, :) = ucov(1, :, :)
168	guez	3
169	guez	65	call start_inter_3d('V', rlonu, rlatu(:jjm), pls(:, :jjm, :), vcov)
170			forall (l = 1: llm) vcov(:iim, :, l) = vcov(:iim, :, l) * cv_2d(:iim, :)
171			vcov(iim + 1, :, :) = vcov(1, :, :)
172	guez	3
173	guez	23	call start_inter_3d('TEMP', rlonu, rlatv, pls, t3d)
174	guez	78	PRINT *, 'minval(t3d) = ', minval(t3d)
175	guez	36	print *, "maxval(t3d) = ", maxval(t3d)
176	guez	3
177	guez	73	teta(:iim, :, :) = t3d(:iim, :, :) * cpp / pk(:iim, :, :)
178			teta(iim + 1, :, :) = teta(1, :, :)
179			DO l = 1, llm
180			teta(:, 1, l) = SUM(aire_2d(:, 1) * teta(:, 1, l)) / apoln
181			teta(:, jjm + 1, l) = SUM(aire_2d(:, jjm + 1) * teta(:, jjm + 1, l)) &
182	guez	3	/ apols
183			ENDDO
184
185	guez	90	! Calcul de l'humidit\'e \`a saturation :
186	guez	36	qsat = q_sat(t3d, pls)
187			PRINT *, "minval(qsat) = ", minval(qsat)
188			print *, "maxval(qsat) = ", maxval(qsat)
189	guez	3	IF (MINVAL(qsat) == MAXVAL(qsat)) stop '"qsat" should not be constant'
190
191			! Water vapor:
192	guez	68	call start_inter_3d('R', rlonu, rlatv, pls, q(:, :, :, 1))
193			q(:, :, :, 1) = 0.01 * q(:, :, :, 1) * qsat
194			WHERE (q(:, :, :, 1) < 0.) q(:, :, :, 1) = 1E-10
195	guez	3	DO l = 1, llm
196	guez	68	q(:, 1, l, 1) = SUM(aire_2d(:, 1) * q(:, 1, l, 1)) / apoln
197			q(:, jjm + 1, l, 1) &
198			= SUM(aire_2d(:, jjm + 1) * q(:, jjm + 1, l, 1)) / apols
199	guez	3	ENDDO
200
201	guez	68	q(:, :, :, 2:4) = 0. ! liquid water, radon and lead
202	guez	3
203	guez	5	if (nqmx >= 5) then
204			! Ozone:
205	guez	7	call regr_lat_time_coefoz
206	guez	97	call regr_pr_o3(p3d, q(:, :, :, 5))
207	guez	7	! Convert from mole fraction to mass fraction:
208	guez	78	q(:, :, :, 5) = q(:, :, :, 5) * 48. / 29.
209	guez	5	end if
210	guez	3
211	guez	13	tsol = pack(tsol_2d, dyn_phy)
212			qsol = pack(qsol_2d, dyn_phy)
213			sn = 0. ! snow
214			radsol = 0.
215			tslab = 0. ! IM "slab" ocean
216			seaice = 0.
217			rugmer = 0.001
218	guez	78	zmea = pack(zmea_2d, dyn_phy)
219	guez	13	zstd = pack(zstd_2d, dyn_phy)
220			zsig = pack(zsig_2d, dyn_phy)
221			zgam = pack(zgam_2d, dyn_phy)
222			zthe = pack(zthe_2d, dyn_phy)
223			zpic = pack(zpic_2d, dyn_phy)
224			zval = pack(zval_2d, dyn_phy)
225	guez	3
226	guez	49	! On initialise les sous-surfaces.
227	guez	3	! Lecture du fichier glace de terre pour fixer la fraction de terre
228	guez	49	! et de glace de terre :
229	guez	68
230	guez	48	call nf95_open("landiceref.nc", nf90_nowrite, ncid)
231	guez	68
232			call nf95_inq_varid(ncid, 'longitude', varid)
233			call nf95_gw_var(ncid, varid, dlon_lic)
234			iml_lic = size(dlon_lic)
235
236			call nf95_inq_varid(ncid, 'latitude', varid)
237			call nf95_gw_var(ncid, varid, dlat_lic)
238			jml_lic = size(dlat_lic)
239
240	guez	48	call nf95_inq_varid(ncid, 'landice', varid)
241	guez	68	ALLOCATE(fraclic(iml_lic, jml_lic))
242	guez	48	call nf95_get_var(ncid, varid, fraclic)
243	guez	68
244	guez	48	call nf95_close(ncid)
245	guez	3
246	guez	90	! Interpolation sur la grille T du mod\`ele :
247	guez	68	PRINT *, 'Dimensions de "landiceref.nc"'
248	guez	3	print *, "iml_lic = ", iml_lic
249			print *, "jml_lic = ", jml_lic
250
251	guez	90	! Si les coordonn\'ees sont en degr\'es, on les transforme :
252	guez	78	IF (MAXVAL(dlon_lic) > pi) THEN
253	guez	68	dlon_lic = dlon_lic * pi / 180.
254	guez	3	ENDIF
255	guez	73	IF (maxval(dlat_lic) > pi) THEN
256	guez	68	dlat_lic = dlat_lic * pi/ 180.
257	guez	3	ENDIF
258
259			flic_tmp(:iim, :) = grille_m(dlon_lic, dlat_lic, fraclic, rlonv(:iim), &
260			rlatu)
261			flic_tmp(iim + 1, :) = flic_tmp(1, :)
262
263	guez	68	deallocate(dlon_lic, dlat_lic) ! pointers
264
265	guez	3	! Passage sur la grille physique
266	guez	15	pctsrf = 0.
267	guez	3	pctsrf(:, is_lic) = pack(flic_tmp, dyn_phy)
268	guez	90	! Ad\'equation avec le maque terre/mer
269	guez	73	WHERE (pctsrf(:, is_lic) < EPSFRA) pctsrf(:, is_lic) = 0.
270	guez	13	WHERE (zmasq < EPSFRA) pctsrf(:, is_lic) = 0.
271			pctsrf(:, is_ter) = zmasq
272			where (zmasq > EPSFRA)
273			where (pctsrf(:, is_lic) >= zmasq)
274			pctsrf(:, is_lic) = zmasq
275	guez	3	pctsrf(:, is_ter) = 0.
276			elsewhere
277	guez	13	pctsrf(:, is_ter) = zmasq - pctsrf(:, is_lic)
278	guez	3	where (pctsrf(:, is_ter) < EPSFRA)
279			pctsrf(:, is_ter) = 0.
280	guez	13	pctsrf(:, is_lic) = zmasq
281	guez	3	end where
282			end where
283			end where
284
285	guez	90	! Sous-surface oc\'ean et glace de mer (pour d\'emarrer on met glace
286			! de mer \`a 0) :
287	guez	13	pctsrf(:, is_oce) = 1. - zmasq
288	guez	3	WHERE (pctsrf(:, is_oce) < EPSFRA) pctsrf(:, is_oce) = 0.
289
290	guez	90	! V\'erification que somme des sous-surfaces vaut 1 :
291	guez	15	ji = count(abs(sum(pctsrf, dim = 2) - 1.) > EPSFRA)
292			IF (ji /= 0) then
293	guez	91	PRINT *, 'Bad surface percentages for ', ji, 'points'
294	guez	15	end IF
295	guez	3
296	guez	90	! Calcul interm\'ediaire :
297	guez	3	CALL massdair(p3d, masse)
298
299			print *, 'ALPHAX = ', alphax
300
301	guez	78	forall (l = 1:llm)
302	guez	3	masse(:, 1, l) = SUM(aire_2d(:iim, 1) * masse(:iim, 1, l)) / apoln
303			masse(:, jjm + 1, l) = &
304			SUM(aire_2d(:iim, jjm + 1) * masse(:iim, jjm + 1, l)) / apols
305			END forall
306
307			! Initialisation pour traceurs:
308	guez	5	call iniadvtrac
309	guez	3	itau_phy = 0
310			day_ref = dayref
311			annee_ref = anneeref
312
313	guez	78	CALL geopot(teta, pk , pks, phis, phi)
314	guez	73	CALL caldyn0(ucov, vcov, teta, ps, masse, pk, phis, phi, w, pbaru, &
315	guez	23	pbarv)
316	guez	5	CALL dynredem0("start.nc", dayref, phis)
317	guez	73	CALL dynredem1("start.nc", vcov, ucov, teta, q, masse, ps, itau=0)
318	guez	3
319			! Initialisations :
320			tsolsrf(:, is_ter) = tsol
321			tsolsrf(:, is_lic) = tsol
322			tsolsrf(:, is_oce) = tsol
323			tsolsrf(:, is_sic) = tsol
324			snsrf(:, is_ter) = sn
325			snsrf(:, is_lic) = sn
326			snsrf(:, is_oce) = sn
327			snsrf(:, is_sic) = sn
328			albe(:, is_ter) = 0.08
329			albe(:, is_lic) = 0.6
330			albe(:, is_oce) = 0.5
331			albe(:, is_sic) = 0.6
332			alblw = albe
333	guez	15	evap = 0.
334	guez	3	qsolsrf(:, is_ter) = 150.
335			qsolsrf(:, is_lic) = 150.
336			qsolsrf(:, is_oce) = 150.
337			qsolsrf(:, is_sic) = 150.
338			tsoil = spread(spread(tsol, 2, nsoilmx), 3, nbsrf)
339			rain_fall = 0.
340			snow_fall = 0.
341			solsw = 165.
342			sollw = -53.
343			t_ancien = 273.15
344			q_ancien = 0.
345			agesno = 0.
346			!IM "slab" ocean
347	guez	13	tslab = tsolsrf(:, is_oce)
348	guez	3	seaice = 0.
349
350	guez	13	frugs(:, is_oce) = rugmer
351	guez	78	frugs(:, is_ter) = MAX(1e-5, zstd * zsig / 2)
352			frugs(:, is_lic) = MAX(1e-5, zstd * zsig / 2)
353	guez	3	frugs(:, is_sic) = 0.001
354			fder = 0.
355			clwcon = 0.
356			rnebcon = 0.
357			ratqs = 0.
358			run_off_lic_0 = 0.
359	guez	72	sig1 = 0.
360			w01 = 0.
361	guez	3
362	guez	13	call phyredem("startphy.nc", latfi, lonfi, pctsrf, &
363	guez	3	tsolsrf, tsoil, tslab, seaice, qsolsrf, qsol, snsrf, albe, alblw, &
364			evap, rain_fall, snow_fall, solsw, sollw, fder, radsol, frugs, &
365	guez	13	agesno, zmea, zstd, zsig, zgam, zthe, zpic, zval, &
366	guez	72	t_ancien, q_ancien, rnebcon, ratqs, clwcon, run_off_lic_0, sig1, w01)
367	guez	3	CALL histclo
368
369			END SUBROUTINE etat0
370
371			end module etat0_mod