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 exner_hyb_m, ONLY: exner_hyb
    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
    use writefield_m, only: writefield

    INTEGER, INTENT(IN):: itau
    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 réanalyses

    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 ncid, dimid
    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)

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

    !!PRINT *, 'Call sequence information: guide'

    first_call: IF (itau == 0) THEN
       CALL conf_guide

       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)
       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
       ncid = -99

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

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