GCC Code Coverage Report


Directory: ./
File: dyn/guide_mod.f90
Date: 2022-01-11 19:19:34
Exec Total Coverage
Lines: 0 940 0.0%
Branches: 0 1572 0.0%

Line Branch Exec Source
1 !
2 ! $Id$
3 !
4 MODULE guide_mod
5
6 !=======================================================================
7 ! Auteur: F.Hourdin
8 ! F. Codron 01/09
9 !=======================================================================
10
11 USE getparam, only: ini_getparam, fin_getparam, getpar
12 USE Write_Field
13 use netcdf, only: nf90_nowrite, nf90_open, nf90_inq_varid, nf90_close, &
14 nf90_inq_dimid, nf90_inquire_dimension
15 use pres2lev_mod, only: pres2lev
16
17 IMPLICIT NONE
18
19 ! ---------------------------------------------
20 ! Declarations des cles logiques et parametres
21 ! ---------------------------------------------
22 INTEGER, PRIVATE, SAVE :: iguide_read,iguide_int,iguide_sav
23 INTEGER, PRIVATE, SAVE :: nlevnc, guide_plevs
24 LOGICAL, PRIVATE, SAVE :: guide_u,guide_v,guide_T,guide_Q,guide_P
25 LOGICAL, PRIVATE, SAVE :: guide_hr,guide_teta
26 LOGICAL, PRIVATE, SAVE :: guide_BL,guide_reg,guide_add,gamma4,guide_zon
27 LOGICAL, PRIVATE, SAVE :: invert_p,invert_y,ini_anal
28 LOGICAL, PRIVATE, SAVE :: guide_2D,guide_sav,guide_modele
29 !FC
30 LOGICAL, PRIVATE, SAVE :: convert_Pa
31
32 REAL, PRIVATE, SAVE :: tau_min_u,tau_max_u
33 REAL, PRIVATE, SAVE :: tau_min_v,tau_max_v
34 REAL, PRIVATE, SAVE :: tau_min_T,tau_max_T
35 REAL, PRIVATE, SAVE :: tau_min_Q,tau_max_Q
36 REAL, PRIVATE, SAVE :: tau_min_P,tau_max_P
37
38 REAL, PRIVATE, SAVE :: lat_min_g,lat_max_g
39 REAL, PRIVATE, SAVE :: lon_min_g,lon_max_g
40 REAL, PRIVATE, SAVE :: tau_lon,tau_lat
41
42 REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE :: alpha_u,alpha_v
43 REAL, ALLOCATABLE, DIMENSION(:, :), PRIVATE, SAVE :: alpha_T,alpha_Q
44 REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE :: alpha_P,alpha_pcor
45
46 ! ---------------------------------------------
47 ! Variables de guidage
48 ! ---------------------------------------------
49 ! Variables des fichiers de guidage
50 REAL, ALLOCATABLE, DIMENSION(:,:,:), PRIVATE, SAVE :: unat1,unat2
51 REAL, ALLOCATABLE, DIMENSION(:,:,:), PRIVATE, SAVE :: vnat1,vnat2
52 REAL, ALLOCATABLE, DIMENSION(:,:,:), PRIVATE, SAVE :: tnat1,tnat2
53 REAL, ALLOCATABLE, DIMENSION(:,:,:), PRIVATE, SAVE :: qnat1,qnat2
54 REAL, ALLOCATABLE, DIMENSION(:,:,:), PRIVATE, SAVE :: pnat1,pnat2
55 REAL, ALLOCATABLE, DIMENSION(:,:), PRIVATE, SAVE :: psnat1,psnat2
56 REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE :: apnc,bpnc
57 ! Variables aux dimensions du modele
58 REAL, ALLOCATABLE, DIMENSION(:,:), PRIVATE, SAVE :: ugui1,ugui2
59 REAL, ALLOCATABLE, DIMENSION(:,:), PRIVATE, SAVE :: vgui1,vgui2
60 REAL, ALLOCATABLE, DIMENSION(:,:), PRIVATE, SAVE :: tgui1,tgui2
61 REAL, ALLOCATABLE, DIMENSION(:,:), PRIVATE, SAVE :: qgui1,qgui2
62 REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE :: psgui1,psgui2
63
64 CONTAINS
65 !=======================================================================
66
67 ✗ SUBROUTINE guide_init
68
69 USE control_mod, ONLY: day_step
70 USE serre_mod, ONLY: grossismx
71
72 IMPLICIT NONE
73
74 INCLUDE "dimensions.h"
75 INCLUDE "paramet.h"
76 INCLUDE "netcdf.inc"
77
78 INTEGER :: error,ncidpl,rid,rcod
79 CHARACTER (len = 80) :: abort_message
80 CHARACTER (len = 20) :: modname = 'guide_init'
81 CHARACTER (len = 20) :: namedim
82
83 ! ---------------------------------------------
84 ! Lecture des parametres:
85 ! ---------------------------------------------
86 ✗ call ini_getparam("nudging_parameters_out.txt")
87 ! Variables guidees
88 ✗ CALL getpar('guide_u',.true.,guide_u,'guidage de u')
89 ✗ CALL getpar('guide_v',.true.,guide_v,'guidage de v')
90 ✗ CALL getpar('guide_T',.true.,guide_T,'guidage de T')
91 ✗ CALL getpar('guide_P',.true.,guide_P,'guidage de P')
92 ✗ CALL getpar('guide_Q',.true.,guide_Q,'guidage de Q')
93 ✗ CALL getpar('guide_hr',.true.,guide_hr,'guidage de Q par H.R')
94 ✗ CALL getpar('guide_teta',.false.,guide_teta,'guidage de T par Teta')
95
96 ✗ CALL getpar('guide_add',.false.,guide_add,'for�age constant?')
97 ✗ CALL getpar('guide_zon',.false.,guide_zon,'guidage moy zonale')
98 ✗ if (guide_zon .and. abs(grossismx - 1.) > 0.01) &
99 call abort_gcm("guide_init", &
100 ✗ "zonal nudging requires grid regular in longitude", 1)
101
102 ! Constantes de rappel. Unite : fraction de jour
103 ✗ CALL getpar('tau_min_u',0.02,tau_min_u,'Cste de rappel min, u')
104 ✗ CALL getpar('tau_max_u', 10.,tau_max_u,'Cste de rappel max, u')
105 ✗ CALL getpar('tau_min_v',0.02,tau_min_v,'Cste de rappel min, v')
106 ✗ CALL getpar('tau_max_v', 10.,tau_max_v,'Cste de rappel max, v')
107 ✗ CALL getpar('tau_min_T',0.02,tau_min_T,'Cste de rappel min, T')
108 ✗ CALL getpar('tau_max_T', 10.,tau_max_T,'Cste de rappel max, T')
109 ✗ CALL getpar('tau_min_Q',0.02,tau_min_Q,'Cste de rappel min, Q')
110 ✗ CALL getpar('tau_max_Q', 10.,tau_max_Q,'Cste de rappel max, Q')
111 ✗ CALL getpar('tau_min_P',0.02,tau_min_P,'Cste de rappel min, P')
112 ✗ CALL getpar('tau_max_P', 10.,tau_max_P,'Cste de rappel max, P')
113 ✗ CALL getpar('gamma4',.false.,gamma4,'Zone sans rappel elargie')
114 ✗ CALL getpar('guide_BL',.true.,guide_BL,'guidage dans C.Lim')
115
116 ! Sauvegarde du for�age
117 ✗ CALL getpar('guide_sav',.false.,guide_sav,'sauvegarde guidage')
118 ✗ CALL getpar('iguide_sav',4,iguide_sav,'freq. sauvegarde guidage')
119 ! frequences f>0: fx/jour; f<0: tous les f jours; f=0: 1 seule fois.
120 ✗ IF (iguide_sav.GT.0) THEN
121 ✗ iguide_sav=day_step/iguide_sav
122 ✗ ELSE if (iguide_sav == 0) then
123 ✗ iguide_sav = huge(0)
124 ELSE
125 ✗ iguide_sav=day_step*iguide_sav
126 ENDIF
127
128 ! Guidage regional seulement (sinon constant ou suivant le zoom)
129 ✗ CALL getpar('guide_reg',.false.,guide_reg,'guidage regional')
130 ✗ CALL getpar('lat_min_g',-90.,lat_min_g,'Latitude mini guidage ')
131 ✗ CALL getpar('lat_max_g', 90.,lat_max_g,'Latitude maxi guidage ')
132 ✗ CALL getpar('lon_min_g',-180.,lon_min_g,'longitude mini guidage ')
133 ✗ CALL getpar('lon_max_g', 180.,lon_max_g,'longitude maxi guidage ')
134 ✗ CALL getpar('tau_lat', 5.,tau_lat,'raideur lat guide regional ')
135 ✗ CALL getpar('tau_lon', 5.,tau_lon,'raideur lon guide regional ')
136
137 ! Parametres pour lecture des fichiers
138 ✗ CALL getpar('iguide_read',4,iguide_read,'freq. lecture guidage')
139 ✗ CALL getpar('iguide_int',4,iguide_int,'freq. interpolation vert')
140 ✗ IF (iguide_int.EQ.0) THEN
141 ✗ iguide_int=1
142 ✗ ELSEIF (iguide_int.GT.0) THEN
143 ✗ iguide_int=day_step/iguide_int
144 ELSE
145 ✗ iguide_int=day_step*iguide_int
146 ENDIF
147 ✗ CALL getpar('guide_plevs',0,guide_plevs,'niveaux pression fichiers guidage')
148 ! Pour compatibilite avec ancienne version avec guide_modele
149 ✗ CALL getpar('guide_modele',.false.,guide_modele,'niveaux pression ap+bp*psol')
150 ✗ IF (guide_modele) THEN
151 ✗ guide_plevs=1
152 ENDIF
153 !FC
154 ✗ CALL getpar('convert_Pa',.true.,convert_Pa,'Convert Pressure levels in Pa')
155 ! Fin raccord
156 ✗ CALL getpar('ini_anal',.false.,ini_anal,'Etat initial = analyse')
157 ✗ CALL getpar('guide_invertp',.true.,invert_p,'niveaux p inverses')
158 ✗ CALL getpar('guide_inverty',.true.,invert_y,'inversion N-S')
159 ✗ CALL getpar('guide_2D',.false.,guide_2D,'fichier guidage lat-P')
160
161 ✗ call fin_getparam
162
163 ! ---------------------------------------------
164 ! Determination du nombre de niveaux verticaux
165 ! des fichiers guidage
166 ! ---------------------------------------------
167 ✗ ncidpl=-99
168 ✗ if (guide_plevs.EQ.1) then
169 if (ncidpl.eq.-99) then
170 ✗ rcod=nf90_open('apbp.nc',Nf90_NOWRITe, ncidpl)
171 ✗ if (rcod.NE.NF_NOERR) THEN
172 ✗ abort_message=' Nudging error -> no file apbp.nc'
173 ✗ CALL abort_gcm(modname,abort_message,1)
174 endif
175 endif
176 ✗ elseif (guide_plevs.EQ.2) then
177 if (ncidpl.EQ.-99) then
178 ✗ rcod=nf90_open('P.nc',Nf90_NOWRITe,ncidpl)
179 ✗ if (rcod.NE.NF_NOERR) THEN
180 ✗ abort_message=' Nudging error -> no file P.nc'
181 ✗ CALL abort_gcm(modname,abort_message,1)
182 endif
183 endif
184
185 ✗ elseif (guide_u) then
186 if (ncidpl.eq.-99) then
187 ✗ rcod=nf90_open('u.nc',Nf90_NOWRITe,ncidpl)
188 ✗ if (rcod.NE.NF_NOERR) THEN
189 CALL abort_gcm(modname, &
190 ✗ ' Nudging error -> no file u.nc',1)
191 endif
192 endif
193
194 ✗ elseif (guide_v) then
195 if (ncidpl.eq.-99) then
196 ✗ rcod=nf90_open('v.nc',nf90_nowrite,ncidpl)
197 ✗ if (rcod.NE.NF_NOERR) THEN
198 CALL abort_gcm(modname, &
199 ✗ ' Nudging error -> no file v.nc',1)
200 endif
201 endif
202 ✗ elseif (guide_T) then
203 if (ncidpl.eq.-99) then
204 ✗ rcod=nf90_open('T.nc',nf90_nowrite,ncidpl)
205 ✗ if (rcod.NE.NF_NOERR) THEN
206 CALL abort_gcm(modname, &
207 ✗ ' Nudging error -> no file T.nc',1)
208 endif
209 endif
210 ✗ elseif (guide_Q) then
211 if (ncidpl.eq.-99) then
212 ✗ rcod=nf90_open('hur.nc',nf90_nowrite, ncidpl)
213 ✗ if (rcod.NE.NF_NOERR) THEN
214 CALL abort_gcm(modname, &
215 ✗ ' Nudging error -> no file hur.nc',1)
216 endif
217 endif
218
219
220 endif
221 ✗ error=NF_INQ_DIMID(ncidpl,'LEVEL',rid)
222 ✗ IF (error.NE.NF_NOERR) error=NF_INQ_DIMID(ncidpl,'PRESSURE',rid)
223 ✗ IF (error.NE.NF_NOERR) THEN
224 ✗ CALL abort_gcm(modname,'Nudging: error reading pressure levels',1)
225 ENDIF
226 ✗ error=NF_INQ_DIMLEN(ncidpl,rid,nlevnc)
227 ✗ write(*,*)trim(modname)//' : number of vertical levels nlevnc', nlevnc
228 ✗ rcod = nf90_close(ncidpl)
229
230 ! ---------------------------------------------
231 ! Allocation des variables
232 ! ---------------------------------------------
233 ✗ abort_message='nudging allocation error'
234
235 ✗ ALLOCATE(apnc(nlevnc), stat = error)
236 ✗ IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
237 ✗ ALLOCATE(bpnc(nlevnc), stat = error)
238 ✗ IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
239 ✗ apnc=0.;bpnc=0.
240
241 ✗ ALLOCATE(alpha_pcor(llm), stat = error)
242 ✗ IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
243 ✗ ALLOCATE(alpha_u(ip1jmp1), stat = error)
244 ✗ IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
245 ✗ ALLOCATE(alpha_v(ip1jm), stat = error)
246 ✗ IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
247 ✗ ALLOCATE(alpha_T(iip1, jjp1), stat = error)
248 ✗ IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
249 ✗ ALLOCATE(alpha_Q(iip1, jjp1), stat = error)
250 ✗ IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
251 ✗ ALLOCATE(alpha_P(ip1jmp1), stat = error)
252 ✗ IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
253 ✗ alpha_u=0.;alpha_v=0;alpha_T=0;alpha_Q=0;alpha_P=0
254
255 ✗ IF (guide_u) THEN
256 ✗ ALLOCATE(unat1(iip1,jjp1,nlevnc), stat = error)
257 ✗ IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
258 ✗ ALLOCATE(ugui1(ip1jmp1,llm), stat = error)
259 ✗ IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
260 ✗ ALLOCATE(unat2(iip1,jjp1,nlevnc), stat = error)
261 ✗ IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
262 ✗ ALLOCATE(ugui2(ip1jmp1,llm), stat = error)
263 ✗ IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
264 ✗ unat1=0.;unat2=0.;ugui1=0.;ugui2=0.
265 ENDIF
266
267 ✗ IF (guide_T) THEN
268 ✗ ALLOCATE(tnat1(iip1,jjp1,nlevnc), stat = error)
269 ✗ IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
270 ✗ ALLOCATE(tgui1(ip1jmp1,llm), stat = error)
271 ✗ IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
272 ✗ ALLOCATE(tnat2(iip1,jjp1,nlevnc), stat = error)
273 ✗ IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
274 ✗ ALLOCATE(tgui2(ip1jmp1,llm), stat = error)
275 ✗ IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
276 ✗ tnat1=0.;tnat2=0.;tgui1=0.;tgui2=0.
277 ENDIF
278
279 ✗ IF (guide_Q) THEN
280 ✗ ALLOCATE(qnat1(iip1,jjp1,nlevnc), stat = error)
281 ✗ IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
282 ✗ ALLOCATE(qgui1(ip1jmp1,llm), stat = error)
283 ✗ IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
284 ✗ ALLOCATE(qnat2(iip1,jjp1,nlevnc), stat = error)
285 ✗ IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
286 ✗ ALLOCATE(qgui2(ip1jmp1,llm), stat = error)
287 ✗ IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
288 ✗ qnat1=0.;qnat2=0.;qgui1=0.;qgui2=0.
289 ENDIF
290
291 ✗ IF (guide_v) THEN
292 ✗ ALLOCATE(vnat1(iip1,jjm,nlevnc), stat = error)
293 ✗ IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
294 ✗ ALLOCATE(vgui1(ip1jm,llm), stat = error)
295 ✗ IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
296 ✗ ALLOCATE(vnat2(iip1,jjm,nlevnc), stat = error)
297 ✗ IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
298 ✗ ALLOCATE(vgui2(ip1jm,llm), stat = error)
299 ✗ IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
300 ✗ vnat1=0.;vnat2=0.;vgui1=0.;vgui2=0.
301 ENDIF
302
303 ✗ IF (guide_plevs.EQ.2) THEN
304 ✗ ALLOCATE(pnat1(iip1,jjp1,nlevnc), stat = error)
305 ✗ IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
306 ✗ ALLOCATE(pnat2(iip1,jjp1,nlevnc), stat = error)
307 ✗ IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
308 ✗ pnat1=0.;pnat2=0.;
309 ENDIF
310
311 ✗ IF (guide_P.OR.guide_plevs.EQ.1) THEN
312 ✗ ALLOCATE(psnat1(iip1,jjp1), stat = error)
313 ✗ IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
314 ✗ ALLOCATE(psnat2(iip1,jjp1), stat = error)
315 ✗ IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
316 ✗ psnat1=0.;psnat2=0.;
317 ENDIF
318 ✗ IF (guide_P) THEN
319 ✗ ALLOCATE(psgui2(ip1jmp1), stat = error)
320 ✗ IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
321 ✗ ALLOCATE(psgui1(ip1jmp1), stat = error)
322 ✗ IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
323 ✗ psgui1=0.;psgui2=0.
324 ENDIF
325
326 ! ---------------------------------------------
327 ! Lecture du premier etat de guidage.
328 ! ---------------------------------------------
329 ✗ IF (guide_2D) THEN
330 ✗ CALL guide_read2D(1)
331 ELSE
332 ✗ CALL guide_read(1)
333 ENDIF
334 ✗ IF (guide_v) vnat1=vnat2
335 ✗ IF (guide_u) unat1=unat2
336 ✗ IF (guide_T) tnat1=tnat2
337 ✗ IF (guide_Q) qnat1=qnat2
338 ✗ IF (guide_plevs.EQ.2) pnat1=pnat2
339 ✗ IF (guide_P.OR.guide_plevs.EQ.1) psnat1=psnat2
340
341 ✗ END SUBROUTINE guide_init
342
343 !=======================================================================
344 ✗ SUBROUTINE guide_main(itau,ucov,vcov,teta,q,masse,ps)
345
346 USE exner_hyb_m, ONLY: exner_hyb
347 USE exner_milieu_m, ONLY: exner_milieu
348 USE control_mod, ONLY: day_step, iperiod
349 USE comconst_mod, ONLY: cpp, dtvr, daysec,kappa
350 USE comvert_mod, ONLY: ap, bp, preff, presnivs, pressure_exner
351
352 IMPLICIT NONE
353
354 INCLUDE "dimensions.h"
355 INCLUDE "paramet.h"
356
357 ! Variables entree
358 INTEGER, INTENT(IN) :: itau !pas de temps
359 REAL, DIMENSION (ip1jmp1,llm), INTENT(INOUT) :: ucov,teta,q,masse
360 REAL, DIMENSION (ip1jm,llm), INTENT(INOUT) :: vcov
361 REAL, DIMENSION (ip1jmp1), INTENT(INOUT) :: ps
362
363 ! Variables locales
364 LOGICAL, SAVE :: first=.TRUE.
365 LOGICAL :: f_out ! sortie guidage
366 REAL, DIMENSION (ip1jmp1,llm) :: f_add ! var aux: champ de guidage
367 REAL :: pk(ip1jmp1,llm) ! Exner at mid-layers
368 REAL :: pks(ip1jmp1) ! Exner at the surface
369 REAL :: unskap ! 1./kappa
370 REAL, DIMENSION (ip1jmp1,llmp1) :: p ! Pressure at inter-layers
371 ! Compteurs temps:
372 INTEGER, SAVE :: step_rea,count_no_rea,itau_test ! lecture guidage
373 REAL :: ditau, dday_step
374 REAL :: tau,reste ! position entre 2 etats de guidage
375 REAL, SAVE :: factt ! pas de temps en fraction de jour
376
377 INTEGER :: l
378 CHARACTER(LEN=20) :: modname="guide_main"
379
380 !-----------------------------------------------------------------------
381 ! Initialisations au premier passage
382 !-----------------------------------------------------------------------
383 ✗ IF (first) THEN
384 ✗ first=.FALSE.
385 ✗ CALL guide_init
386 ✗ itau_test=1001
387 ✗ step_rea=1
388 ✗ count_no_rea=0
389 ! Calcul des constantes de rappel
390 ✗ factt=dtvr*iperiod/daysec
391 ✗ call tau2alpha(3,iip1,jjm ,factt,tau_min_v,tau_max_v,alpha_v)
392 ✗ call tau2alpha(2,iip1,jjp1,factt,tau_min_u,tau_max_u,alpha_u)
393 ✗ call tau2alpha(1,iip1,jjp1,factt,tau_min_T,tau_max_T,alpha_T)
394 ✗ call tau2alpha(1,iip1,jjp1,factt,tau_min_P,tau_max_P,alpha_P)
395 ✗ call tau2alpha(1,iip1,jjp1,factt,tau_min_Q,tau_max_Q,alpha_Q)
396 ! correction de rappel dans couche limite
397 ✗ if (guide_BL) then
398 ✗ alpha_pcor(:)=1.
399 else
400 ✗ do l=1,llm
401 ✗ alpha_pcor(l)=(1.+tanh((0.85-presnivs(l)/preff)/0.05))/2.
402 enddo
403 endif
404 ! ini_anal: etat initial egal au guidage
405 ✗ IF (ini_anal) THEN
406 ✗ CALL guide_interp(ps,teta)
407 ✗ IF (guide_u) ucov=ugui2
408 ✗ IF (guide_v) vcov=ugui2
409 ✗ IF (guide_T) teta=tgui2
410 ✗ IF (guide_Q) q=qgui2
411 ✗ IF (guide_P) THEN
412 ✗ ps=psgui2
413 ✗ CALL pression(ip1jmp1,ap,bp,ps,p)
414 ✗ CALL massdair(p,masse)
415 ENDIF
416 ✗ RETURN
417 ENDIF
418 ! Verification structure guidage
419 ! IF (guide_u) THEN
420 ! CALL writefield('unat',unat1)
421 ! CALL writefield('ucov',RESHAPE(ucov,(/iip1,jjp1,llm/)))
422 ! ENDIF
423 ! IF (guide_T) THEN
424 ! CALL writefield('tnat',tnat1)
425 ! CALL writefield('teta',RESHAPE(teta,(/iip1,jjp1,llm/)))
426 ! ENDIF
427
428 ENDIF !first
429
430 !-----------------------------------------------------------------------
431 ! Lecture des fichiers de guidage ?
432 !-----------------------------------------------------------------------
433 ✗ IF (iguide_read.NE.0) THEN
434 ✗ ditau=real(itau)
435 ✗ dday_step=real(day_step)
436 ✗ IF (iguide_read.LT.0) THEN
437 ✗ tau=ditau/dday_step/REAL(iguide_read)
438 ELSE
439 ✗ tau=REAL(iguide_read)*ditau/dday_step
440 ENDIF
441 ✗ reste=tau-AINT(tau)
442 ✗ IF (reste.EQ.0.) THEN
443 ✗ IF (itau_test.EQ.itau) THEN
444 ✗ write(*,*)trim(modname)//' second pass in advreel at itau=',&
445 ✗ itau
446 ✗ stop
447 ELSE
448 ✗ IF (guide_v) vnat1=vnat2
449 ✗ IF (guide_u) unat1=unat2
450 ✗ IF (guide_T) tnat1=tnat2
451 ✗ IF (guide_Q) qnat1=qnat2
452 ✗ IF (guide_plevs.EQ.2) pnat1=pnat2
453 ✗ IF (guide_P.OR.guide_plevs.EQ.1) psnat1=psnat2
454 ✗ step_rea=step_rea+1
455 ✗ itau_test=itau
456 ✗ write(*,*)trim(modname)//' Reading nudging files, step ',&
457 ✗ step_rea,'after ',count_no_rea,' skips'
458 ✗ IF (guide_2D) THEN
459 ✗ CALL guide_read2D(step_rea)
460 ELSE
461 ✗ CALL guide_read(step_rea)
462 ENDIF
463 ✗ count_no_rea=0
464 ENDIF
465 ELSE
466 ✗ count_no_rea=count_no_rea+1
467
468 ENDIF
469 ENDIF !iguide_read=0
470
471 !-----------------------------------------------------------------------
472 ! Interpolation et conversion des champs de guidage
473 !-----------------------------------------------------------------------
474 ✗ IF (MOD(itau,iguide_int).EQ.0) THEN
475 ✗ CALL guide_interp(ps,teta)
476 ENDIF
477 ! Repartition entre 2 etats de guidage
478 ✗ IF (iguide_read.NE.0) THEN
479 tau=reste
480 ELSE
481 tau=1.
482 ENDIF
483
484 !-----------------------------------------------------------------------
485 ! Ajout des champs de guidage
486 !-----------------------------------------------------------------------
487 ! Sauvegarde du guidage?
488 ✗ f_out=((MOD(itau,iguide_sav).EQ.0).AND.guide_sav)
489 IF (f_out) THEN
490 ! compute pressures at layer interfaces
491 ✗ CALL pression(ip1jmp1,ap,bp,ps,p)
492 ✗ if (pressure_exner) then
493 ✗ call exner_hyb(ip1jmp1,ps,p,pks,pk)
494 else
495 ✗ call exner_milieu(ip1jmp1,ps,p,pks,pk)
496 endif
497 ✗ unskap=1./kappa
498 ! Now compute pressures at mid-layer
499 ✗ do l=1,llm
500 ✗ p(:,l)=preff*(pk(:,l)/cpp)**unskap
501 enddo
502 ✗ CALL guide_out("SP",jjp1,llm,p(:,1:llm))
503 ENDIF
504
505 ✗ if (guide_u) then
506 ✗ if (guide_add) then
507 ✗ f_add=(1.-tau)*ugui1+tau*ugui2
508 else
509 ✗ f_add=(1.-tau)*ugui1+tau*ugui2-ucov
510 endif
511 ✗ if (guide_zon) CALL guide_zonave(1,jjp1,llm,f_add)
512 ✗ CALL guide_addfield(ip1jmp1,llm,f_add,alpha_u)
513 ✗ IF (f_out) CALL guide_out("ua",jjp1,llm,(1.-tau)*ugui1+tau*ugui2)
514 ✗ IF (f_out) CALL guide_out("u",jjp1,llm,ucov)
515 ✗ IF (f_out) CALL guide_out("ucov",jjp1,llm,f_add/factt)
516 ✗ ucov=ucov+f_add
517 endif
518
519 ✗ if (guide_T) then
520 ✗ if (guide_add) then
521 ✗ f_add=(1.-tau)*tgui1+tau*tgui2
522 else
523 ✗ f_add=(1.-tau)*tgui1+tau*tgui2-teta
524 endif
525 ✗ if (guide_zon) CALL guide_zonave(2,jjp1,llm,f_add)
526 ✗ CALL guide_addfield(ip1jmp1,llm,f_add,alpha_T)
527 ✗ IF (f_out) CALL guide_out("teta",jjp1,llm,f_add/factt)
528 ✗ teta=teta+f_add
529 endif
530
531 ✗ if (guide_P) then
532 ✗ if (guide_add) then
533 ✗ f_add(1:ip1jmp1,1)=(1.-tau)*psgui1+tau*psgui2
534 else
535 ✗ f_add(1:ip1jmp1,1)=(1.-tau)*psgui1+tau*psgui2-ps
536 endif
537 ✗ if (guide_zon) CALL guide_zonave(2,jjp1,1,f_add(1:ip1jmp1,1))
538 ✗ CALL guide_addfield(ip1jmp1,1,f_add(1:ip1jmp1,1),alpha_P)
539 ! IF (f_out) CALL guide_out("ps",jjp1,1,f_add(1:ip1jmp1,1)/factt)
540 ✗ ps=ps+f_add(1:ip1jmp1,1)
541 ✗ CALL pression(ip1jmp1,ap,bp,ps,p)
542 ✗ CALL massdair(p,masse)
543 endif
544
545 ✗ if (guide_Q) then
546 ✗ if (guide_add) then
547 ✗ f_add=(1.-tau)*qgui1+tau*qgui2
548 else
549 ✗ f_add=(1.-tau)*qgui1+tau*qgui2-q
550 endif
551 ✗ if (guide_zon) CALL guide_zonave(2,jjp1,llm,f_add)
552 ✗ CALL guide_addfield(ip1jmp1,llm,f_add,alpha_Q)
553 ✗ IF (f_out) CALL guide_out("q",jjp1,llm,f_add/factt)
554 ✗ q=q+f_add
555 endif
556
557 ✗ if (guide_v) then
558 ✗ if (guide_add) then
559 ✗ f_add(1:ip1jm,:)=(1.-tau)*vgui1+tau*vgui2
560 else
561 ✗ f_add(1:ip1jm,:)=(1.-tau)*vgui1+tau*vgui2-vcov
562 endif
563 ✗ if (guide_zon) CALL guide_zonave(2,jjm,llm,f_add(1:ip1jm,:))
564 ✗ CALL guide_addfield(ip1jm,llm,f_add(1:ip1jm,:),alpha_v)
565 ✗ IF (f_out) CALL guide_out("v",jjm,llm,vcov)
566 ✗ IF (f_out) CALL guide_out("va",jjm,llm,(1.-tau)*vgui1+tau*vgui2)
567 ✗ IF (f_out) CALL guide_out("vcov",jjm,llm,f_add(1:ip1jm,:)/factt)
568 ✗ vcov=vcov+f_add(1:ip1jm,:)
569 endif
570 END SUBROUTINE guide_main
571
572 !=======================================================================
573 ✗ SUBROUTINE guide_addfield(hsize,vsize,field,alpha)
574 ! field1=a*field1+alpha*field2
575
576 IMPLICIT NONE
577
578 ! input variables
579 INTEGER, INTENT(IN) :: hsize
580 INTEGER, INTENT(IN) :: vsize
581 REAL, DIMENSION(hsize), INTENT(IN) :: alpha
582 REAL, DIMENSION(hsize,vsize), INTENT(INOUT) :: field
583
584 ! Local variables
585 INTEGER :: l
586
587 ✗ do l=1,vsize
588 ✗ field(:,l)=alpha*field(:,l)*alpha_pcor(l)
589 enddo
590
591 ✗ END SUBROUTINE guide_addfield
592
593 !=======================================================================
594 ✗ SUBROUTINE guide_zonave(typ,hsize,vsize,field)
595
596 USE comconst_mod, ONLY: pi
597
598 IMPLICIT NONE
599
600 INCLUDE "dimensions.h"
601 INCLUDE "paramet.h"
602 INCLUDE "comgeom.h"
603
604 ! input/output variables
605 INTEGER, INTENT(IN) :: typ
606 INTEGER, INTENT(IN) :: vsize
607 INTEGER, INTENT(IN) :: hsize
608 REAL, DIMENSION(hsize*iip1,vsize), INTENT(INOUT) :: field
609
610 ! Local variables
611 LOGICAL, SAVE :: first=.TRUE.
612 INTEGER, DIMENSION (2), SAVE :: imin, imax ! averaging domain
613 INTEGER :: i,j,l,ij
614 REAL, DIMENSION (iip1) :: lond ! longitude in Deg.
615 ✗ REAL, DIMENSION (hsize,vsize):: fieldm ! zon-averaged field
616
617 ✗ IF (first) THEN
618 ✗ first=.FALSE.
619 !Compute domain for averaging
620 ✗ lond=rlonu*180./pi
621 ✗ imin(1)=1;imax(1)=iip1;
622 ✗ imin(2)=1;imax(2)=iip1;
623 ✗ IF (guide_reg) THEN
624 ✗ DO i=1,iim
625 ✗ IF (lond(i).LT.lon_min_g) imin(1)=i
626 ✗ IF (lond(i).LE.lon_max_g) imax(1)=i
627 ENDDO
628 ✗ lond=rlonv*180./pi
629 ✗ DO i=1,iim
630 ✗ IF (lond(i).LT.lon_min_g) imin(2)=i
631 ✗ IF (lond(i).LE.lon_max_g) imax(2)=i
632 ENDDO
633 ENDIF
634 ENDIF
635
636 ✗ fieldm=0.
637 ✗ DO l=1,vsize
638 ! Compute zonal average
639 ✗ DO j=1,hsize
640 ✗ DO i=imin(typ),imax(typ)
641 ✗ ij=(j-1)*iip1+i
642 ✗ fieldm(j,l)=fieldm(j,l)+field(ij,l)
643 ENDDO
644 ENDDO
645 ✗ fieldm(:,l)=fieldm(:,l)/REAL(imax(typ)-imin(typ)+1)
646 ! Compute forcing
647 ✗ DO j=1,hsize
648 ✗ DO i=1,iip1
649 ✗ ij=(j-1)*iip1+i
650 ✗ field(ij,l)=fieldm(j,l)
651 ENDDO
652 ENDDO
653 ENDDO
654
655 ✗ END SUBROUTINE guide_zonave
656
657 !=======================================================================
658 ✗ SUBROUTINE guide_interp(psi,teta)
659
660 use exner_hyb_m, only: exner_hyb
661 use exner_milieu_m, only: exner_milieu
662 use comconst_mod, only: kappa, cpp
663 use comvert_mod, only: preff, pressure_exner, bp, ap
664 IMPLICIT NONE
665
666 include "dimensions.h"
667 include "paramet.h"
668 include "comgeom2.h"
669
670 REAL, DIMENSION (iip1,jjp1), INTENT(IN) :: psi ! Psol gcm
671 REAL, DIMENSION (iip1,jjp1,llm), INTENT(IN) :: teta ! Temp. Pot. gcm
672
673 LOGICAL, SAVE :: first=.TRUE.
674 ! Variables pour niveaux pression:
675 ✗ REAL, DIMENSION (iip1,jjp1,nlevnc) :: plnc1,plnc2 !niveaux pression guidage
676 REAL, DIMENSION (iip1,jjp1,llm) :: plunc,plsnc !niveaux pression modele
677 REAL, DIMENSION (iip1,jjm,llm) :: plvnc !niveaux pression modele
678 REAL, DIMENSION (iip1,jjp1,llmp1) :: p ! pression intercouches
679 REAL, DIMENSION (iip1,jjp1,llm) :: pls, pext ! var intermediaire
680 REAL, DIMENSION (iip1,jjp1,llm) :: pbarx
681 REAL, DIMENSION (iip1,jjm,llm) :: pbary
682 ! Variables pour fonction Exner (P milieu couche)
683 REAL, DIMENSION (iip1,jjp1,llm) :: pk
684 REAL, DIMENSION (iip1,jjp1) :: pks
685 REAL :: prefkap,unskap
686 ! Pression de vapeur saturante
687 REAL, DIMENSION (ip1jmp1,llm) :: qsat
688 !Variables intermediaires interpolation
689 REAL, DIMENSION (iip1,jjp1,llm) :: zu1,zu2
690 REAL, DIMENSION (iip1,jjm,llm) :: zv1,zv2
691
692 INTEGER :: i,j,l,ij
693 CHARACTER(LEN=20),PARAMETER :: modname="guide_interp"
694
695 ✗ write(*,*)trim(modname)//': interpolate nudging variables'
696 ! -----------------------------------------------------------------
697 ! Calcul des niveaux de pression champs guidage
698 ! -----------------------------------------------------------------
699 ✗ if (guide_modele) then
700 ✗ do i=1,iip1
701 ✗ do j=1,jjp1
702 ✗ do l=1,nlevnc
703 ✗ plnc2(i,j,l)=apnc(l)+bpnc(l)*psnat2(i,j)
704 ✗ plnc1(i,j,l)=apnc(l)+bpnc(l)*psnat1(i,j)
705 enddo
706 enddo
707 enddo
708 else
709 ✗ do i=1,iip1
710 ✗ do j=1,jjp1
711 ✗ do l=1,nlevnc
712 ✗ plnc2(i,j,l)=apnc(l)
713 ✗ plnc1(i,j,l)=apnc(l)
714 enddo
715 enddo
716 enddo
717
718 endif
719 ✗ if (first) then
720 ✗ first=.FALSE.
721 ✗ write(*,*)trim(modname)//' : check vertical level order'
722 ✗ write(*,*)trim(modname)//' LMDZ :'
723 ✗ do l=1,llm
724 ✗ write(*,*)trim(modname)//' PL(',l,')=',(ap(l)+ap(l+1))/2. &
725 ✗ +psi(1,jjp1)*(bp(l)+bp(l+1))/2.
726 enddo
727 ✗ write(*,*)trim(modname)//' nudging file :'
728 ✗ do l=1,nlevnc
729 ✗ write(*,*)trim(modname)//' PL(',l,')=',plnc2(1,1,l)
730 enddo
731 ✗ write(*,*)trim(modname)//' invert ordering: invert_p=',invert_p
732 ✗ if (guide_u) then
733 ✗ do l=1,nlevnc
734 ✗ write(*,*)trim(modname)//' U(',l,')=',unat2(1,1,l)
735 enddo
736 endif
737 ✗ if (guide_T) then
738 ✗ do l=1,nlevnc
739 ✗ write(*,*)trim(modname)//' T(',l,')=',tnat2(1,1,l)
740 enddo
741 endif
742 endif
743
744 ! -----------------------------------------------------------------
745 ! Calcul niveaux pression modele
746 ! -----------------------------------------------------------------
747 ✗ CALL pression( ip1jmp1, ap, bp, psi, p )
748 ✗ if (pressure_exner) then
749 ✗ CALL exner_hyb(ip1jmp1,psi,p,pks,pk)
750 else
751 ✗ CALL exner_milieu(ip1jmp1,psi,p,pks,pk)
752 endif
753 ! .... Calcul de pls , pression au milieu des couches ,en Pascals
754 ✗ unskap=1./kappa
755 ✗ prefkap = preff ** kappa
756 ✗ DO l = 1, llm
757 ✗ DO j=1,jjp1
758 ✗ DO i =1, iip1
759 ✗ pls(i,j,l) = preff * ( pk(i,j,l)/cpp) ** unskap
760 ENDDO
761 ENDDO
762 ENDDO
763
764 ! calcul des pressions pour les grilles u et v
765 ✗ do l=1,llm
766 ✗ do j=1,jjp1
767 ✗ do i=1,iip1
768 ✗ pext(i,j,l)=pls(i,j,l)*aire(i,j)
769 enddo
770 enddo
771 enddo
772 ✗ call massbar(pext, pbarx, pbary )
773 ✗ do l=1,llm
774 ✗ do j=1,jjp1
775 ✗ do i=1,iip1
776 ✗ plunc(i,j,l)=pbarx(i,j,l)/aireu(i,j)
777 ✗ plsnc(i,j,l)=pls(i,j,l)
778 enddo
779 enddo
780 enddo
781 ✗ do l=1,llm
782 ✗ do j=1,jjm
783 ✗ do i=1,iip1
784 ✗ plvnc(i,j,l)=pbary(i,j,l)/airev(i,j)
785 enddo
786 enddo
787 enddo
788
789 ! -----------------------------------------------------------------
790 ! Interpolation champs guidage sur niveaux modele (+inversion N/S)
791 ! Conversion en variables gcm (ucov, vcov...)
792 ! -----------------------------------------------------------------
793 ✗ if (guide_P) then
794 ✗ do j=1,jjp1
795 ✗ do i=1,iim
796 ✗ ij=(j-1)*iip1+i
797 ✗ psgui1(ij)=psnat1(i,j)
798 ✗ psgui2(ij)=psnat2(i,j)
799 enddo
800 ✗ psgui1(iip1*j)=psnat1(1,j)
801 ✗ psgui2(iip1*j)=psnat2(1,j)
802 enddo
803 endif
804
805 ✗ IF (guide_u) THEN
806 ✗ CALL pres2lev(unat1,zu1,nlevnc,llm,plnc1,plunc,iip1,jjp1,invert_p)
807 ✗ CALL pres2lev(unat2,zu2,nlevnc,llm,plnc2,plunc,iip1,jjp1,invert_p)
808 ✗ do l=1,llm
809 ✗ do j=1,jjp1
810 ✗ do i=1,iim
811 ✗ ij=(j-1)*iip1+i
812 ✗ ugui1(ij,l)=zu1(i,j,l)*cu(i,j)
813 ✗ ugui2(ij,l)=zu2(i,j,l)*cu(i,j)
814 enddo
815 ✗ ugui1(j*iip1,l)=ugui1((j-1)*iip1+1,l)
816 ✗ ugui2(j*iip1,l)=ugui2((j-1)*iip1+1,l)
817 enddo
818 ✗ do i=1,iip1
819 ✗ ugui1(i,l)=0.
820 ✗ ugui1(ip1jm+i,l)=0.
821 ✗ ugui2(i,l)=0.
822 ✗ ugui2(ip1jm+i,l)=0.
823 enddo
824 enddo
825 ENDIF
826
827 ✗ IF (guide_T) THEN
828 ✗ CALL pres2lev(tnat1,zu1,nlevnc,llm,plnc1,plsnc,iip1,jjp1,invert_p)
829 ✗ CALL pres2lev(tnat2,zu2,nlevnc,llm,plnc2,plsnc,iip1,jjp1,invert_p)
830 ✗ do l=1,llm
831 ✗ do j=1,jjp1
832 ✗ IF (guide_teta) THEN
833 ✗ do i=1,iim
834 ✗ ij=(j-1)*iip1+i
835 ✗ tgui1(ij,l)=zu1(i,j,l)
836 ✗ tgui2(ij,l)=zu2(i,j,l)
837 enddo
838 ELSE
839 ✗ do i=1,iim
840 ✗ ij=(j-1)*iip1+i
841 ✗ tgui1(ij,l)=zu1(i,j,l)*cpp/pk(i,j,l)
842 ✗ tgui2(ij,l)=zu2(i,j,l)*cpp/pk(i,j,l)
843 enddo
844 ENDIF
845 ✗ tgui1(j*iip1,l)=tgui1((j-1)*iip1+1,l)
846 ✗ tgui2(j*iip1,l)=tgui2((j-1)*iip1+1,l)
847 enddo
848 ✗ do i=1,iip1
849 ✗ tgui1(i,l)=tgui1(1,l)
850 ✗ tgui1(ip1jm+i,l)=tgui1(ip1jm+1,l)
851 ✗ tgui2(i,l)=tgui2(1,l)
852 ✗ tgui2(ip1jm+i,l)=tgui2(ip1jm+1,l)
853 enddo
854 enddo
855 ENDIF
856
857 ✗ IF (guide_v) THEN
858
859 ✗ CALL pres2lev(vnat1,zv1,nlevnc,llm,plnc1(:,:jjm,:),plvnc,iip1,jjm,invert_p)
860 ✗ CALL pres2lev(vnat2,zv2,nlevnc,llm,plnc2(:,:jjm,:),plvnc,iip1,jjm,invert_p)
861
862 ✗ do l=1,llm
863 ✗ do j=1,jjm
864 ✗ do i=1,iim
865 ✗ ij=(j-1)*iip1+i
866 ✗ vgui1(ij,l)=zv1(i,j,l)*cv(i,j)
867 ✗ vgui2(ij,l)=zv2(i,j,l)*cv(i,j)
868 enddo
869 ✗ vgui1(j*iip1,l)=vgui1((j-1)*iip1+1,l)
870 ✗ vgui2(j*iip1,l)=vgui2((j-1)*iip1+1,l)
871 enddo
872 enddo
873 ENDIF
874
875 ✗ IF (guide_Q) THEN
876 ! On suppose qu'on a la bonne variable dans le fichier de guidage:
877 ! Hum.Rel si guide_hr, Hum.Spec. sinon.
878 ✗ CALL pres2lev(qnat1,zu1,nlevnc,llm,plnc1,plsnc,iip1,jjp1,invert_p)
879 ✗ CALL pres2lev(qnat2,zu2,nlevnc,llm,plnc2,plsnc,iip1,jjp1,invert_p)
880 ✗ do l=1,llm
881 ✗ do j=1,jjp1
882 ✗ do i=1,iim
883 ✗ ij=(j-1)*iip1+i
884 ✗ qgui1(ij,l)=zu1(i,j,l)
885 ✗ qgui2(ij,l)=zu2(i,j,l)
886 enddo
887 ✗ qgui1(j*iip1,l)=qgui1((j-1)*iip1+1,l)
888 ✗ qgui2(j*iip1,l)=qgui2((j-1)*iip1+1,l)
889 enddo
890 ✗ do i=1,iip1
891 ✗ qgui1(i,l)=qgui1(1,l)
892 ✗ qgui1(ip1jm+i,l)=qgui1(ip1jm+1,l)
893 ✗ qgui2(i,l)=qgui2(1,l)
894 ✗ qgui2(ip1jm+i,l)=qgui2(ip1jm+1,l)
895 enddo
896 enddo
897 ✗ IF (guide_hr) THEN
898 ✗ CALL q_sat(iip1*jjp1*llm,teta*pk/cpp,plsnc,qsat)
899 ✗ qgui1=qgui1*qsat*0.01 !hum. rel. en %
900 ✗ qgui2=qgui2*qsat*0.01
901 ENDIF
902 ENDIF
903
904 ✗ END SUBROUTINE guide_interp
905
906 !=======================================================================
907 ✗ SUBROUTINE tau2alpha(typ,pim,pjm,factt,taumin,taumax,alpha)
908
909 ! Calcul des constantes de rappel alpha (=1/tau)
910
911 use comconst_mod, only: pi
912 use serre_mod, only: clon, clat, grossismx, grossismy
913
914 implicit none
915
916 include "dimensions.h"
917 include "paramet.h"
918 include "comgeom2.h"
919
920 ! input arguments :
921 INTEGER, INTENT(IN) :: typ ! u(2),v(3), ou scalaire(1)
922 INTEGER, INTENT(IN) :: pim,pjm ! dimensions en lat, lon
923 REAL, INTENT(IN) :: factt ! pas de temps en fraction de jour
924 REAL, INTENT(IN) :: taumin,taumax
925 ! output arguments:
926 REAL, DIMENSION(pim,pjm), INTENT(OUT) :: alpha
927
928 ! local variables:
929 LOGICAL, SAVE :: first=.TRUE.
930 REAL, SAVE :: gamma,dxdy_min,dxdy_max
931 REAL, DIMENSION (iip1,jjp1) :: zdx,zdy
932 REAL, DIMENSION (iip1,jjp1) :: dxdys,dxdyu
933 REAL, DIMENSION (iip1,jjm) :: dxdyv
934 real dxdy_
935 real zlat,zlon
936 real alphamin,alphamax,xi
937 integer i,j,ilon,ilat
938 character(len=20),parameter :: modname="tau2alpha"
939
940
941 ✗ alphamin=factt/taumax
942 ✗ alphamax=factt/taumin
943 ✗ IF (guide_reg.OR.guide_add) THEN
944 ✗ alpha=alphamax
945 !-----------------------------------------------------------------------
946 ! guide_reg: alpha=alpha_min dans region, 0. sinon.
947 !-----------------------------------------------------------------------
948 ✗ IF (guide_reg) THEN
949 ✗ do j=1,pjm
950 ✗ do i=1,pim
951 ✗ if (typ.eq.2) then
952 ✗ zlat=rlatu(j)*180./pi
953 ✗ zlon=rlonu(i)*180./pi
954 ✗ elseif (typ.eq.1) then
955 ✗ zlat=rlatu(j)*180./pi
956 ✗ zlon=rlonv(i)*180./pi
957 ✗ elseif (typ.eq.3) then
958 ✗ zlat=rlatv(j)*180./pi
959 ✗ zlon=rlonv(i)*180./pi
960 endif
961 alpha(i,j)=alphamax/16.* &
962 (1.+tanh((zlat-lat_min_g)/tau_lat))* &
963 (1.+tanh((lat_max_g-zlat)/tau_lat))* &
964 (1.+tanh((zlon-lon_min_g)/tau_lon))* &
965 ✗ (1.+tanh((lon_max_g-zlon)/tau_lon))
966 enddo
967 enddo
968 ENDIF
969 ELSE
970 !-----------------------------------------------------------------------
971 ! Sinon, alpha varie entre alpha_min et alpha_max suivant le zoom.
972 !-----------------------------------------------------------------------
973 !Calcul de l'aire des mailles
974 ✗ do j=2,jjm
975 ✗ do i=2,iip1
976 ✗ zdx(i,j)=0.5*(cu(i-1,j)+cu(i,j))/cos(rlatu(j))
977 enddo
978 ✗ zdx(1,j)=zdx(iip1,j)
979 enddo
980 ✗ do j=2,jjm
981 ✗ do i=1,iip1
982 ✗ zdy(i,j)=0.5*(cv(i,j-1)+cv(i,j))
983 enddo
984 enddo
985 ✗ do i=1,iip1
986 ✗ zdx(i,1)=zdx(i,2)
987 ✗ zdx(i,jjp1)=zdx(i,jjm)
988 ✗ zdy(i,1)=zdy(i,2)
989 ✗ zdy(i,jjp1)=zdy(i,jjm)
990 enddo
991 ✗ do j=1,jjp1
992 ✗ do i=1,iip1
993 ✗ dxdys(i,j)=sqrt(zdx(i,j)*zdx(i,j)+zdy(i,j)*zdy(i,j))
994 enddo
995 enddo
996 ✗ IF (typ.EQ.2) THEN
997 ✗ do j=1,jjp1
998 ✗ do i=1,iim
999 ✗ dxdyu(i,j)=0.5*(dxdys(i,j)+dxdys(i+1,j))
1000 enddo
1001 ✗ dxdyu(iip1,j)=dxdyu(1,j)
1002 enddo
1003 ENDIF
1004 ✗ IF (typ.EQ.3) THEN
1005 ✗ do j=1,jjm
1006 ✗ do i=1,iip1
1007 ✗ dxdyv(i,j)=0.5*(dxdys(i,j)+dxdys(i,j+1))
1008 enddo
1009 enddo
1010 ENDIF
1011 ! Premier appel: calcul des aires min et max et de gamma.
1012 ✗ IF (first) THEN
1013 ✗ first=.FALSE.
1014 ! coordonnees du centre du zoom
1015 ✗ CALL coordij(clon,clat,ilon,ilat)
1016 ! aire de la maille au centre du zoom
1017 ✗ dxdy_min=dxdys(ilon,ilat)
1018 ! dxdy maximale de la maille
1019 ✗ dxdy_max=0.
1020 ✗ do j=1,jjp1
1021 ✗ do i=1,iip1
1022 ✗ dxdy_max=max(dxdy_max,dxdys(i,j))
1023 enddo
1024 enddo
1025 ! Calcul de gamma
1026 ✗ if (abs(grossismx-1.).lt.0.1.or.abs(grossismy-1.).lt.0.1) then
1027 ✗ write(*,*)trim(modname)//' ATTENTION modele peu zoome'
1028 ✗ write(*,*)trim(modname)//' ATTENTION on prend une constante de guidage cste'
1029 ✗ gamma=0.
1030 else
1031 ✗ gamma=(dxdy_max-2.*dxdy_min)/(dxdy_max-dxdy_min)
1032 ✗ write(*,*)trim(modname)//' gamma=',gamma
1033 ✗ if (gamma.lt.1.e-5) then
1034 ✗ write(*,*)trim(modname)//' gamma =',gamma,'<1e-5'
1035 ✗ stop
1036 endif
1037 ✗ gamma=log(0.5)/log(gamma)
1038 ✗ if (gamma4) then
1039 ✗ gamma=min(gamma,4.)
1040 endif
1041 ✗ write(*,*)trim(modname)//' gamma=',gamma
1042 endif
1043 ENDIF !first
1044
1045 ✗ do j=1,pjm
1046 ✗ do i=1,pim
1047 ✗ if (typ.eq.1) then
1048 ✗ dxdy_=dxdys(i,j)
1049 ✗ zlat=rlatu(j)*180./pi
1050 ✗ elseif (typ.eq.2) then
1051 ✗ dxdy_=dxdyu(i,j)
1052 ✗ zlat=rlatu(j)*180./pi
1053 ✗ elseif (typ.eq.3) then
1054 ✗ dxdy_=dxdyv(i,j)
1055 ✗ zlat=rlatv(j)*180./pi
1056 endif
1057 ✗ if (abs(grossismx-1.).lt.0.1.or.abs(grossismy-1.).lt.0.1) then
1058 ! pour une grille reguliere, xi=xxx**0=1 -> alpha=alphamin
1059 ✗ alpha(i,j)=alphamin
1060 else
1061 ✗ xi=((dxdy_max-dxdy_)/(dxdy_max-dxdy_min))**gamma
1062 ✗ xi=min(xi,1.)
1063 ✗ if(lat_min_g.le.zlat .and. zlat.le.lat_max_g) then
1064 ✗ alpha(i,j)=xi*alphamin+(1.-xi)*alphamax
1065 else
1066 ✗ alpha(i,j)=0.
1067 endif
1068 endif
1069 enddo
1070 enddo
1071 ENDIF ! guide_reg
1072
1073 ✗ if (.not. guide_add) alpha = 1. - exp(- alpha)
1074
1075 ✗ END SUBROUTINE tau2alpha
1076
1077 !=======================================================================
1078 ✗ SUBROUTINE guide_read(timestep)
1079
1080 IMPLICIT NONE
1081
1082 include "netcdf.inc"
1083 include "dimensions.h"
1084 include "paramet.h"
1085
1086 INTEGER, INTENT(IN) :: timestep
1087
1088 LOGICAL, SAVE :: first=.TRUE.
1089 ! Identification fichiers et variables NetCDF:
1090 INTEGER, SAVE :: ncidu,varidu,ncidv,varidv,ncidp,varidp
1091 INTEGER, SAVE :: ncidQ,varidQ,ncidt,varidt,ncidps,varidps
1092 INTEGER :: ncidpl,varidpl,varidap,varidbp,dimid,lendim
1093 ! Variables auxiliaires NetCDF:
1094 INTEGER, DIMENSION(4) :: start,count
1095 INTEGER :: status,rcode
1096 CHARACTER (len = 80) :: abort_message
1097 CHARACTER (len = 20) :: modname = 'guide_read'
1098 CHARACTER (len = 20) :: namedim
1099
1100 ! -----------------------------------------------------------------
1101 ! Premier appel: initialisation de la lecture des fichiers
1102 ! -----------------------------------------------------------------
1103 ✗ if (first) then
1104 ✗ ncidpl=-99
1105 ✗ write(*,*),trim(modname)//': opening nudging files '
1106 ! Niveaux de pression si non constants
1107 ✗ if (guide_plevs.EQ.1) then
1108 ✗ write(*,*),trim(modname)//' Reading nudging on model levels'
1109 ✗ rcode = nf90_open('apbp.nc', nf90_nowrite, ncidpl)
1110 ✗ IF (rcode.NE.NF_NOERR) THEN
1111 ✗ abort_message='Nudging: error -> no file apbp.nc'
1112 ✗ CALL abort_gcm(modname,abort_message,1)
1113 ENDIF
1114 ✗ rcode = nf90_inq_varid(ncidpl, 'AP', varidap)
1115 ✗ IF (rcode.NE.NF_NOERR) THEN
1116 ✗ abort_message='Nudging: error -> no AP variable in file apbp.nc'
1117 ✗ CALL abort_gcm(modname,abort_message,1)
1118 ENDIF
1119 ✗ rcode = nf90_inq_varid(ncidpl, 'BP', varidbp)
1120 ✗ IF (rcode.NE.NF_NOERR) THEN
1121 ✗ abort_message='Nudging: error -> no BP variable in file apbp.nc'
1122 ✗ CALL abort_gcm(modname,abort_message,1)
1123 ENDIF
1124 ✗ write(*,*),trim(modname)//' ncidpl,varidap',ncidpl,varidap
1125 endif
1126
1127 ! Pression si guidage sur niveaux P variables
1128 ✗ if (guide_plevs.EQ.2) then
1129 ✗ rcode = nf90_open('P.nc', nf90_nowrite, ncidp)
1130 ✗ IF (rcode.NE.NF_NOERR) THEN
1131 ✗ abort_message='Nudging: error -> no file P.nc'
1132 ✗ CALL abort_gcm(modname,abort_message,1)
1133 ENDIF
1134 ✗ rcode = nf90_inq_varid(ncidp, 'PRES', varidp)
1135 ✗ IF (rcode.NE.NF_NOERR) THEN
1136 ✗ abort_message='Nudging: error -> no PRES variable in file P.nc'
1137 ✗ CALL abort_gcm(modname,abort_message,1)
1138 ENDIF
1139 ✗ write(*,*),trim(modname)//' ncidp,varidp',ncidp,varidp
1140 ✗ if (ncidpl.eq.-99) ncidpl=ncidp
1141 endif
1142
1143 ! Vent zonal
1144 ✗ if (guide_u) then
1145 ✗ rcode = nf90_open('u.nc', nf90_nowrite, ncidu)
1146 ✗ IF (rcode.NE.NF_NOERR) THEN
1147 ✗ abort_message='Nudging: error -> no file u.nc'
1148 ✗ CALL abort_gcm(modname,abort_message,1)
1149 ENDIF
1150 ✗ rcode = nf90_inq_varid(ncidu, 'UWND', varidu)
1151 ✗ IF (rcode.NE.NF_NOERR) THEN
1152 ✗ abort_message='Nudging: error -> no UWND variable in file u.nc'
1153 ✗ CALL abort_gcm(modname,abort_message,1)
1154 ENDIF
1155 ✗ write(*,*),trim(modname)//' ncidu,varidu',ncidu,varidu
1156 ✗ if (ncidpl.eq.-99) ncidpl=ncidu
1157
1158 ✗ status=NF90_INQ_DIMID(ncidu, "LONU", dimid)
1159 ✗ status=NF90_INQUIRE_DIMENSION(ncidu,dimid,namedim,lendim)
1160 ✗ IF (lendim .NE. iip1) THEN
1161 ✗ abort_message='dimension LONU different from iip1 in u.nc'
1162 ✗ CALL abort_gcm(modname,abort_message,1)
1163 ENDIF
1164
1165 ✗ status=NF90_INQ_DIMID(ncidu, "LATU", dimid)
1166 ✗ status=NF90_INQUIRE_DIMENSION(ncidu,dimid,namedim,lendim)
1167 ✗ IF (lendim .NE. jjp1) THEN
1168 ✗ abort_message='dimension LATU different from jjp1 in u.nc'
1169 ✗ CALL abort_gcm(modname,abort_message,1)
1170 ENDIF
1171
1172 endif
1173
1174 ! Vent meridien
1175 ✗ if (guide_v) then
1176 ✗ rcode = nf90_open('v.nc', nf90_nowrite, ncidv)
1177 ✗ IF (rcode.NE.NF_NOERR) THEN
1178 ✗ abort_message='Nudging: error -> no file v.nc'
1179 ✗ CALL abort_gcm(modname,abort_message,1)
1180 ENDIF
1181 ✗ rcode = nf90_inq_varid(ncidv, 'VWND', varidv)
1182 ✗ IF (rcode.NE.NF_NOERR) THEN
1183 ✗ abort_message='Nudging: error -> no VWND variable in file v.nc'
1184 ✗ CALL abort_gcm(modname,abort_message,1)
1185 ENDIF
1186 ✗ write(*,*),trim(modname)//' ncidv,varidv',ncidv,varidv
1187 ✗ if (ncidpl.eq.-99) ncidpl=ncidv
1188
1189 ✗ status=NF90_INQ_DIMID(ncidv, "LONV", dimid)
1190 ✗ status=NF90_INQUIRE_DIMENSION(ncidv,dimid,namedim,lendim)
1191
1192 ✗ IF (lendim .NE. iip1) THEN
1193 ✗ abort_message='dimension LONV different from iip1 in v.nc'
1194 ✗ CALL abort_gcm(modname,abort_message,1)
1195 ENDIF
1196
1197
1198 ✗ status=NF90_INQ_DIMID(ncidv, "LATV", dimid)
1199 ✗ status=NF90_INQUIRE_DIMENSION(ncidv,dimid,namedim,lendim)
1200 ✗ IF (lendim .NE. jjm) THEN
1201 ✗ abort_message='dimension LATV different from jjm in v.nc'
1202 ✗ CALL abort_gcm(modname,abort_message,1)
1203 ENDIF
1204
1205 endif
1206
1207 ! Temperature
1208 ✗ if (guide_T) then
1209 ✗ rcode = nf90_open('T.nc', nf90_nowrite, ncidt)
1210 ✗ IF (rcode.NE.NF_NOERR) THEN
1211 ✗ abort_message='Nudging: error -> no file T.nc'
1212 ✗ CALL abort_gcm(modname,abort_message,1)
1213 ENDIF
1214 ✗ rcode = nf90_inq_varid(ncidt, 'AIR', varidt)
1215 ✗ IF (rcode.NE.NF_NOERR) THEN
1216 ✗ abort_message='Nudging: error -> no AIR variable in file T.nc'
1217 ✗ CALL abort_gcm(modname,abort_message,1)
1218 ENDIF
1219 ✗ write(*,*),trim(modname)//' ncidT,varidT',ncidt,varidt
1220 ✗ if (ncidpl.eq.-99) ncidpl=ncidt
1221
1222 ✗ status=NF90_INQ_DIMID(ncidt, "LONV", dimid)
1223 ✗ status=NF90_INQUIRE_DIMENSION(ncidt,dimid,namedim,lendim)
1224 ✗ IF (lendim .NE. iip1) THEN
1225 ✗ abort_message='dimension LONV different from iip1 in T.nc'
1226 ✗ CALL abort_gcm(modname,abort_message,1)
1227 ENDIF
1228
1229 ✗ status=NF90_INQ_DIMID(ncidt, "LATU", dimid)
1230 ✗ status=NF90_INQUIRE_DIMENSION(ncidt,dimid,namedim,lendim)
1231 ✗ IF (lendim .NE. jjp1) THEN
1232 ✗ abort_message='dimension LATU different from jjp1 in T.nc'
1233 ✗ CALL abort_gcm(modname,abort_message,1)
1234 ENDIF
1235
1236 endif
1237
1238 ! Humidite
1239 ✗ if (guide_Q) then
1240 ✗ rcode = nf90_open('hur.nc', nf90_nowrite, ncidQ)
1241 ✗ IF (rcode.NE.NF_NOERR) THEN
1242 ✗ abort_message='Nudging: error -> no file hur.nc'
1243 ✗ CALL abort_gcm(modname,abort_message,1)
1244 ENDIF
1245 ✗ rcode = nf90_inq_varid(ncidQ, 'RH', varidQ)
1246 ✗ IF (rcode.NE.NF_NOERR) THEN
1247 ✗ abort_message='Nudging: error -> no RH variable in file hur.nc'
1248 ✗ CALL abort_gcm(modname,abort_message,1)
1249 ENDIF
1250 ✗ write(*,*),trim(modname)//' ncidQ,varidQ',ncidQ,varidQ
1251 ✗ if (ncidpl.eq.-99) ncidpl=ncidQ
1252
1253 ✗ status=NF90_INQ_DIMID(ncidQ, "LONV", dimid)
1254 ✗ status=NF90_INQUIRE_DIMENSION(ncidQ,dimid,namedim,lendim)
1255 ✗ IF (lendim .NE. iip1) THEN
1256 ✗ abort_message='dimension LONV different from iip1 in hur.nc'
1257 ✗ CALL abort_gcm(modname,abort_message,1)
1258 ENDIF
1259
1260 ✗ status=NF90_INQ_DIMID(ncidQ, "LATU", dimid)
1261 ✗ status=NF90_INQUIRE_DIMENSION(ncidQ,dimid,namedim,lendim)
1262 ✗ IF (lendim .NE. jjp1) THEN
1263 ✗ abort_message='dimension LATU different from jjp1 in hur.nc'
1264 ✗ CALL abort_gcm(modname,abort_message,1)
1265 ENDIF
1266
1267 endif
1268
1269 ! Pression de surface
1270 ✗ if ((guide_P).OR.(guide_modele)) then
1271 ✗ rcode = nf90_open('ps.nc', nf90_nowrite, ncidps)
1272 ✗ IF (rcode.NE.NF_NOERR) THEN
1273 ✗ abort_message='Nudging: error -> no file ps.nc'
1274 ✗ CALL abort_gcm(modname,abort_message,1)
1275 ENDIF
1276 ✗ rcode = nf90_inq_varid(ncidps, 'SP', varidps)
1277 ✗ IF (rcode.NE.NF_NOERR) THEN
1278 ✗ abort_message='Nudging: error -> no SP variable in file ps.nc'
1279 ✗ CALL abort_gcm(modname,abort_message,1)
1280 ENDIF
1281 ✗ write(*,*),trim(modname)//' ncidps,varidps',ncidps,varidps
1282 endif
1283 ! Coordonnee verticale
1284 ✗ if (guide_plevs.EQ.0) then
1285 ✗ rcode = nf90_inq_varid(ncidpl, 'LEVEL', varidpl)
1286 ✗ IF (rcode.NE.0) rcode = nf90_inq_varid(ncidpl, 'PRESSURE', varidpl)
1287 ✗ write(*,*),trim(modname)//' ncidpl,varidpl',ncidpl,varidpl
1288 endif
1289 ! Coefs ap, bp pour calcul de la pression aux differents niveaux
1290 ✗ if (guide_plevs.EQ.1) then
1291 ✗ status=NF_GET_VARA_DOUBLE(ncidpl,varidap,1,nlevnc,apnc)
1292 ✗ status=NF_GET_VARA_DOUBLE(ncidpl,varidbp,1,nlevnc,bpnc)
1293 ✗ ELSEIF (guide_plevs.EQ.0) THEN
1294 ✗ status=NF_GET_VARA_DOUBLE(ncidpl,varidpl,1,nlevnc,apnc)
1295 !FC Pour les corrections la pression est deja en Pascals on commente la ligne ci-dessous
1296 ✗ IF(convert_Pa) apnc=apnc*100.! conversion en Pascals
1297 ✗ bpnc(:)=0.
1298 endif
1299 ✗ first=.FALSE.
1300 endif ! (first)
1301
1302 ! -----------------------------------------------------------------
1303 ! lecture des champs u, v, T, Q, ps
1304 ! -----------------------------------------------------------------
1305
1306 ! dimensions pour les champs scalaires et le vent zonal
1307 ✗ start(1)=1
1308 ✗ start(2)=1
1309 ✗ start(3)=1
1310 ✗ start(4)=timestep
1311
1312 ✗ count(1)=iip1
1313 ✗ count(2)=jjp1
1314 ✗ count(3)=nlevnc
1315 ✗ count(4)=1
1316
1317 ! Pression
1318 ✗ if (guide_plevs.EQ.2) then
1319 ✗ status=NF_GET_VARA_DOUBLE(ncidp,varidp,start,count,pnat2)
1320 ✗ IF (invert_y) THEN
1321 ! PRINT*,"Invertion impossible actuellement"
1322 ! CALL abort_gcm(modname,abort_message,1)
1323 ✗ CALL invert_lat(iip1,jjp1,nlevnc,pnat2)
1324 ENDIF
1325 endif
1326
1327 ! Vent zonal
1328 ✗ if (guide_u) then
1329 ✗ status=NF_GET_VARA_DOUBLE(ncidu,varidu,start,count,unat2)
1330 ✗ IF (invert_y) THEN
1331 ✗ CALL invert_lat(iip1,jjp1,nlevnc,unat2)
1332 ENDIF
1333 endif
1334
1335 ! Temperature
1336 ✗ if (guide_T) then
1337 ✗ status=NF_GET_VARA_DOUBLE(ncidt,varidt,start,count,tnat2)
1338 ✗ IF (invert_y) THEN
1339 ✗ CALL invert_lat(iip1,jjp1,nlevnc,tnat2)
1340 ENDIF
1341 endif
1342
1343 ! Humidite
1344 ✗ if (guide_Q) then
1345 ✗ status=NF_GET_VARA_DOUBLE(ncidQ,varidQ,start,count,qnat2)
1346 ✗ IF (invert_y) THEN
1347 ✗ CALL invert_lat(iip1,jjp1,nlevnc,qnat2)
1348 ENDIF
1349
1350 endif
1351
1352 ! Vent meridien
1353 ✗ if (guide_v) then
1354 ✗ count(2)=jjm
1355 ✗ status=NF_GET_VARA_DOUBLE(ncidv,varidv,start,count,vnat2)
1356 ✗ IF (invert_y) THEN
1357 ✗ CALL invert_lat(iip1,jjm,nlevnc,vnat2)
1358 ENDIF
1359 endif
1360
1361 ! Pression de surface
1362 ✗ if ((guide_P).OR.(guide_modele)) then
1363 ✗ start(3)=timestep
1364 ✗ start(4)=0
1365 ✗ count(2)=jjp1
1366 ✗ count(3)=1
1367 ✗ count(4)=0
1368 ✗ status=NF_GET_VARA_DOUBLE(ncidps,varidps,start,count,psnat2)
1369 ✗ IF (invert_y) THEN
1370 ✗ CALL invert_lat(iip1,jjp1,1,psnat2)
1371 ENDIF
1372 endif
1373
1374 ✗ END SUBROUTINE guide_read
1375
1376 !=======================================================================
1377 ✗ SUBROUTINE guide_read2D(timestep)
1378
1379 IMPLICIT NONE
1380
1381 include "netcdf.inc"
1382 include "dimensions.h"
1383 include "paramet.h"
1384
1385 INTEGER, INTENT(IN) :: timestep
1386
1387 LOGICAL, SAVE :: first=.TRUE.
1388 ! Identification fichiers et variables NetCDF:
1389 INTEGER, SAVE :: ncidu,varidu,ncidv,varidv,ncidp,varidp
1390 INTEGER, SAVE :: ncidQ,varidQ,ncidt,varidt,ncidps,varidps
1391 INTEGER :: ncidpl,varidpl,varidap,varidbp
1392 ! Variables auxiliaires NetCDF:
1393 INTEGER, DIMENSION(4) :: start,count
1394 INTEGER :: status,rcode
1395 ! Variables for 3D extension:
1396 REAL, DIMENSION (jjp1,llm) :: zu
1397 REAL, DIMENSION (jjm,llm) :: zv
1398 INTEGER :: i
1399
1400 CHARACTER (len = 80) :: abort_message
1401 CHARACTER (len = 20) :: modname = 'guide_read2D'
1402 ! -----------------------------------------------------------------
1403 ! Premier appel: initialisation de la lecture des fichiers
1404 ! -----------------------------------------------------------------
1405 ✗ if (first) then
1406 ✗ ncidpl=-99
1407 ✗ write(*,*)trim(modname)//' : opening nudging files '
1408 ! Ap et Bp si niveaux de pression hybrides
1409 ✗ if (guide_plevs.EQ.1) then
1410 ✗ write(*,*)trim(modname)//' Reading nudging on model levels'
1411 ✗ rcode = nf90_open('apbp.nc', nf90_nowrite, ncidpl)
1412 ✗ IF (rcode.NE.NF_NOERR) THEN
1413 ✗ abort_message='Nudging: error -> no file apbp.nc'
1414 ✗ CALL abort_gcm(modname,abort_message,1)
1415 ENDIF
1416 ✗ rcode = nf90_inq_varid(ncidpl, 'AP', varidap)
1417 ✗ IF (rcode.NE.NF_NOERR) THEN
1418 ✗ abort_message='Nudging: error -> no AP variable in file apbp.nc'
1419 ✗ CALL abort_gcm(modname,abort_message,1)
1420 ENDIF
1421 ✗ rcode = nf90_inq_varid(ncidpl, 'BP', varidbp)
1422 ✗ IF (rcode.NE.NF_NOERR) THEN
1423 ✗ abort_message='Nudging: error -> no BP variable in file apbp.nc'
1424 ✗ CALL abort_gcm(modname,abort_message,1)
1425 ENDIF
1426 ✗ write(*,*)trim(modname)//'ncidpl,varidap',ncidpl,varidap
1427 endif
1428 ! Pression
1429 ✗ if (guide_plevs.EQ.2) then
1430 ✗ rcode = nf90_open('P.nc', nf90_nowrite, ncidp)
1431 ✗ IF (rcode.NE.NF_NOERR) THEN
1432 ✗ abort_message='Nudging: error -> no file P.nc'
1433 ✗ CALL abort_gcm(modname,abort_message,1)
1434 ENDIF
1435 ✗ rcode = nf90_inq_varid(ncidp, 'PRES', varidp)
1436 ✗ IF (rcode.NE.NF_NOERR) THEN
1437 ✗ abort_message='Nudging: error -> no PRES variable in file P.nc'
1438 ✗ CALL abort_gcm(modname,abort_message,1)
1439 ENDIF
1440 ✗ write(*,*)trim(modname)//' ncidp,varidp',ncidp,varidp
1441 ✗ if (ncidpl.eq.-99) ncidpl=ncidp
1442 endif
1443 ! Vent zonal
1444 ✗ if (guide_u) then
1445 ✗ rcode = nf90_open('u.nc', nf90_nowrite, ncidu)
1446 ✗ IF (rcode.NE.NF_NOERR) THEN
1447 ✗ abort_message='Nudging: error -> no file u.nc'
1448 ✗ CALL abort_gcm(modname,abort_message,1)
1449 ENDIF
1450 ✗ rcode = nf90_inq_varid(ncidu, 'UWND', varidu)
1451 ✗ IF (rcode.NE.NF_NOERR) THEN
1452 ✗ abort_message='Nudging: error -> no UWND variable in file u.nc'
1453 ✗ CALL abort_gcm(modname,abort_message,1)
1454 ENDIF
1455 ✗ write(*,*)trim(modname)//' ncidu,varidu',ncidu,varidu
1456 ✗ if (ncidpl.eq.-99) ncidpl=ncidu
1457 endif
1458 ! Vent meridien
1459 ✗ if (guide_v) then
1460 ✗ rcode = nf90_open('v.nc', nf90_nowrite, ncidv)
1461 ✗ IF (rcode.NE.NF_NOERR) THEN
1462 ✗ abort_message='Nudging: error -> no file v.nc'
1463 ✗ CALL abort_gcm(modname,abort_message,1)
1464 ENDIF
1465 ✗ rcode = nf90_inq_varid(ncidv, 'VWND', varidv)
1466 ✗ IF (rcode.NE.NF_NOERR) THEN
1467 ✗ abort_message='Nudging: error -> no VWND variable in file v.nc'
1468 ✗ CALL abort_gcm(modname,abort_message,1)
1469 ENDIF
1470 ✗ write(*,*)trim(modname)//' ncidv,varidv',ncidv,varidv
1471 ✗ if (ncidpl.eq.-99) ncidpl=ncidv
1472 endif
1473 ! Temperature
1474 ✗ if (guide_T) then
1475 ✗ rcode = nf90_open('T.nc', nf90_nowrite, ncidt)
1476 ✗ IF (rcode.NE.NF_NOERR) THEN
1477 ✗ abort_message='Nudging: error -> no file T.nc'
1478 ✗ CALL abort_gcm(modname,abort_message,1)
1479 ENDIF
1480 ✗ rcode = nf90_inq_varid(ncidt, 'AIR', varidt)
1481 ✗ IF (rcode.NE.NF_NOERR) THEN
1482 ✗ abort_message='Nudging: error -> no AIR variable in file T.nc'
1483 ✗ CALL abort_gcm(modname,abort_message,1)
1484 ENDIF
1485 ✗ write(*,*)trim(modname)//' ncidT,varidT',ncidt,varidt
1486 ✗ if (ncidpl.eq.-99) ncidpl=ncidt
1487 endif
1488 ! Humidite
1489 ✗ if (guide_Q) then
1490 ✗ rcode = nf90_open('hur.nc', nf90_nowrite, ncidQ)
1491 ✗ IF (rcode.NE.NF_NOERR) THEN
1492 ✗ abort_message='Nudging: error -> no file hur.nc'
1493 ✗ CALL abort_gcm(modname,abort_message,1)
1494 ENDIF
1495 ✗ rcode = nf90_inq_varid(ncidQ, 'RH', varidQ)
1496 ✗ IF (rcode.NE.NF_NOERR) THEN
1497 ✗ abort_message='Nudging: error -> no RH,variable in file hur.nc'
1498 ✗ CALL abort_gcm(modname,abort_message,1)
1499 ENDIF
1500 ✗ write(*,*)trim(modname)//' ncidQ,varidQ',ncidQ,varidQ
1501 ✗ if (ncidpl.eq.-99) ncidpl=ncidQ
1502 endif
1503 ! Pression de surface
1504 ✗ if ((guide_P).OR.(guide_modele)) then
1505 ✗ rcode = nf90_open('ps.nc', nf90_nowrite, ncidps)
1506 ✗ IF (rcode.NE.NF_NOERR) THEN
1507 ✗ abort_message='Nudging: error -> no file ps.nc'
1508 ✗ CALL abort_gcm(modname,abort_message,1)
1509 ENDIF
1510 ✗ rcode = nf90_inq_varid(ncidps, 'SP', varidps)
1511 ✗ IF (rcode.NE.NF_NOERR) THEN
1512 ✗ abort_message='Nudging: error -> no SP variable in file ps.nc'
1513 ✗ CALL abort_gcm(modname,abort_message,1)
1514 ENDIF
1515 ✗ write(*,*)trim(modname)//' ncidps,varidps',ncidps,varidps
1516 endif
1517 ! Coordonnee verticale
1518 ✗ if (guide_plevs.EQ.0) then
1519 ✗ rcode = nf90_inq_varid(ncidpl, 'LEVEL', varidpl)
1520 ✗ IF (rcode.NE.0) rcode = nf90_inq_varid(ncidpl, 'PRESSURE', varidpl)
1521 ✗ write(*,*)trim(modname)//' ncidpl,varidpl',ncidpl,varidpl
1522 endif
1523 ! Coefs ap, bp pour calcul de la pression aux differents niveaux
1524 ✗ if (guide_plevs.EQ.1) then
1525 ✗ status=NF_GET_VARA_DOUBLE(ncidpl,varidap,1,nlevnc,apnc)
1526 ✗ status=NF_GET_VARA_DOUBLE(ncidpl,varidbp,1,nlevnc,bpnc)
1527 ✗ elseif (guide_plevs.EQ.0) THEN
1528 ✗ status=NF_GET_VARA_DOUBLE(ncidpl,varidpl,1,nlevnc,apnc)
1529 ✗ apnc=apnc*100.! conversion en Pascals
1530 ✗ bpnc(:)=0.
1531 endif
1532 ✗ first=.FALSE.
1533 endif ! (first)
1534
1535 ! -----------------------------------------------------------------
1536 ! lecture des champs u, v, T, Q, ps
1537 ! -----------------------------------------------------------------
1538
1539 ! dimensions pour les champs scalaires et le vent zonal
1540 ✗ start(1)=1
1541 ✗ start(2)=1
1542 ✗ start(3)=1
1543 ✗ start(4)=timestep
1544
1545 ✗ count(1)=1
1546 ✗ count(2)=jjp1
1547 ✗ count(3)=nlevnc
1548 ✗ count(4)=1
1549
1550 ! Pression
1551 ✗ if (guide_plevs.EQ.2) then
1552 ✗ status=NF_GET_VARA_DOUBLE(ncidp,varidp,start,count,zu)
1553 ✗ DO i=1,iip1
1554 ✗ pnat2(i,:,:)=zu(:,:)
1555 ENDDO
1556
1557 ✗ IF (invert_y) THEN
1558 ! PRINT*,"Invertion impossible actuellement"
1559 ! CALL abort_gcm(modname,abort_message,1)
1560 ✗ CALL invert_lat(iip1,jjp1,nlevnc,pnat2)
1561 ENDIF
1562 endif
1563 ! Vent zonal
1564 ✗ if (guide_u) then
1565 ✗ status=NF_GET_VARA_DOUBLE(ncidu,varidu,start,count,zu)
1566 ✗ DO i=1,iip1
1567 ✗ unat2(i,:,:)=zu(:,:)
1568 ENDDO
1569
1570 ✗ IF (invert_y) THEN
1571 ✗ CALL invert_lat(iip1,jjp1,nlevnc,unat2)
1572 ENDIF
1573
1574 endif
1575
1576 ! Temperature
1577 ✗ if (guide_T) then
1578 ✗ status=NF_GET_VARA_DOUBLE(ncidt,varidt,start,count,zu)
1579 ✗ DO i=1,iip1
1580 ✗ tnat2(i,:,:)=zu(:,:)
1581 ENDDO
1582
1583 ✗ IF (invert_y) THEN
1584 ✗ CALL invert_lat(iip1,jjp1,nlevnc,tnat2)
1585 ENDIF
1586
1587 endif
1588
1589 ! Humidite
1590 ✗ if (guide_Q) then
1591 ✗ status=NF_GET_VARA_DOUBLE(ncidQ,varidQ,start,count,zu)
1592 ✗ DO i=1,iip1
1593 ✗ qnat2(i,:,:)=zu(:,:)
1594 ENDDO
1595
1596 ✗ IF (invert_y) THEN
1597 ✗ CALL invert_lat(iip1,jjp1,nlevnc,qnat2)
1598 ENDIF
1599
1600 endif
1601
1602 ! Vent meridien
1603 ✗ if (guide_v) then
1604 ✗ count(2)=jjm
1605 ✗ status=NF_GET_VARA_DOUBLE(ncidv,varidv,start,count,zv)
1606 ✗ DO i=1,iip1
1607 ✗ vnat2(i,:,:)=zv(:,:)
1608 ENDDO
1609
1610 ✗ IF (invert_y) THEN
1611 ✗ CALL invert_lat(iip1,jjm,nlevnc,vnat2)
1612 ENDIF
1613
1614 endif
1615
1616 ! Pression de surface
1617 ✗ if ((guide_P).OR.(guide_plevs.EQ.1)) then
1618 ✗ start(3)=timestep
1619 ✗ start(4)=0
1620 ✗ count(2)=jjp1
1621 ✗ count(3)=1
1622 ✗ count(4)=0
1623 ✗ status=NF_GET_VARA_DOUBLE(ncidps,varidps,start,count,zu(:,1))
1624 ✗ DO i=1,iip1
1625 ✗ psnat2(i,:)=zu(:,1)
1626 ENDDO
1627
1628 ✗ IF (invert_y) THEN
1629 ✗ CALL invert_lat(iip1,jjp1,1,psnat2)
1630 ENDIF
1631
1632 endif
1633
1634 ✗ END SUBROUTINE guide_read2D
1635
1636 !=======================================================================
1637 ✗ SUBROUTINE guide_out(varname,hsize,vsize,field)
1638
1639 USE comconst_mod, ONLY: pi
1640 USE comvert_mod, ONLY: presnivs
1641 use netcdf95, only: nf95_def_var, nf95_put_var
1642 use netcdf, only: nf90_float
1643
1644 IMPLICIT NONE
1645
1646 INCLUDE "dimensions.h"
1647 INCLUDE "paramet.h"
1648 INCLUDE "netcdf.inc"
1649 INCLUDE "comgeom2.h"
1650
1651 ! Variables entree
1652 CHARACTER*(*), INTENT(IN) :: varname
1653 INTEGER, INTENT (IN) :: hsize,vsize
1654 REAL, DIMENSION (iip1,hsize,vsize), INTENT(IN) :: field
1655
1656 ! Variables locales
1657 INTEGER, SAVE :: timestep=0
1658 ! Identites fichier netcdf
1659 INTEGER :: nid, id_lonu, id_lonv, id_latu, id_latv, id_tim, id_lev
1660 INTEGER :: vid_lonu,vid_lonv,vid_latu,vid_latv,vid_cu,vid_cv,vid_lev
1661 INTEGER :: vid_au,vid_av, varid_alpha_t, varid_alpha_q
1662 INTEGER, DIMENSION (3) :: dim3
1663 INTEGER, DIMENSION (4) :: dim4,count,start
1664 INTEGER :: ierr, varid,l
1665 ✗ REAL, DIMENSION (iip1,hsize,vsize) :: field2
1666 CHARACTER(LEN=20),PARAMETER :: modname="guide_out"
1667
1668 ✗ write(*,*)trim(modname)//': output timestep',timestep,'var ',varname
1669 ✗ IF (timestep.EQ.0) THEN
1670 ! ----------------------------------------------
1671 ! initialisation fichier de sortie
1672 ! ----------------------------------------------
1673 ! Ouverture du fichier
1674 ✗ ierr=NF_CREATE("guide_ins.nc",IOR(NF_CLOBBER,NF_64BIT_OFFSET),nid)
1675 ! Definition des dimensions
1676 ✗ ierr=NF_DEF_DIM(nid,"LONU",iip1,id_lonu)
1677 ✗ ierr=NF_DEF_DIM(nid,"LONV",iip1,id_lonv)
1678 ✗ ierr=NF_DEF_DIM(nid,"LATU",jjp1,id_latu)
1679 ✗ ierr=NF_DEF_DIM(nid,"LATV",jjm,id_latv)
1680 ✗ ierr=NF_DEF_DIM(nid,"LEVEL",llm,id_lev)
1681 ✗ ierr=NF_DEF_DIM(nid,"TIME",NF_UNLIMITED,id_tim)
1682
1683 ! Creation des variables dimensions
1684 ✗ ierr=NF_DEF_VAR(nid,"LONU",NF_FLOAT,1,id_lonu,vid_lonu)
1685 ✗ ierr=NF_DEF_VAR(nid,"LONV",NF_FLOAT,1,id_lonv,vid_lonv)
1686 ✗ ierr=NF_DEF_VAR(nid,"LATU",NF_FLOAT,1,id_latu,vid_latu)
1687 ✗ ierr=NF_DEF_VAR(nid,"LATV",NF_FLOAT,1,id_latv,vid_latv)
1688 ✗ ierr=NF_DEF_VAR(nid,"LEVEL",NF_FLOAT,1,id_lev,vid_lev)
1689 ✗ ierr=NF_DEF_VAR(nid,"cu",NF_FLOAT,2,(/id_lonu,id_latu/),vid_cu)
1690 ✗ ierr=NF_DEF_VAR(nid,"cv",NF_FLOAT,2,(/id_lonv,id_latv/),vid_cv)
1691 ✗ ierr=NF_DEF_VAR(nid,"au",NF_FLOAT,2,(/id_lonu,id_latu/),vid_au)
1692 ✗ ierr=NF_DEF_VAR(nid,"av",NF_FLOAT,2,(/id_lonv,id_latv/),vid_av)
1693 call nf95_def_var(nid, "alpha_T", nf90_float, (/id_lonv, id_latu/), &
1694 ✗ varid_alpha_t)
1695 call nf95_def_var(nid, "alpha_q", nf90_float, (/id_lonv, id_latu/), &
1696 ✗ varid_alpha_q)
1697
1698 ✗ ierr=NF_ENDDEF(nid)
1699
1700 ! Enregistrement des variables dimensions
1701 ✗ ierr = NF_PUT_VAR_DOUBLE(nid,vid_lonu,rlonu*180./pi)
1702 ✗ ierr = NF_PUT_VAR_DOUBLE(nid,vid_lonv,rlonv*180./pi)
1703 ✗ ierr = NF_PUT_VAR_DOUBLE(nid,vid_latu,rlatu*180./pi)
1704 ✗ ierr = NF_PUT_VAR_DOUBLE(nid,vid_latv,rlatv*180./pi)
1705 ✗ ierr = NF_PUT_VAR_DOUBLE(nid,vid_lev,presnivs)
1706 ✗ ierr = NF_PUT_VAR_DOUBLE(nid,vid_cu,cu)
1707 ✗ ierr = NF_PUT_VAR_DOUBLE(nid,vid_cv,cv)
1708 ✗ ierr = NF_PUT_VAR_DOUBLE(nid,vid_au,alpha_u)
1709 ✗ ierr = NF_PUT_VAR_DOUBLE(nid,vid_av,alpha_v)
1710 ✗ call nf95_put_var(nid, varid_alpha_t, alpha_t)
1711 ✗ call nf95_put_var(nid, varid_alpha_q, alpha_q)
1712 ! --------------------------------------------------------------------
1713 ! Cr�ation des variables sauvegard�es
1714 ! --------------------------------------------------------------------
1715 ✗ ierr = NF_REDEF(nid)
1716 ! Pressure (GCM)
1717 ✗ dim4=(/id_lonv,id_latu,id_lev,id_tim/)
1718 ✗ ierr = NF_DEF_VAR(nid,"SP",NF_FLOAT,4,dim4,varid)
1719 ! Surface pressure (guidage)
1720 ✗ IF (guide_P) THEN
1721 ✗ dim3=(/id_lonv,id_latu,id_tim/)
1722 ✗ ierr = NF_DEF_VAR(nid,"ps",NF_FLOAT,3,dim3,varid)
1723 ENDIF
1724 ! Zonal wind
1725 ✗ IF (guide_u) THEN
1726 ✗ dim4=(/id_lonu,id_latu,id_lev,id_tim/)
1727 ✗ ierr = NF_DEF_VAR(nid,"u",NF_FLOAT,4,dim4,varid)
1728 ✗ ierr = NF_DEF_VAR(nid,"ua",NF_FLOAT,4,dim4,varid)
1729 ✗ ierr = NF_DEF_VAR(nid,"ucov",NF_FLOAT,4,dim4,varid)
1730 ENDIF
1731 ! Merid. wind
1732 ✗ IF (guide_v) THEN
1733 ✗ dim4=(/id_lonv,id_latv,id_lev,id_tim/)
1734 ✗ ierr = NF_DEF_VAR(nid,"v",NF_FLOAT,4,dim4,varid)
1735 ✗ ierr = NF_DEF_VAR(nid,"va",NF_FLOAT,4,dim4,varid)
1736 ✗ ierr = NF_DEF_VAR(nid,"vcov",NF_FLOAT,4,dim4,varid)
1737 ENDIF
1738 ! Pot. Temperature
1739 ✗ IF (guide_T) THEN
1740 ✗ dim4=(/id_lonv,id_latu,id_lev,id_tim/)
1741 ✗ ierr = NF_DEF_VAR(nid,"teta",NF_FLOAT,4,dim4,varid)
1742 ENDIF
1743 ! Specific Humidity
1744 ✗ IF (guide_Q) THEN
1745 ✗ dim4=(/id_lonv,id_latu,id_lev,id_tim/)
1746 ✗ ierr = NF_DEF_VAR(nid,"q",NF_FLOAT,4,dim4,varid)
1747 ENDIF
1748
1749 ✗ ierr = NF_ENDDEF(nid)
1750 ✗ ierr = NF_CLOSE(nid)
1751 ENDIF ! timestep=0
1752
1753 ! --------------------------------------------------------------------
1754 ! Enregistrement du champ
1755 ! --------------------------------------------------------------------
1756 ✗ ierr=NF_OPEN("guide_ins.nc",NF_WRITE,nid)
1757
1758 ✗ IF (varname=="SP") timestep=timestep+1
1759
1760 ✗ ierr = NF_INQ_VARID(nid,varname,varid)
1761 ✗ SELECT CASE (varname)
1762 CASE ("SP","ps")
1763 ✗ start=(/1,1,1,timestep/)
1764 ✗ count=(/iip1,jjp1,llm,1/)
1765 CASE ("v","va","vcov")
1766 ✗ start=(/1,1,1,timestep/)
1767 ✗ count=(/iip1,jjm,llm,1/)
1768 CASE DEFAULT
1769 ✗ start=(/1,1,1,timestep/)
1770 ✗ count=(/iip1,jjp1,llm,1/)
1771 END SELECT
1772
1773 SELECT CASE (varname)
1774 CASE("u","ua")
1775 ✗ DO l=1,llm ; field2(:,2:jjm,l)=field(:,2:jjm,l)/cu(:,2:jjm) ; ENDDO
1776 ✗ field2(:,1,:)=0. ; field2(:,jjp1,:)=0.
1777 CASE("v","va")
1778 ✗ DO l=1,llm ; field2(:,:,l)=field(:,:,l)/cv(:,:) ; ENDDO
1779 CASE DEFAULT
1780 ✗ field2=field
1781 END SELECT
1782
1783
1784 ✗ ierr = NF_PUT_VARA_DOUBLE(nid,varid,start,count,field2)
1785
1786 ✗ ierr = NF_CLOSE(nid)
1787
1788 ✗ END SUBROUTINE guide_out
1789
1790
1791 !===========================================================================
1792 ✗ subroutine correctbid(iim,nl,x)
1793 integer iim,nl
1794 real x(iim+1,nl)
1795 integer i,l
1796 real zz
1797
1798 ✗ do l=1,nl
1799 ✗ do i=2,iim-1
1800 ✗ if(abs(x(i,l)).gt.1.e10) then
1801 ✗ zz=0.5*(x(i-1,l)+x(i+1,l))
1802 ✗ print*,'correction ',i,l,x(i,l),zz
1803 ✗ x(i,l)=zz
1804 endif
1805 enddo
1806 enddo
1807 ✗ return
1808 end subroutine correctbid
1809
1810 !===========================================================================
1811 END MODULE guide_mod
1812