libf/dyn3d/etat0.f90

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 inidissip_m, only: inidissip
    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
    CALL inidissip
    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	27	use inidissip_m, only: inidissip
39	guez	54	use inifilr_m, only: inifilr
40	guez	27	use inigeom_m, only: inigeom
41	guez	67	use massdair_m, only: massdair
42	guez	48	use netcdf, only: nf90_nowrite
43	guez	68	use netcdf95, only: nf95_close, nf95_get_var, nf95_gw_var, &
44			nf95_inq_varid, nf95_open
45	guez	39	use nr_util, only: pi
46	guez	27	use paramet_m, only: ip1jm, ip1jmp1
47			use phyredem_m, only: phyredem
48	guez	10	use pressure_var, only: pls, p3d
49	guez	27	use q_sat_m, only: q_sat
50	guez	7	use regr_lat_time_coefoz_m, only: regr_lat_time_coefoz
51			use regr_pr_o3_m, only: regr_pr_o3
52	guez	27	use serre, only: alphax
53	guez	68	use startdyn, only: start_init_dyn
54	guez	27	USE start_init_orog_m, only: start_init_orog, mask, phis
55	guez	43	use start_init_phys_m, only: start_init_phys
56	guez	54	use start_inter_3d_m, only: start_inter_3d
57	guez	28	use temps, only: itau_phy, annee_ref, day_ref
58	guez	3
59			! Variables local to the procedure:
60
61			REAL latfi(klon), lonfi(klon)
62			! (latitude and longitude of a point of the scalar grid identified
63			! by a simple index, in °)
64
65	guez	65	REAL, dimension(iim + 1, jjm + 1, llm):: ucov, t3d, tpot
66			REAL vcov(iim + 1, jjm, llm)
67	guez	3
68	guez	68	REAL q(iim + 1, jjm + 1, llm, nqmx)
69	guez	3	! (mass fractions of trace species
70	guez	68	! "q(i, j, l)" is at longitude "rlonv(i)", latitude "rlatu(j)"
71	guez	3	! and pressure level "pls(i, j, l)".)
72
73			real qsat(iim + 1, jjm + 1, llm) ! mass fraction of saturating water vapor
74			REAL tsol(klon), qsol(klon), sn(klon)
75			REAL tsolsrf(klon, nbsrf), qsolsrf(klon, nbsrf), snsrf(klon, nbsrf)
76			REAL albe(klon, nbsrf), evap(klon, nbsrf)
77			REAL alblw(klon, nbsrf)
78			REAL tsoil(klon, nsoilmx, nbsrf)
79			REAL radsol(klon), rain_fall(klon), snow_fall(klon)
80			REAL solsw(klon), sollw(klon), fder(klon)
81			!IM "slab" ocean
82			REAL tslab(klon)
83			real seaice(klon) ! kg m-2
84			REAL frugs(klon, nbsrf), agesno(klon, nbsrf)
85			REAL rugmer(klon)
86			real, dimension(iim + 1, jjm + 1):: relief, zstd_2d, zsig_2d, zgam_2d
87			real, dimension(iim + 1, jjm + 1):: zthe_2d, zpic_2d, zval_2d
88	guez	43	real, dimension(iim + 1, jjm + 1):: tsol_2d, qsol_2d, psol
89	guez	3	REAL zmea(klon), zstd(klon)
90			REAL zsig(klon), zgam(klon)
91			REAL zthe(klon)
92			REAL zpic(klon), zval(klon)
93			REAL t_ancien(klon, llm), q_ancien(klon, llm) !
94			REAL run_off_lic_0(klon)
95			real clwcon(klon, llm), rnebcon(klon, llm), ratqs(klon, llm)
96	guez	68
97			! Déclarations pour lecture glace de mer :
98			INTEGER iml_lic, jml_lic
99			INTEGER ncid, varid
100			REAL, pointer:: dlon_lic(:), dlat_lic(:)
101	guez	3	REAL, ALLOCATABLE:: fraclic(:, :) ! fraction land ice
102	guez	68	REAL flic_tmp(iim + 1, jjm + 1) ! fraction land ice temporary
103	guez	3
104			INTEGER l, ji
105
106			REAL pk(iim + 1, jjm + 1, llm) ! fonction d'Exner aux milieux des couches
107			real pks(iim + 1, jjm + 1)
108
109			REAL masse(iim + 1, jjm + 1, llm)
110	guez	70	REAL phi(iim + 1, jjm + 1, llm)
111	guez	3	REAL pbaru(ip1jmp1, llm), pbarv(ip1jm, llm)
112			REAL w(ip1jmp1, llm)
113			REAL phystep
114
115			!---------------------------------
116
117			print *, "Call sequence information: etat0"
118
119			dtvr = daysec / real(day_step)
120			print *, 'dtvr = ', dtvr
121
122	guez	36	! Construct a grid:
123
124	guez	3	pa = 5e4
125			CALL iniconst
126			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			call start_init_orog(relief, zstd_2d, zsig_2d, zgam_2d, zthe_2d, zpic_2d, &
140	guez	15	zval_2d) ! also compute "mask" and "phis"
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			CALL start_init_dyn(tsol_2d, psol)
147	guez	3
148			! Compute pressure on intermediate levels:
149	guez	15	forall(l = 1: llm + 1) p3d(:, :, l) = ap(l) + bp(l) * psol
150	guez	3	CALL exner_hyb(psol, 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	36	PRINT *, 'minval(t3d) = ', minval(t3d)
170			print *, "maxval(t3d) = ", maxval(t3d)
171	guez	3
172			tpot(:iim, :, :) = t3d(:iim, :, :) * cpp / pk(:iim, :, :)
173			tpot(iim + 1, :, :) = tpot(1, :, :)
174			DO l=1, llm
175			tpot(:, 1, l) = SUM(aire_2d(:, 1) * tpot(:, 1, l)) / apoln
176			tpot(:, jjm + 1, l) = SUM(aire_2d(:, jjm + 1) * tpot(:, jjm + 1, l)) &
177			/ 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	68	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			zmea = pack(relief, dyn_phy)
214			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	68	IF (MAXVAL( dlon_lic ) > pi) THEN
248			dlon_lic = dlon_lic * pi / 180.
249	guez	3	ENDIF
250	guez	68	IF (maxval( dlat_lic ) > pi) THEN
251			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			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			! 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			! Calcul intermédiaire:
292			CALL massdair(p3d, masse)
293
294			print *, 'ALPHAX = ', alphax
295
296			forall (l = 1:llm)
297			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	27	CALL inidissip
305	guez	3	itau_phy = 0
306			day_ref = dayref
307			annee_ref = anneeref
308
309	guez	70	CALL geopot(tpot, pk , pks, phis, phi)
310	guez	65	CALL caldyn0(ucov, vcov, tpot, psol, masse, pk, phis, phi, w, pbaru, &
311	guez	23	pbarv)
312	guez	5	CALL dynredem0("start.nc", dayref, phis)
313	guez	68	CALL dynredem1("start.nc", vcov, ucov, tpot, q, masse, psol, itau=0)
314	guez	3
315			! Ecriture état initial physique:
316			print *, "iphysiq = ", iphysiq
317			phystep = dtvr * REAL(iphysiq)
318			print *, 'phystep = ', phystep
319
320			! Initialisations :
321			tsolsrf(:, is_ter) = tsol
322			tsolsrf(:, is_lic) = tsol
323			tsolsrf(:, is_oce) = tsol
324			tsolsrf(:, is_sic) = tsol
325			snsrf(:, is_ter) = sn
326			snsrf(:, is_lic) = sn
327			snsrf(:, is_oce) = sn
328			snsrf(:, is_sic) = sn
329			albe(:, is_ter) = 0.08
330			albe(:, is_lic) = 0.6
331			albe(:, is_oce) = 0.5
332			albe(:, is_sic) = 0.6
333			alblw = albe
334	guez	15	evap = 0.
335	guez	3	qsolsrf(:, is_ter) = 150.
336			qsolsrf(:, is_lic) = 150.
337			qsolsrf(:, is_oce) = 150.
338			qsolsrf(:, is_sic) = 150.
339			tsoil = spread(spread(tsol, 2, nsoilmx), 3, nbsrf)
340			rain_fall = 0.
341			snow_fall = 0.
342			solsw = 165.
343			sollw = -53.
344			t_ancien = 273.15
345			q_ancien = 0.
346			agesno = 0.
347			!IM "slab" ocean
348	guez	13	tslab = tsolsrf(:, is_oce)
349	guez	3	seaice = 0.
350
351	guez	13	frugs(:, is_oce) = rugmer
352			frugs(:, is_ter) = MAX(1.e-05, zstd * zsig / 2)
353			frugs(:, is_lic) = MAX(1.e-05, zstd * zsig / 2)
354	guez	3	frugs(:, is_sic) = 0.001
355			fder = 0.
356			clwcon = 0.
357			rnebcon = 0.
358			ratqs = 0.
359			run_off_lic_0 = 0.
360
361	guez	13	call phyredem("startphy.nc", latfi, lonfi, pctsrf, &
362	guez	3	tsolsrf, tsoil, tslab, seaice, qsolsrf, qsol, snsrf, albe, alblw, &
363			evap, rain_fall, snow_fall, solsw, sollw, fder, radsol, frugs, &
364	guez	13	agesno, zmea, zstd, zsig, zgam, zthe, zpic, zval, &
365	guez	3	t_ancien, q_ancien, rnebcon, ratqs, clwcon, run_off_lic_0)
366			CALL histclo
367
368			END SUBROUTINE etat0
369
370			end module etat0_mod