source: LMDZ4/trunk/libf/dyn3dpar/guide_p_mod.F90 @ 5353

Last change on this file since 5353 was 1418, checked in by jghattas, 14 years ago

Following bug corrections are already done at branch LMDZ4_AR5 rev
1417,1416 :

Bug corrections for nudged run (pres2lev.F90, guide_p_mod.F90) :

  • now the results are the same for sequentiel and parallel mode(if adjust=n and compiled with vsafe at mercure).
  • the results are the same as the sequential mode in previous revision.
  • test done only with guide_u=y,guide_v=y
  • copied optimized version of pres2lev.F90 from dyn3dpar to dyn3d


Added condition read_climoz for the variable O3daySTD(calcul_STDlev.h,
undefSTD.F)

Corrected bug in calculation of the diagnostic variable ec550aer
"Extinction at 550nm" (aeropt_5wv.F90) (Maria Raffaella Vuolo, LSCE)

Added stop if use_filtre_fft=y with dyn3d (conf_gcm.F) : this option is
not implemented in dyn3d.

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