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