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=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

    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)

       CALL writefield("alpha_u", alpha_u)
       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
       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	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			! 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	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	103	CALL tau2alpha(3, factt, tau_min_v, tau_max_v, alpha_v)
108			CALL tau2alpha(2, factt, tau_min_u, tau_max_u, alpha_u)
109			CALL tau2alpha(1, factt, tau_min_t, tau_max_t, alpha_t)
110			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	102	alpha_t = 0.
114			alpha_u = 0.
115			alpha_v = 0.
116	guez	103	alpha_q = 0.
117	guez	102	END IF
118	guez	3
119	guez	102	step_rea = 1
120			count_no_rea = 0
121	guez	108	ncid = -99
122	guez	3
123	guez	102	! lecture d'un fichier netcdf pour determiner le nombre de niveaux
124	guez	108	if (guide_u) call nf95_open('u.nc',Nf90_NOWRITe,ncid)
125			if (guide_v) call nf95_open('v.nc',nf90_nowrite,ncid)
126			if (guide_T) call nf95_open('T.nc',nf90_nowrite,ncid)
127			if (guide_Q) call nf95_open('hur.nc',nf90_nowrite, ncid)
128	guez	3
129	guez	102	IF (ncep) THEN
130	guez	108	call nf95_inq_dimid(ncid, 'LEVEL', dimid)
131	guez	102	ELSE
132	guez	108	call nf95_inq_dimid(ncid, 'PRESSURE', dimid)
133	guez	102	END IF
134	guez	108	call nf95_inquire_dimension(ncid, dimid, nclen=nlev)
135	guez	102	PRINT *, 'nlev', nlev
136	guez	108	call nf95_close(ncid)
137	guez	102	! Lecture du premier etat des reanalyses.
138			CALL read_reanalyse(1, ps, ucovrea2, vcovrea2, tetarea2, qrea2, &
139			masserea2, nlev)
140			qrea2 = max(qrea2, 0.1)
141	guez	109
142			CALL writefield("alpha_u", alpha_u)
143			CALL writefield("alpha_t", alpha_t)
144	guez	102	END IF first_call
145	guez	3
146	guez	102	! IMPORTATION DES VENTS, PRESSION ET TEMPERATURE REELS:
147	guez	3
148	guez	102	! Nudging fields are given 4 times per day:
149			IF (mod(itau, day_step / 4) == 0) THEN
150			vcovrea1 = vcovrea2
151			ucovrea1 = ucovrea2
152			tetarea1 = tetarea2
153			qrea1 = qrea2
154	guez	3
155	guez	102	PRINT *, 'LECTURE REANALYSES, pas ', step_rea, 'apres ', &
156			count_no_rea, ' non lectures'
157			step_rea = step_rea + 1
158			CALL read_reanalyse(step_rea, ps, ucovrea2, vcovrea2, tetarea2, qrea2, &
159			masserea2, nlev)
160			qrea2 = max(qrea2, 0.1)
161			factt = dtvr * iperiod / daysec
162	guez	108	CALL writefield("ucov", ucov)
163			CALL writefield("ucovrea2", ucovrea2)
164			CALL writefield("teta", teta)
165			CALL writefield("tetarea2", tetarea2)
166			CALL writefield("qrea2", qrea2)
167			CALL writefield("q", q)
168	guez	102	ELSE
169			count_no_rea = count_no_rea + 1
170			END IF
171	guez	3
172	guez	102	! Guidage
173	guez	3
174	guez	102	tau = mod(real(itau) / real(day_step / 4), 1.)
175	guez	3
176	guez	102	! x_gcm = a * x_gcm + (1 - a) * x_reanalyses
177	guez	3
178	guez	102	IF (guide_u) THEN
179			IF (itau == 0 .AND. ini_anal) then
180			ucov = ucovrea1
181			else
182			forall (l = 1: llm) ucov(:, :, l) = (1. - alpha_u) * ucov(:, :, l) &
183			+ alpha_u * ((1. - tau) * ucovrea1(:, :, l) &
184			+ tau * ucovrea2(:, :, l))
185			end IF
186			END IF
187	guez	3
188	guez	102	IF (guide_t) THEN
189			IF (itau == 0 .AND. ini_anal) then
190			teta = tetarea1
191			else
192			forall (l = 1: llm) teta(:, :, l) = (1. - alpha_t) * teta(:, :, l) &
193			+ alpha_t * ((1. - tau) * tetarea1(:, :, l) &
194			+ tau * tetarea2(:, :, l))
195			end IF
196			END IF
197	guez	3
198	guez	102	IF (guide_q) THEN
199			! Calcul de l'humidité saturante :
200			forall (l = 1: llm + 1) p(:, :, l) = ap(l) + bp(l) * ps
201			CALL exner_hyb(ps, p, pks, pk)
202			qsat = q_sat(pk * teta / cpp, preff * (pk / cpp)**(1. / kappa))
203	guez	3
204	guez	102	! humidité relative en % -> humidité spécifique
205			IF (itau == 0 .AND. ini_anal) then
206			q = qsat * qrea1 * 0.01
207			else
208			forall (l = 1: llm) q(:, :, l) = (1. - alpha_q) * q(:, :, l) &
209			+ alpha_q * (qsat(:, :, l) * ((1. - tau) * qrea1(:, :, l) &
210			+ tau * qrea2(:, :, l)) * 0.01)
211			end IF
212			END IF
213	guez	3
214	guez	102	IF (guide_v) THEN
215			IF (itau == 0 .AND. ini_anal) then
216			vcov = vcovrea1
217			else
218			forall (l = 1: llm) vcov(:, :, l) = (1. - alpha_v) * vcov(:, :, l) &
219			+ alpha_v * ((1. - tau) * vcovrea1(:, :, l) &
220			+ tau * vcovrea2(:, :, l))
221			end IF
222			END IF
223	guez	3
224	guez	29	END SUBROUTINE guide
225	guez	20
226			END MODULE guide_m