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
    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	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
28	use netcdf95, only: nf95_close, nf95_inq_dimid, nf95_inquire_dimension, &
29	nf95_open
30	use nr_util, only: pi
31	USE paramet_m, ONLY: iip1, ip1jmp1, jjp1, llmp1
32	USE q_sat_m, ONLY: q_sat
33	use read_reanalyse_m, only: read_reanalyse
34	USE serre, ONLY: clat, clon
35	use tau2alpha_m, only: tau2alpha, dxdys
36
37	INTEGER, INTENT(IN):: itau
38
39	! variables dynamiques
40
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	REAL, intent(inout):: teta(iim + 1, jjm + 1, llm) ! température potentielle
45	REAL, intent(inout):: q(iim + 1, jjm + 1, llm)
46	REAL, intent(in):: ps(:, :) ! (iim + 1, jjm + 1) pression au sol
47
48	! Local:
49
50	! variables dynamiques pour les reanalyses.
51
52	REAL, save:: ucovrea1(iim + 1, jjm + 1, llm), vcovrea1(iim + 1, jjm, llm)
53	! vents covariants reanalyses
54
55	REAL, save:: tetarea1(iim + 1, jjm + 1, llm) ! temp pot reales
56	REAL, save:: qrea1(iim + 1, jjm + 1, llm) ! temp pot reales
57
58	REAL, save:: ucovrea2(iim + 1, jjm + 1, llm), vcovrea2(iim + 1, jjm, llm)
59	! vents covariants reanalyses
60
61	REAL, save:: tetarea2(iim + 1, jjm + 1, llm) ! temp pot reales
62	REAL, save:: qrea2(iim + 1, jjm + 1, llm) ! temp pot reales
63	REAL, save:: masserea2(ip1jmp1, llm) ! masse
64
65	! 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	REAL, save:: alpha_q(iim + 1, jjm + 1)
69	REAL, save:: alpha_t(iim + 1, jjm + 1)
70	REAL, save:: alpha_u(iim + 1, jjm + 1), alpha_v(iim + 1, jjm)
71
72	INTEGER, save:: step_rea, count_no_rea
73
74	INTEGER ilon, ilat
75	REAL factt ! pas de temps entre deux appels au guidage, en fraction de jour
76	real ztau(iim + 1, jjm + 1)
77
78	INTEGER ij, l
79	INTEGER ncidpl
80	INTEGER rid
81	REAL tau
82	INTEGER, SAVE:: nlev
83
84	! TEST SUR QSAT
85	REAL p(iim + 1, jjm + 1, llmp1)
86	real pk(iim + 1, jjm + 1, llm), pks(iim + 1, jjm + 1)
87
88	REAL qsat(iim + 1, jjm + 1, llm)
89
90	INTEGER, parameter:: igrads = 2
91	REAL:: dtgrads = 100.
92
93	!-----------------------------------------------------------------------
94
95	PRINT *, 'Call sequence information: guide'
96
97	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
102	IF (online) THEN
103	! Constantes de temps de rappel en jour
104
105	! 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
115	factt = dtvr * iperiod / daysec
116
117	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
122	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	! Cas où on force exactement par les variables analysées
127	alpha_t = 0.
128	alpha_u = 0.
129	alpha_v = 0.
130	alpha_q = 0.
131	END IF
132
133	step_rea = 1
134	count_no_rea = 0
135	ncidpl = -99
136
137	! 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
143	IF (ncep) THEN
144	call nf95_inq_dimid(ncidpl, 'LEVEL', rid)
145	ELSE
146	call nf95_inq_dimid(ncidpl, 'PRESSURE', rid)
147	END IF
148	call nf95_inquire_dimension(ncidpl, rid, nclen=nlev)
149	PRINT *, 'nlev', nlev
150	call nf95_close(ncidpl)
151	! 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
157	! IMPORTATION DES VENTS, PRESSION ET TEMPERATURE REELS:
158
159	! 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
166	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
189	! Guidage
190
191	tau = mod(real(itau) / real(day_step / 4), 1.)
192
193	! x_gcm = a * x_gcm + (1 - a) * x_reanalyses
194
195	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
205	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
215	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
221	! 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
231	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
241	END SUBROUTINE guide
242
243	END MODULE guide_m