dyn3d/Guide/guide.f

MODULE guide_m

  ! From dyn3d/guide.F, version 1.3 2005/05/25 13:10:09
  ! and dyn3d/guide.h, version 1.1.1.1 2004/05/19 12:53:06

  IMPLICIT NONE

  REAL aire_min, aire_max

CONTAINS

  SUBROUTINE guide(itau, ucov, vcov, teta, q, ps)

    ! Author: F.Hourdin

    USE comconst, ONLY: cpp, daysec, dtvr, kappa
    USE comgeom, ONLY: aire, rlatu, rlonv
    USE conf_gcm_m, ONLY: day_step, iperiod
    use conf_guide_m, only: conf_guide, guide_u, guide_v, guide_t, guide_q, &
         ncep, ini_anal, tau_min_u, tau_max_u, tau_min_v, tau_max_v, &
         tau_min_t, tau_max_t, tau_min_q, tau_max_q, tau_min_p, tau_max_p, &
         online
    USE dimens_m, ONLY: iim, jjm, llm
    USE disvert_m, ONLY: ap, bp, preff, presnivs
    use dump2d_m, only: dump2d
    USE exner_hyb_m, ONLY: exner_hyb
    USE inigrads_m, ONLY: inigrads
    use massdair_m, only: massdair
    use netcdf, only: nf90_nowrite, nf90_close, nf90_inq_dimid
    use netcdf95, only: nf95_inquire_dimension, nf95_open
    use nr_util, only: pi
    USE paramet_m, ONLY: iip1, ip1jm, ip1jmp1, jjp1, llmp1
    USE q_sat_m, ONLY: q_sat
    use read_reanalyse_m, only: read_reanalyse
    USE serre, ONLY: clat, clon
    use tau2alpha_m, only: tau2alpha, dxdys

    INTEGER, INTENT(IN):: itau

    ! variables dynamiques

    REAL, intent(inout):: ucov(:, :, :) ! (iim + 1, jjm + 1, llm) vent covariant
    REAL, intent(inout):: vcov(:, :, :) ! (iim + 1, jjm, llm) ! vent covariant

    REAL, intent(inout):: teta(iim + 1, jjm + 1, llm) ! température potentielle 
    REAL, intent(inout):: q(iim + 1, jjm + 1, llm)
    REAL, intent(in):: ps(:, :) ! (iim + 1, jjm + 1) pression au sol

    ! Local:

    ! variables dynamiques pour les reanalyses.

    REAL, save:: ucovrea1(iim + 1, jjm + 1, llm), vcovrea1(iim + 1, jjm, llm)
    ! vents covariants reanalyses

    REAL, save:: tetarea1(iim + 1, jjm + 1, llm) ! temp pot reales
    REAL, save:: qrea1(iim + 1, jjm + 1, llm) ! temp pot reales

    REAL, save:: ucovrea2(iim + 1, jjm + 1, llm), vcovrea2(iim + 1, jjm, llm)
    ! vents covariants reanalyses

    REAL, save:: tetarea2(iim + 1, jjm + 1, llm) ! temp pot reales
    REAL, save:: qrea2(iim + 1, jjm + 1, llm) ! temp pot reales
    REAL, save:: masserea2(ip1jmp1, llm) ! masse

    ! alpha determine la part des injections de donnees a chaque etape
    ! alpha=1 signifie pas d'injection
    ! alpha=0 signifie injection totale
    REAL, save:: alpha_q(iim + 1, jjm + 1)
    REAL, save:: alpha_t(iim + 1, jjm + 1), alpha_p(ip1jmp1)
    REAL, save:: alpha_u(iim + 1, jjm + 1), alpha_v(iim + 1, jjm)

    INTEGER, save:: step_rea, count_no_rea

    INTEGER ilon, ilat
    REAL factt ! pas de temps entre deux appels au guidage, en fraction de jour
    real ztau(iim + 1, jjm + 1)

    INTEGER ij, l
    INTEGER ncidpl, status
    INTEGER rcod, rid
    REAL tau
    INTEGER, SAVE:: nlev

    ! TEST SUR QSAT
    REAL p(iim + 1, jjm + 1, llmp1)
    real pk(iim + 1, jjm + 1, llm), pks(iim + 1, jjm + 1)

    REAL qsat(iim + 1, jjm + 1, llm)

    INTEGER, parameter:: igrads = 2
    REAL:: dtgrads = 100.

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

    PRINT *, 'Call sequence information: guide'

    first_call: IF (itau == 0) THEN
       CALL conf_guide
       CALL inigrads(igrads, rlonv, 180. / pi, -180., 180., rlatu, -90., &
            90., 180. / pi, presnivs, 1., dtgrads, 'guide', 'dyn_zon ')

       IF (online) THEN
          ! Constantes de temps de rappel en jour

          ! coordonnees du centre du zoom
          CALL coordij(clon, clat, ilon, ilat)
          ! aire de la maille au centre du zoom
          aire_min = aire(ilon+(ilat - 1) * iip1)
          ! aire maximale de la maille
          aire_max = 0.
          DO ij = 1, ip1jmp1
             aire_max = max(aire_max, aire(ij))
          END DO

          factt = dtvr * iperiod / daysec

          CALL tau2alpha(3, iip1, jjm, factt, tau_min_v, tau_max_v, alpha_v)
          CALL tau2alpha(2, iip1, jjp1, factt, tau_min_u, tau_max_u, alpha_u)
          CALL tau2alpha(1, iip1, jjp1, factt, tau_min_t, tau_max_t, alpha_t)
          CALL tau2alpha(1, iip1, jjp1, factt, tau_min_p, tau_max_p, alpha_p)
          CALL tau2alpha(1, iip1, jjp1, factt, tau_min_q, tau_max_q, alpha_q)

          CALL dump2d(iip1, jjp1, aire, 'AIRE MAILLe ')
          CALL dump2d(iip1, jjp1, alpha_u, 'COEFF U ')
          CALL dump2d(iip1, jjp1, alpha_t, 'COEFF T ')
       ELSE
          ! Cas ou on force exactement par les variables analysees
          alpha_t = 0.
          alpha_u = 0.
          alpha_v = 0.
          alpha_p = 0.
       END IF

       step_rea = 1
       count_no_rea = 0
       ncidpl = -99

       ! lecture d'un fichier netcdf pour determiner le nombre de niveaux
       if (guide_u) call nf95_open('u.nc',Nf90_NOWRITe,ncidpl)
       if (guide_v) call nf95_open('v.nc',nf90_nowrite,ncidpl)
       if (guide_T) call nf95_open('T.nc',nf90_nowrite,ncidpl)
       if (guide_Q) call nf95_open('hur.nc',nf90_nowrite, ncidpl)

       IF (ncep) THEN
          status = nf90_inq_dimid(ncidpl, 'LEVEL', rid)
       ELSE
          status = nf90_inq_dimid(ncidpl, 'PRESSURE', rid)
       END IF
       call nf95_inquire_dimension(ncidpl, rid, nclen=nlev)
       PRINT *, 'nlev', nlev
       rcod = nf90_close(ncidpl)
       ! Lecture du premier etat des reanalyses.
       CALL read_reanalyse(1, ps, ucovrea2, vcovrea2, tetarea2, qrea2, &
            masserea2, nlev)
       qrea2 = max(qrea2, 0.1)
    END IF first_call

    ! IMPORTATION DES VENTS, PRESSION ET TEMPERATURE REELS:

    ! Nudging fields are given 4 times per day:
    IF (mod(itau, day_step / 4) == 0) THEN
       vcovrea1 = vcovrea2
       ucovrea1 = ucovrea2
       tetarea1 = tetarea2
       qrea1 = qrea2

       PRINT *, 'LECTURE REANALYSES, pas ', step_rea, 'apres ', &
            count_no_rea, ' non lectures'
       step_rea = step_rea + 1
       CALL read_reanalyse(step_rea, ps, ucovrea2, vcovrea2, tetarea2, qrea2, &
            masserea2, nlev)
       qrea2 = max(qrea2, 0.1)
       factt = dtvr * iperiod / daysec
       ztau = factt / max(alpha_t, 1E-10)
       CALL wrgrads(igrads, 1, aire, 'aire ', 'aire ')
       CALL wrgrads(igrads, 1, dxdys, 'dxdy ', 'dxdy ')
       CALL wrgrads(igrads, 1, alpha_u, 'au ', 'au ')
       CALL wrgrads(igrads, 1, alpha_t, 'at ', 'at ')
       CALL wrgrads(igrads, 1, ztau, 'taut ', 'taut ')
       CALL wrgrads(igrads, llm, ucov, 'u ', 'u ')
       CALL wrgrads(igrads, llm, ucovrea2, 'ua ', 'ua ')
       CALL wrgrads(igrads, llm, teta, 'T ', 'T ')
       CALL wrgrads(igrads, llm, tetarea2, 'Ta ', 'Ta ')
       CALL wrgrads(igrads, llm, qrea2, 'Qa ', 'Qa ')
       CALL wrgrads(igrads, llm, q, 'Q ', 'Q ')
    ELSE
       count_no_rea = count_no_rea + 1
    END IF

    ! Guidage

    tau = mod(real(itau) / real(day_step / 4), 1.)

    ! x_gcm = a * x_gcm + (1 - a) * x_reanalyses

    IF (guide_u) THEN
       IF (itau == 0 .AND. ini_anal) then
          ucov = ucovrea1
       else
          forall (l = 1: llm) ucov(:, :, l) = (1. - alpha_u) * ucov(:, :, l) &
               + alpha_u * ((1. - tau) * ucovrea1(:, :, l) &
               + tau * ucovrea2(:, :, l))
       end IF
    END IF

    IF (guide_t) THEN
       IF (itau == 0 .AND. ini_anal) then
          teta = tetarea1
       else
          forall (l = 1: llm) teta(:, :, l) = (1. - alpha_t) * teta(:, :, l) &
               + alpha_t * ((1. - tau) * tetarea1(:, :, l) &
               + tau * tetarea2(:, :, l))
       end IF
    END IF

    IF (guide_q) THEN
       ! Calcul de l'humidité saturante :
       forall (l = 1: llm + 1) p(:, :, l) = ap(l) + bp(l) * ps
       CALL exner_hyb(ps, p, pks, pk)
       qsat = q_sat(pk * teta / cpp, preff * (pk / cpp)**(1. / kappa))

       ! humidité relative en % -> humidité spécifique
       IF (itau == 0 .AND. ini_anal) then
          q = qsat * qrea1 * 0.01
       else
          forall (l = 1: llm) q(:, :, l) = (1. - alpha_q) * q(:, :, l) &
               + alpha_q * (qsat(:, :, l) * ((1. - tau) * qrea1(:, :, l) &
               + tau * qrea2(:, :, l)) * 0.01)
       end IF
    END IF

    IF (guide_v) THEN
       IF (itau == 0 .AND. ini_anal) then
          vcov = vcovrea1
       else
          forall (l = 1: llm) vcov(:, :, l) = (1. - alpha_v) * vcov(:, :, l) &
               + alpha_v * ((1. - tau) * vcovrea1(:, :, l) &
               + tau * vcovrea2(:, :, l))
       end IF
    END IF

  END SUBROUTINE guide

END MODULE guide_m
1	guez	20	MODULE guide_m
2	guez	3
3	guez	29	! From dyn3d/guide.F, version 1.3 2005/05/25 13:10:09
4			! and dyn3d/guide.h, version 1.1.1.1 2004/05/19 12:53:06
5	guez	3
6	guez	37	IMPLICIT NONE
7
8	guez	36	REAL aire_min, aire_max
9	guez	3
10	guez	20	CONTAINS
11	guez	3
12	guez	102	SUBROUTINE guide(itau, ucov, vcov, teta, q, ps)
13	guez	3
14	guez	29	! Author: F.Hourdin
15	guez	3
16	guez	83	USE comconst, ONLY: cpp, daysec, dtvr, kappa
17			USE comgeom, ONLY: aire, rlatu, rlonv
18			USE conf_gcm_m, ONLY: day_step, iperiod
19	guez	44	use conf_guide_m, only: conf_guide, guide_u, guide_v, guide_t, guide_q, &
20	guez	85	ncep, ini_anal, tau_min_u, tau_max_u, tau_min_v, tau_max_v, &
21	guez	44	tau_min_t, tau_max_t, tau_min_q, tau_max_q, tau_min_p, tau_max_p, &
22			online
23	guez	88	USE dimens_m, ONLY: iim, jjm, llm
24	guez	83	USE disvert_m, ONLY: ap, bp, preff, presnivs
25	guez	102	use dump2d_m, only: dump2d
26	guez	83	USE exner_hyb_m, ONLY: exner_hyb
27			USE inigrads_m, ONLY: inigrads
28	guez	67	use massdair_m, only: massdair
29	guez	102	use netcdf, only: nf90_nowrite, nf90_close, nf90_inq_dimid
30			use netcdf95, only: nf95_inquire_dimension, nf95_open
31	guez	39	use nr_util, only: pi
32	guez	83	USE paramet_m, ONLY: iip1, ip1jm, ip1jmp1, jjp1, llmp1
33			USE q_sat_m, ONLY: q_sat
34	guez	88	use read_reanalyse_m, only: read_reanalyse
35	guez	83	USE serre, ONLY: clat, clon
36	guez	44	use tau2alpha_m, only: tau2alpha, dxdys
37	guez	3
38	guez	83	INTEGER, INTENT(IN):: itau
39
40	guez	44	! variables dynamiques
41	guez	102
42			REAL, intent(inout):: ucov(:, :, :) ! (iim + 1, jjm + 1, llm) vent covariant
43			REAL, intent(inout):: vcov(:, :, :) ! (iim + 1, jjm, llm) ! vent covariant
44
45	guez	90	REAL, intent(inout):: teta(iim + 1, jjm + 1, llm) ! température potentielle
46	guez	102	REAL, intent(inout):: q(iim + 1, jjm + 1, llm)
47	guez	88	REAL, intent(in):: ps(:, :) ! (iim + 1, jjm + 1) pression au sol
48	guez	3
49	guez	83	! Local:
50
51	guez	44	! variables dynamiques pour les reanalyses.
52	guez	102
53			REAL, save:: ucovrea1(iim + 1, jjm + 1, llm), vcovrea1(iim + 1, jjm, llm)
54			! vents covariants reanalyses
55
56	guez	90	REAL, save:: tetarea1(iim + 1, jjm + 1, llm) ! temp pot reales
57			REAL, save:: qrea1(iim + 1, jjm + 1, llm) ! temp pot reales
58	guez	102
59			REAL, save:: ucovrea2(iim + 1, jjm + 1, llm), vcovrea2(iim + 1, jjm, llm)
60			! vents covariants reanalyses
61
62	guez	90	REAL, save:: tetarea2(iim + 1, jjm + 1, llm) ! temp pot reales
63			REAL, save:: qrea2(iim + 1, jjm + 1, llm) ! temp pot reales
64	guez	44	REAL, save:: masserea2(ip1jmp1, llm) ! masse
65	guez	3
66	guez	102	! alpha determine la part des injections de donnees a chaque etape
67			! alpha=1 signifie pas d'injection
68			! alpha=0 signifie injection totale
69	guez	90	REAL, save:: alpha_q(iim + 1, jjm + 1)
70			REAL, save:: alpha_t(iim + 1, jjm + 1), alpha_p(ip1jmp1)
71	guez	102	REAL, save:: alpha_u(iim + 1, jjm + 1), alpha_v(iim + 1, jjm)
72
73	guez	44	INTEGER, save:: step_rea, count_no_rea
74	guez	3
75	guez	36	INTEGER ilon, ilat
76	guez	102	REAL factt ! pas de temps entre deux appels au guidage, en fraction de jour
77			real ztau(iim + 1, jjm + 1)
78	guez	3
79	guez	102	INTEGER ij, l
80	guez	90	INTEGER ncidpl, status
81	guez	36	INTEGER rcod, rid
82	guez	102	REAL tau
83	guez	44	INTEGER, SAVE:: nlev
84	guez	3
85	guez	44	! TEST SUR QSAT
86	guez	90	REAL p(iim + 1, jjm + 1, llmp1)
87			real pk(iim + 1, jjm + 1, llm), pks(iim + 1, jjm + 1)
88	guez	3
89	guez	90	REAL qsat(iim + 1, jjm + 1, llm)
90	guez	3
91	guez	102	INTEGER, parameter:: igrads = 2
92	guez	44	REAL:: dtgrads = 100.
93	guez	3
94	guez	29	!-----------------------------------------------------------------------
95	guez	3
96	guez	29	PRINT *, 'Call sequence information: guide'
97	guez	3
98	guez	102	first_call: IF (itau == 0) THEN
99			CALL conf_guide
100			CALL inigrads(igrads, rlonv, 180. / pi, -180., 180., rlatu, -90., &
101			90., 180. / pi, presnivs, 1., dtgrads, 'guide', 'dyn_zon ')
102	guez	3
103	guez	102	IF (online) THEN
104			! Constantes de temps de rappel en jour
105	guez	3
106	guez	102	! coordonnees du centre du zoom
107			CALL coordij(clon, clat, ilon, ilat)
108			! aire de la maille au centre du zoom
109			aire_min = aire(ilon+(ilat - 1) * iip1)
110			! aire maximale de la maille
111			aire_max = 0.
112			DO ij = 1, ip1jmp1
113			aire_max = max(aire_max, aire(ij))
114			END DO
115	guez	3
116	guez	102	factt = dtvr * iperiod / daysec
117	guez	3
118	guez	102	CALL tau2alpha(3, iip1, jjm, factt, tau_min_v, tau_max_v, alpha_v)
119			CALL tau2alpha(2, iip1, jjp1, factt, tau_min_u, tau_max_u, alpha_u)
120			CALL tau2alpha(1, iip1, jjp1, factt, tau_min_t, tau_max_t, alpha_t)
121			CALL tau2alpha(1, iip1, jjp1, factt, tau_min_p, tau_max_p, alpha_p)
122			CALL tau2alpha(1, iip1, jjp1, factt, tau_min_q, tau_max_q, alpha_q)
123	guez	3
124	guez	102	CALL dump2d(iip1, jjp1, aire, 'AIRE MAILLe ')
125			CALL dump2d(iip1, jjp1, alpha_u, 'COEFF U ')
126			CALL dump2d(iip1, jjp1, alpha_t, 'COEFF T ')
127			ELSE
128			! Cas ou on force exactement par les variables analysees
129			alpha_t = 0.
130			alpha_u = 0.
131			alpha_v = 0.
132			alpha_p = 0.
133			END IF
134	guez	3
135	guez	102	step_rea = 1
136			count_no_rea = 0
137			ncidpl = -99
138	guez	3
139	guez	102	! lecture d'un fichier netcdf pour determiner le nombre de niveaux
140			if (guide_u) call nf95_open('u.nc',Nf90_NOWRITe,ncidpl)
141			if (guide_v) call nf95_open('v.nc',nf90_nowrite,ncidpl)
142			if (guide_T) call nf95_open('T.nc',nf90_nowrite,ncidpl)
143			if (guide_Q) call nf95_open('hur.nc',nf90_nowrite, ncidpl)
144	guez	3
145	guez	102	IF (ncep) THEN
146			status = nf90_inq_dimid(ncidpl, 'LEVEL', rid)
147			ELSE
148			status = nf90_inq_dimid(ncidpl, 'PRESSURE', rid)
149			END IF
150			call nf95_inquire_dimension(ncidpl, rid, nclen=nlev)
151			PRINT *, 'nlev', nlev
152			rcod = nf90_close(ncidpl)
153			! Lecture du premier etat des reanalyses.
154			CALL read_reanalyse(1, ps, ucovrea2, vcovrea2, tetarea2, qrea2, &
155			masserea2, nlev)
156			qrea2 = max(qrea2, 0.1)
157			END IF first_call
158	guez	3
159	guez	102	! IMPORTATION DES VENTS, PRESSION ET TEMPERATURE REELS:
160	guez	3
161	guez	102	! Nudging fields are given 4 times per day:
162			IF (mod(itau, day_step / 4) == 0) THEN
163			vcovrea1 = vcovrea2
164			ucovrea1 = ucovrea2
165			tetarea1 = tetarea2
166			qrea1 = qrea2
167	guez	3
168	guez	102	PRINT *, 'LECTURE REANALYSES, pas ', step_rea, 'apres ', &
169			count_no_rea, ' non lectures'
170			step_rea = step_rea + 1
171			CALL read_reanalyse(step_rea, ps, ucovrea2, vcovrea2, tetarea2, qrea2, &
172			masserea2, nlev)
173			qrea2 = max(qrea2, 0.1)
174			factt = dtvr * iperiod / daysec
175			ztau = factt / max(alpha_t, 1E-10)
176			CALL wrgrads(igrads, 1, aire, 'aire ', 'aire ')
177			CALL wrgrads(igrads, 1, dxdys, 'dxdy ', 'dxdy ')
178			CALL wrgrads(igrads, 1, alpha_u, 'au ', 'au ')
179			CALL wrgrads(igrads, 1, alpha_t, 'at ', 'at ')
180			CALL wrgrads(igrads, 1, ztau, 'taut ', 'taut ')
181			CALL wrgrads(igrads, llm, ucov, 'u ', 'u ')
182			CALL wrgrads(igrads, llm, ucovrea2, 'ua ', 'ua ')
183			CALL wrgrads(igrads, llm, teta, 'T ', 'T ')
184			CALL wrgrads(igrads, llm, tetarea2, 'Ta ', 'Ta ')
185			CALL wrgrads(igrads, llm, qrea2, 'Qa ', 'Qa ')
186			CALL wrgrads(igrads, llm, q, 'Q ', 'Q ')
187			ELSE
188			count_no_rea = count_no_rea + 1
189			END IF
190	guez	3
191	guez	102	! Guidage
192	guez	3
193	guez	102	tau = mod(real(itau) / real(day_step / 4), 1.)
194	guez	3
195	guez	102	! x_gcm = a * x_gcm + (1 - a) * x_reanalyses
196	guez	3
197	guez	102	IF (guide_u) THEN
198			IF (itau == 0 .AND. ini_anal) then
199			ucov = ucovrea1
200			else
201			forall (l = 1: llm) ucov(:, :, l) = (1. - alpha_u) * ucov(:, :, l) &
202			+ alpha_u * ((1. - tau) * ucovrea1(:, :, l) &
203			+ tau * ucovrea2(:, :, l))
204			end IF
205			END IF
206	guez	3
207	guez	102	IF (guide_t) THEN
208			IF (itau == 0 .AND. ini_anal) then
209			teta = tetarea1
210			else
211			forall (l = 1: llm) teta(:, :, l) = (1. - alpha_t) * teta(:, :, l) &
212			+ alpha_t * ((1. - tau) * tetarea1(:, :, l) &
213			+ tau * tetarea2(:, :, l))
214			end IF
215			END IF
216	guez	3
217	guez	102	IF (guide_q) THEN
218			! Calcul de l'humidité saturante :
219			forall (l = 1: llm + 1) p(:, :, l) = ap(l) + bp(l) * ps
220			CALL exner_hyb(ps, p, pks, pk)
221			qsat = q_sat(pk * teta / cpp, preff * (pk / cpp)**(1. / kappa))
222	guez	3
223	guez	102	! humidité relative en % -> humidité spécifique
224			IF (itau == 0 .AND. ini_anal) then
225			q = qsat * qrea1 * 0.01
226			else
227			forall (l = 1: llm) q(:, :, l) = (1. - alpha_q) * q(:, :, l) &
228			+ alpha_q * (qsat(:, :, l) * ((1. - tau) * qrea1(:, :, l) &
229			+ tau * qrea2(:, :, l)) * 0.01)
230			end IF
231			END IF
232	guez	3
233	guez	102	IF (guide_v) THEN
234			IF (itau == 0 .AND. ini_anal) then
235			vcov = vcovrea1
236			else
237			forall (l = 1: llm) vcov(:, :, l) = (1. - alpha_v) * vcov(:, :, l) &
238			+ alpha_v * ((1. - tau) * vcovrea1(:, :, l) &
239			+ tau * vcovrea2(:, :, l))
240			end IF
241			END IF
242	guez	3
243	guez	29	END SUBROUTINE guide
244	guez	20
245			END MODULE guide_m