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
    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=0 signifie pas d'injection
    ! alpha=1 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

    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

          if (guide_u) CALL tau2alpha(3, factt, tau_min_v, tau_max_v, alpha_v)
          if (guide_v) CALL tau2alpha(2, factt, tau_min_u, tau_max_u, alpha_u)
          if (guide_t) CALL tau2alpha(1, factt, tau_min_t, tau_max_t, alpha_t)
          if (guide_q) CALL tau2alpha(1, factt, tau_min_q, tau_max_q, alpha_q)
       ELSE
          ! Cas où on force exactement par les variables analysées
          if (guide_u) alpha_t = 1.
          if (guide_v) alpha_u = 1.
          if (guide_t) alpha_v = 1.
          if (guide_q) alpha_q = 1.
       END IF

       step_rea = 1
       count_no_rea = 0

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

       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)

       if (guide_u) CALL writefield("alpha_u", alpha_u)
       if (guide_t) CALL writefield("alpha_t", alpha_t)
    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

       if (guide_u) then
          CALL writefield("ucov", ucov)
          CALL writefield("ucovrea2", ucovrea2)
       end if

       if (guide_t) then
          CALL writefield("teta", teta)
          CALL writefield("tetarea2", tetarea2)
       end if

       if (guide_q) then
          CALL writefield("qrea2", qrea2)
          CALL writefield("q", q)
       end if
    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	109	use tau2alpha_m, only: tau2alpha
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	guez	112	! alpha=0 signifie pas d'injection
65			! alpha=1 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	guez	3
75	guez	102	INTEGER ij, l
76	guez	108	INTEGER ncid, dimid
77	guez	102	REAL tau
78	guez	44	INTEGER, SAVE:: nlev
79	guez	3
80	guez	44	! TEST SUR QSAT
81	guez	90	REAL p(iim + 1, jjm + 1, llmp1)
82			real pk(iim + 1, jjm + 1, llm), pks(iim + 1, jjm + 1)
83			REAL qsat(iim + 1, jjm + 1, llm)
84	guez	3
85	guez	29	!-----------------------------------------------------------------------
86	guez	3
87	guez	108	!!PRINT *, 'Call sequence information: guide'
88	guez	3
89	guez	102	first_call: IF (itau == 0) THEN
90			CALL conf_guide
91	guez	3
92	guez	102	IF (online) THEN
93			! Constantes de temps de rappel en jour
94	guez	3
95	guez	102	! coordonnees du centre du zoom
96			CALL coordij(clon, clat, ilon, ilat)
97			! aire de la maille au centre du zoom
98			aire_min = aire(ilon+(ilat - 1) * iip1)
99			! aire maximale de la maille
100			aire_max = 0.
101			DO ij = 1, ip1jmp1
102			aire_max = max(aire_max, aire(ij))
103			END DO
104	guez	3
105	guez	102	factt = dtvr * iperiod / daysec
106	guez	3
107	guez	112	if (guide_u) CALL tau2alpha(3, factt, tau_min_v, tau_max_v, alpha_v)
108			if (guide_v) CALL tau2alpha(2, factt, tau_min_u, tau_max_u, alpha_u)
109			if (guide_t) CALL tau2alpha(1, factt, tau_min_t, tau_max_t, alpha_t)
110			if (guide_q) CALL tau2alpha(1, factt, tau_min_q, tau_max_q, alpha_q)
111	guez	102	ELSE
112	guez	107	! Cas où on force exactement par les variables analysées
113	guez	112	if (guide_u) alpha_t = 1.
114			if (guide_v) alpha_u = 1.
115			if (guide_t) alpha_v = 1.
116			if (guide_q) alpha_q = 1.
117	guez	102	END IF
118	guez	3
119	guez	102	step_rea = 1
120			count_no_rea = 0
121	guez	3
122	guez	102	! lecture d'un fichier netcdf pour determiner le nombre de niveaux
123	guez	112	if (guide_u) then
124			call nf95_open('u.nc',Nf90_NOWRITe,ncid)
125			else if (guide_v) then
126			call nf95_open('v.nc',nf90_nowrite,ncid)
127			else if (guide_T) then
128			call nf95_open('T.nc',nf90_nowrite,ncid)
129			else
130			call nf95_open('hur.nc',nf90_nowrite, ncid)
131			end if
132	guez	3
133	guez	102	IF (ncep) THEN
134	guez	108	call nf95_inq_dimid(ncid, 'LEVEL', dimid)
135	guez	102	ELSE
136	guez	108	call nf95_inq_dimid(ncid, 'PRESSURE', dimid)
137	guez	102	END IF
138	guez	108	call nf95_inquire_dimension(ncid, dimid, nclen=nlev)
139	guez	102	PRINT *, 'nlev', nlev
140	guez	108	call nf95_close(ncid)
141	guez	102	! Lecture du premier etat des reanalyses.
142			CALL read_reanalyse(1, ps, ucovrea2, vcovrea2, tetarea2, qrea2, &
143			masserea2, nlev)
144			qrea2 = max(qrea2, 0.1)
145	guez	109
146	guez	112	if (guide_u) CALL writefield("alpha_u", alpha_u)
147			if (guide_t) CALL writefield("alpha_t", alpha_t)
148	guez	102	END IF first_call
149	guez	3
150	guez	102	! IMPORTATION DES VENTS, PRESSION ET TEMPERATURE REELS:
151	guez	3
152	guez	102	! Nudging fields are given 4 times per day:
153			IF (mod(itau, day_step / 4) == 0) THEN
154			vcovrea1 = vcovrea2
155			ucovrea1 = ucovrea2
156			tetarea1 = tetarea2
157			qrea1 = qrea2
158	guez	3
159	guez	102	PRINT *, 'LECTURE REANALYSES, pas ', step_rea, 'apres ', &
160			count_no_rea, ' non lectures'
161			step_rea = step_rea + 1
162			CALL read_reanalyse(step_rea, ps, ucovrea2, vcovrea2, tetarea2, qrea2, &
163			masserea2, nlev)
164			qrea2 = max(qrea2, 0.1)
165			factt = dtvr * iperiod / daysec
166	guez	112
167			if (guide_u) then
168			CALL writefield("ucov", ucov)
169			CALL writefield("ucovrea2", ucovrea2)
170			end if
171
172			if (guide_t) then
173			CALL writefield("teta", teta)
174			CALL writefield("tetarea2", tetarea2)
175			end if
176
177			if (guide_q) then
178			CALL writefield("qrea2", qrea2)
179			CALL writefield("q", q)
180			end if
181	guez	102	ELSE
182			count_no_rea = count_no_rea + 1
183			END IF
184	guez	3
185	guez	102	! Guidage
186	guez	3
187	guez	102	tau = mod(real(itau) / real(day_step / 4), 1.)
188	guez	3
189	guez	102	! x_gcm = a * x_gcm + (1 - a) * x_reanalyses
190	guez	3
191	guez	102	IF (guide_u) THEN
192			IF (itau == 0 .AND. ini_anal) then
193			ucov = ucovrea1
194			else
195			forall (l = 1: llm) ucov(:, :, l) = (1. - alpha_u) * ucov(:, :, l) &
196			+ alpha_u * ((1. - tau) * ucovrea1(:, :, l) &
197			+ tau * ucovrea2(:, :, l))
198			end IF
199			END IF
200	guez	3
201	guez	102	IF (guide_t) THEN
202			IF (itau == 0 .AND. ini_anal) then
203			teta = tetarea1
204			else
205			forall (l = 1: llm) teta(:, :, l) = (1. - alpha_t) * teta(:, :, l) &
206			+ alpha_t * ((1. - tau) * tetarea1(:, :, l) &
207			+ tau * tetarea2(:, :, l))
208			end IF
209			END IF
210	guez	3
211	guez	102	IF (guide_q) THEN
212			! Calcul de l'humidité saturante :
213			forall (l = 1: llm + 1) p(:, :, l) = ap(l) + bp(l) * ps
214			CALL exner_hyb(ps, p, pks, pk)
215			qsat = q_sat(pk * teta / cpp, preff * (pk / cpp)**(1. / kappa))
216	guez	3
217	guez	102	! humidité relative en % -> humidité spécifique
218			IF (itau == 0 .AND. ini_anal) then
219			q = qsat * qrea1 * 0.01
220			else
221			forall (l = 1: llm) q(:, :, l) = (1. - alpha_q) * q(:, :, l) &
222			+ alpha_q * (qsat(:, :, l) * ((1. - tau) * qrea1(:, :, l) &
223			+ tau * qrea2(:, :, l)) * 0.01)
224			end IF
225			END IF
226	guez	3
227	guez	102	IF (guide_v) THEN
228			IF (itau == 0 .AND. ini_anal) then
229			vcov = vcovrea1
230			else
231			forall (l = 1: llm) vcov(:, :, l) = (1. - alpha_v) * vcov(:, :, l) &
232			+ alpha_v * ((1. - tau) * vcovrea1(:, :, l) &
233			+ tau * vcovrea2(:, :, l))
234			end IF
235			END IF
236	guez	3
237	guez	29	END SUBROUTINE guide
238	guez	20
239			END MODULE guide_m