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