source: LMDZ5/branches/LF-private/libf/dyn3dpar/guide_p_mod.F90 @ 3642

Last change on this file since 3642 was 1823, checked in by Ehouarn Millour, 11 years ago

Remplacement de parallel.F90 (en conflit avec orchidée) par parallel_lmdz.F90.
UG
.........................................
Renaming parallel.F90 (conflicting with orchidée) into parallel_lmdz.F90.
UG

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