source: LMDZ4/branches/LMDZ4-dev/libf/dyn3d/guide_mod.F90 @ 1270

Last change on this file since 1270 was 1188, checked in by lguez, 15 years ago

Translated calls using NetCDF 2.4 interface to calls using NetCDF 3.6
Fortran 90 interface.
Corrected a few comments in "output.def".

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