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 comvert, only: ap, bp, preff, pa
    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 dynredem0_m, only: dynredem0
    use dynredem1_m, only: dynredem1
    use exner_hyb_m, only: exner_hyb
    USE flincom, only: flinclo, flinopen_nozoom, flininfo
    use geopot_m, only: geopot
    use grid_atob, only: grille_m
    use grid_change, only: init_dyn_phy, dyn_phy
    use histcom, 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 inigeom_m, only: inigeom
    use netcdf, only: nf90_nowrite
    use netcdf95, only: nf95_open, nf95_close, nf95_get_var, nf95_inq_varid
    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 start_init_orog_m, only: start_init_orog, mask, phis
    use start_init_phys_m, only: start_init_phys
    use startdyn, only: start_inter_3d, start_init_dyn
    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):: uvent, t3d, tpot
    REAL vvent(iim + 1, jjm, llm)

    REAL q3d(iim + 1, jjm + 1, llm, nqmx)
    ! (mass fractions of trace species
    ! "q3d(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, llm_tmp, ttm_tmp
    INTEGER itaul(1), fid, ncid, varid
    REAL lev(1), date
    REAL, ALLOCATABLE:: lon_lic(:, :), lat_lic(:, :)
    REAL, ALLOCATABLE:: 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(ip1jmp1, llm)
    REAL pbaru(ip1jmp1, llm), pbarv(ip1jm, llm)
    REAL w(ip1jmp1, llm)
    REAL phystep
    real trash

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

    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, uvent)
    forall (l = 1: llm) uvent(:iim, :, l) = uvent(:iim, :, l) * cu_2d(:iim, :)
    uvent(iim+1, :, :) = uvent(1, :, :)

    call start_inter_3d('V', rlonu, rlatu(:jjm), pls(:, :jjm, :), vvent)
    forall (l = 1: llm) vvent(:iim, :, l) = vvent(:iim, :, l) * cv_2d(:iim, :)
    vvent(iim + 1, :, :) = vvent(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, q3d(:, :, :, 1))
    q3d(:, :, :, 1) = 0.01 * q3d(:, :, :, 1) * qsat
    WHERE (q3d(:, :, :, 1) < 0.) q3d(:, :, :, 1) = 1E-10
    DO l = 1, llm
       q3d(:, 1, l, 1) = SUM(aire_2d(:, 1) * q3d(:, 1, l, 1)) / apoln
       q3d(:, jjm + 1, l, 1) &
            = SUM(aire_2d(:, jjm + 1) * q3d(:, jjm + 1, l, 1)) / apols
    ENDDO

    q3d(:, :, :, 2:4) = 0. ! liquid water, radon and lead

    if (nqmx >= 5) then
       ! Ozone:
       call regr_lat_time_coefoz
       call regr_pr_o3(q3d(:, :, :, 5))
       ! Convert from mole fraction to mass fraction:
       q3d(:, :, :, 5) = q3d(:, :, :, 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 flininfo("landiceref.nc", iml_lic, jml_lic, llm_tmp, &
         ttm_tmp, fid)
    ALLOCATE(lat_lic(iml_lic, jml_lic))
    ALLOCATE(lon_lic(iml_lic, jml_lic))
    ALLOCATE(dlon_lic(iml_lic))
    ALLOCATE(dlat_lic(jml_lic))
    ALLOCATE(fraclic(iml_lic, jml_lic))
    CALL flinopen_nozoom(iml_lic, jml_lic, &
         llm_tmp, lon_lic, lat_lic, lev, ttm_tmp, itaul, date, trash,  &
         fid)
    CALL flinclo(fid)
    call nf95_open("landiceref.nc", nf90_nowrite, ncid)
    call nf95_inq_varid(ncid, 'landice', varid)
    call nf95_get_var(ncid, varid, fraclic)
    call nf95_close(ncid)

    ! Interpolation sur la grille T du modèle :
    PRINT *, 'Dimensions de "landice"'
    print *, "iml_lic = ", iml_lic
    print *, "jml_lic = ", jml_lic

    ! Si les coordonnées sont en degrés, on les transforme :
    IF (MAXVAL( lon_lic ) > pi)  THEN
       lon_lic = lon_lic * pi / 180.
    ENDIF
    IF (maxval( lat_lic ) > pi) THEN 
       lat_lic = lat_lic * pi/ 180.
    ENDIF

    dlon_lic = lon_lic(:, 1)
    dlat_lic = lat_lic(1, :) 

    flic_tmp(:iim, :) = grille_m(dlon_lic, dlat_lic, fraclic, rlonv(:iim), &
         rlatu)
    flic_tmp(iim + 1, :) = flic_tmp(1, :)

    ! 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(ip1jmp1, tpot, pk , pks,  phis, phi)
    CALL caldyn0(uvent, vvent, tpot, psol, masse, pk, phis, phi, w, pbaru, &
         pbarv)
    CALL dynredem0("start.nc", dayref, phis)
    CALL dynredem1("start.nc", vvent, uvent, tpot, q3d, 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	guez	3	use comvert, only: ap, bp, preff, pa
25	guez	27	use conf_gcm_m, only: day_step, iphysiq, dayref, anneeref
26			use dimens_m, only: iim, jjm, llm, nqmx
27	guez	3	use dimphy, only: zmasq
28			use dimsoil, only: nsoilmx
29	guez	27	use dynredem0_m, only: dynredem0
30			use dynredem1_m, only: dynredem1
31			use exner_hyb_m, only: exner_hyb
32	guez	39	USE flincom, only: flinclo, flinopen_nozoom, flininfo
33	guez	43	use geopot_m, only: geopot
34	guez	3	use grid_atob, only: grille_m
35			use grid_change, only: init_dyn_phy, dyn_phy
36	guez	39	use histcom, only: histclo
37	guez	27	use indicesol, only: is_oce, is_sic, is_ter, is_lic, epsfra
38	guez	18	use iniadvtrac_m, only: iniadvtrac
39	guez	27	use inidissip_m, only: inidissip
40			use inigeom_m, only: inigeom
41	guez	48	use netcdf, only: nf90_nowrite
42			use netcdf95, only: nf95_open, nf95_close, nf95_get_var, nf95_inq_varid
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			USE start_init_orog_m, only: start_init_orog, mask, phis
52	guez	43	use start_init_phys_m, only: start_init_phys
53	guez	27	use startdyn, only: start_inter_3d, start_init_dyn
54	guez	28	use temps, only: itau_phy, annee_ref, day_ref
55	guez	3
56			! Variables local to the procedure:
57
58			REAL latfi(klon), lonfi(klon)
59			! (latitude and longitude of a point of the scalar grid identified
60			! by a simple index, in °)
61
62			REAL, dimension(iim + 1, jjm + 1, llm):: uvent, t3d, tpot
63			REAL vvent(iim + 1, jjm, llm)
64
65			REAL q3d(iim + 1, jjm + 1, llm, nqmx)
66			! (mass fractions of trace species
67			! "q3d(i, j, l)" is at longitude "rlonv(i)", latitude "rlatu(j)"
68			! and pressure level "pls(i, j, l)".)
69
70			real qsat(iim + 1, jjm + 1, llm) ! mass fraction of saturating water vapor
71			REAL tsol(klon), qsol(klon), sn(klon)
72			REAL tsolsrf(klon, nbsrf), qsolsrf(klon, nbsrf), snsrf(klon, nbsrf)
73			REAL albe(klon, nbsrf), evap(klon, nbsrf)
74			REAL alblw(klon, nbsrf)
75			REAL tsoil(klon, nsoilmx, nbsrf)
76			REAL radsol(klon), rain_fall(klon), snow_fall(klon)
77			REAL solsw(klon), sollw(klon), fder(klon)
78			!IM "slab" ocean
79			REAL tslab(klon)
80			real seaice(klon) ! kg m-2
81			REAL frugs(klon, nbsrf), agesno(klon, nbsrf)
82			REAL rugmer(klon)
83			real, dimension(iim + 1, jjm + 1):: relief, zstd_2d, zsig_2d, zgam_2d
84			real, dimension(iim + 1, jjm + 1):: zthe_2d, zpic_2d, zval_2d
85	guez	43	real, dimension(iim + 1, jjm + 1):: tsol_2d, qsol_2d, psol
86	guez	3	REAL zmea(klon), zstd(klon)
87			REAL zsig(klon), zgam(klon)
88			REAL zthe(klon)
89			REAL zpic(klon), zval(klon)
90			REAL t_ancien(klon, llm), q_ancien(klon, llm) !
91			REAL run_off_lic_0(klon)
92			real clwcon(klon, llm), rnebcon(klon, llm), ratqs(klon, llm)
93			! déclarations pour lecture glace de mer
94			INTEGER iml_lic, jml_lic, llm_tmp, ttm_tmp
95	guez	48	INTEGER itaul(1), fid, ncid, varid
96	guez	3	REAL lev(1), date
97			REAL, ALLOCATABLE:: lon_lic(:, :), lat_lic(:, :)
98			REAL, ALLOCATABLE:: dlon_lic(:), dlat_lic(:)
99			REAL, ALLOCATABLE:: fraclic(:, :) ! fraction land ice
100			REAL flic_tmp(iim + 1, jjm + 1) !fraction land ice temporary
101
102			INTEGER l, ji
103
104			REAL pk(iim + 1, jjm + 1, llm) ! fonction d'Exner aux milieux des couches
105			real pks(iim + 1, jjm + 1)
106
107			REAL masse(iim + 1, jjm + 1, llm)
108			REAL phi(ip1jmp1, llm)
109			REAL pbaru(ip1jmp1, llm), pbarv(ip1jm, llm)
110			REAL w(ip1jmp1, llm)
111			REAL phystep
112	guez	28	real trash
113	guez	3
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	23	call start_inter_3d('U', rlonv, rlatv, pls, uvent)
160	guez	3	forall (l = 1: llm) uvent(:iim, :, l) = uvent(:iim, :, l) * cu_2d(:iim, :)
161			uvent(iim+1, :, :) = uvent(1, :, :)
162
163	guez	23	call start_inter_3d('V', rlonu, rlatu(:jjm), pls(:, :jjm, :), vvent)
164	guez	3	forall (l = 1: llm) vvent(:iim, :, l) = vvent(:iim, :, l) * cv_2d(:iim, :)
165			vvent(iim + 1, :, :) = vvent(1, :, :)
166
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	23	call start_inter_3d('R', rlonu, rlatv, pls, q3d(:, :, :, 1))
187			q3d(:, :, :, 1) = 0.01 * q3d(:, :, :, 1) * qsat
188	guez	3	WHERE (q3d(:, :, :, 1) < 0.) q3d(:, :, :, 1) = 1E-10
189			DO l = 1, llm
190			q3d(:, 1, l, 1) = SUM(aire_2d(:, 1) * q3d(:, 1, l, 1)) / apoln
191			q3d(:, jjm + 1, l, 1) &
192			= SUM(aire_2d(:, jjm + 1) * q3d(:, jjm + 1, l, 1)) / apols
193			ENDDO
194
195			q3d(:, :, :, 2:4) = 0. ! liquid water, radon and lead
196
197	guez	5	if (nqmx >= 5) then
198			! Ozone:
199	guez	7	call regr_lat_time_coefoz
200			call regr_pr_o3(q3d(:, :, :, 5))
201			! Convert from mole fraction to mass fraction:
202			q3d(:, :, :, 5) = q3d(:, :, :, 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	3	CALL flininfo("landiceref.nc", iml_lic, jml_lic, llm_tmp, &
224			ttm_tmp, fid)
225			ALLOCATE(lat_lic(iml_lic, jml_lic))
226			ALLOCATE(lon_lic(iml_lic, jml_lic))
227			ALLOCATE(dlon_lic(iml_lic))
228			ALLOCATE(dlat_lic(jml_lic))
229			ALLOCATE(fraclic(iml_lic, jml_lic))
230	guez	32	CALL flinopen_nozoom(iml_lic, jml_lic, &
231	guez	28	llm_tmp, lon_lic, lat_lic, lev, ttm_tmp, itaul, date, trash, &
232	guez	3	fid)
233	guez	49	CALL flinclo(fid)
234	guez	48	call nf95_open("landiceref.nc", nf90_nowrite, ncid)
235			call nf95_inq_varid(ncid, 'landice', varid)
236			call nf95_get_var(ncid, varid, fraclic)
237			call nf95_close(ncid)
238	guez	3
239			! Interpolation sur la grille T du modèle :
240			PRINT *, 'Dimensions de "landice"'
241			print *, "iml_lic = ", iml_lic
242			print *, "jml_lic = ", jml_lic
243
244			! Si les coordonnées sont en degrés, on les transforme :
245	guez	15	IF (MAXVAL( lon_lic ) > pi) THEN
246			lon_lic = lon_lic * pi / 180.
247	guez	3	ENDIF
248	guez	15	IF (maxval( lat_lic ) > pi) THEN
249			lat_lic = lat_lic * pi/ 180.
250	guez	3	ENDIF
251
252	guez	13	dlon_lic = lon_lic(:, 1)
253			dlat_lic = lat_lic(1, :)
254	guez	3
255			flic_tmp(:iim, :) = grille_m(dlon_lic, dlat_lic, fraclic, rlonv(:iim), &
256			rlatu)
257			flic_tmp(iim + 1, :) = flic_tmp(1, :)
258
259			! 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	27	CALL inidissip
304	guez	3	itau_phy = 0
305			day_ref = dayref
306			annee_ref = anneeref
307
308	guez	22	CALL geopot(ip1jmp1, tpot, pk , pks, phis, phi)
309	guez	23	CALL caldyn0(uvent, vvent, tpot, psol, masse, pk, phis, phi, w, pbaru, &
310			pbarv)
311	guez	5	CALL dynredem0("start.nc", dayref, phis)
312	guez	27	CALL dynredem1("start.nc", vvent, uvent, tpot, q3d, masse, psol, itau=0)
313	guez	3
314			! Ecriture état initial physique:
315			print *, "iphysiq = ", iphysiq
316			phystep = dtvr * REAL(iphysiq)
317			print *, 'phystep = ', phystep
318
319			! Initialisations :
320			tsolsrf(:, is_ter) = tsol
321			tsolsrf(:, is_lic) = tsol
322			tsolsrf(:, is_oce) = tsol
323			tsolsrf(:, is_sic) = tsol
324			snsrf(:, is_ter) = sn
325			snsrf(:, is_lic) = sn
326			snsrf(:, is_oce) = sn
327			snsrf(:, is_sic) = sn
328			albe(:, is_ter) = 0.08
329			albe(:, is_lic) = 0.6
330			albe(:, is_oce) = 0.5
331			albe(:, is_sic) = 0.6
332			alblw = albe
333	guez	15	evap = 0.
334	guez	3	qsolsrf(:, is_ter) = 150.
335			qsolsrf(:, is_lic) = 150.
336			qsolsrf(:, is_oce) = 150.
337			qsolsrf(:, is_sic) = 150.
338			tsoil = spread(spread(tsol, 2, nsoilmx), 3, nbsrf)
339			rain_fall = 0.
340			snow_fall = 0.
341			solsw = 165.
342			sollw = -53.
343			t_ancien = 273.15
344			q_ancien = 0.
345			agesno = 0.
346			!IM "slab" ocean
347	guez	13	tslab = tsolsrf(:, is_oce)
348	guez	3	seaice = 0.
349
350	guez	13	frugs(:, is_oce) = rugmer
351			frugs(:, is_ter) = MAX(1.e-05, zstd * zsig / 2)
352			frugs(:, is_lic) = MAX(1.e-05, zstd * zsig / 2)
353	guez	3	frugs(:, is_sic) = 0.001
354			fder = 0.
355			clwcon = 0.
356			rnebcon = 0.
357			ratqs = 0.
358			run_off_lic_0 = 0.
359
360	guez	13	call phyredem("startphy.nc", latfi, lonfi, pctsrf, &
361	guez	3	tsolsrf, tsoil, tslab, seaice, qsolsrf, qsol, snsrf, albe, alblw, &
362			evap, rain_fall, snow_fall, solsw, sollw, fder, radsol, frugs, &
363	guez	13	agesno, zmea, zstd, zsig, zgam, zthe, zpic, zval, &
364	guez	3	t_ancien, q_ancien, rnebcon, ratqs, clwcon, run_off_lic_0)
365			CALL histclo
366
367			END SUBROUTINE etat0
368
369			end module etat0_mod