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