source: LMDZ5/branches/LMDZ5_AR5/libf/dyn3d/guide_mod.F90 @ 3110

Last change on this file since 3110 was 1520, checked in by Ehouarn Millour, 13 years ago

Implementation of a different vertical discretization (from/for planets, but
can in principle also be used for Earth).
Choice of vertical discretization is set by flag 'disvert_type';
'disvert_type=1' is Earth standard (default; ie set to 1 if
planet_type=="earth") case.
With 'disvert_type=2', approximate altitude of layers and reference atmospheric
scale height must be given using an input file ("z2sig.def", first line
should give scale height, in km, following lines must specify the altitude,
in km above surface, of mid-layers, one per line; see disvert_noterre.F).

Checked that these changes do not impact on 'bench' results, on Vargas.

EM.

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