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, assert
    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 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(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)
    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(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 *, 'Probl\`eme r\'epartition sous maille pour ', 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	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	guez	90	! by a simple index, in degrees)
62	guez	3
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	guez	90	! D\'eclarations pour lecture glace de mer :
97	guez	68	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	90	call assert(MINVAL(pk) /= MAXVAL(pk), '"pk" should not be constant')
152	guez	3
153	guez	36	pls = preff * (pk / cpp)**(1. / kappa)
154			PRINT *, "minval(pls) = ", minval(pls)
155			print *, "maxval(pls) = ", maxval(pls)
156	guez	3
157	guez	65	call start_inter_3d('U', rlonv, rlatv, pls, ucov)
158			forall (l = 1: llm) ucov(:iim, :, l) = ucov(:iim, :, l) * cu_2d(:iim, :)
159			ucov(iim+1, :, :) = ucov(1, :, :)
160	guez	3
161	guez	65	call start_inter_3d('V', rlonu, rlatu(:jjm), pls(:, :jjm, :), vcov)
162			forall (l = 1: llm) vcov(:iim, :, l) = vcov(:iim, :, l) * cv_2d(:iim, :)
163			vcov(iim + 1, :, :) = vcov(1, :, :)
164	guez	3
165	guez	23	call start_inter_3d('TEMP', rlonu, rlatv, pls, t3d)
166	guez	78	PRINT *, 'minval(t3d) = ', minval(t3d)
167	guez	36	print *, "maxval(t3d) = ", maxval(t3d)
168	guez	3
169	guez	73	teta(:iim, :, :) = t3d(:iim, :, :) * cpp / pk(:iim, :, :)
170			teta(iim + 1, :, :) = teta(1, :, :)
171			DO l = 1, llm
172			teta(:, 1, l) = SUM(aire_2d(:, 1) * teta(:, 1, l)) / apoln
173			teta(:, jjm + 1, l) = SUM(aire_2d(:, jjm + 1) * teta(:, jjm + 1, l)) &
174	guez	3	/ apols
175			ENDDO
176
177	guez	90	! Calcul de l'humidit\'e \`a saturation :
178	guez	36	qsat = q_sat(t3d, pls)
179			PRINT *, "minval(qsat) = ", minval(qsat)
180			print *, "maxval(qsat) = ", maxval(qsat)
181	guez	3	IF (MINVAL(qsat) == MAXVAL(qsat)) stop '"qsat" should not be constant'
182
183			! Water vapor:
184	guez	68	call start_inter_3d('R', rlonu, rlatv, pls, q(:, :, :, 1))
185			q(:, :, :, 1) = 0.01 * q(:, :, :, 1) * qsat
186			WHERE (q(:, :, :, 1) < 0.) q(:, :, :, 1) = 1E-10
187	guez	3	DO l = 1, llm
188	guez	68	q(:, 1, l, 1) = SUM(aire_2d(:, 1) * q(:, 1, l, 1)) / apoln
189			q(:, jjm + 1, l, 1) &
190			= SUM(aire_2d(:, jjm + 1) * q(:, jjm + 1, l, 1)) / apols
191	guez	3	ENDDO
192
193	guez	68	q(:, :, :, 2:4) = 0. ! liquid water, radon and lead
194	guez	3
195	guez	5	if (nqmx >= 5) then
196			! Ozone:
197	guez	7	call regr_lat_time_coefoz
198	guez	68	call regr_pr_o3(q(:, :, :, 5))
199	guez	7	! Convert from mole fraction to mass fraction:
200	guez	78	q(:, :, :, 5) = q(:, :, :, 5) * 48. / 29.
201	guez	5	end if
202	guez	3
203	guez	13	tsol = pack(tsol_2d, dyn_phy)
204			qsol = pack(qsol_2d, dyn_phy)
205			sn = 0. ! snow
206			radsol = 0.
207			tslab = 0. ! IM "slab" ocean
208			seaice = 0.
209			rugmer = 0.001
210	guez	78	zmea = pack(zmea_2d, dyn_phy)
211	guez	13	zstd = pack(zstd_2d, dyn_phy)
212			zsig = pack(zsig_2d, dyn_phy)
213			zgam = pack(zgam_2d, dyn_phy)
214			zthe = pack(zthe_2d, dyn_phy)
215			zpic = pack(zpic_2d, dyn_phy)
216			zval = pack(zval_2d, dyn_phy)
217	guez	3
218	guez	49	! On initialise les sous-surfaces.
219	guez	3	! Lecture du fichier glace de terre pour fixer la fraction de terre
220	guez	49	! et de glace de terre :
221	guez	68
222	guez	48	call nf95_open("landiceref.nc", nf90_nowrite, ncid)
223	guez	68
224			call nf95_inq_varid(ncid, 'longitude', varid)
225			call nf95_gw_var(ncid, varid, dlon_lic)
226			iml_lic = size(dlon_lic)
227
228			call nf95_inq_varid(ncid, 'latitude', varid)
229			call nf95_gw_var(ncid, varid, dlat_lic)
230			jml_lic = size(dlat_lic)
231
232	guez	48	call nf95_inq_varid(ncid, 'landice', varid)
233	guez	68	ALLOCATE(fraclic(iml_lic, jml_lic))
234	guez	48	call nf95_get_var(ncid, varid, fraclic)
235	guez	68
236	guez	48	call nf95_close(ncid)
237	guez	3
238	guez	90	! Interpolation sur la grille T du mod\`ele :
239	guez	68	PRINT *, 'Dimensions de "landiceref.nc"'
240	guez	3	print *, "iml_lic = ", iml_lic
241			print *, "jml_lic = ", jml_lic
242
243	guez	90	! Si les coordonn\'ees sont en degr\'es, on les transforme :
244	guez	78	IF (MAXVAL(dlon_lic) > pi) THEN
245	guez	68	dlon_lic = dlon_lic * pi / 180.
246	guez	3	ENDIF
247	guez	73	IF (maxval(dlat_lic) > pi) THEN
248	guez	68	dlat_lic = dlat_lic * pi/ 180.
249	guez	3	ENDIF
250
251			flic_tmp(:iim, :) = grille_m(dlon_lic, dlat_lic, fraclic, rlonv(:iim), &
252			rlatu)
253			flic_tmp(iim + 1, :) = flic_tmp(1, :)
254
255	guez	68	deallocate(dlon_lic, dlat_lic) ! pointers
256
257	guez	3	! Passage sur la grille physique
258	guez	15	pctsrf = 0.
259	guez	3	pctsrf(:, is_lic) = pack(flic_tmp, dyn_phy)
260	guez	90	! Ad\'equation avec le maque terre/mer
261	guez	73	WHERE (pctsrf(:, is_lic) < EPSFRA) pctsrf(:, is_lic) = 0.
262	guez	13	WHERE (zmasq < EPSFRA) pctsrf(:, is_lic) = 0.
263			pctsrf(:, is_ter) = zmasq
264			where (zmasq > EPSFRA)
265			where (pctsrf(:, is_lic) >= zmasq)
266			pctsrf(:, is_lic) = zmasq
267	guez	3	pctsrf(:, is_ter) = 0.
268			elsewhere
269	guez	13	pctsrf(:, is_ter) = zmasq - pctsrf(:, is_lic)
270	guez	3	where (pctsrf(:, is_ter) < EPSFRA)
271			pctsrf(:, is_ter) = 0.
272	guez	13	pctsrf(:, is_lic) = zmasq
273	guez	3	end where
274			end where
275			end where
276
277	guez	90	! Sous-surface oc\'ean et glace de mer (pour d\'emarrer on met glace
278			! de mer \`a 0) :
279	guez	13	pctsrf(:, is_oce) = 1. - zmasq
280	guez	3	WHERE (pctsrf(:, is_oce) < EPSFRA) pctsrf(:, is_oce) = 0.
281
282	guez	90	! V\'erification que somme des sous-surfaces vaut 1 :
283	guez	15	ji = count(abs(sum(pctsrf, dim = 2) - 1.) > EPSFRA)
284			IF (ji /= 0) then
285	guez	90	PRINT *, 'Probl\`eme r\'epartition sous maille pour ', ji, 'points'
286	guez	15	end IF
287	guez	3
288	guez	90	! Calcul interm\'ediaire :
289	guez	3	CALL massdair(p3d, masse)
290
291			print *, 'ALPHAX = ', alphax
292
293	guez	78	forall (l = 1:llm)
294	guez	3	masse(:, 1, l) = SUM(aire_2d(:iim, 1) * masse(:iim, 1, l)) / apoln
295			masse(:, jjm + 1, l) = &
296			SUM(aire_2d(:iim, jjm + 1) * masse(:iim, jjm + 1, l)) / apols
297			END forall
298
299			! Initialisation pour traceurs:
300	guez	5	call iniadvtrac
301	guez	3	itau_phy = 0
302			day_ref = dayref
303			annee_ref = anneeref
304
305	guez	78	CALL geopot(teta, pk , pks, phis, phi)
306	guez	73	CALL caldyn0(ucov, vcov, teta, ps, masse, pk, phis, phi, w, pbaru, &
307	guez	23	pbarv)
308	guez	5	CALL dynredem0("start.nc", dayref, phis)
309	guez	73	CALL dynredem1("start.nc", vcov, ucov, teta, q, masse, ps, itau=0)
310	guez	3
311			! Initialisations :
312			tsolsrf(:, is_ter) = tsol
313			tsolsrf(:, is_lic) = tsol
314			tsolsrf(:, is_oce) = tsol
315			tsolsrf(:, is_sic) = tsol
316			snsrf(:, is_ter) = sn
317			snsrf(:, is_lic) = sn
318			snsrf(:, is_oce) = sn
319			snsrf(:, is_sic) = sn
320			albe(:, is_ter) = 0.08
321			albe(:, is_lic) = 0.6
322			albe(:, is_oce) = 0.5
323			albe(:, is_sic) = 0.6
324			alblw = albe
325	guez	15	evap = 0.
326	guez	3	qsolsrf(:, is_ter) = 150.
327			qsolsrf(:, is_lic) = 150.
328			qsolsrf(:, is_oce) = 150.
329			qsolsrf(:, is_sic) = 150.
330			tsoil = spread(spread(tsol, 2, nsoilmx), 3, nbsrf)
331			rain_fall = 0.
332			snow_fall = 0.
333			solsw = 165.
334			sollw = -53.
335			t_ancien = 273.15
336			q_ancien = 0.
337			agesno = 0.
338			!IM "slab" ocean
339	guez	13	tslab = tsolsrf(:, is_oce)
340	guez	3	seaice = 0.
341
342	guez	13	frugs(:, is_oce) = rugmer
343	guez	78	frugs(:, is_ter) = MAX(1e-5, zstd * zsig / 2)
344			frugs(:, is_lic) = MAX(1e-5, zstd * zsig / 2)
345	guez	3	frugs(:, is_sic) = 0.001
346			fder = 0.
347			clwcon = 0.
348			rnebcon = 0.
349			ratqs = 0.
350			run_off_lic_0 = 0.
351	guez	72	sig1 = 0.
352			w01 = 0.
353	guez	3
354	guez	13	call phyredem("startphy.nc", latfi, lonfi, pctsrf, &
355	guez	3	tsolsrf, tsoil, tslab, seaice, qsolsrf, qsol, snsrf, albe, alblw, &
356			evap, rain_fall, snow_fall, solsw, sollw, fder, radsol, frugs, &
357	guez	13	agesno, zmea, zstd, zsig, zgam, zthe, zpic, zval, &
358	guez	72	t_ancien, q_ancien, rnebcon, ratqs, clwcon, run_off_lic_0, sig1, w01)
359	guez	3	CALL histclo
360
361			END SUBROUTINE etat0
362
363			end module etat0_mod