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