trunk/dyn3d/etat0.f90

module etat0_m

  IMPLICIT NONE

contains

  SUBROUTINE etat0(phis, pctsrf)

    ! From "etat0_netcdf.F", version 1.3, 2005/05/25 13:10:09

    use caldyn0_m, only: caldyn0
    use comconst, only: cpp, kappa, iniconst
    use comgeom, only: aire_2d, apoln, apols, cu_2d, cv_2d, inigeom
    use conf_gcm_m, only: nday
    use dimensions, only: iim, jjm, llm, nqmx
    use dimphy, only: klon
    use dimsoil, only: nsoilmx
    use disvert_m, only: ap, bp, preff, disvert
    use dynetat0_m, only: rlatu, rlatv, rlonu, rlonv, fyhyp, fxhyp
    use dynetat0_chosen_m, only: day_ref
    use dynredem0_m, only: dynredem0
    use dynredem1_m, only: dynredem1
    use exner_hyb_m, only: exner_hyb
    use geopot_m, only: geopot
    use grille_m_m, only: grille_m
    use grid_change, only: init_dyn_phy, dyn_phy
    use indicesol, only: is_oce, is_sic, is_ter, is_lic, epsfra, nbsrf
    use infotrac_init_m, only: infotrac_init
    use inifilr_m, only: inifilr
    use massdair_m, only: massdair
    use netcdf, only: nf90_nowrite
    use netcdf95, only: nf95_close, nf95_get_var, nf95_gw_var, nf95_put_var, &
         nf95_inq_varid, nf95_open
    use nr_util, only: pi, assert
    use phyetat0_m, only: masque, phyetat0_new
    use phyredem0_m, only: phyredem0, ncid_restartphy
    use phyredem_m, only: phyredem
    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 startdyn, only: start_init_dyn
    USE start_init_orog_m, only: start_init_orog
    use start_init_phys_m, only: start_init_phys
    use start_inter_3d_m, only: start_inter_3d
    use test_disvert_m, only: test_disvert

    REAL, intent(out):: phis(:, :) ! (iim + 1, jjm + 1)
    ! surface geopotential, in m2 s-2

    REAL, intent(out):: pctsrf(:, :) ! (klon, nbsrf)
    ! "pctsrf(i, :)" is the composition of the surface at horizontal
    ! position "i".

    ! Local:

    REAL, dimension(iim + 1, jjm + 1, llm):: ucov, t3d, teta
    REAL vcov(iim + 1, jjm, llm)

    REAL q(iim + 1, jjm + 1, llm, nqmx)
    ! (mass fractions of trace species
    ! "q(i, j, l)" is at longitude "rlonv(i)", latitude "rlatu(j)"
    ! and pressure level "pls(i, j, l)".)

    real qsat(iim + 1, jjm + 1, llm) ! mass fraction of saturating water vapor
    REAL qsolsrf(klon, nbsrf), fsnow(klon, nbsrf) 
    REAL falbe(klon, nbsrf)
    REAL ftsoil(klon, nsoilmx, nbsrf) 
    REAL null_array(klon)
    REAL solsw(klon), sollw(klon)
    !IM "slab" ocean
    REAL frugs(klon, nbsrf), agesno(klon, nbsrf)
    REAL rugmer(klon)
    real, dimension(iim + 1, jjm + 1):: zmea_2d, zstd_2d, zsig_2d, zgam_2d
    real, dimension(iim + 1, jjm + 1):: zthe_2d, zpic_2d, zval_2d

    real tsol_2d(iim + 1, jjm + 1)
    ! both soil temperature and surface temperature, in K

    real, dimension(iim + 1, jjm + 1):: qsol_2d, ps
    REAL zmea(klon), zstd(klon)
    REAL zsig(klon), zgam(klon)
    REAL zthe(klon)
    REAL zpic(klon), zval(klon)
    REAL t_ancien(klon, llm), q_ancien(klon, llm)
    real clwcon(klon, llm), rnebcon(klon, llm), ratqs(klon, llm)

    ! D\'eclarations pour lecture glace de mer :
    INTEGER iml_lic, jml_lic
    INTEGER ncid, varid
    REAL, ALLOCATABLE:: dlon_lic(:), dlat_lic(:)
    REAL, ALLOCATABLE:: landice(:, :) ! 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 sig1(klon, llm) ! section adiabatic updraft
    real w01(klon, llm) ! vertical velocity within adiabatic updraft

    real pls(iim + 1, jjm + 1, llm)
    ! (pressure at mid-layer of LMDZ grid, in Pa)
    ! "pls(i, j, l)" is at longitude "rlonv(i)", latitude "rlatu(j)",
    ! for layer "l")

    REAL p3d(iim + 1, jjm + 1, llm+1) ! pressure at layer interfaces, in Pa
    ! ("p3d(i, j, l)" is at longitude "rlonv(i)", latitude "rlatu(j)",
    ! for interface "l")

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

    print *, "Call sequence information: etat0"

    CALL iniconst

    ! Construct a grid:

    CALL disvert
    call test_disvert
    CALL fyhyp
    CALL fxhyp
    CALL inigeom
    CALL inifilr
    call start_init_orog(phis, zmea_2d, zstd_2d, zsig_2d, zgam_2d, zthe_2d, &
         zpic_2d, zval_2d) ! also compute "mask"
    call init_dyn_phy ! define the mask "dyn_phy" for distinct grid points
    call phyetat0_new

    call start_init_phys(tsol_2d, qsol_2d)
    CALL start_init_dyn(tsol_2d, phis, ps)

    ! Compute pressure on intermediate levels:
    forall(l = 1: llm + 1) p3d(:, :, l) = ap(l) + bp(l) * ps
    CALL exner_hyb(ps, p3d, pks, pk)
    call assert(MINVAL(pk) /= MAXVAL(pk), '"pk" should not be constant')

    pls = preff * (pk / cpp)**(1. / kappa)
    PRINT *, "minval(pls) = ", minval(pls)
    print *, "maxval(pls) = ", maxval(pls)

    call start_inter_3d('U', rlonv, rlatv, pls, ucov)
    forall (l = 1: llm) ucov(:iim, :, l) = ucov(:iim, :, l) * cu_2d(:iim, :)
    ucov(iim+1, :, :) = ucov(1, :, :)

    call start_inter_3d('V', rlonu, rlatu(:jjm), pls(:, :jjm, :), vcov)
    forall (l = 1: llm) vcov(:iim, :, l) = vcov(:iim, :, l) * cv_2d(:iim, :)
    vcov(iim + 1, :, :) = vcov(1, :, :)

    call start_inter_3d('TEMP', rlonu, rlatv, pls, t3d)
    PRINT *, 'minval(t3d) = ', minval(t3d)
    print *, "maxval(t3d) = ", maxval(t3d)

    teta(:iim, :, :) = t3d(:iim, :, :) * cpp / pk(:iim, :, :)
    teta(iim + 1, :, :) = teta(1, :, :)
    DO l = 1, llm
       teta(:, 1, l) = SUM(aire_2d(:, 1) * teta(:, 1, l)) / apoln
       teta(:, jjm + 1, l) = SUM(aire_2d(:, jjm + 1) * teta(:, jjm + 1, l)) &
            / apols
    ENDDO

    ! Calcul de l'humidit\'e \`a saturation :
    qsat = q_sat(t3d, pls)
    PRINT *, "minval(qsat) = ", minval(qsat)
    print *, "maxval(qsat) = ", maxval(qsat)
    IF (MINVAL(qsat) == MAXVAL(qsat)) stop '"qsat" should not be constant'

    ! Water vapor:
    call start_inter_3d('R', rlonu, rlatv, pls, q(:, :, :, 1))
    q(:, :, :, 1) = 0.01 * q(:, :, :, 1) * qsat
    WHERE (q(:, :, :, 1) < 0.) q(:, :, :, 1) = 1E-10
    DO l = 1, llm
       q(:, 1, l, 1) = SUM(aire_2d(:, 1) * q(:, 1, l, 1)) / apoln
       q(:, jjm + 1, l, 1) &
            = SUM(aire_2d(:, jjm + 1) * q(:, jjm + 1, l, 1)) / apols
    ENDDO

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

    if (nqmx >= 5) then
       ! Ozone:
       call regr_lat_time_coefoz
       call regr_pr_o3(p3d, q(:, :, :, 5))
       ! Convert from mole fraction to mass fraction:
       q(:, :, :, 5) = q(:, :, :, 5) * 48. / 29.
    end if

    null_array = 0.
    rugmer = 0.001
    zmea = pack(zmea_2d, dyn_phy)
    zstd = pack(zstd_2d, dyn_phy)
    zsig = pack(zsig_2d, dyn_phy)
    zgam = pack(zgam_2d, dyn_phy)
    zthe = pack(zthe_2d, dyn_phy)
    zpic = pack(zpic_2d, dyn_phy)
    zval = pack(zval_2d, dyn_phy)

    ! On initialise les sous-surfaces.
    ! Lecture du fichier glace de terre pour fixer la fraction de terre 
    ! et de glace de terre :

    call nf95_open("landiceref.nc", nf90_nowrite, ncid)

    call nf95_inq_varid(ncid, 'longitude', varid)
    call nf95_gw_var(ncid, varid, dlon_lic)
    iml_lic = size(dlon_lic)

    call nf95_inq_varid(ncid, 'latitude', varid)
    call nf95_gw_var(ncid, varid, dlat_lic)
    jml_lic = size(dlat_lic)

    call nf95_inq_varid(ncid, 'landice', varid)
    ALLOCATE(landice(iml_lic, jml_lic))
    call nf95_get_var(ncid, varid, landice)

    call nf95_close(ncid)

    ! Interpolation sur la grille T du mod\`ele :
    PRINT *, 'Dimensions de "landiceref.nc"'
    print *, "iml_lic = ", iml_lic
    print *, "jml_lic = ", jml_lic

    ! Si les coordonn\'ees sont en degr\'es, on les transforme :
    IF (MAXVAL(dlon_lic) > pi) THEN
       dlon_lic = dlon_lic * pi / 180.
    ENDIF
    IF (maxval(dlat_lic) > pi) THEN 
       dlat_lic = dlat_lic * pi/ 180.
    ENDIF

    flic_tmp(:iim, :) = grille_m(dlon_lic, dlat_lic, landice, 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\'equation avec le maque terre/mer :
    WHERE (pctsrf(:, is_lic) < EPSFRA) pctsrf(:, is_lic) = 0.
    WHERE (masque < EPSFRA) pctsrf(:, is_lic) = 0.
    where (masque <= EPSFRA) pctsrf(:, is_ter) = masque
    where (masque > EPSFRA)
       where (pctsrf(:, is_lic) >= masque)
          pctsrf(:, is_lic) = masque
          pctsrf(:, is_ter) = 0.
       elsewhere
          pctsrf(:, is_ter) = masque - pctsrf(:, is_lic)
          where (pctsrf(:, is_ter) < EPSFRA)
             pctsrf(:, is_ter) = 0.
             pctsrf(:, is_lic) = masque
          end where
       end where
    end where

    ! Sous-surface oc\'ean et glace de mer (pour d\'emarrer on met glace
    ! de mer \`a 0) :
    pctsrf(:, is_oce) = 1. - masque
    WHERE (pctsrf(:, is_oce) < EPSFRA) pctsrf(:, is_oce) = 0.

    ! V\'erification que la somme des sous-surfaces vaut 1 :
    ji = count(abs(sum(pctsrf, dim = 2) - 1.) > EPSFRA)
    IF (ji /= 0) then
       PRINT *, 'Bad surface percentages for ', ji, 'points'
    end IF

    ! Calcul interm\'ediaire :
    CALL massdair(p3d, masse)

    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

    call infotrac_init
    CALL geopot(teta, pk , pks, phis, phi)
    CALL caldyn0(ucov, vcov, teta, ps, pk, phis, phi)
    CALL dynredem0(day_ref, phis)
    CALL dynredem1(vcov, ucov, teta, q, masse, ps, itau = 0)

    ! Initialisations :
    fsnow = 0.
    falbe(:, is_ter) = 0.08
    falbe(:, is_lic) = 0.6
    falbe(:, is_oce) = 0.5
    falbe(:, is_sic) = 0.6
    qsolsrf = 150.
    ftsoil = spread(spread(pack(tsol_2d, dyn_phy), 2, nsoilmx), 3, nbsrf)
    solsw = 165.
    sollw = -53.
    t_ancien = 273.15
    q_ancien = 0.
    agesno = 0.

    frugs(:, is_oce) = rugmer
    frugs(:, is_ter) = MAX(1e-5, zstd * zsig / 2)
    frugs(:, is_lic) = MAX(1e-5, zstd * zsig / 2)
    frugs(:, is_sic) = 0.001
    clwcon = 0.
    rnebcon = 0.
    ratqs = 0.
    sig1 = 0.
    w01 = 0.

    call phyredem0(0)

    call nf95_inq_varid(ncid_restartphy, "trs", varid)
    call nf95_put_var(ncid_restartphy, varid, null_array)

    call phyredem(pctsrf, ftsoil(:, 1, :), ftsoil, qsolsrf, &
         pack(qsol_2d, dyn_phy), fsnow, falbe, null_array, null_array, solsw, &
         sollw, null_array, null_array, frugs, agesno, zmea, zstd, zsig, zgam, &
         zthe, zpic, zval, t_ancien, q_ancien, rnebcon, ratqs, clwcon, &
         null_array, sig1, w01)

  END SUBROUTINE etat0

end module etat0_m
1	guez	276	module etat0_m
2	guez	3
3			IMPLICIT NONE
4
5			contains
6
7	guez	276	SUBROUTINE etat0(phis, pctsrf)
8	guez	3
9	guez	90	! From "etat0_netcdf.F", version 1.3, 2005/05/25 13:10:09
10	guez	3
11	guez	27	use caldyn0_m, only: caldyn0
12	guez	79	use comconst, only: cpp, kappa, iniconst
13	guez	139	use comgeom, only: aire_2d, apoln, apols, cu_2d, cv_2d, inigeom
14	guez	207	use conf_gcm_m, only: nday
15	guez	265	use dimensions, only: iim, jjm, llm, nqmx
16	guez	276	use dimphy, only: klon
17	guez	3	use dimsoil, only: nsoilmx
18	guez	313	use disvert_m, only: ap, bp, preff, disvert
19			use dynetat0_m, only: rlatu, rlatv, rlonu, rlonv, fyhyp, fxhyp
20			use dynetat0_chosen_m, only: day_ref
21	guez	27	use dynredem0_m, only: dynredem0
22			use dynredem1_m, only: dynredem1
23			use exner_hyb_m, only: exner_hyb
24	guez	43	use geopot_m, only: geopot
25	guez	212	use grille_m_m, only: grille_m
26	guez	3	use grid_change, only: init_dyn_phy, dyn_phy
27	guez	276	use indicesol, only: is_oce, is_sic, is_ter, is_lic, epsfra, nbsrf
28	guez	321	use infotrac_init_m, only: infotrac_init
29	guez	54	use inifilr_m, only: inifilr
30	guez	67	use massdair_m, only: massdair
31	guez	48	use netcdf, only: nf90_nowrite
32	guez	157	use netcdf95, only: nf95_close, nf95_get_var, nf95_gw_var, nf95_put_var, &
33	guez	68	nf95_inq_varid, nf95_open
34	guez	90	use nr_util, only: pi, assert
35	guez	324	use phyetat0_m, only: masque, phyetat0_new
36	guez	157	use phyredem0_m, only: phyredem0, ncid_restartphy
37	guez	27	use phyredem_m, only: phyredem
38			use q_sat_m, only: q_sat
39	guez	7	use regr_lat_time_coefoz_m, only: regr_lat_time_coefoz
40			use regr_pr_o3_m, only: regr_pr_o3
41	guez	68	use startdyn, only: start_init_dyn
42	guez	278	USE start_init_orog_m, only: start_init_orog
43	guez	43	use start_init_phys_m, only: start_init_phys
44	guez	54	use start_inter_3d_m, only: start_inter_3d
45	guez	99	use test_disvert_m, only: test_disvert
46	guez	3
47	guez	107	REAL, intent(out):: phis(:, :) ! (iim + 1, jjm + 1)
48			! surface geopotential, in m2 s-2
49	guez	3
50	guez	276	REAL, intent(out):: pctsrf(:, :) ! (klon, nbsrf)
51			! "pctsrf(i, :)" is the composition of the surface at horizontal
52			! position "i".
53
54	guez	107	! Local:
55
56	guez	73	REAL, dimension(iim + 1, jjm + 1, llm):: ucov, t3d, teta
57	guez	65	REAL vcov(iim + 1, jjm, llm)
58	guez	3
59	guez	68	REAL q(iim + 1, jjm + 1, llm, nqmx)
60	guez	3	! (mass fractions of trace species
61	guez	68	! "q(i, j, l)" is at longitude "rlonv(i)", latitude "rlatu(j)"
62	guez	3	! and pressure level "pls(i, j, l)".)
63
64			real qsat(iim + 1, jjm + 1, llm) ! mass fraction of saturating water vapor
65	guez	330	REAL qsolsrf(klon, nbsrf), fsnow(klon, nbsrf)
66	guez	331	REAL falbe(klon, nbsrf)
67	guez	344	REAL ftsoil(klon, nsoilmx, nbsrf)
68	guez	157	REAL null_array(klon)
69			REAL solsw(klon), sollw(klon)
70	guez	3	!IM "slab" ocean
71			REAL frugs(klon, nbsrf), agesno(klon, nbsrf)
72			REAL rugmer(klon)
73	guez	78	real, dimension(iim + 1, jjm + 1):: zmea_2d, zstd_2d, zsig_2d, zgam_2d
74	guez	3	real, dimension(iim + 1, jjm + 1):: zthe_2d, zpic_2d, zval_2d
75	guez	344
76			real tsol_2d(iim + 1, jjm + 1)
77			! both soil temperature and surface temperature, in K
78
79			real, dimension(iim + 1, jjm + 1):: qsol_2d, ps
80	guez	3	REAL zmea(klon), zstd(klon)
81			REAL zsig(klon), zgam(klon)
82			REAL zthe(klon)
83			REAL zpic(klon), zval(klon)
84	guez	78	REAL t_ancien(klon, llm), q_ancien(klon, llm)
85	guez	3	real clwcon(klon, llm), rnebcon(klon, llm), ratqs(klon, llm)
86	guez	68
87	guez	90	! D\'eclarations pour lecture glace de mer :
88	guez	68	INTEGER iml_lic, jml_lic
89			INTEGER ncid, varid
90	guez	225	REAL, ALLOCATABLE:: dlon_lic(:), dlat_lic(:)
91	guez	324	REAL, ALLOCATABLE:: landice(:, :) ! fraction land ice
92	guez	68	REAL flic_tmp(iim + 1, jjm + 1) ! fraction land ice temporary
93	guez	3
94			INTEGER l, ji
95
96			REAL pk(iim + 1, jjm + 1, llm) ! fonction d'Exner aux milieux des couches
97			real pks(iim + 1, jjm + 1)
98			REAL masse(iim + 1, jjm + 1, llm)
99	guez	70	REAL phi(iim + 1, jjm + 1, llm)
100	guez	72	real sig1(klon, llm) ! section adiabatic updraft
101			real w01(klon, llm) ! vertical velocity within adiabatic updraft
102
103	guez	97	real pls(iim + 1, jjm + 1, llm)
104			! (pressure at mid-layer of LMDZ grid, in Pa)
105			! "pls(i, j, l)" is at longitude "rlonv(i)", latitude "rlatu(j)",
106			! for layer "l")
107
108			REAL p3d(iim + 1, jjm + 1, llm+1) ! pressure at layer interfaces, in Pa
109			! ("p3d(i, j, l)" is at longitude "rlonv(i)", latitude "rlatu(j)",
110			! for interface "l")
111
112	guez	3	!---------------------------------
113
114			print *, "Call sequence information: etat0"
115
116	guez	79	CALL iniconst
117	guez	3
118	guez	36	! Construct a grid:
119
120	guez	79	CALL disvert
121	guez	99	call test_disvert
122	guez	277	CALL fyhyp
123			CALL fxhyp
124	guez	3	CALL inigeom
125			CALL inifilr
126	guez	78	call start_init_orog(phis, zmea_2d, zstd_2d, zsig_2d, zgam_2d, zthe_2d, &
127			zpic_2d, zval_2d) ! also compute "mask"
128	guez	3	call init_dyn_phy ! define the mask "dyn_phy" for distinct grid points
129	guez	278	call phyetat0_new
130	guez	3
131	guez	43	call start_init_phys(tsol_2d, qsol_2d)
132	guez	78	CALL start_init_dyn(tsol_2d, phis, ps)
133	guez	3
134			! Compute pressure on intermediate levels:
135	guez	73	forall(l = 1: llm + 1) p3d(:, :, l) = ap(l) + bp(l) * ps
136			CALL exner_hyb(ps, p3d, pks, pk)
137	guez	90	call assert(MINVAL(pk) /= MAXVAL(pk), '"pk" should not be constant')
138	guez	3
139	guez	36	pls = preff * (pk / cpp)**(1. / kappa)
140			PRINT *, "minval(pls) = ", minval(pls)
141			print *, "maxval(pls) = ", maxval(pls)
142	guez	3
143	guez	65	call start_inter_3d('U', rlonv, rlatv, pls, ucov)
144			forall (l = 1: llm) ucov(:iim, :, l) = ucov(:iim, :, l) * cu_2d(:iim, :)
145			ucov(iim+1, :, :) = ucov(1, :, :)
146	guez	3
147	guez	65	call start_inter_3d('V', rlonu, rlatu(:jjm), pls(:, :jjm, :), vcov)
148			forall (l = 1: llm) vcov(:iim, :, l) = vcov(:iim, :, l) * cv_2d(:iim, :)
149			vcov(iim + 1, :, :) = vcov(1, :, :)
150	guez	3
151	guez	23	call start_inter_3d('TEMP', rlonu, rlatv, pls, t3d)
152	guez	78	PRINT *, 'minval(t3d) = ', minval(t3d)
153	guez	36	print *, "maxval(t3d) = ", maxval(t3d)
154	guez	3
155	guez	73	teta(:iim, :, :) = t3d(:iim, :, :) * cpp / pk(:iim, :, :)
156			teta(iim + 1, :, :) = teta(1, :, :)
157			DO l = 1, llm
158			teta(:, 1, l) = SUM(aire_2d(:, 1) * teta(:, 1, l)) / apoln
159			teta(:, jjm + 1, l) = SUM(aire_2d(:, jjm + 1) * teta(:, jjm + 1, l)) &
160	guez	3	/ apols
161			ENDDO
162
163	guez	90	! Calcul de l'humidit\'e \`a saturation :
164	guez	36	qsat = q_sat(t3d, pls)
165			PRINT *, "minval(qsat) = ", minval(qsat)
166			print *, "maxval(qsat) = ", maxval(qsat)
167	guez	3	IF (MINVAL(qsat) == MAXVAL(qsat)) stop '"qsat" should not be constant'
168
169			! Water vapor:
170	guez	68	call start_inter_3d('R', rlonu, rlatv, pls, q(:, :, :, 1))
171			q(:, :, :, 1) = 0.01 * q(:, :, :, 1) * qsat
172			WHERE (q(:, :, :, 1) < 0.) q(:, :, :, 1) = 1E-10
173	guez	3	DO l = 1, llm
174	guez	68	q(:, 1, l, 1) = SUM(aire_2d(:, 1) * q(:, 1, l, 1)) / apoln
175			q(:, jjm + 1, l, 1) &
176			= SUM(aire_2d(:, jjm + 1) * q(:, jjm + 1, l, 1)) / apols
177	guez	3	ENDDO
178
179	guez	68	q(:, :, :, 2:4) = 0. ! liquid water, radon and lead
180	guez	3
181	guez	5	if (nqmx >= 5) then
182			! Ozone:
183	guez	7	call regr_lat_time_coefoz
184	guez	97	call regr_pr_o3(p3d, q(:, :, :, 5))
185	guez	7	! Convert from mole fraction to mass fraction:
186	guez	78	q(:, :, :, 5) = q(:, :, :, 5) * 48. / 29.
187	guez	5	end if
188	guez	3
189	guez	157	null_array = 0.
190	guez	13	rugmer = 0.001
191	guez	78	zmea = pack(zmea_2d, dyn_phy)
192	guez	13	zstd = pack(zstd_2d, dyn_phy)
193			zsig = pack(zsig_2d, dyn_phy)
194			zgam = pack(zgam_2d, dyn_phy)
195			zthe = pack(zthe_2d, dyn_phy)
196			zpic = pack(zpic_2d, dyn_phy)
197			zval = pack(zval_2d, dyn_phy)
198	guez	3
199	guez	49	! On initialise les sous-surfaces.
200	guez	3	! Lecture du fichier glace de terre pour fixer la fraction de terre
201	guez	49	! et de glace de terre :
202	guez	68
203	guez	48	call nf95_open("landiceref.nc", nf90_nowrite, ncid)
204	guez	68
205			call nf95_inq_varid(ncid, 'longitude', varid)
206			call nf95_gw_var(ncid, varid, dlon_lic)
207			iml_lic = size(dlon_lic)
208
209			call nf95_inq_varid(ncid, 'latitude', varid)
210			call nf95_gw_var(ncid, varid, dlat_lic)
211			jml_lic = size(dlat_lic)
212
213	guez	48	call nf95_inq_varid(ncid, 'landice', varid)
214	guez	324	ALLOCATE(landice(iml_lic, jml_lic))
215			call nf95_get_var(ncid, varid, landice)
216	guez	68
217	guez	48	call nf95_close(ncid)
218	guez	3
219	guez	90	! Interpolation sur la grille T du mod\`ele :
220	guez	68	PRINT *, 'Dimensions de "landiceref.nc"'
221	guez	3	print *, "iml_lic = ", iml_lic
222			print *, "jml_lic = ", jml_lic
223
224	guez	90	! Si les coordonn\'ees sont en degr\'es, on les transforme :
225	guez	78	IF (MAXVAL(dlon_lic) > pi) THEN
226	guez	68	dlon_lic = dlon_lic * pi / 180.
227	guez	3	ENDIF
228	guez	73	IF (maxval(dlat_lic) > pi) THEN
229	guez	68	dlat_lic = dlat_lic * pi/ 180.
230	guez	3	ENDIF
231
232	guez	324	flic_tmp(:iim, :) = grille_m(dlon_lic, dlat_lic, landice, rlonv(:iim), &
233	guez	3	rlatu)
234			flic_tmp(iim + 1, :) = flic_tmp(1, :)
235
236	guez	279	! Passage sur la grille physique :
237	guez	15	pctsrf = 0.
238	guez	3	pctsrf(:, is_lic) = pack(flic_tmp, dyn_phy)
239	guez	279
240			! Ad\'equation avec le maque terre/mer :
241	guez	73	WHERE (pctsrf(:, is_lic) < EPSFRA) pctsrf(:, is_lic) = 0.
242	guez	324	WHERE (masque < EPSFRA) pctsrf(:, is_lic) = 0.
243			where (masque <= EPSFRA) pctsrf(:, is_ter) = masque
244			where (masque > EPSFRA)
245			where (pctsrf(:, is_lic) >= masque)
246			pctsrf(:, is_lic) = masque
247	guez	3	pctsrf(:, is_ter) = 0.
248			elsewhere
249	guez	324	pctsrf(:, is_ter) = masque - pctsrf(:, is_lic)
250	guez	3	where (pctsrf(:, is_ter) < EPSFRA)
251			pctsrf(:, is_ter) = 0.
252	guez	324	pctsrf(:, is_lic) = masque
253	guez	3	end where
254			end where
255			end where
256
257	guez	90	! Sous-surface oc\'ean et glace de mer (pour d\'emarrer on met glace
258			! de mer \`a 0) :
259	guez	324	pctsrf(:, is_oce) = 1. - masque
260	guez	3	WHERE (pctsrf(:, is_oce) < EPSFRA) pctsrf(:, is_oce) = 0.
261
262	guez	163	! V\'erification que la somme des sous-surfaces vaut 1 :
263	guez	15	ji = count(abs(sum(pctsrf, dim = 2) - 1.) > EPSFRA)
264			IF (ji /= 0) then
265	guez	91	PRINT *, 'Bad surface percentages for ', ji, 'points'
266	guez	15	end IF
267	guez	3
268	guez	90	! Calcul interm\'ediaire :
269	guez	3	CALL massdair(p3d, masse)
270
271	guez	78	forall (l = 1:llm)
272	guez	3	masse(:, 1, l) = SUM(aire_2d(:iim, 1) * masse(:iim, 1, l)) / apoln
273			masse(:, jjm + 1, l) = &
274			SUM(aire_2d(:iim, jjm + 1) * masse(:iim, jjm + 1, l)) / apols
275			END forall
276
277	guez	321	call infotrac_init
278	guez	78	CALL geopot(teta, pk , pks, phis, phi)
279	guez	163	CALL caldyn0(ucov, vcov, teta, ps, pk, phis, phi)
280	guez	157	CALL dynredem0(day_ref, phis)
281			CALL dynredem1(vcov, ucov, teta, q, masse, ps, itau = 0)
282	guez	3
283			! Initialisations :
284	guez	330	fsnow = 0.
285	guez	331	falbe(:, is_ter) = 0.08
286			falbe(:, is_lic) = 0.6
287			falbe(:, is_oce) = 0.5
288			falbe(:, is_sic) = 0.6
289	guez	99	qsolsrf = 150.
290	guez	344	ftsoil = spread(spread(pack(tsol_2d, dyn_phy), 2, nsoilmx), 3, nbsrf)
291	guez	3	solsw = 165.
292			sollw = -53.
293			t_ancien = 273.15
294			q_ancien = 0.
295			agesno = 0.
296
297	guez	13	frugs(:, is_oce) = rugmer
298	guez	78	frugs(:, is_ter) = MAX(1e-5, zstd * zsig / 2)
299			frugs(:, is_lic) = MAX(1e-5, zstd * zsig / 2)
300	guez	3	frugs(:, is_sic) = 0.001
301			clwcon = 0.
302			rnebcon = 0.
303			ratqs = 0.
304	guez	72	sig1 = 0.
305			w01 = 0.
306	guez	3
307	guez	350	call phyredem0(0)
308	guez	157
309			call nf95_inq_varid(ncid_restartphy, "trs", varid)
310			call nf95_put_var(ncid_restartphy, varid, null_array)
311
312	guez	344	call phyredem(pctsrf, ftsoil(:, 1, :), ftsoil, qsolsrf, &
313	guez	331	pack(qsol_2d, dyn_phy), fsnow, falbe, null_array, null_array, solsw, &
314	guez	304	sollw, null_array, null_array, frugs, agesno, zmea, zstd, zsig, zgam, &
315			zthe, zpic, zval, t_ancien, q_ancien, rnebcon, ratqs, clwcon, &
316			null_array, sig1, w01)
317	guez	3
318			END SUBROUTINE etat0
319
320	guez	276	end module etat0_m