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 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	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			USE exner_hyb_m, ONLY: exner_hyb
25	guez	107	use netcdf, only: nf90_nowrite
26			use netcdf95, only: nf95_close, nf95_inq_dimid, nf95_inquire_dimension, &
27			nf95_open
28	guez	39	use nr_util, only: pi
29	guez	103	USE paramet_m, ONLY: iip1, ip1jmp1, jjp1, llmp1
30	guez	83	USE q_sat_m, ONLY: q_sat
31	guez	88	use read_reanalyse_m, only: read_reanalyse
32	guez	83	USE serre, ONLY: clat, clon
33	guez	44	use tau2alpha_m, only: tau2alpha, dxdys
34	guez	108	use writefield_m, only: writefield
35	guez	3
36	guez	83	INTEGER, INTENT(IN):: itau
37	guez	102	REAL, intent(inout):: ucov(:, :, :) ! (iim + 1, jjm + 1, llm) vent covariant
38			REAL, intent(inout):: vcov(:, :, :) ! (iim + 1, jjm, llm) ! vent covariant
39
40	guez	108	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	guez	88	REAL, intent(in):: ps(:, :) ! (iim + 1, jjm + 1) pression au sol
45	guez	3
46	guez	83	! Local:
47
48	guez	108	! variables dynamiques pour les réanalyses
49	guez	102
50			REAL, save:: ucovrea1(iim + 1, jjm + 1, llm), vcovrea1(iim + 1, jjm, llm)
51			! vents covariants reanalyses
52
53	guez	90	REAL, save:: tetarea1(iim + 1, jjm + 1, llm) ! temp pot reales
54			REAL, save:: qrea1(iim + 1, jjm + 1, llm) ! temp pot reales
55	guez	102
56			REAL, save:: ucovrea2(iim + 1, jjm + 1, llm), vcovrea2(iim + 1, jjm, llm)
57			! vents covariants reanalyses
58
59	guez	90	REAL, save:: tetarea2(iim + 1, jjm + 1, llm) ! temp pot reales
60			REAL, save:: qrea2(iim + 1, jjm + 1, llm) ! temp pot reales
61	guez	44	REAL, save:: masserea2(ip1jmp1, llm) ! masse
62	guez	3
63	guez	102	! 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	guez	90	REAL, save:: alpha_q(iim + 1, jjm + 1)
67	guez	103	REAL, save:: alpha_t(iim + 1, jjm + 1)
68	guez	102	REAL, save:: alpha_u(iim + 1, jjm + 1), alpha_v(iim + 1, jjm)
69
70	guez	44	INTEGER, save:: step_rea, count_no_rea
71	guez	3
72	guez	36	INTEGER ilon, ilat
73	guez	102	REAL factt ! pas de temps entre deux appels au guidage, en fraction de jour
74			real ztau(iim + 1, jjm + 1)
75	guez	3
76	guez	102	INTEGER ij, l
77	guez	108	INTEGER ncid, dimid
78	guez	102	REAL tau
79	guez	44	INTEGER, SAVE:: nlev
80	guez	3
81	guez	44	! TEST SUR QSAT
82	guez	90	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	guez	3
86	guez	29	!-----------------------------------------------------------------------
87	guez	3
88	guez	108	!!PRINT *, 'Call sequence information: guide'
89	guez	3
90	guez	102	first_call: IF (itau == 0) THEN
91			CALL conf_guide
92	guez	3
93	guez	102	IF (online) THEN
94			! Constantes de temps de rappel en jour
95	guez	3
96	guez	102	! 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	guez	3
106	guez	102	factt = dtvr * iperiod / daysec
107	guez	3
108	guez	103	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	guez	102	ELSE
113	guez	107	! Cas où on force exactement par les variables analysées
114	guez	102	alpha_t = 0.
115			alpha_u = 0.
116			alpha_v = 0.
117	guez	103	alpha_q = 0.
118	guez	102	END IF
119	guez	3
120	guez	102	step_rea = 1
121			count_no_rea = 0
122	guez	108	ncid = -99
123	guez	3
124	guez	102	! lecture d'un fichier netcdf pour determiner le nombre de niveaux
125	guez	108	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	guez	3
130	guez	102	IF (ncep) THEN
131	guez	108	call nf95_inq_dimid(ncid, 'LEVEL', dimid)
132	guez	102	ELSE
133	guez	108	call nf95_inq_dimid(ncid, 'PRESSURE', dimid)
134	guez	102	END IF
135	guez	108	call nf95_inquire_dimension(ncid, dimid, nclen=nlev)
136	guez	102	PRINT *, 'nlev', nlev
137	guez	108	call nf95_close(ncid)
138	guez	102	! 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	guez	3
144	guez	102	! IMPORTATION DES VENTS, PRESSION ET TEMPERATURE REELS:
145	guez	3
146	guez	102	! 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	guez	3
153	guez	102	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	guez	108	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	guez	102	ELSE
173			count_no_rea = count_no_rea + 1
174			END IF
175	guez	3
176	guez	102	! Guidage
177	guez	3
178	guez	102	tau = mod(real(itau) / real(day_step / 4), 1.)
179	guez	3
180	guez	102	! x_gcm = a * x_gcm + (1 - a) * x_reanalyses
181	guez	3
182	guez	102	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	guez	3
192	guez	102	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	guez	3
202	guez	102	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	guez	3
208	guez	102	! 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	guez	3
218	guez	102	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	guez	3
228	guez	29	END SUBROUTINE guide
229	guez	20
230			END MODULE guide_m