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: 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
    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
    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 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é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 sig1(klon, llm) ! section adiabatic updraft
    real w01(klon, llm) ! vertical velocity within adiabatic updraft

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

    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)
    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(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è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)

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