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