dyn3d/Guide/guide.f90

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

CONTAINS

  SUBROUTINE guide(itau, ucov, vcov, teta, q, ps)

    ! Author: F. Hourdin

    USE comconst, ONLY: cpp, kappa
    USE conf_gcm_m, ONLY: day_step
    use conf_guide_m, only: guide_u, guide_v, guide_t, guide_q, 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, factt
    USE dimens_m, ONLY: iim, jjm, llm
    USE disvert_m, ONLY: ap, bp, preff
    use dynetat0_m, only: grossismx, grossismy, rlatu, rlatv
    USE exner_hyb_m, ONLY: exner_hyb
    use init_tau2alpha_m, only: init_tau2alpha
    USE paramet_m, ONLY: iip1, jjp1
    USE q_sat_m, ONLY: q_sat
    use read_reanalyse_m, only: read_reanalyse
    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 r\'eanalyses

    REAL, save:: tetarea1(iim + 1, jjm + 1, llm)
    ! potential temperture from reanalysis
    
    REAL, save:: qrea1(iim + 1, jjm + 1, llm)

    REAL, save:: ucovrea2(iim + 1, jjm + 1, llm), vcovrea2(iim + 1, jjm, llm)
    ! vents covariants reanalyses

    REAL, save:: tetarea2(iim + 1, jjm + 1, llm)
    ! potential temperture from reanalysis
    
    REAL, save:: qrea2(iim + 1, jjm + 1, llm)

    ! alpha détermine la part des injections de données à chaque étape
    ! 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 l
    REAL tau

    ! TEST SUR QSAT
    REAL p(iim + 1, jjm + 1, llm + 1)
    real pk(iim + 1, jjm + 1, llm), pks(iim + 1, jjm + 1)
    REAL qsat(iim + 1, jjm + 1, llm)

    REAL dxdys(iip1, jjp1), dxdyu(iip1, jjp1), dxdyv(iip1, jjm)

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

    IF (itau == 0) THEN
       IF (online) THEN
          IF (abs(grossismx - 1.) < 0.1 .OR. abs(grossismy - 1.) < 0.1) THEN
             ! grille regulière
             if (guide_u) alpha_u = 1. - exp(- factt / tau_max_u)
             if (guide_v) alpha_v = 1. - exp(- factt / tau_max_v)
             if (guide_t) alpha_t = 1. - exp(- factt / tau_max_t)
             if (guide_q) alpha_q = 1. - exp(- factt / tau_max_q)
          else
             call init_tau2alpha(dxdys, dxdyu, dxdyv)

             if (guide_u) then
                CALL tau2alpha(dxdyu, rlatu, tau_min_u, tau_max_u, alpha_u)
                CALL writefield("alpha_u", alpha_u)
             end if

             if (guide_v) then
                CALL tau2alpha(dxdyv, rlatv, tau_min_v, tau_max_v, alpha_v)
                CALL writefield("alpha_v", alpha_v)
             end if

             if (guide_t) then
                CALL tau2alpha(dxdys, rlatu, tau_min_t, tau_max_t, alpha_t)
                CALL writefield("alpha_t", alpha_t)
             end if

             if (guide_q)  then
                CALL tau2alpha(dxdys, rlatu, tau_min_q, tau_max_q, alpha_q)
                CALL writefield("alpha_q", alpha_q)
             end if
          end IF
       ELSE
          ! Cas où on force exactement par les variables analysées
          if (guide_u) alpha_u = 1.
          if (guide_v) alpha_v = 1.
          if (guide_t) alpha_t = 1.
          if (guide_q) alpha_q = 1.
       END IF

       ! Lecture du premier état des réanalyses :
       CALL read_reanalyse(ps, ucovrea2, vcovrea2, tetarea2, qrea2)
       qrea2 = max(qrea2, 0.1)

       if (ini_anal) then
          IF (guide_u) ucov = ucovrea2
          IF (guide_v) vcov = vcovrea2
          IF (guide_t) teta = tetarea2

          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)
             q = q_sat(pk * teta / cpp, preff * (pk / cpp)**(1. / kappa)) &
                  * qrea2 * 0.01
          end IF
       end if
    END IF

    ! Importation des vents, pression et temp\'erature r\'eels :

    ! Nudging fields are given 4 times per day:
    IF (mod(itau, day_step / 4) == 0) THEN
       vcovrea1 = vcovrea2
       ucovrea1 = ucovrea2
       tetarea1 = tetarea2
       qrea1 = qrea2

       CALL read_reanalyse(ps, ucovrea2, vcovrea2, tetarea2, qrea2)
       qrea2 = max(qrea2, 0.1)

       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
    END IF

    ! Guidage

    tau = mod(real(itau) / real(day_step / 4), 1.)

    ! x_gcm = a * x_gcm + (1 - a) * x_reanalyses

    IF (guide_u) forall (l = 1: llm) ucov(:, :, l) = (1. - alpha_u) &
         * ucov(:, :, l) + alpha_u * ((1. - tau) * ucovrea1(:, :, l) + tau &
         * ucovrea2(:, :, l))

    IF (guide_v) forall (l = 1: llm) vcov(:, :, l) = (1. - alpha_v) &
         * vcov(:, :, l) + alpha_v * ((1. - tau) * vcovrea1(:, :, l) + tau &
         * vcovrea2(:, :, l))

    IF (guide_t) forall (l = 1: llm) teta(:, :, l) = (1. - alpha_t) &
         * teta(:, :, l) + alpha_t * ((1. - tau) * tetarea1(:, :, l) + tau &
         * tetarea2(:, :, l))

    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
       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 SUBROUTINE guide

END MODULE guide_m
1	guez	20	MODULE guide_m
2	guez	3
3	guez	210	! 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	20	CONTAINS
9	guez	3
10	guez	102	SUBROUTINE guide(itau, ucov, vcov, teta, q, ps)
11	guez	3
12	guez	210	! Author: F. Hourdin
13	guez	3
14	guez	115	USE comconst, ONLY: cpp, kappa
15			USE conf_gcm_m, ONLY: day_step
16	guez	173	use conf_guide_m, only: guide_u, guide_v, guide_t, guide_q, ini_anal, &
17			tau_min_u, tau_max_u, tau_min_v, tau_max_v, tau_min_t, tau_max_t, &
18			tau_min_q, tau_max_q, online, factt
19	guez	88	USE dimens_m, ONLY: iim, jjm, llm
20	guez	135	USE disvert_m, ONLY: ap, bp, preff
21	guez	139	use dynetat0_m, only: grossismx, grossismy, rlatu, rlatv
22	guez	83	USE exner_hyb_m, ONLY: exner_hyb
23	guez	115	use init_tau2alpha_m, only: init_tau2alpha
24	guez	210	USE paramet_m, ONLY: iip1, jjp1
25	guez	83	USE q_sat_m, ONLY: q_sat
26	guez	88	use read_reanalyse_m, only: read_reanalyse
27	guez	109	use tau2alpha_m, only: tau2alpha
28	guez	108	use writefield_m, only: writefield
29	guez	3
30	guez	83	INTEGER, INTENT(IN):: itau
31	guez	102	REAL, intent(inout):: ucov(:, :, :) ! (iim + 1, jjm + 1, llm) vent covariant
32			REAL, intent(inout):: vcov(:, :, :) ! (iim + 1, jjm, llm) ! vent covariant
33
34	guez	108	REAL, intent(inout):: teta(:, :, :) ! (iim + 1, jjm + 1, llm)
35			! température potentielle
36
37			REAL, intent(inout):: q(:, :, :) ! (iim + 1, jjm + 1, llm)
38	guez	88	REAL, intent(in):: ps(:, :) ! (iim + 1, jjm + 1) pression au sol
39	guez	3
40	guez	83	! Local:
41
42	guez	210	! Variables dynamiques pour les réanalyses
43	guez	102
44			REAL, save:: ucovrea1(iim + 1, jjm + 1, llm), vcovrea1(iim + 1, jjm, llm)
45	guez	210	! vents covariants r\'eanalyses
46	guez	102
47	guez	210	REAL, save:: tetarea1(iim + 1, jjm + 1, llm)
48			! potential temperture from reanalysis
49
50			REAL, save:: qrea1(iim + 1, jjm + 1, llm)
51	guez	102
52			REAL, save:: ucovrea2(iim + 1, jjm + 1, llm), vcovrea2(iim + 1, jjm, llm)
53			! vents covariants reanalyses
54
55	guez	210	REAL, save:: tetarea2(iim + 1, jjm + 1, llm)
56			! potential temperture from reanalysis
57
58			REAL, save:: qrea2(iim + 1, jjm + 1, llm)
59	guez	3
60	guez	116	! alpha détermine la part des injections de données à chaque étape
61	guez	112	! alpha=0 signifie pas d'injection
62			! alpha=1 signifie injection totale
63	guez	90	REAL, save:: alpha_q(iim + 1, jjm + 1)
64	guez	103	REAL, save:: alpha_t(iim + 1, jjm + 1)
65	guez	102	REAL, save:: alpha_u(iim + 1, jjm + 1), alpha_v(iim + 1, jjm)
66
67	guez	115	INTEGER l
68	guez	102	REAL tau
69	guez	3
70	guez	44	! TEST SUR QSAT
71	guez	210	REAL p(iim + 1, jjm + 1, llm + 1)
72	guez	90	real pk(iim + 1, jjm + 1, llm), pks(iim + 1, jjm + 1)
73			REAL qsat(iim + 1, jjm + 1, llm)
74	guez	3
75	guez	115	REAL dxdys(iip1, jjp1), dxdyu(iip1, jjp1), dxdyv(iip1, jjm)
76
77	guez	29	!-----------------------------------------------------------------------
78	guez	3
79	guez	172	IF (itau == 0) THEN
80	guez	102	IF (online) THEN
81	guez	115	IF (abs(grossismx - 1.) < 0.1 .OR. abs(grossismy - 1.) < 0.1) THEN
82			! grille regulière
83	guez	210	if (guide_u) alpha_u = 1. - exp(- factt / tau_max_u)
84			if (guide_v) alpha_v = 1. - exp(- factt / tau_max_v)
85			if (guide_t) alpha_t = 1. - exp(- factt / tau_max_t)
86			if (guide_q) alpha_q = 1. - exp(- factt / tau_max_q)
87	guez	115	else
88			call init_tau2alpha(dxdys, dxdyu, dxdyv)
89	guez	3
90	guez	115	if (guide_u) then
91			CALL tau2alpha(dxdyu, rlatu, tau_min_u, tau_max_u, alpha_u)
92			CALL writefield("alpha_u", alpha_u)
93			end if
94	guez	3
95	guez	115	if (guide_v) then
96			CALL tau2alpha(dxdyv, rlatv, tau_min_v, tau_max_v, alpha_v)
97			CALL writefield("alpha_v", alpha_v)
98			end if
99	guez	3
100	guez	115	if (guide_t) then
101			CALL tau2alpha(dxdys, rlatu, tau_min_t, tau_max_t, alpha_t)
102			CALL writefield("alpha_t", alpha_t)
103			end if
104
105			if (guide_q) then
106			CALL tau2alpha(dxdys, rlatu, tau_min_q, tau_max_q, alpha_q)
107			CALL writefield("alpha_q", alpha_q)
108			end if
109			end IF
110	guez	102	ELSE
111	guez	107	! Cas où on force exactement par les variables analysées
112	guez	115	if (guide_u) alpha_u = 1.
113			if (guide_v) alpha_v = 1.
114			if (guide_t) alpha_t = 1.
115	guez	112	if (guide_q) alpha_q = 1.
116	guez	102	END IF
117	guez	3
118	guez	172	! Lecture du premier état des réanalyses :
119			CALL read_reanalyse(ps, ucovrea2, vcovrea2, tetarea2, qrea2)
120			qrea2 = max(qrea2, 0.1)
121	guez	3
122	guez	172	if (ini_anal) then
123			IF (guide_u) ucov = ucovrea2
124			IF (guide_v) vcov = vcovrea2
125			IF (guide_t) teta = tetarea2
126	guez	115
127	guez	172	IF (guide_q) then
128			! Calcul de l'humidité saturante :
129			forall (l = 1: llm + 1) p(:, :, l) = ap(l) + bp(l) * ps
130			CALL exner_hyb(ps, p, pks, pk)
131			q = q_sat(pk * teta / cpp, preff * (pk / cpp)**(1. / kappa)) &
132			* qrea2 * 0.01
133			end IF
134	guez	112	end if
135	guez	172	END IF
136	guez	3
137	guez	172	! Importation des vents, pression et temp\'erature r\'eels :
138	guez	115
139	guez	102	! Nudging fields are given 4 times per day:
140			IF (mod(itau, day_step / 4) == 0) THEN
141			vcovrea1 = vcovrea2
142			ucovrea1 = ucovrea2
143			tetarea1 = tetarea2
144			qrea1 = qrea2
145	guez	3
146	guez	172	CALL read_reanalyse(ps, ucovrea2, vcovrea2, tetarea2, qrea2)
147	guez	102	qrea2 = max(qrea2, 0.1)
148	guez	112
149			if (guide_u) then
150			CALL writefield("ucov", ucov)
151			CALL writefield("ucovrea2", ucovrea2)
152			end if
153
154			if (guide_t) then
155			CALL writefield("teta", teta)
156			CALL writefield("tetarea2", tetarea2)
157			end if
158
159			if (guide_q) then
160			CALL writefield("qrea2", qrea2)
161			CALL writefield("q", q)
162			end if
163	guez	102	END IF
164	guez	3
165	guez	102	! Guidage
166	guez	3
167	guez	102	tau = mod(real(itau) / real(day_step / 4), 1.)
168	guez	3
169	guez	102	! x_gcm = a * x_gcm + (1 - a) * x_reanalyses
170	guez	3
171	guez	172	IF (guide_u) forall (l = 1: llm) ucov(:, :, l) = (1. - alpha_u) &
172			* ucov(:, :, l) + alpha_u * ((1. - tau) * ucovrea1(:, :, l) + tau &
173			* ucovrea2(:, :, l))
174	guez	3
175	guez	172	IF (guide_v) forall (l = 1: llm) vcov(:, :, l) = (1. - alpha_v) &
176			* vcov(:, :, l) + alpha_v * ((1. - tau) * vcovrea1(:, :, l) + tau &
177			* vcovrea2(:, :, l))
178	guez	3
179	guez	172	IF (guide_t) forall (l = 1: llm) teta(:, :, l) = (1. - alpha_t) &
180			* teta(:, :, l) + alpha_t * ((1. - tau) * tetarea1(:, :, l) + tau &
181			* tetarea2(:, :, l))
182
183	guez	102	IF (guide_q) THEN
184			! Calcul de l'humidité saturante :
185			forall (l = 1: llm + 1) p(:, :, l) = ap(l) + bp(l) * ps
186			CALL exner_hyb(ps, p, pks, pk)
187			qsat = q_sat(pk * teta / cpp, preff * (pk / cpp)**(1. / kappa))
188	guez	3
189	guez	102	! humidité relative en % -> humidité spécifique
190	guez	172	forall (l = 1: llm) q(:, :, l) = (1. - alpha_q) * q(:, :, l) &
191			+ alpha_q * (qsat(:, :, l) * ((1. - tau) * qrea1(:, :, l) &
192			+ tau * qrea2(:, :, l)) * 0.01)
193	guez	102	END IF
194	guez	3
195	guez	29	END SUBROUTINE guide
196	guez	20
197			END MODULE guide_m