Sources/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: dtvr, daysec, cpp, kappa
    use comgeom, only: rlatu, rlonv, rlonu, rlatv, aire_2d, apoln, apols, &
         cu_2d, cv_2d
    use conf_gcm_m, only: day_step, iphysiq, 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
    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 inigeom_m, only: inigeom
    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
    use paramet_m, only: ip1jm, ip1jmp1
    use phyredem_m, only: phyredem
    use pressure_var, only: pls, p3d
    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, phis
    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 °)

    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, dimension(iim + 1, jjm + 1):: relief, 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éclarations 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(ip1jmp1, llm)
    REAL phystep

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

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

    print *, "Call sequence information: etat0"

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

    ! Construct a grid:

    pa = 5e4
    CALL iniconst
    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(relief, zstd_2d, zsig_2d, zgam_2d, zthe_2d, zpic_2d, &
         zval_2d) ! also compute "mask" and "phis"
    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, 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)
    IF (MINVAL(pk) == MAXVAL(pk)) then
       print *, '"pk" should not be constant'
       stop 1
    end IF

    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é à 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(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(relief, 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èle :
    PRINT *, 'Dimensions de "landiceref.nc"'
    print *, "iml_lic = ", iml_lic
    print *, "jml_lic = ", jml_lic

    ! Si les coordonnées sont en degrés, 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équation 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éan et glace de mer (pour démarrer on met glace
    ! de mer à 0) :
    pctsrf(:, is_oce) = 1. - zmasq
    WHERE (pctsrf(:, is_oce) < EPSFRA) pctsrf(:, is_oce) = 0.

    ! Vérification que somme des sous-surfaces vaut 1 :
    ji = count(abs(sum(pctsrf, dim = 2) - 1.) > EPSFRA)
    IF (ji /= 0) then
       PRINT *, 'Problème répartition sous maille pour ', ji, 'points'
    end IF

    ! Calcul intermédiaire :
    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)

    ! Ecriture état initial physique:
    print *, "iphysiq = ", iphysiq
    phystep   = dtvr * REAL(iphysiq)
    print *, 'phystep = ', phystep

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