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, tpot
    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, psol
    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

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

    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, psol)

    ! Compute pressure on intermediate levels:
    forall(l = 1: llm + 1) p3d(:, :, l) = ap(l) + bp(l) * psol
    CALL exner_hyb(psol, 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)

    tpot(:iim, :, :) = t3d(:iim, :, :) * cpp / pk(:iim, :, :)
    tpot(iim + 1, :, :) = tpot(1, :, :)
    DO l=1, llm
       tpot(:, 1, l) = SUM(aire_2d(:, 1) * tpot(:, 1, l)) / apoln
       tpot(:, jjm + 1, l) = SUM(aire_2d(:, jjm + 1) * tpot(:, 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(tpot, pk , pks,  phis, phi)
    CALL caldyn0(ucov, vcov, tpot, psol, masse, pk, phis, phi, w, pbaru, &
         pbarv)
    CALL dynredem0("start.nc", dayref, phis)
    CALL dynredem1("start.nc", vcov, ucov, tpot, q, masse, psol, 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.

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