source: LMDZ6/branches/Ocean_skin/libf/dyn3d/guide_mod.F90

Last change on this file was 3103, checked in by oboucher, 7 years ago

Removing x permission from these files

  • 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
File size: 59.2 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, ONLY: day_step
66    USE serre_mod, ONLY: grossismx
67
68    IMPLICIT NONE
69 
70    INCLUDE "dimensions.h"
71    INCLUDE "paramet.h"
72    INCLUDE "netcdf.inc"
73
74    INTEGER                :: error,ncidpl,rid,rcod
75    CHARACTER (len = 80)   :: abort_message
76    CHARACTER (len = 20)   :: modname = 'guide_init'
77
78! ---------------------------------------------
79! Lecture des parametres: 
80! ---------------------------------------------
81    call ini_getparam("nudging_parameters_out.txt")
82! Variables guidees
83    CALL getpar('guide_u',.true.,guide_u,'guidage de u')
84    CALL getpar('guide_v',.true.,guide_v,'guidage de v')
85    CALL getpar('guide_T',.true.,guide_T,'guidage de T')
86    CALL getpar('guide_P',.true.,guide_P,'guidage de P')
87    CALL getpar('guide_Q',.true.,guide_Q,'guidage de Q')
88    CALL getpar('guide_hr',.true.,guide_hr,'guidage de Q par H.R')
89    CALL getpar('guide_teta',.false.,guide_teta,'guidage de T par Teta')
90
91    CALL getpar('guide_add',.false.,guide_add,'for�age constant?')
92    CALL getpar('guide_zon',.false.,guide_zon,'guidage moy zonale')
93    if (guide_zon .and. abs(grossismx - 1.) > 0.01) &
94         call abort_gcm("guide_init", &
95         "zonal nudging requires grid regular in longitude", 1)
96
97!   Constantes de rappel. Unite : fraction de jour
98    CALL getpar('tau_min_u',0.02,tau_min_u,'Cste de rappel min, u')
99    CALL getpar('tau_max_u', 10.,tau_max_u,'Cste de rappel max, u')
100    CALL getpar('tau_min_v',0.02,tau_min_v,'Cste de rappel min, v')
101    CALL getpar('tau_max_v', 10.,tau_max_v,'Cste de rappel max, v')
102    CALL getpar('tau_min_T',0.02,tau_min_T,'Cste de rappel min, T')
103    CALL getpar('tau_max_T', 10.,tau_max_T,'Cste de rappel max, T')
104    CALL getpar('tau_min_Q',0.02,tau_min_Q,'Cste de rappel min, Q')
105    CALL getpar('tau_max_Q', 10.,tau_max_Q,'Cste de rappel max, Q')
106    CALL getpar('tau_min_P',0.02,tau_min_P,'Cste de rappel min, P')
107    CALL getpar('tau_max_P', 10.,tau_max_P,'Cste de rappel max, P')
108    CALL getpar('gamma4',.false.,gamma4,'Zone sans rappel elargie')
109    CALL getpar('guide_BL',.true.,guide_BL,'guidage dans C.Lim')
110
111! Sauvegarde du for�age
112    CALL getpar('guide_sav',.false.,guide_sav,'sauvegarde guidage')
113    CALL getpar('iguide_sav',4,iguide_sav,'freq. sauvegarde guidage')
114    ! frequences f>0: fx/jour; f<0: tous les f jours; f=0: 1 seule fois.
115    IF (iguide_sav.GT.0) THEN
116       iguide_sav=day_step/iguide_sav
117    ELSE if (iguide_sav == 0) then
118       iguide_sav = huge(0)
119    ELSE
120       iguide_sav=day_step*iguide_sav
121    ENDIF
122
123! Guidage regional seulement (sinon constant ou suivant le zoom)
124    CALL getpar('guide_reg',.false.,guide_reg,'guidage regional')
125    CALL getpar('lat_min_g',-90.,lat_min_g,'Latitude mini guidage ')
126    CALL getpar('lat_max_g', 90.,lat_max_g,'Latitude maxi guidage ')
127    CALL getpar('lon_min_g',-180.,lon_min_g,'longitude mini guidage ')
128    CALL getpar('lon_max_g', 180.,lon_max_g,'longitude maxi guidage ')
129    CALL getpar('tau_lat', 5.,tau_lat,'raideur lat guide regional ')
130    CALL getpar('tau_lon', 5.,tau_lon,'raideur lon guide regional ')
131
132! Parametres pour lecture des fichiers
133    CALL getpar('iguide_read',4,iguide_read,'freq. lecture guidage')
134    CALL getpar('iguide_int',4,iguide_int,'freq. interpolation vert')
135    IF (iguide_int.EQ.0) THEN
136        iguide_int=1
137    ELSEIF (iguide_int.GT.0) THEN
138        iguide_int=day_step/iguide_int
139    ELSE
140        iguide_int=day_step*iguide_int
141    ENDIF
142    CALL getpar('guide_modele',.false.,guide_modele,'guidage niveaux modele')
143    CALL getpar('ini_anal',.false.,ini_anal,'Etat initial = analyse')
144    CALL getpar('guide_invertp',.true.,invert_p,'niveaux p inverses')
145    CALL getpar('guide_inverty',.true.,invert_y,'inversion N-S')
146    CALL getpar('guide_2D',.false.,guide_2D,'fichier guidage lat-P')
147
148    call fin_getparam
149   
150! ---------------------------------------------
151! Determination du nombre de niveaux verticaux
152! des fichiers guidage
153! ---------------------------------------------
154    ncidpl=-99
155    if (guide_modele) then
156       if (ncidpl.eq.-99) then
157          rcod=nf90_open('apbp.nc',Nf90_NOWRITe, ncidpl)
158          if (rcod.NE.NF_NOERR) THEN
159             CALL abort_gcm(modname, &
160                  'Guide: probleme -> pas de fichier apbp.nc',1)
161          endif
162       endif
163    else
164         if (guide_u) then
165           if (ncidpl.eq.-99) then
166               rcod=nf90_open('u.nc',Nf90_NOWRITe,ncidpl)
167               if (rcod.NE.NF_NOERR) THEN
168                  CALL abort_gcm(modname, &
169                       'Guide: probleme -> pas de fichier u.nc',1)
170               endif
171           endif
172         elseif (guide_v) then
173           if (ncidpl.eq.-99) then
174               rcod=nf90_open('v.nc',nf90_nowrite,ncidpl)
175               if (rcod.NE.NF_NOERR) THEN
176                  CALL abort_gcm(modname, &
177                       'Guide: probleme -> pas de fichier v.nc',1)
178               endif
179           endif
180         elseif (guide_T) then
181           if (ncidpl.eq.-99) then
182               rcod=nf90_open('T.nc',nf90_nowrite,ncidpl)
183               if (rcod.NE.NF_NOERR) THEN
184                  CALL abort_gcm(modname, &
185                       'Guide: probleme -> pas de fichier T.nc',1)
186               endif
187           endif
188         elseif (guide_Q) then
189           if (ncidpl.eq.-99) then
190               rcod=nf90_open('hur.nc',nf90_nowrite, ncidpl)
191               if (rcod.NE.NF_NOERR) THEN
192                  CALL abort_gcm(modname, &
193                       'Guide: probleme -> pas de fichier hur.nc',1)
194               endif
195           endif
196         endif
197    endif
198    error=NF_INQ_DIMID(ncidpl,'LEVEL',rid)
199    IF (error.NE.NF_NOERR) error=NF_INQ_DIMID(ncidpl,'PRESSURE',rid)
200    IF (error.NE.NF_NOERR) THEN
201        CALL abort_gcm(modname,'Guide: probleme lecture niveaux pression',1)
202    ENDIF
203    error=NF_INQ_DIMLEN(ncidpl,rid,nlevnc)
204    print *,'Guide: nombre niveaux vert. nlevnc', nlevnc
205    rcod = nf90_close(ncidpl)
206
207! ---------------------------------------------
208! Allocation des variables
209! ---------------------------------------------
210    abort_message='pb in allocation guide'
211
212    ALLOCATE(apnc(nlevnc), stat = error)
213    IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
214    ALLOCATE(bpnc(nlevnc), stat = error)
215    IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
216    apnc=0.;bpnc=0.
217
218    ALLOCATE(alpha_pcor(llm), stat = error)
219    IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
220    ALLOCATE(alpha_u(ip1jmp1), stat = error)
221    IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
222    ALLOCATE(alpha_v(ip1jm), stat = error)
223    IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
224    ALLOCATE(alpha_T(iip1, jjp1), stat = error)
225    IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
226    ALLOCATE(alpha_Q(iip1, jjp1), stat = error)
227    IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
228    ALLOCATE(alpha_P(ip1jmp1), stat = error)
229    IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
230    alpha_u=0.;alpha_v=0;alpha_T=0;alpha_Q=0;alpha_P=0
231   
232    IF (guide_u) THEN
233        ALLOCATE(unat1(iip1,jjp1,nlevnc), stat = error)
234        IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
235        ALLOCATE(ugui1(ip1jmp1,llm), stat = error)
236        IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
237        ALLOCATE(unat2(iip1,jjp1,nlevnc), stat = error)
238        IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
239        ALLOCATE(ugui2(ip1jmp1,llm), stat = error)
240        IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
241        unat1=0.;unat2=0.;ugui1=0.;ugui2=0.
242    ENDIF
243
244    IF (guide_T) THEN
245        ALLOCATE(tnat1(iip1,jjp1,nlevnc), stat = error)
246        IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
247        ALLOCATE(tgui1(ip1jmp1,llm), stat = error)
248        IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
249        ALLOCATE(tnat2(iip1,jjp1,nlevnc), stat = error)
250        IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
251        ALLOCATE(tgui2(ip1jmp1,llm), stat = error)
252        IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
253        tnat1=0.;tnat2=0.;tgui1=0.;tgui2=0.
254    ENDIF
255     
256    IF (guide_Q) THEN
257        ALLOCATE(qnat1(iip1,jjp1,nlevnc), stat = error)
258        IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
259        ALLOCATE(qgui1(ip1jmp1,llm), stat = error)
260        IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
261        ALLOCATE(qnat2(iip1,jjp1,nlevnc), stat = error)
262        IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
263        ALLOCATE(qgui2(ip1jmp1,llm), stat = error)
264        IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
265        qnat1=0.;qnat2=0.;qgui1=0.;qgui2=0.
266    ENDIF
267
268    IF (guide_v) THEN
269        ALLOCATE(vnat1(iip1,jjm,nlevnc), stat = error)
270        IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
271        ALLOCATE(vgui1(ip1jm,llm), stat = error)
272        IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
273        ALLOCATE(vnat2(iip1,jjm,nlevnc), stat = error)
274        IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
275        ALLOCATE(vgui2(ip1jm,llm), stat = error)
276        IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
277        vnat1=0.;vnat2=0.;vgui1=0.;vgui2=0.
278    ENDIF
279
280    IF (guide_P.OR.guide_modele) THEN
281        ALLOCATE(psnat1(iip1,jjp1), stat = error)
282        IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
283        ALLOCATE(psnat2(iip1,jjp1), stat = error)
284        IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
285        psnat1=0.;psnat2=0.;
286    ENDIF
287    IF (guide_P) THEN
288        ALLOCATE(psgui2(ip1jmp1), stat = error)
289        IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
290        ALLOCATE(psgui1(ip1jmp1), stat = error)
291        IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
292        psgui1=0.;psgui2=0.
293    ENDIF
294
295! ---------------------------------------------
296!   Lecture du premier etat de guidage.
297! ---------------------------------------------
298    IF (guide_2D) THEN
299        CALL guide_read2D(1)
300    ELSE
301        CALL guide_read(1)
302    ENDIF
303    IF (guide_v) vnat1=vnat2
304    IF (guide_u) unat1=unat2
305    IF (guide_T) tnat1=tnat2
306    IF (guide_Q) qnat1=qnat2
307    IF (guide_P.OR.guide_modele) psnat1=psnat2
308
309  END SUBROUTINE guide_init
310
311!=======================================================================
312  SUBROUTINE guide_main(itau,ucov,vcov,teta,q,masse,ps)
313
314    USE control_mod, ONLY: day_step, iperiod
315    USE comconst_mod, ONLY: dtvr, daysec
316    USE comvert_mod, ONLY: ap, bp, preff, presnivs
317 
318    IMPLICIT NONE
319 
320    INCLUDE "dimensions.h"
321    INCLUDE "paramet.h"
322
323    ! Variables entree
324    INTEGER,                       INTENT(IN)    :: itau !pas de temps
325    REAL, DIMENSION (ip1jmp1,llm), INTENT(INOUT) :: ucov,teta,q,masse
326    REAL, DIMENSION (ip1jm,llm),   INTENT(INOUT) :: vcov
327    REAL, DIMENSION (ip1jmp1),     INTENT(INOUT) :: ps
328
329    ! Variables locales
330    LOGICAL, SAVE :: first=.TRUE.
331    LOGICAL       :: f_out ! sortie guidage
332    REAL, DIMENSION (ip1jmp1,llm) :: f_add ! var aux: champ de guidage
333    REAL, DIMENSION (ip1jmp1,llm) :: p ! besoin si guide_P
334    ! Compteurs temps:
335    INTEGER, SAVE :: step_rea,count_no_rea,itau_test ! lecture guidage
336    REAL          :: ditau, dday_step
337    REAL          :: tau,reste ! position entre 2 etats de guidage
338    REAL, SAVE    :: factt ! pas de temps en fraction de jour
339   
340    INTEGER       :: l
341
342!-----------------------------------------------------------------------
343! Initialisations au premier passage
344!-----------------------------------------------------------------------
345    IF (first) THEN
346        first=.FALSE.
347        CALL guide_init
348        itau_test=1001
349        step_rea=1
350        count_no_rea=0
351! Calcul des constantes de rappel
352        factt=dtvr*iperiod/daysec
353        call tau2alpha(3,iip1,jjm ,factt,tau_min_v,tau_max_v,alpha_v)
354        call tau2alpha(2,iip1,jjp1,factt,tau_min_u,tau_max_u,alpha_u)
355        call tau2alpha(1,iip1,jjp1,factt,tau_min_T,tau_max_T,alpha_T)
356        call tau2alpha(1,iip1,jjp1,factt,tau_min_P,tau_max_P,alpha_P)
357        call tau2alpha(1,iip1,jjp1,factt,tau_min_Q,tau_max_Q,alpha_Q)
358! correction de rappel dans couche limite
359        if (guide_BL) then
360             alpha_pcor(:)=1.
361        else
362            do l=1,llm
363                alpha_pcor(l)=(1.+tanh((0.85-presnivs(l)/preff)/0.05))/2.
364            enddo
365        endif
366! ini_anal: etat initial egal au guidage       
367        IF (ini_anal) THEN
368            CALL guide_interp(ps,teta)
369            IF (guide_u) ucov=ugui2
370            IF (guide_v) vcov=ugui2
371            IF (guide_T) teta=tgui2
372            IF (guide_Q) q=qgui2
373            IF (guide_P) THEN
374                ps=psgui2
375                CALL pression(ip1jmp1,ap,bp,ps,p)
376                CALL massdair(p,masse)
377            ENDIF
378            RETURN
379        ENDIF
380! Verification structure guidage
381        IF (guide_u) THEN
382            CALL writefield('unat',unat1)
383            CALL writefield('ucov',RESHAPE(ucov,(/iip1,jjp1,llm/)))
384        ENDIF
385        IF (guide_T) THEN
386            CALL writefield('tnat',tnat1)
387            CALL writefield('teta',RESHAPE(teta,(/iip1,jjp1,llm/)))
388        ENDIF
389
390    ENDIF !first
391
392!-----------------------------------------------------------------------
393! Lecture des fichiers de guidage ?
394!-----------------------------------------------------------------------
395    IF (iguide_read.NE.0) THEN
396      ditau=real(itau)
397      dday_step=real(day_step)
398      IF (iguide_read.LT.0) THEN
399          tau=ditau/dday_step/REAL(iguide_read)
400      ELSE
401          tau=REAL(iguide_read)*ditau/dday_step
402      ENDIF
403      reste=tau-AINT(tau)
404      IF (reste.EQ.0.) THEN
405          IF (itau_test.EQ.itau) THEN
406              write(*,*)'deuxieme passage de advreel a itau=',itau
407              stop
408          ELSE
409              IF (guide_v) vnat1=vnat2
410              IF (guide_u) unat1=unat2
411              IF (guide_T) tnat1=tnat2
412              IF (guide_Q) qnat1=qnat2
413              IF (guide_P.OR.guide_modele) psnat1=psnat2
414              step_rea=step_rea+1
415              itau_test=itau
416              print*,'Lecture fichiers guidage, pas ',step_rea, &
417                    'apres ',count_no_rea,' non lectures'
418              IF (guide_2D) THEN
419                  CALL guide_read2D(step_rea)
420              ELSE
421                  CALL guide_read(step_rea)
422              ENDIF
423              count_no_rea=0
424          ENDIF
425      ELSE
426        count_no_rea=count_no_rea+1
427
428      ENDIF
429    ENDIF !iguide_read=0
430
431!-----------------------------------------------------------------------
432! Interpolation et conversion des champs de guidage
433!-----------------------------------------------------------------------
434    IF (MOD(itau,iguide_int).EQ.0) THEN
435        CALL guide_interp(ps,teta)
436    ENDIF
437! Repartition entre 2 etats de guidage
438    IF (iguide_read.NE.0) THEN
439        tau=reste
440    ELSE
441        tau=1.
442    ENDIF
443
444!-----------------------------------------------------------------------
445!   Ajout des champs de guidage
446!-----------------------------------------------------------------------
447! Sauvegarde du guidage?
448    f_out=((MOD(itau,iguide_sav).EQ.0).AND.guide_sav) 
449    IF (f_out) CALL guide_out("SP",jjp1,1,ps)
450   
451    if (guide_u) then
452        if (guide_add) then
453           f_add=(1.-tau)*ugui1+tau*ugui2
454        else
455           f_add=(1.-tau)*ugui1+tau*ugui2-ucov
456        endif
457        if (guide_zon) CALL guide_zonave(1,jjp1,llm,f_add)
458        CALL guide_addfield(ip1jmp1,llm,f_add,alpha_u)
459        IF (f_out) CALL guide_out("ua",jjp1,llm,(1.-tau)*ugui1+tau*ugui2)
460        IF (f_out) CALL guide_out("u",jjp1,llm,ucov)
461        IF (f_out) CALL guide_out("ucov",jjp1,llm,f_add/factt)
462        ucov=ucov+f_add
463    endif
464
465    if (guide_T) then
466        if (guide_add) then
467           f_add=(1.-tau)*tgui1+tau*tgui2
468        else
469           f_add=(1.-tau)*tgui1+tau*tgui2-teta
470        endif
471        if (guide_zon) CALL guide_zonave(2,jjp1,llm,f_add)
472        CALL guide_addfield(ip1jmp1,llm,f_add,alpha_T)
473        IF (f_out) CALL guide_out("teta",jjp1,llm,f_add/factt)
474        teta=teta+f_add
475    endif
476
477    if (guide_P) then
478        if (guide_add) then
479           f_add(1:ip1jmp1,1)=(1.-tau)*psgui1+tau*psgui2
480        else
481           f_add(1:ip1jmp1,1)=(1.-tau)*psgui1+tau*psgui2-ps
482        endif
483        if (guide_zon) CALL guide_zonave(2,jjp1,1,f_add(1:ip1jmp1,1))
484        CALL guide_addfield(ip1jmp1,1,f_add(1:ip1jmp1,1),alpha_P)
485        IF (f_out) CALL guide_out("ps",jjp1,1,f_add(1:ip1jmp1,1)/factt)
486        ps=ps+f_add(1:ip1jmp1,1)
487        CALL pression(ip1jmp1,ap,bp,ps,p)
488        CALL massdair(p,masse)
489    endif
490
491    if (guide_Q) then
492        if (guide_add) then
493           f_add=(1.-tau)*qgui1+tau*qgui2
494        else
495           f_add=(1.-tau)*qgui1+tau*qgui2-q
496        endif
497        if (guide_zon) CALL guide_zonave(2,jjp1,llm,f_add)
498        CALL guide_addfield(ip1jmp1,llm,f_add,alpha_Q)
499        IF (f_out) CALL guide_out("q",jjp1,llm,f_add/factt)
500        q=q+f_add
501    endif
502
503    if (guide_v) then
504        if (guide_add) then
505           f_add(1:ip1jm,:)=(1.-tau)*vgui1+tau*vgui2
506        else
507           f_add(1:ip1jm,:)=(1.-tau)*vgui1+tau*vgui2-vcov
508        endif
509        if (guide_zon) CALL guide_zonave(2,jjm,llm,f_add(1:ip1jm,:))
510        CALL guide_addfield(ip1jm,llm,f_add(1:ip1jm,:),alpha_v)
511        IF (f_out) CALL guide_out("v",jjm,llm,vcov)
512        IF (f_out) CALL guide_out("va",jjm,llm,(1.-tau)*vgui1+tau*vgui2)
513        IF (f_out) CALL guide_out("vcov",jjm,llm,f_add(1:ip1jm,:)/factt)
514        vcov=vcov+f_add(1:ip1jm,:)
515    endif
516  END SUBROUTINE guide_main
517
518!=======================================================================
519  SUBROUTINE guide_addfield(hsize,vsize,field,alpha)
520! field1=a*field1+alpha*field2
521
522    IMPLICIT NONE
523
524    ! input variables
525    INTEGER,                      INTENT(IN)    :: hsize
526    INTEGER,                      INTENT(IN)    :: vsize
527    REAL, DIMENSION(hsize),       INTENT(IN)    :: alpha
528    REAL, DIMENSION(hsize,vsize), INTENT(INOUT) :: field
529
530    ! Local variables
531    INTEGER :: l
532
533    do l=1,vsize
534        field(:,l)=alpha*field(:,l)*alpha_pcor(l)
535    enddo
536
537  END SUBROUTINE guide_addfield
538
539!=======================================================================
540  SUBROUTINE guide_zonave(typ,hsize,vsize,field)
541
542    USE comconst_mod, ONLY: pi
543   
544    IMPLICIT NONE
545
546    INCLUDE "dimensions.h"
547    INCLUDE "paramet.h"
548    INCLUDE "comgeom.h"
549   
550    ! input/output variables
551    INTEGER,                           INTENT(IN)    :: typ
552    INTEGER,                           INTENT(IN)    :: vsize
553    INTEGER,                           INTENT(IN)    :: hsize
554    REAL, DIMENSION(hsize*iip1,vsize), INTENT(INOUT) :: field
555
556    ! Local variables
557    LOGICAL, SAVE                :: first=.TRUE.
558    INTEGER, DIMENSION (2), SAVE :: imin, imax ! averaging domain
559    INTEGER                      :: i,j,l,ij
560    REAL, DIMENSION (iip1)       :: lond       ! longitude in Deg.
561    REAL, DIMENSION (hsize,vsize):: fieldm     ! zon-averaged field
562
563    IF (first) THEN
564        first=.FALSE.
565!Compute domain for averaging
566        lond=rlonu*180./pi
567        imin(1)=1;imax(1)=iip1;
568        imin(2)=1;imax(2)=iip1;
569        IF (guide_reg) THEN
570            DO i=1,iim
571                IF (lond(i).LT.lon_min_g) imin(1)=i
572                IF (lond(i).LE.lon_max_g) imax(1)=i
573            ENDDO
574            lond=rlonv*180./pi
575            DO i=1,iim
576                IF (lond(i).LT.lon_min_g) imin(2)=i
577                IF (lond(i).LE.lon_max_g) imax(2)=i
578            ENDDO
579        ENDIF
580    ENDIF
581
582    fieldm=0.
583    DO l=1,vsize
584    ! Compute zonal average
585        DO j=1,hsize
586            DO i=imin(typ),imax(typ)
587                ij=(j-1)*iip1+i
588                fieldm(j,l)=fieldm(j,l)+field(ij,l)
589            ENDDO
590        ENDDO
591        fieldm(:,l)=fieldm(:,l)/REAL(imax(typ)-imin(typ)+1)
592    ! Compute forcing
593        DO j=1,hsize
594            DO i=1,iip1
595                ij=(j-1)*iip1+i
596                field(ij,l)=fieldm(j,l)
597            ENDDO
598        ENDDO
599    ENDDO
600
601  END SUBROUTINE guide_zonave
602
603!=======================================================================
604  SUBROUTINE guide_interp(psi,teta)
605 
606  use exner_hyb_m, only: exner_hyb
607  use exner_milieu_m, only: exner_milieu
608  use comconst_mod, only: kappa, cpp
609  use comvert_mod, only: preff, pressure_exner, bp, ap
610  IMPLICIT NONE
611
612  include "dimensions.h"
613  include "paramet.h"
614  include "comgeom2.h"
615
616  REAL, DIMENSION (iip1,jjp1),     INTENT(IN) :: psi ! Psol gcm
617  REAL, DIMENSION (iip1,jjp1,llm), INTENT(IN) :: teta ! Temp. Pot. gcm
618
619  LOGICAL, SAVE                      :: first=.TRUE.
620  ! Variables pour niveaux pression:
621  REAL, DIMENSION (iip1,jjp1,nlevnc) :: plnc1,plnc2 !niveaux pression guidage
622  REAL, DIMENSION (iip1,jjp1,llm)    :: plunc,plsnc !niveaux pression modele
623  REAL, DIMENSION (iip1,jjm,llm)     :: plvnc       !niveaux pression modele
624  REAL, DIMENSION (iip1,jjp1,llmp1)  :: p           ! pression intercouches
625  REAL, DIMENSION (iip1,jjp1,llm)    :: pls, pext   ! var intermediaire
626  REAL, DIMENSION (iip1,jjp1,llm)    :: pbarx
627  REAL, DIMENSION (iip1,jjm,llm)     :: pbary
628  ! Variables pour fonction Exner (P milieu couche)
629  REAL, DIMENSION (iip1,jjp1,llm)    :: pk
630  REAL, DIMENSION (iip1,jjp1)        :: pks   
631  REAL                               :: prefkap,unskap
632  ! Pression de vapeur saturante
633  REAL, DIMENSION (ip1jmp1,llm)      :: qsat
634  !Variables intermediaires interpolation
635  REAL, DIMENSION (iip1,jjp1,llm)    :: zu1,zu2
636  REAL, DIMENSION (iip1,jjm,llm)     :: zv1,zv2
637 
638  INTEGER                            :: i,j,l,ij
639 
640    print *,'Guide: conversion variables guidage'
641! -----------------------------------------------------------------
642! Calcul des niveaux de pression champs guidage
643! -----------------------------------------------------------------
644if (guide_modele) then
645    do i=1,iip1
646        do j=1,jjp1
647            do l=1,nlevnc
648                plnc2(i,j,l)=apnc(l)+bpnc(l)*psnat2(i,j)
649                plnc1(i,j,l)=apnc(l)+bpnc(l)*psnat1(i,j)
650            enddo
651        enddo
652    enddo
653else
654    do i=1,iip1
655        do j=1,jjp1
656            do l=1,nlevnc
657                plnc2(i,j,l)=apnc(l)
658                plnc1(i,j,l)=apnc(l)
659           enddo
660        enddo
661    enddo
662
663endif
664    if (first) then
665        first=.FALSE.
666        print*,'Guide: verification ordre niveaux verticaux'
667        print*,'LMDZ :'
668        do l=1,llm
669            print*,'PL(',l,')=',(ap(l)+ap(l+1))/2. &
670                  +psi(1,jjp1)*(bp(l)+bp(l+1))/2.
671        enddo
672        print*,'Fichiers guidage'
673        do l=1,nlevnc
674             print*,'PL(',l,')=',plnc2(1,1,l)
675        enddo
676        print *,'inversion de l''ordre: invert_p=',invert_p
677        if (guide_u) then
678            do l=1,nlevnc
679                print*,'U(',l,')=',unat2(1,1,l)
680            enddo
681        endif
682        if (guide_T) then
683            do l=1,nlevnc
684                print*,'T(',l,')=',tnat2(1,1,l)
685            enddo
686        endif
687    endif
688   
689! -----------------------------------------------------------------
690! Calcul niveaux pression modele
691! -----------------------------------------------------------------
692    CALL pression( ip1jmp1, ap, bp, psi, p )
693    if (pressure_exner) then
694      CALL exner_hyb(ip1jmp1,psi,p,pks,pk)
695    else
696      CALL exner_milieu(ip1jmp1,psi,p,pks,pk)
697    endif
698!    ....  Calcul de pls , pression au milieu des couches ,en Pascals
699    unskap=1./kappa
700    prefkap =  preff  ** kappa
701    DO l = 1, llm
702        DO j=1,jjp1
703            DO i =1, iip1
704                pls(i,j,l) = preff * ( pk(i,j,l)/cpp) ** unskap
705            ENDDO
706        ENDDO
707    ENDDO
708
709!   calcul des pressions pour les grilles u et v
710    do l=1,llm
711        do j=1,jjp1
712            do i=1,iip1
713                pext(i,j,l)=pls(i,j,l)*aire(i,j)
714            enddo
715        enddo
716    enddo
717    call massbar(pext, pbarx, pbary )
718    do l=1,llm
719        do j=1,jjp1
720            do i=1,iip1
721                plunc(i,j,l)=pbarx(i,j,l)/aireu(i,j)
722                plsnc(i,j,l)=pls(i,j,l)
723            enddo
724        enddo
725    enddo
726    do l=1,llm
727        do j=1,jjm
728            do i=1,iip1
729                plvnc(i,j,l)=pbary(i,j,l)/airev(i,j)
730            enddo
731        enddo
732    enddo
733
734! -----------------------------------------------------------------
735! Interpolation champs guidage sur niveaux modele (+inversion N/S)
736! Conversion en variables gcm (ucov, vcov...)
737! -----------------------------------------------------------------
738    if (guide_P) then
739        do j=1,jjp1
740            do i=1,iim
741                ij=(j-1)*iip1+i
742                psgui1(ij)=psnat1(i,j)
743                psgui2(ij)=psnat2(i,j)
744            enddo
745            psgui1(iip1*j)=psnat1(1,j)
746            psgui2(iip1*j)=psnat2(1,j)
747        enddo
748    endif
749
750    IF (guide_u) THEN
751        CALL pres2lev(unat1,zu1,nlevnc,llm,plnc1,plunc,iip1,jjp1,invert_p)
752        CALL pres2lev(unat2,zu2,nlevnc,llm,plnc2,plunc,iip1,jjp1,invert_p)
753        do l=1,llm
754            do j=1,jjp1
755                do i=1,iim
756                    ij=(j-1)*iip1+i
757                    ugui1(ij,l)=zu1(i,j,l)*cu(i,j)
758                    ugui2(ij,l)=zu2(i,j,l)*cu(i,j)
759                enddo
760                ugui1(j*iip1,l)=ugui1((j-1)*iip1+1,l)   
761                ugui2(j*iip1,l)=ugui2((j-1)*iip1+1,l)   
762            enddo
763            do i=1,iip1
764                ugui1(i,l)=0.
765                ugui1(ip1jm+i,l)=0.
766                ugui2(i,l)=0.
767                ugui2(ip1jm+i,l)=0.
768            enddo
769        enddo
770    ENDIF
771   
772    IF (guide_T) THEN
773        CALL pres2lev(tnat1,zu1,nlevnc,llm,plnc1,plsnc,iip1,jjp1,invert_p)
774        CALL pres2lev(tnat2,zu2,nlevnc,llm,plnc2,plsnc,iip1,jjp1,invert_p)
775        do l=1,llm
776            do j=1,jjp1
777                IF (guide_teta) THEN
778                    do i=1,iim
779                        ij=(j-1)*iip1+i
780                        tgui1(ij,l)=zu1(i,j,l)
781                        tgui2(ij,l)=zu2(i,j,l)
782                    enddo
783                ELSE
784                    do i=1,iim
785                        ij=(j-1)*iip1+i
786                        tgui1(ij,l)=zu1(i,j,l)*cpp/pk(i,j,l)
787                        tgui2(ij,l)=zu2(i,j,l)*cpp/pk(i,j,l)
788                    enddo
789                ENDIF
790                tgui1(j*iip1,l)=tgui1((j-1)*iip1+1,l)   
791                tgui2(j*iip1,l)=tgui2((j-1)*iip1+1,l)   
792            enddo
793            do i=1,iip1
794                tgui1(i,l)=tgui1(1,l)
795                tgui1(ip1jm+i,l)=tgui1(ip1jm+1,l)
796                tgui2(i,l)=tgui2(1,l)
797                tgui2(ip1jm+i,l)=tgui2(ip1jm+1,l)
798            enddo
799        enddo
800    ENDIF
801
802    IF (guide_v) THEN
803
804        CALL pres2lev(vnat1,zv1,nlevnc,llm,plnc1(:,:jjm,:),plvnc,iip1,jjm,invert_p)
805        CALL pres2lev(vnat2,zv2,nlevnc,llm,plnc2(:,:jjm,:),plvnc,iip1,jjm,invert_p)
806
807        do l=1,llm
808            do j=1,jjm
809                do i=1,iim
810                    ij=(j-1)*iip1+i
811                    vgui1(ij,l)=zv1(i,j,l)*cv(i,j)
812                    vgui2(ij,l)=zv2(i,j,l)*cv(i,j)
813                enddo
814                vgui1(j*iip1,l)=vgui1((j-1)*iip1+1,l)   
815                vgui2(j*iip1,l)=vgui2((j-1)*iip1+1,l)   
816            enddo
817        enddo
818    ENDIF
819   
820    IF (guide_Q) THEN
821        ! On suppose qu'on a la bonne variable dans le fichier de guidage:
822        ! Hum.Rel si guide_hr, Hum.Spec. sinon.
823        CALL pres2lev(qnat1,zu1,nlevnc,llm,plnc1,plsnc,iip1,jjp1,invert_p)
824        CALL pres2lev(qnat2,zu2,nlevnc,llm,plnc2,plsnc,iip1,jjp1,invert_p)
825        do l=1,llm
826            do j=1,jjp1
827                do i=1,iim
828                    ij=(j-1)*iip1+i
829                    qgui1(ij,l)=zu1(i,j,l)
830                    qgui2(ij,l)=zu2(i,j,l)
831                enddo
832                qgui1(j*iip1,l)=qgui1((j-1)*iip1+1,l)   
833                qgui2(j*iip1,l)=qgui2((j-1)*iip1+1,l)   
834            enddo
835            do i=1,iip1
836                qgui1(i,l)=qgui1(1,l)
837                qgui1(ip1jm+i,l)=qgui1(ip1jm+1,l)
838                qgui2(i,l)=qgui2(1,l)
839                qgui2(ip1jm+i,l)=qgui2(ip1jm+1,l)
840            enddo
841        enddo
842        IF (guide_hr) THEN
843            CALL q_sat(iip1*jjp1*llm,teta*pk/cpp,plsnc,qsat)
844            qgui1=qgui1*qsat*0.01 !hum. rel. en %
845            qgui2=qgui2*qsat*0.01
846        ENDIF
847    ENDIF
848
849  END SUBROUTINE guide_interp
850
851!=======================================================================
852  SUBROUTINE tau2alpha(typ,pim,pjm,factt,taumin,taumax,alpha)
853
854! Calcul des constantes de rappel alpha (=1/tau)
855
856    use comconst_mod, only: pi
857    use serre_mod, only: clon, clat, grossismx, grossismy
858   
859    implicit none
860
861    include "dimensions.h"
862    include "paramet.h"
863    include "comgeom2.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    USE comconst_mod, ONLY: pi
1519    USE comvert_mod, ONLY: presnivs
1520    use netcdf95, only: nf95_def_var, nf95_put_var
1521    use netcdf, only: nf90_float
1522   
1523    IMPLICIT NONE
1524
1525    INCLUDE "dimensions.h"
1526    INCLUDE "paramet.h"
1527    INCLUDE "netcdf.inc"
1528    INCLUDE "comgeom2.h"
1529   
1530    ! Variables entree
1531    CHARACTER*(*), INTENT(IN)                          :: varname
1532    INTEGER,   INTENT (IN)                         :: hsize,vsize
1533    REAL, DIMENSION (iip1,hsize,vsize), INTENT(IN) :: field
1534
1535    ! Variables locales
1536    INTEGER, SAVE :: timestep=0
1537    ! Identites fichier netcdf
1538    INTEGER       :: nid, id_lonu, id_lonv, id_latu, id_latv, id_tim, id_lev
1539    INTEGER       :: vid_lonu,vid_lonv,vid_latu,vid_latv,vid_cu,vid_cv,vid_lev
1540    INTEGER       :: vid_au,vid_av, varid_alpha_t, varid_alpha_q
1541    INTEGER, DIMENSION (3) :: dim3
1542    INTEGER, DIMENSION (4) :: dim4,count,start
1543    INTEGER                :: ierr, varid,l
1544    REAL, DIMENSION (iip1,hsize,vsize) :: field2
1545
1546    print *,'Guide: output timestep',timestep,'var ',varname
1547    IF (timestep.EQ.0) THEN
1548! ----------------------------------------------
1549! initialisation fichier de sortie
1550! ----------------------------------------------
1551! Ouverture du fichier
1552        ierr=NF_CREATE("guide_ins.nc",NF_CLOBBER,nid)
1553! Definition des dimensions
1554        ierr=NF_DEF_DIM(nid,"LONU",iip1,id_lonu)
1555        ierr=NF_DEF_DIM(nid,"LONV",iip1,id_lonv)
1556        ierr=NF_DEF_DIM(nid,"LATU",jjp1,id_latu)
1557        ierr=NF_DEF_DIM(nid,"LATV",jjm,id_latv)
1558        ierr=NF_DEF_DIM(nid,"LEVEL",llm,id_lev)
1559        ierr=NF_DEF_DIM(nid,"TIME",NF_UNLIMITED,id_tim)
1560
1561! Creation des variables dimensions
1562        ierr=NF_DEF_VAR(nid,"LONU",NF_FLOAT,1,id_lonu,vid_lonu)
1563        ierr=NF_DEF_VAR(nid,"LONV",NF_FLOAT,1,id_lonv,vid_lonv)
1564        ierr=NF_DEF_VAR(nid,"LATU",NF_FLOAT,1,id_latu,vid_latu)
1565        ierr=NF_DEF_VAR(nid,"LATV",NF_FLOAT,1,id_latv,vid_latv)
1566        ierr=NF_DEF_VAR(nid,"LEVEL",NF_FLOAT,1,id_lev,vid_lev)
1567        ierr=NF_DEF_VAR(nid,"cu",NF_FLOAT,2,(/id_lonu,id_latu/),vid_cu)
1568        ierr=NF_DEF_VAR(nid,"au",NF_FLOAT,2,(/id_lonu,id_latu/),vid_au)
1569        ierr=NF_DEF_VAR(nid,"cv",NF_FLOAT,2,(/id_lonv,id_latv/),vid_cv)
1570        ierr=NF_DEF_VAR(nid,"av",NF_FLOAT,2,(/id_lonv,id_latv/),vid_av)
1571        call nf95_def_var(nid, "alpha_T", nf90_float, (/id_lonv, id_latu/), &
1572             varid_alpha_t)
1573        call nf95_def_var(nid, "alpha_q", nf90_float, (/id_lonv, id_latu/), &
1574             varid_alpha_q)
1575       
1576        ierr=NF_ENDDEF(nid)
1577
1578! Enregistrement des variables dimensions
1579#ifdef NC_DOUBLE
1580        ierr = NF_PUT_VAR_DOUBLE(nid,vid_lonu,rlonu*180./pi)
1581        ierr = NF_PUT_VAR_DOUBLE(nid,vid_lonv,rlonv*180./pi)
1582        ierr = NF_PUT_VAR_DOUBLE(nid,vid_latu,rlatu*180./pi)
1583        ierr = NF_PUT_VAR_DOUBLE(nid,vid_latv,rlatv*180./pi)
1584        ierr = NF_PUT_VAR_DOUBLE(nid,vid_lev,presnivs)
1585        ierr = NF_PUT_VAR_DOUBLE(nid,vid_cu,cu)
1586        ierr = NF_PUT_VAR_DOUBLE(nid,vid_cv,cv)
1587        ierr = NF_PUT_VAR_DOUBLE(nid,vid_au,alpha_u)
1588        ierr = NF_PUT_VAR_DOUBLE(nid,vid_av,alpha_v)
1589#else
1590        ierr = NF_PUT_VAR_REAL(nid,vid_lonu,rlonu*180./pi)
1591        ierr = NF_PUT_VAR_REAL(nid,vid_lonv,rlonv*180./pi)
1592        ierr = NF_PUT_VAR_REAL(nid,vid_latu,rlatu*180./pi)
1593        ierr = NF_PUT_VAR_REAL(nid,vid_latv,rlatv*180./pi)
1594        ierr = NF_PUT_VAR_REAL(nid,vid_lev,presnivs)
1595        ierr = NF_PUT_VAR_REAL(nid,vid_cu,cu)
1596        ierr = NF_PUT_VAR_REAL(nid,vid_cv,cv)
1597        ierr = NF_PUT_VAR_REAL(nid,vid_au,alpha_u)
1598        ierr = NF_PUT_VAR_REAL(nid,vid_av,alpha_v)
1599#endif
1600        call nf95_put_var(nid, varid_alpha_t, alpha_t)
1601        call nf95_put_var(nid, varid_alpha_q, alpha_q)
1602! --------------------------------------------------------------------
1603! Cr�ation des variables sauvegard�es
1604! --------------------------------------------------------------------
1605        ierr = NF_REDEF(nid)
1606! Surface pressure (GCM)
1607        dim3=(/id_lonv,id_latu,id_tim/)
1608        ierr = NF_DEF_VAR(nid,"SP",NF_FLOAT,3,dim3,varid)
1609! Surface pressure (guidage)
1610        IF (guide_P) THEN
1611            dim3=(/id_lonv,id_latu,id_tim/)
1612            ierr = NF_DEF_VAR(nid,"ps",NF_FLOAT,3,dim3,varid)
1613        ENDIF
1614! Zonal wind
1615        IF (guide_u) THEN
1616            dim4=(/id_lonu,id_latu,id_lev,id_tim/)
1617            ierr = NF_DEF_VAR(nid,"u",NF_FLOAT,4,dim4,varid)
1618            ierr = NF_DEF_VAR(nid,"ua",NF_FLOAT,4,dim4,varid)
1619            ierr = NF_DEF_VAR(nid,"ucov",NF_FLOAT,4,dim4,varid)
1620        ENDIF
1621! Merid. wind
1622        IF (guide_v) THEN
1623            dim4=(/id_lonv,id_latv,id_lev,id_tim/)
1624            ierr = NF_DEF_VAR(nid,"v",NF_FLOAT,4,dim4,varid)
1625            ierr = NF_DEF_VAR(nid,"va",NF_FLOAT,4,dim4,varid)
1626            ierr = NF_DEF_VAR(nid,"vcov",NF_FLOAT,4,dim4,varid)
1627        ENDIF
1628! Pot. Temperature
1629        IF (guide_T) THEN
1630            dim4=(/id_lonv,id_latu,id_lev,id_tim/)
1631            ierr = NF_DEF_VAR(nid,"teta",NF_FLOAT,4,dim4,varid)
1632        ENDIF
1633! Specific Humidity
1634        IF (guide_Q) THEN
1635            dim4=(/id_lonv,id_latu,id_lev,id_tim/)
1636            ierr = NF_DEF_VAR(nid,"q",NF_FLOAT,4,dim4,varid)
1637        ENDIF
1638       
1639        ierr = NF_ENDDEF(nid)
1640        ierr = NF_CLOSE(nid)
1641    ENDIF ! timestep=0
1642
1643! --------------------------------------------------------------------
1644! Enregistrement du champ
1645! --------------------------------------------------------------------
1646    ierr=NF_OPEN("guide_ins.nc",NF_WRITE,nid)
1647
1648    IF (varname=="SP") timestep=timestep+1
1649
1650    ierr = NF_INQ_VARID(nid,varname,varid)
1651    SELECT CASE (varname)
1652    CASE ("SP","ps")
1653        start=(/1,1,timestep,0/)
1654        count=(/iip1,jjp1,1,0/)
1655    CASE ("v","va","vcov")
1656        start=(/1,1,1,timestep/)
1657        count=(/iip1,jjm,llm,1/)
1658    CASE DEFAULT
1659        start=(/1,1,1,timestep/)
1660        count=(/iip1,jjp1,llm,1/)
1661    END SELECT
1662
1663    SELECT CASE (varname)
1664    CASE("u","ua")
1665        DO l=1,llm ; field2(:,2:jjm,l)=field(:,2:jjm,l)/cu(:,2:jjm) ; ENDDO
1666        field2(:,1,:)=0. ; field2(:,jjp1,:)=0.
1667    CASE("v","va")
1668        DO l=1,llm ; field2(:,:,l)=field(:,:,l)/cv(:,:) ; ENDDO
1669    CASE DEFAULT
1670        field2=field
1671    END SELECT
1672
1673
1674#ifdef NC_DOUBLE
1675    ierr = NF_PUT_VARA_DOUBLE(nid,varid,start,count,field2)
1676#else
1677    ierr = NF_PUT_VARA_REAL(nid,varid,start,count,field2)
1678#endif
1679
1680    ierr = NF_CLOSE(nid)
1681
1682  END SUBROUTINE guide_out
1683   
1684 
1685!===========================================================================
1686  subroutine correctbid(iim,nl,x)
1687    integer iim,nl
1688    real x(iim+1,nl)
1689    integer i,l
1690    real zz
1691
1692    do l=1,nl
1693        do i=2,iim-1
1694            if(abs(x(i,l)).gt.1.e10) then
1695               zz=0.5*(x(i-1,l)+x(i+1,l))
1696              print*,'correction ',i,l,x(i,l),zz
1697               x(i,l)=zz
1698            endif
1699         enddo
1700     enddo
1701     return
1702  end subroutine correctbid
1703
1704!===========================================================================
1705END MODULE guide_mod
Note: See TracBrowser for help on using the repository browser.