source: LMDZ5/branches/testing/libf/dyn3d/guide_mod.F90 @ 2435

Last change on this file since 2435 was 2298, checked in by Laurent Fairhead, 10 years ago

Merged trunk changes -r2237:2291 into testing branch

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