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

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, 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, 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 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 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 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, save:: step_rea, count_no_rea

    INTEGER 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)

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

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

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

    first_call: 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 = factt / tau_max_u
             if (guide_v) alpha_v = factt / tau_max_v
             if (guide_t) alpha_t = factt / tau_max_t
             if (guide_q) alpha_q = 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

       step_rea = 1
       count_no_rea = 0

       ! Lecture d'un fichier NetCDF pour d\'eterminer 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 état des réanalyses :
       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 fichiers guidage, 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)

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