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

Last change on this file since 1348 was 1304, checked in by musat, 15 years ago

Bug correction for the vertical dimension in the calling routine for the latitude inversion of nudging fields
FC/IM

  • Property svn:executable set to *
File size: 56.0 KB
Line 
1!
2! $Header: /home/cvsroot/LMDZ4/libf/dyn3d/guide.F,v 1.3.4.1 2006/11/06 15:51:16 fairhead Exp $
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    IMPLICIT NONE
69 
70    INCLUDE "dimensions.h"
71    INCLUDE "paramet.h"
72    INCLUDE "netcdf.inc"
73    INCLUDE "control.h"
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   
277    IMPLICIT NONE
278 
279    INCLUDE "dimensions.h"
280    INCLUDE "paramet.h"
281    INCLUDE "control.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/FLOAT(iguide_read)
383      ELSE
384          tau=FLOAT(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)/FLOAT(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)/FLOAT(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.