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