source: LMDZ6/branches/DYNAMICO-conv/libf/dyn3dpar/guide_p_mod.F90 @ 3356

Last change on this file since 3356 was 3356, checked in by Laurent Fairhead, 6 years ago

First attempt at merging with trunk

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