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	MODULE guide_m
2
3	! 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
6	IMPLICIT NONE
7
8	REAL aire_min, aire_max
9
10	CONTAINS
11
12	SUBROUTINE guide(itau, ucov, vcov, teta, q, ps)
13
14	! Author: F.Hourdin
15
16	USE comconst, ONLY: cpp, daysec, dtvr, kappa
17	USE comgeom, ONLY: aire, rlatu, rlonv
18	USE conf_gcm_m, ONLY: day_step, iperiod
19	use conf_guide_m, only: conf_guide, guide_u, guide_v, guide_t, guide_q, &
20	ncep, ini_anal, tau_min_u, tau_max_u, tau_min_v, tau_max_v, &
21	tau_min_t, tau_max_t, tau_min_q, tau_max_q, online
22	USE dimens_m, ONLY: iim, jjm, llm
23	USE disvert_m, ONLY: ap, bp, preff, presnivs
24	use dump2d_m, only: dump2d
25	USE exner_hyb_m, ONLY: exner_hyb
26	USE inigrads_m, ONLY: inigrads
27	use netcdf, only: nf90_nowrite, nf90_close, nf90_inq_dimid
28	use netcdf95, only: nf95_inquire_dimension, nf95_open
29	use nr_util, only: pi
30	USE paramet_m, ONLY: iip1, ip1jmp1, jjp1, llmp1
31	USE q_sat_m, ONLY: q_sat
32	use read_reanalyse_m, only: read_reanalyse
33	USE serre, ONLY: clat, clon
34	use tau2alpha_m, only: tau2alpha, dxdys
35
36	INTEGER, INTENT(IN):: itau
37
38	! variables dynamiques
39
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	REAL, intent(inout):: teta(iim + 1, jjm + 1, llm) ! température potentielle
44	REAL, intent(inout):: q(iim + 1, jjm + 1, llm)
45	REAL, intent(in):: ps(:, :) ! (iim + 1, jjm + 1) pression au sol
46
47	! Local:
48
49	! variables dynamiques pour les reanalyses.
50
51	REAL, save:: ucovrea1(iim + 1, jjm + 1, llm), vcovrea1(iim + 1, jjm, llm)
52	! vents covariants reanalyses
53
54	REAL, save:: tetarea1(iim + 1, jjm + 1, llm) ! temp pot reales
55	REAL, save:: qrea1(iim + 1, jjm + 1, llm) ! temp pot reales
56
57	REAL, save:: ucovrea2(iim + 1, jjm + 1, llm), vcovrea2(iim + 1, jjm, llm)
58	! vents covariants reanalyses
59
60	REAL, save:: tetarea2(iim + 1, jjm + 1, llm) ! temp pot reales
61	REAL, save:: qrea2(iim + 1, jjm + 1, llm) ! temp pot reales
62	REAL, save:: masserea2(ip1jmp1, llm) ! masse
63
64	! 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	REAL, save:: alpha_q(iim + 1, jjm + 1)
68	REAL, save:: alpha_t(iim + 1, jjm + 1)
69	REAL, save:: alpha_u(iim + 1, jjm + 1), alpha_v(iim + 1, jjm)
70
71	INTEGER, save:: step_rea, count_no_rea
72
73	INTEGER ilon, ilat
74	REAL factt ! pas de temps entre deux appels au guidage, en fraction de jour
75	real ztau(iim + 1, jjm + 1)
76
77	INTEGER ij, l
78	INTEGER ncidpl, status
79	INTEGER rcod, rid
80	REAL tau
81	INTEGER, SAVE:: nlev
82
83	! TEST SUR QSAT
84	REAL p(iim + 1, jjm + 1, llmp1)
85	real pk(iim + 1, jjm + 1, llm), pks(iim + 1, jjm + 1)
86
87	REAL qsat(iim + 1, jjm + 1, llm)
88
89	INTEGER, parameter:: igrads = 2
90	REAL:: dtgrads = 100.
91
92	!-----------------------------------------------------------------------
93
94	PRINT *, 'Call sequence information: guide'
95
96	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
101	IF (online) THEN
102	! Constantes de temps de rappel en jour
103
104	! 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
114	factt = dtvr * iperiod / daysec
115
116	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
121	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	alpha_q = 0.
130	END IF
131
132	step_rea = 1
133	count_no_rea = 0
134	ncidpl = -99
135
136	! 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
142	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
156	! IMPORTATION DES VENTS, PRESSION ET TEMPERATURE REELS:
157
158	! 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
165	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
188	! Guidage
189
190	tau = mod(real(itau) / real(day_step / 4), 1.)
191
192	! x_gcm = a * x_gcm + (1 - a) * x_reanalyses
193
194	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
204	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
214	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
220	! 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
230	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
240	END SUBROUTINE guide
241
242	END MODULE guide_m