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