source: LMDZ5/trunk/libf/dyn3dpar/guide_p_mod.F90 @ 2182

Last change on this file since 2182 was 2134, checked in by lguez, 10 years ago

In nudging procedures, replaced explicit Euler integration of nudged
fields by exact integration. This does not change anything if
guide_add is true, but it changes the value of alpha if guide_add is
false. We could have taken into account the variation of the nudging
field during a nudging time step. This would be a small correction. We
choose not to take it into account for the time being. Also, we add a
restriction on zonal nudging: we allow it only for a grid which is
regular in longitude. It does not seem to make sense otherwise and the
exact integration would take more programming for an irregular grid.

In the sequential version, copying the parallel versions, set
iguide_int to 1 when the input value is 0.

  • Property copyright set to
    Name of program: LMDZ
    Creation date: 1984
    Version: LMDZ5
    License: CeCILL version 2
    Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
    See the license file in the root directory
  • Property svn:executable set to *
File size: 71.9 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
71
72    IMPLICIT NONE
73 
74    INCLUDE "dimensions.h"
75    INCLUDE "paramet.h"
76    INCLUDE "netcdf.inc"
77
78    ! For grossismx:
79    include "serre.h"
80
81    INTEGER                :: error,ncidpl,rid,rcod
82    CHARACTER (len = 80)   :: abort_message
83    CHARACTER (len = 20)   :: modname = 'guide_init'
84
85! ---------------------------------------------
86! Lecture des parametres: 
87! ---------------------------------------------
88! Variables guidees
89    CALL getpar('guide_u',.true.,guide_u,'guidage de u')
90    CALL getpar('guide_v',.true.,guide_v,'guidage de v')
91    CALL getpar('guide_T',.true.,guide_T,'guidage de T')
92    CALL getpar('guide_P',.true.,guide_P,'guidage de P')
93    CALL getpar('guide_Q',.true.,guide_Q,'guidage de Q')
94    CALL getpar('guide_hr',.true.,guide_hr,'guidage de Q par H.R')
95    CALL getpar('guide_teta',.false.,guide_teta,'guidage de T par Teta')
96
97    CALL getpar('guide_add',.false.,guide_add,'for�age constant?')
98    CALL getpar('guide_zon',.false.,guide_zon,'guidage moy zonale')
99    if (guide_zon .and. abs(grossismx - 1.) > 0.01) &
100         call abort_gcm("guide_init", &
101         "zonal nudging requires grid regular in longitude", 1)
102
103!   Constantes de rappel. Unite : fraction de jour
104    CALL getpar('tau_min_u',0.02,tau_min_u,'Cste de rappel min, u')
105    CALL getpar('tau_max_u', 10.,tau_max_u,'Cste de rappel max, u')
106    CALL getpar('tau_min_v',0.02,tau_min_v,'Cste de rappel min, v')
107    CALL getpar('tau_max_v', 10.,tau_max_v,'Cste de rappel max, v')
108    CALL getpar('tau_min_T',0.02,tau_min_T,'Cste de rappel min, T')
109    CALL getpar('tau_max_T', 10.,tau_max_T,'Cste de rappel max, T')
110    CALL getpar('tau_min_Q',0.02,tau_min_Q,'Cste de rappel min, Q')
111    CALL getpar('tau_max_Q', 10.,tau_max_Q,'Cste de rappel max, Q')
112    CALL getpar('tau_min_P',0.02,tau_min_P,'Cste de rappel min, P')
113    CALL getpar('tau_max_P', 10.,tau_max_P,'Cste de rappel max, P')
114    CALL getpar('gamma4',.false.,gamma4,'Zone sans rappel elargie')
115    CALL getpar('guide_BL',.true.,guide_BL,'guidage dans C.Lim')
116   
117! Sauvegarde du for�age
118    CALL getpar('guide_sav',.false.,guide_sav,'sauvegarde guidage')
119    CALL getpar('iguide_sav',4,iguide_sav,'freq. sauvegarde guidage')
120    ! frequences f>0: fx/jour; f<0: tous les f jours; f=0: 1 seule fois.
121    IF (iguide_sav.GT.0) THEN
122       iguide_sav=day_step/iguide_sav
123    ELSE if (iguide_sav == 0) then
124       iguide_sav = huge(0)
125    ELSE
126       iguide_sav=day_step*iguide_sav
127    ENDIF
128
129! Guidage regional seulement (sinon constant ou suivant le zoom)
130    CALL getpar('guide_reg',.false.,guide_reg,'guidage regional')
131    CALL getpar('lat_min_g',-90.,lat_min_g,'Latitude mini guidage ')
132    CALL getpar('lat_max_g', 90.,lat_max_g,'Latitude maxi guidage ')
133    CALL getpar('lon_min_g',-180.,lon_min_g,'longitude mini guidage ')
134    CALL getpar('lon_max_g', 180.,lon_max_g,'longitude maxi guidage ')
135    CALL getpar('tau_lat', 5.,tau_lat,'raideur lat guide regional ')
136    CALL getpar('tau_lon', 5.,tau_lon,'raideur lon guide regional ')
137
138! Parametres pour lecture des fichiers
139    CALL getpar('iguide_read',4,iguide_read,'freq. lecture guidage')
140    CALL getpar('iguide_int',4,iguide_int,'freq. interpolation vert')
141    IF (iguide_int.EQ.0) THEN
142        iguide_int=1
143    ELSEIF (iguide_int.GT.0) THEN
144        iguide_int=day_step/iguide_int
145    ELSE
146        iguide_int=day_step*iguide_int
147    ENDIF
148    CALL getpar('guide_plevs',0,guide_plevs,'niveaux pression fichiers guidage')
149    ! Pour compatibilite avec ancienne version avec guide_modele
150    CALL getpar('guide_modele',.false.,guide_modele,'niveaux pression ap+bp*psol')
151    IF (guide_modele) THEN
152        guide_plevs=1
153    ENDIF
154    ! Fin raccord
155    CALL getpar('ini_anal',.false.,ini_anal,'Etat initial = analyse')
156    CALL getpar('guide_invertp',.true.,invert_p,'niveaux p inverses')
157    CALL getpar('guide_inverty',.true.,invert_y,'inversion N-S')
158    CALL getpar('guide_2D',.false.,guide_2D,'fichier guidage lat-P')
159
160! ---------------------------------------------
161! Determination du nombre de niveaux verticaux
162! des fichiers guidage
163! ---------------------------------------------
164    ncidpl=-99
165    if (guide_plevs.EQ.1) then
166       if (ncidpl.eq.-99) then
167          rcod=nf90_open('apbp.nc',Nf90_NOWRITe, ncidpl)
168          if (rcod.NE.NF_NOERR) THEN
169             print *,'Guide: probleme -> pas de fichier apbp.nc'
170             CALL abort_gcm(modname,abort_message,1)
171          endif
172       endif
173    elseif (guide_plevs.EQ.2) then
174       if (ncidpl.EQ.-99) then
175          rcod=nf90_open('P.nc',Nf90_NOWRITe,ncidpl)
176          if (rcod.NE.NF_NOERR) THEN
177             print *,'Guide: probleme -> pas de fichier P.nc'
178             CALL abort_gcm(modname,abort_message,1)
179          endif
180       endif
181    elseif (guide_u) then
182       if (ncidpl.eq.-99) then
183          rcod=nf90_open('u.nc',Nf90_NOWRITe,ncidpl)
184          if (rcod.NE.NF_NOERR) THEN
185             print *,'Guide: probleme -> pas de fichier u.nc'
186             CALL abort_gcm(modname,abort_message,1)
187          endif
188       endif
189    elseif (guide_v) then
190       if (ncidpl.eq.-99) then
191          rcod=nf90_open('v.nc',nf90_nowrite,ncidpl)
192          if (rcod.NE.NF_NOERR) THEN
193             print *,'Guide: probleme -> pas de fichier v.nc'
194             CALL abort_gcm(modname,abort_message,1)
195          endif
196       endif
197    elseif (guide_T) then
198       if (ncidpl.eq.-99) then
199          rcod=nf90_open('T.nc',nf90_nowrite,ncidpl)
200          if (rcod.NE.NF_NOERR) THEN
201             print *,'Guide: probleme -> pas de fichier T.nc'
202             CALL abort_gcm(modname,abort_message,1)
203          endif
204       endif
205    elseif (guide_Q) then
206       if (ncidpl.eq.-99) then
207          rcod=nf90_open('hur.nc',nf90_nowrite, ncidpl)
208          if (rcod.NE.NF_NOERR) THEN
209             print *,'Guide: probleme -> pas de fichier hur.nc'
210             CALL abort_gcm(modname,abort_message,1)
211          endif
212       endif
213    endif
214    error=NF_INQ_DIMID(ncidpl,'LEVEL',rid)
215    IF (error.NE.NF_NOERR) error=NF_INQ_DIMID(ncidpl,'PRESSURE',rid)
216    IF (error.NE.NF_NOERR) THEN
217        print *,'Guide: probleme lecture niveaux pression'
218        CALL abort_gcm(modname,abort_message,1)
219    ENDIF
220    error=NF_INQ_DIMLEN(ncidpl,rid,nlevnc)
221    print *,'Guide: nombre niveaux vert. nlevnc', nlevnc
222    rcod = nf90_close(ncidpl)
223
224! ---------------------------------------------
225! Allocation des variables
226! ---------------------------------------------
227    abort_message='pb in allocation guide'
228
229    ALLOCATE(apnc(nlevnc), stat = error)
230    IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
231    ALLOCATE(bpnc(nlevnc), stat = error)
232    IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
233    apnc=0.;bpnc=0.
234
235    ALLOCATE(alpha_pcor(llm), stat = error)
236    IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
237    ALLOCATE(alpha_u(ip1jmp1), stat = error)
238    IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
239    ALLOCATE(alpha_v(ip1jm), stat = error)
240    IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
241    ALLOCATE(alpha_T(ip1jmp1), stat = error)
242    IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
243    ALLOCATE(alpha_Q(ip1jmp1), stat = error)
244    IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
245    ALLOCATE(alpha_P(ip1jmp1), stat = error)
246    IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
247    alpha_u=0.;alpha_v=0;alpha_T=0;alpha_Q=0;alpha_P=0
248   
249    IF (guide_u) THEN
250        ALLOCATE(unat1(iip1,jjp1,nlevnc), stat = error)
251        IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
252        ALLOCATE(ugui1(ip1jmp1,llm), stat = error)
253        IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
254        ALLOCATE(unat2(iip1,jjp1,nlevnc), stat = error)
255        IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
256        ALLOCATE(ugui2(ip1jmp1,llm), stat = error)
257        IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
258        unat1=0.;unat2=0.;ugui1=0.;ugui2=0.
259    ENDIF
260
261    IF (guide_T) THEN
262        ALLOCATE(tnat1(iip1,jjp1,nlevnc), stat = error)
263        IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
264        ALLOCATE(tgui1(ip1jmp1,llm), stat = error)
265        IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
266        ALLOCATE(tnat2(iip1,jjp1,nlevnc), stat = error)
267        IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
268        ALLOCATE(tgui2(ip1jmp1,llm), stat = error)
269        IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
270        tnat1=0.;tnat2=0.;tgui1=0.;tgui2=0.
271    ENDIF
272     
273    IF (guide_Q) THEN
274        ALLOCATE(qnat1(iip1,jjp1,nlevnc), stat = error)
275        IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
276        ALLOCATE(qgui1(ip1jmp1,llm), stat = error)
277        IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
278        ALLOCATE(qnat2(iip1,jjp1,nlevnc), stat = error)
279        IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
280        ALLOCATE(qgui2(ip1jmp1,llm), stat = error)
281        IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
282        qnat1=0.;qnat2=0.;qgui1=0.;qgui2=0.
283    ENDIF
284
285    IF (guide_v) THEN
286        ALLOCATE(vnat1(iip1,jjm,nlevnc), stat = error)
287        IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
288        ALLOCATE(vgui1(ip1jm,llm), stat = error)
289        IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
290        ALLOCATE(vnat2(iip1,jjm,nlevnc), stat = error)
291        IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
292        ALLOCATE(vgui2(ip1jm,llm), stat = error)
293        IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
294        vnat1=0.;vnat2=0.;vgui1=0.;vgui2=0.
295    ENDIF
296
297    IF (guide_plevs.EQ.2) THEN
298        ALLOCATE(pnat1(iip1,jjp1,nlevnc), stat = error)
299        IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
300        ALLOCATE(pnat2(iip1,jjp1,nlevnc), stat = error)
301        IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
302        pnat1=0.;pnat2=0.;
303    ENDIF
304
305    IF (guide_P.OR.guide_plevs.EQ.1) THEN
306        ALLOCATE(psnat1(iip1,jjp1), stat = error)
307        IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
308        ALLOCATE(psnat2(iip1,jjp1), stat = error)
309        IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
310        psnat1=0.;psnat2=0.;
311    ENDIF
312    IF (guide_P) THEN
313        ALLOCATE(psgui2(ip1jmp1), stat = error)
314        IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
315        ALLOCATE(psgui1(ip1jmp1), stat = error)
316        IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
317        psgui1=0.;psgui2=0.
318    ENDIF
319
320! ---------------------------------------------
321!   Lecture du premier etat de guidage.
322! ---------------------------------------------
323    IF (guide_2D) THEN
324        CALL guide_read2D(1)
325    ELSE
326        CALL guide_read(1)
327    ENDIF
328    IF (guide_v) vnat1=vnat2
329    IF (guide_u) unat1=unat2
330    IF (guide_T) tnat1=tnat2
331    IF (guide_Q) qnat1=qnat2
332    IF (guide_plevs.EQ.2) pnat1=pnat2
333    IF (guide_P.OR.guide_plevs.EQ.1) psnat1=psnat2
334
335  END SUBROUTINE guide_init
336
337!=======================================================================
338  SUBROUTINE guide_main(itau,ucov,vcov,teta,q,masse,ps)
339    use exner_hyb_p_m, only: exner_hyb_p
340    use exner_milieu_p_m, only: exner_milieu_p
341    USE parallel_lmdz
342    USE control_mod
343   
344    IMPLICIT NONE
345 
346    INCLUDE "dimensions.h"
347    INCLUDE "paramet.h"
348    INCLUDE "comconst.h"
349    INCLUDE "comvert.h"
350
351    ! Variables entree
352    INTEGER,                       INTENT(IN)    :: itau !pas de temps
353    REAL, DIMENSION (ip1jmp1,llm), INTENT(INOUT) :: ucov,teta,q,masse
354    REAL, DIMENSION (ip1jm,llm),   INTENT(INOUT) :: vcov
355    REAL, DIMENSION (ip1jmp1),     INTENT(INOUT) :: ps
356
357    ! Variables locales
358    LOGICAL, SAVE :: first=.TRUE.
359    LOGICAL       :: f_out ! sortie guidage
360    REAL, DIMENSION (ip1jmp1,llm) :: f_add ! var aux: champ de guidage
361    ! Variables pour fonction Exner (P milieu couche)
362    REAL, DIMENSION (iip1,jjp1,llm)    :: pk
363    REAL, DIMENSION (iip1,jjp1)        :: pks   
364    REAL                               :: unskap
365    REAL, DIMENSION (ip1jmp1,llmp1)    :: p ! besoin si guide_P
366    ! Compteurs temps:
367    INTEGER, SAVE :: step_rea,count_no_rea,itau_test ! lecture guidage
368    REAL          :: ditau, dday_step
369    REAL          :: tau,reste ! position entre 2 etats de guidage
370    REAL, SAVE    :: factt ! pas de temps en fraction de jour
371   
372    INTEGER       :: i,j,l
373   
374    ijb_u=ij_begin ; ije_u=ij_end ; ijn_u=ije_u-ijb_u+1 
375    jjb_u=jj_begin ; jje_u=jj_end ; jjn_u=jje_u-jjb_u+1
376    ijb_v=ij_begin ; ije_v=ij_end ; ijn_v=ije_v-ijb_v+1   
377    jjb_v=jj_begin ; jje_v=jj_end ; jjn_v=jje_v-jjb_v+1
378    IF (pole_sud) THEN
379      ije_v=ij_end-iip1
380      jje_v=jj_end-1
381      ijn_v=ije_v-ijb_v+1
382      jjn_v=jje_v-jjb_v+1
383    ENDIF
384     
385!    CALL AllGather_Field(ucov,ip1jmp1,llm)
386!    CALL AllGather_Field(vcov,ip1jm,llm)
387!    CALL AllGather_Field(teta,ip1jmp1,llm)
388!    CALL AllGather_Field(ps,ip1jmp1,1)
389!    CALL AllGather_Field(q,ip1jmp1,llm)
390   
391!-----------------------------------------------------------------------
392! Initialisations au premier passage
393!-----------------------------------------------------------------------
394
395    IF (first) THEN
396        first=.FALSE.
397        CALL guide_init
398        itau_test=1001
399        step_rea=1
400        count_no_rea=0
401! Calcul des constantes de rappel
402        factt=dtvr*iperiod/daysec
403        call tau2alpha(3,iip1,jjm ,factt,tau_min_v,tau_max_v,alpha_v)
404        call tau2alpha(2,iip1,jjp1,factt,tau_min_u,tau_max_u,alpha_u)
405        call tau2alpha(1,iip1,jjp1,factt,tau_min_T,tau_max_T,alpha_T)
406        call tau2alpha(1,iip1,jjp1,factt,tau_min_P,tau_max_P,alpha_P)
407        call tau2alpha(1,iip1,jjp1,factt,tau_min_Q,tau_max_Q,alpha_Q)
408! correction de rappel dans couche limite
409        if (guide_BL) then
410             alpha_pcor(:)=1.
411        else
412            do l=1,llm
413                alpha_pcor(l)=(1.+tanh((0.85-presnivs(l)/preff)/0.05))/2.
414            enddo
415        endif
416! ini_anal: etat initial egal au guidage       
417        IF (ini_anal) THEN
418            CALL guide_interp(ps,teta)
419            IF (guide_u) ucov(ijb_u:ije_u,:)=ugui2(ijb_u:ije_u,:)
420            IF (guide_v) vcov(ijb_v:ije_v,:)=ugui2(ijb_v:ije_v,:)
421            IF (guide_T) teta(ijb_u:ije_u,:)=tgui2(ijb_u:ije_u,:)
422            IF (guide_Q) q(ijb_u:ije_u,:)=qgui2(ijb_u:ije_u,:)
423            IF (guide_P) THEN
424                ps(ijb_u:ije_u)=psgui2(ijb_u:ije_u)
425                CALL pression_p(ip1jmp1,ap,bp,ps,p)
426                CALL massdair_p(p,masse)
427            ENDIF
428            RETURN
429        ENDIF
430! Verification structure guidage
431        IF (guide_u) THEN
432            CALL writefield_p('unat',unat1)
433            CALL writefield_p('ucov',RESHAPE(ucov,(/iip1,jjp1,llm/)))
434        ENDIF
435        IF (guide_T) THEN
436            CALL writefield_p('tnat',tnat1)
437            CALL writefield_p('teta',RESHAPE(teta,(/iip1,jjp1,llm/)))
438        ENDIF
439
440    ENDIF !first
441
442!-----------------------------------------------------------------------
443! Lecture des fichiers de guidage ?
444!-----------------------------------------------------------------------
445    IF (iguide_read.NE.0) THEN
446      ditau=real(itau)
447      dday_step=real(day_step)
448      IF (iguide_read.LT.0) THEN
449          tau=ditau/dday_step/REAL(iguide_read)
450      ELSE
451          tau=REAL(iguide_read)*ditau/dday_step
452      ENDIF
453      reste=tau-AINT(tau)
454      IF (reste.EQ.0.) THEN
455          IF (itau_test.EQ.itau) THEN
456              write(*,*)'deuxieme passage de advreel a itau=',itau
457              stop
458          ELSE
459              IF (guide_v) vnat1(:,jjb_v:jje_v,:)=vnat2(:,jjb_v:jje_v,:)
460              IF (guide_u) unat1(:,jjb_u:jje_u,:)=unat2(:,jjb_u:jje_u,:)
461              IF (guide_T) tnat1(:,jjb_u:jje_u,:)=tnat2(:,jjb_u:jje_u,:)
462              IF (guide_Q) qnat1(:,jjb_u:jje_u,:)=qnat2(:,jjb_u:jje_u,:)
463              IF (guide_plevs.EQ.2) pnat1(:,jjb_u:jje_u,:)=pnat2(:,jjb_u:jje_u,:)
464              IF (guide_P.OR.guide_plevs.EQ.1) psnat1(:,jjb_u:jje_u)=psnat2(:,jjb_u:jje_u)
465              step_rea=step_rea+1
466              itau_test=itau
467              print*,'Lecture fichiers guidage, pas ',step_rea, &
468                    'apres ',count_no_rea,' non lectures'
469              IF (guide_2D) THEN
470                  CALL guide_read2D(step_rea)
471              ELSE
472                  CALL guide_read(step_rea)
473              ENDIF
474              count_no_rea=0
475          ENDIF
476      ELSE
477        count_no_rea=count_no_rea+1
478
479      ENDIF
480    ENDIF !iguide_read=0
481
482!-----------------------------------------------------------------------
483! Interpolation et conversion des champs de guidage
484!-----------------------------------------------------------------------
485    IF (MOD(itau,iguide_int).EQ.0) THEN
486        CALL guide_interp(ps,teta)
487    ENDIF
488! Repartition entre 2 etats de guidage
489    IF (iguide_read.NE.0) THEN
490        tau=reste
491    ELSE
492        tau=1.
493    ENDIF
494
495!-----------------------------------------------------------------------
496!   Ajout des champs de guidage
497!-----------------------------------------------------------------------
498! Sauvegarde du guidage?
499    f_out=((MOD(itau,iguide_sav).EQ.0).AND.guide_sav) 
500    IF (f_out) THEN
501!       Calcul niveaux pression milieu de couches
502        CALL pression_p( ip1jmp1, ap, bp, ps, p )
503        if (pressure_exner) then
504          CALL exner_hyb_p(ip1jmp1,ps,p,pks,pk)
505        else
506          CALL exner_milieu_p(ip1jmp1,ps,p,pks,pk)
507        endif
508        unskap=1./kappa
509        DO l = 1, llm
510           DO j=jjb_u,jje_u
511              DO i =1, iip1
512                 p(i+(j-1)*iip1,l) = preff * ( pk(i,j,l)/cpp) ** unskap
513              ENDDO
514           ENDDO
515        ENDDO
516        CALL guide_out("SP",jjp1,llm,p(1:ip1jmp1,1:llm),1.)
517    ENDIF
518   
519    if (guide_u) then
520        if (guide_add) then
521           f_add(ijb_u:ije_u,:)=(1.-tau)*ugui1(ijb_u:ije_u,:)+tau*ugui2(ijb_u:ije_u,:)
522        else
523           f_add(ijb_u:ije_u,:)=(1.-tau)*ugui1(ijb_u:ije_u,:)+tau*ugui2(ijb_u:ije_u,:)-ucov(ijb_u:ije_u,:)
524        endif
525
526        if (guide_zon) CALL guide_zonave(1,jjp1,llm,f_add)
527        CALL guide_addfield(ip1jmp1,llm,f_add,alpha_u)
528        IF (f_out) CALL guide_out("u",jjp1,llm,ucov,factt)
529        IF (f_out) CALL guide_out("ua",jjp1,llm,(1.-tau)*ugui1(:,:)+tau*ugui2(:,:),factt)
530        IF (f_out) CALL guide_out("ucov",jjp1,llm,f_add(:,:)/factt,factt)
531        ucov(ijb_u:ije_u,:)=ucov(ijb_u:ije_u,:)+f_add(ijb_u:ije_u,:)
532    endif
533
534    if (guide_T) then
535        if (guide_add) then
536           f_add(ijb_u:ije_u,:)=(1.-tau)*tgui1(ijb_u:ije_u,:)+tau*tgui2(ijb_u:ije_u,:)
537        else
538           f_add(ijb_u:ije_u,:)=(1.-tau)*tgui1(ijb_u:ije_u,:)+tau*tgui2(ijb_u:ije_u,:)-teta(ijb_u:ije_u,:)
539        endif
540        if (guide_zon) CALL guide_zonave(2,jjp1,llm,f_add)
541        CALL guide_addfield(ip1jmp1,llm,f_add,alpha_T)
542        IF (f_out) CALL guide_out("teta",jjp1,llm,f_add(:,:)/factt,factt)
543        teta(ijb_u:ije_u,:)=teta(ijb_u:ije_u,:)+f_add(ijb_u:ije_u,:)
544    endif
545
546    if (guide_P) then
547        if (guide_add) then
548           f_add(ijb_u:ije_u,1)=(1.-tau)*psgui1(ijb_u:ije_u)+tau*psgui2(ijb_u:ije_u)
549        else
550           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)
551        endif
552        if (guide_zon) CALL guide_zonave(2,jjp1,1,f_add(1:ip1jmp1,1))
553        CALL guide_addfield(ip1jmp1,1,f_add(1:ip1jmp1,1),alpha_P)
554        IF (f_out) CALL guide_out("ps",jjp1,1,f_add(1:ip1jmp1,1)/factt,factt)
555        ps(ijb_u:ije_u)=ps(ijb_u:ije_u)+f_add(ijb_u:ije_u,1)
556        CALL pression_p(ip1jmp1,ap,bp,ps,p)
557        CALL massdair_p(p,masse)
558    endif
559
560    if (guide_Q) then
561        if (guide_add) then
562           f_add(ijb_u:ije_u,:)=(1.-tau)*qgui1(ijb_u:ije_u,:)+tau*qgui2(ijb_u:ije_u,:)
563        else
564           f_add(ijb_u:ije_u,:)=(1.-tau)*qgui1(ijb_u:ije_u,:)+tau*qgui2(ijb_u:ije_u,:)-q(ijb_u:ije_u,:)
565        endif
566        if (guide_zon) CALL guide_zonave(2,jjp1,llm,f_add)
567        CALL guide_addfield(ip1jmp1,llm,f_add,alpha_Q)
568        IF (f_out) CALL guide_out("q",jjp1,llm,f_add(:,:)/factt,factt)
569        q(ijb_u:ije_u,:)=q(ijb_u:ije_u,:)+f_add(ijb_u:ije_u,:)
570    endif
571
572    if (guide_v) then
573        if (guide_add) then
574           f_add(ijb_v:ije_v,:)=(1.-tau)*vgui1(ijb_v:ije_v,:)+tau*vgui2(ijb_v:ije_v,:)
575        else
576           f_add(ijb_v:ije_v,:)=(1.-tau)*vgui1(ijb_v:ije_v,:)+tau*vgui2(ijb_v:ije_v,:)-vcov(ijb_v:ije_v,:)
577        endif
578       
579        if (guide_zon) CALL guide_zonave(2,jjm,llm,f_add(1:ip1jm,:))
580        CALL guide_addfield(ip1jm,llm,f_add(1:ip1jm,:),alpha_v)
581        IF (f_out) CALL guide_out("v",jjm,llm,vcov(1:ip1jm,:),factt)
582        IF (f_out) CALL guide_out("va",jjm,llm,(1.-tau)*vgui1(:,:)+tau*vgui2(:,:),factt)
583        IF (f_out) CALL guide_out("vcov",jjm,llm,f_add(1:ip1jm,:)/factt,factt)
584        vcov(ijb_v:ije_v,:)=vcov(ijb_v:ije_v,:)+f_add(ijb_v:ije_v,:)
585    endif
586
587  END SUBROUTINE guide_main
588
589!=======================================================================
590  SUBROUTINE guide_addfield(hsize,vsize,field,alpha)
591! field1=a*field1+alpha*field2
592
593    IMPLICIT NONE
594    INCLUDE "dimensions.h"
595    INCLUDE "paramet.h"
596
597    ! input variables
598    INTEGER,                      INTENT(IN)    :: hsize
599    INTEGER,                      INTENT(IN)    :: vsize
600    REAL, DIMENSION(hsize),       INTENT(IN)    :: alpha
601    REAL, DIMENSION(hsize,vsize), INTENT(INOUT) :: field
602
603    ! Local variables
604    INTEGER :: l
605
606    IF (hsize==ip1jm) THEN
607      do l=1,vsize
608        field(ijb_v:ije_v,l)=alpha(ijb_v:ije_v)*field(ijb_v:ije_v,l)*alpha_pcor(l)
609      enddo
610    ELSE
611      do l=1,vsize
612        field(ijb_u:ije_u,l)=alpha(ijb_u:ije_u)*field(ijb_u:ije_u,l)*alpha_pcor(l)
613      enddo
614    ENDIF   
615
616  END SUBROUTINE guide_addfield
617
618!=======================================================================
619  SUBROUTINE guide_zonave(typ,hsize,vsize,field)
620
621    IMPLICIT NONE
622
623    INCLUDE "dimensions.h"
624    INCLUDE "paramet.h"
625    INCLUDE "comgeom.h"
626    INCLUDE "comconst.h"
627   
628    ! input/output variables
629    INTEGER,                           INTENT(IN)    :: typ
630    INTEGER,                           INTENT(IN)    :: vsize
631    INTEGER,                           INTENT(IN)    :: hsize
632    REAL, DIMENSION(hsize*iip1,vsize), INTENT(INOUT) :: field
633
634    ! Local variables
635    LOGICAL, SAVE                :: first=.TRUE.
636    INTEGER, DIMENSION (2), SAVE :: imin, imax ! averaging domain
637    INTEGER                      :: i,j,l,ij
638    REAL, DIMENSION (iip1)       :: lond       ! longitude in Deg.
639    REAL, DIMENSION (hsize,vsize):: fieldm     ! zon-averaged field
640
641    IF (first) THEN
642        first=.FALSE.
643!Compute domain for averaging
644        lond=rlonu*180./pi
645        imin(1)=1;imax(1)=iip1;
646        imin(2)=1;imax(2)=iip1;
647        IF (guide_reg) THEN
648            DO i=1,iim
649                IF (lond(i).LT.lon_min_g) imin(1)=i
650                IF (lond(i).LE.lon_max_g) imax(1)=i
651            ENDDO
652            lond=rlonv*180./pi
653            DO i=1,iim
654                IF (lond(i).LT.lon_min_g) imin(2)=i
655                IF (lond(i).LE.lon_max_g) imax(2)=i
656            ENDDO
657        ENDIF
658    ENDIF
659
660    fieldm=0.
661   
662    IF (hsize==jjm) THEN
663      DO l=1,vsize
664      ! Compute zonal average
665          DO j=jjb_v,jje_v
666              DO i=imin(typ),imax(typ)
667                  ij=(j-1)*iip1+i
668                  fieldm(j,l)=fieldm(j,l)+field(ij,l)
669              ENDDO
670          ENDDO
671          fieldm(:,l)=fieldm(:,l)/REAL(imax(typ)-imin(typ)+1)
672    ! Compute forcing
673          DO j=jjb_v,jje_v
674              DO i=1,iip1
675                  ij=(j-1)*iip1+i
676                  field(ij,l)=fieldm(j,l)
677              ENDDO
678          ENDDO
679      ENDDO
680    ELSE
681      DO l=1,vsize
682      ! Compute zonal average
683          DO j=jjb_v,jje_v
684              DO i=imin(typ),imax(typ)
685                  ij=(j-1)*iip1+i
686                  fieldm(j,l)=fieldm(j,l)+field(ij,l)
687              ENDDO
688          ENDDO
689          fieldm(:,l)=fieldm(:,l)/REAL(imax(typ)-imin(typ)+1)
690    ! Compute forcing
691          DO j=jjb_u,jje_u
692              DO i=1,iip1
693                  ij=(j-1)*iip1+i
694                  field(ij,l)=fieldm(j,l)
695              ENDDO
696          ENDDO
697      ENDDO
698    ENDIF   
699
700  END SUBROUTINE guide_zonave
701
702!=======================================================================
703  SUBROUTINE guide_interp(psi,teta)
704    use exner_hyb_p_m, only: exner_hyb_p
705    use exner_milieu_p_m, only: exner_milieu_p
706  USE parallel_lmdz
707  USE mod_hallo
708  USE Bands
709  IMPLICIT NONE
710
711  include "dimensions.h"
712  include "paramet.h"
713  include "comvert.h"
714  include "comgeom2.h"
715  include "comconst.h"
716
717  REAL, DIMENSION (iip1,jjp1),     INTENT(IN) :: psi ! Psol gcm
718  REAL, DIMENSION (iip1,jjp1,llm), INTENT(IN) :: teta ! Temp. Pot. gcm
719
720  LOGICAL, SAVE                      :: first=.TRUE.
721  ! Variables pour niveaux pression:
722  REAL, DIMENSION (iip1,jjp1,nlevnc) :: plnc1,plnc2 !niveaux pression guidage
723  REAL, DIMENSION (iip1,jjp1,llm)    :: plunc,plsnc !niveaux pression modele
724  REAL, DIMENSION (iip1,jjm,llm)     :: plvnc       !niveaux pression modele
725  REAL, DIMENSION (iip1,jjp1,llmp1)  :: p           ! pression intercouches
726  REAL, DIMENSION (iip1,jjp1,llm)    :: pls, pext   ! var intermediaire
727  REAL, DIMENSION (iip1,jjp1,llm)    :: pbarx
728  REAL, DIMENSION (iip1,jjm,llm)     :: pbary
729  ! Variables pour fonction Exner (P milieu couche)
730  REAL, DIMENSION (iip1,jjp1,llm)    :: pk
731  REAL, DIMENSION (iip1,jjp1)        :: pks   
732  REAL                               :: unskap
733  ! Pression de vapeur saturante
734  REAL, DIMENSION (ip1jmp1,llm)      :: qsat
735  !Variables intermediaires interpolation
736  REAL, DIMENSION (iip1,jjp1,llm)    :: zu1,zu2
737  REAL, DIMENSION (iip1,jjm,llm)     :: zv1,zv2
738 
739  INTEGER                            :: i,j,l,ij
740  TYPE(Request) :: Req 
741
742    print *,'Guide: conversion variables guidage'
743! -----------------------------------------------------------------
744! Calcul des niveaux de pression champs guidage (pour T et Q)
745! -----------------------------------------------------------------
746    IF (guide_plevs.EQ.0) THEN
747        DO l=1,nlevnc
748            DO j=jjb_u,jje_u
749                DO i=1,iip1
750                    plnc2(i,j,l)=apnc(l)
751                    plnc1(i,j,l)=apnc(l)
752               ENDDO
753            ENDDO
754        ENDDO
755    ENDIF   
756
757    if (first) then
758        first=.FALSE.
759        print*,'Guide: verification ordre niveaux verticaux'
760        print*,'LMDZ :'
761        do l=1,llm
762            print*,'PL(',l,')=',(ap(l)+ap(l+1))/2. &
763                  +psi(1,jje_u)*(bp(l)+bp(l+1))/2.
764        enddo
765        print*,'Fichiers guidage'
766        SELECT CASE (guide_plevs)
767        CASE (0)
768            do l=1,nlevnc
769                 print*,'PL(',l,')=',plnc2(1,jjb_u,l)
770            enddo
771        CASE (1)
772            DO l=1,nlevnc
773                 print*,'PL(',l,')=',apnc(l)+bpnc(l)*psnat2(i,jjb_u)
774             ENDDO
775        CASE (2)
776            do l=1,nlevnc
777                 print*,'PL(',l,')=',pnat2(1,jjb_u,l)
778            enddo
779        END SELECT
780        print *,'inversion de l''ordre: invert_p=',invert_p
781        if (guide_u) then
782            do l=1,nlevnc
783                print*,'U(',l,')=',unat2(1,jjb_u,l)
784            enddo
785        endif
786        if (guide_T) then
787            do l=1,nlevnc
788                print*,'T(',l,')=',tnat2(1,jjb_u,l)
789            enddo
790        endif
791    endif
792   
793! -----------------------------------------------------------------
794! Calcul niveaux pression modele
795! -----------------------------------------------------------------
796
797!    ....  Calcul de pls , pression au milieu des couches ,en Pascals
798    IF (guide_plevs.EQ.1) THEN
799        DO l=1,llm
800            DO j=jjb_u,jje_u
801                DO i =1, iip1
802                    pls(i,j,l)=(ap(l)+ap(l+1))/2.+psi(i,j)*(bp(l)+bp(l+1))/2.
803                ENDDO
804            ENDDO
805        ENDDO
806    ELSE
807        CALL pression_p( ip1jmp1, ap, bp, psi, p )
808        if (pressure_exner) then
809          CALL exner_hyb_p(ip1jmp1,psi,p,pks,pk)
810        else
811          CALL exner_milieu_p(ip1jmp1,psi,p,pks,pk)
812        endif
813        unskap=1./kappa
814        DO l = 1, llm
815            DO j=jjb_u,jje_u
816                DO i =1, iip1
817                    pls(i,j,l) = preff * ( pk(i,j,l)/cpp) ** unskap
818                ENDDO
819            ENDDO
820        ENDDO
821    ENDIF
822
823!   calcul des pressions pour les grilles u et v
824    do l=1,llm
825        do j=jjb_u,jje_u
826            do i=1,iip1
827                pext(i,j,l)=pls(i,j,l)*aire(i,j)
828            enddo
829        enddo
830    enddo
831
832     CALL Register_SwapFieldHallo(pext,pext,ip1jmp1,llm,jj_Nb_caldyn,1,2,Req)
833     CALL SendRequest(Req)
834     CALL WaitRequest(Req)
835
836     call massbar_p(pext, pbarx, pbary )
837    do l=1,llm
838        do j=jjb_u,jje_u
839            do i=1,iip1
840                plunc(i,j,l)=pbarx(i,j,l)/aireu(i,j)
841                plsnc(i,j,l)=pls(i,j,l)
842            enddo
843        enddo
844    enddo
845    do l=1,llm
846        do j=jjb_v,jje_v
847            do i=1,iip1
848                plvnc(i,j,l)=pbary(i,j,l)/airev(i,j)
849            enddo
850        enddo
851    enddo
852
853! -----------------------------------------------------------------
854! Interpolation verticale champs guidage sur niveaux modele
855! Conversion en variables gcm (ucov, vcov...)
856! -----------------------------------------------------------------
857    if (guide_P) then
858        do j=jjb_u,jje_u
859            do i=1,iim
860                ij=(j-1)*iip1+i
861                psgui1(ij)=psnat1(i,j)
862                psgui2(ij)=psnat2(i,j)
863            enddo
864            psgui1(iip1*j)=psnat1(1,j)
865            psgui2(iip1*j)=psnat2(1,j)
866        enddo
867    endif
868
869    IF (guide_T) THEN
870        ! Calcul des nouvelles valeurs des niveaux de pression du guidage
871        IF (guide_plevs.EQ.1) THEN
872            DO l=1,nlevnc
873                DO j=jjb_u,jje_u
874                    DO i=1,iip1
875                        plnc2(i,j,l)=apnc(l)+bpnc(l)*psnat2(i,j)
876                        plnc1(i,j,l)=apnc(l)+bpnc(l)*psnat1(i,j)
877                    ENDDO
878                ENDDO
879            ENDDO
880        ELSE IF (guide_plevs.EQ.2) THEN
881            DO l=1,nlevnc
882                DO j=jjb_u,jje_u
883                    DO i=1,iip1
884                        plnc2(i,j,l)=pnat2(i,j,l)
885                        plnc1(i,j,l)=pnat1(i,j,l)
886                    ENDDO
887                ENDDO
888            ENDDO
889        ENDIF
890
891        ! Interpolation verticale
892        CALL pres2lev(tnat1(:,jjb_u:jje_u,:),zu1(:,jjb_u:jje_u,:),nlevnc,llm,           &
893                    plnc1(:,jjb_u:jje_u,:),plsnc(:,jjb_u:jje_u,:),iip1,jjn_u,invert_p)
894        CALL pres2lev(tnat2(:,jjb_u:jje_u,:),zu2(:,jjb_u:jje_u,:),nlevnc,llm,           &
895                    plnc2(:,jjb_u:jje_u,:),plsnc(:,jjb_u:jje_u,:),iip1,jjn_u,invert_p)
896
897        ! Conversion en variables GCM
898        do l=1,llm
899            do j=jjb_u,jje_u
900                IF (guide_teta) THEN
901                    do i=1,iim
902                        ij=(j-1)*iip1+i
903                        tgui1(ij,l)=zu1(i,j,l)
904                        tgui2(ij,l)=zu2(i,j,l)
905                    enddo
906                ELSE
907                    do i=1,iim
908                        ij=(j-1)*iip1+i
909                        tgui1(ij,l)=zu1(i,j,l)*cpp/pk(i,j,l)
910                        tgui2(ij,l)=zu2(i,j,l)*cpp/pk(i,j,l)
911                    enddo
912                ENDIF
913                tgui1(j*iip1,l)=tgui1((j-1)*iip1+1,l)   
914                tgui2(j*iip1,l)=tgui2((j-1)*iip1+1,l)   
915            enddo
916            do i=1,iip1
917                tgui1(i,l)=tgui1(1,l)
918                tgui1(ip1jm+i,l)=tgui1(ip1jm+1,l)
919                tgui2(i,l)=tgui2(1,l)
920                tgui2(ip1jm+i,l)=tgui2(ip1jm+1,l)
921            enddo
922        enddo
923    ENDIF
924
925    IF (guide_Q) THEN
926        ! Calcul des nouvelles valeurs des niveaux de pression du guidage
927        IF (guide_plevs.EQ.1) THEN
928            DO l=1,nlevnc
929                DO j=jjb_u,jje_u
930                    DO i=1,iip1
931                        plnc2(i,j,l)=apnc(l)+bpnc(l)*psnat2(i,j)
932                        plnc1(i,j,l)=apnc(l)+bpnc(l)*psnat1(i,j)
933                    ENDDO
934                ENDDO
935            ENDDO
936        ELSE IF (guide_plevs.EQ.2) THEN
937            DO l=1,nlevnc
938                DO j=jjb_u,jje_u
939                    DO i=1,iip1
940                        plnc2(i,j,l)=pnat2(i,j,l)
941                        plnc1(i,j,l)=pnat1(i,j,l)
942                    ENDDO
943                ENDDO
944            ENDDO
945        ENDIF
946
947        ! Interpolation verticale
948        CALL pres2lev(qnat1(:,jjb_u:jje_u,:),zu1(:,jjb_u:jje_u,:),nlevnc,llm,             &
949                      plnc1(:,jjb_u:jje_u,:),plsnc(:,jjb_u:jje_u,:),iip1,jjn_u,invert_p)
950        CALL pres2lev(qnat2(:,jjb_u:jje_u,:),zu2(:,jjb_u:jje_u,:),nlevnc,llm,             &
951                      plnc2(:,jjb_u:jje_u,:),plsnc(:,jjb_u:jje_u,:),iip1,jjn_u,invert_p)
952
953        ! Conversion en variables GCM
954        ! On suppose qu'on a la bonne variable dans le fichier de guidage:
955        ! Hum.Rel si guide_hr, Hum.Spec. sinon.
956        do l=1,llm
957            do j=jjb_u,jje_u
958                do i=1,iim
959                    ij=(j-1)*iip1+i
960                    qgui1(ij,l)=zu1(i,j,l)
961                    qgui2(ij,l)=zu2(i,j,l)
962                enddo
963                qgui1(j*iip1,l)=qgui1((j-1)*iip1+1,l)   
964                qgui2(j*iip1,l)=qgui2((j-1)*iip1+1,l)   
965            enddo
966            do i=1,iip1
967                qgui1(i,l)=qgui1(1,l)
968                qgui1(ip1jm+i,l)=qgui1(ip1jm+1,l)
969                qgui2(i,l)=qgui2(1,l)
970                qgui2(ip1jm+i,l)=qgui2(ip1jm+1,l)
971            enddo
972        enddo
973        IF (guide_hr) THEN
974            CALL q_sat(iip1*jjn_u*llm,teta(:,jjb_u:jje_u,:)*pk(:,jjb_u:jje_u,:)/cpp,       &
975                       plsnc(:,jjb_u:jje_u,:),qsat(ijb_u:ije_u,:))
976            qgui1(ijb_u:ije_u,:)=qgui1(ijb_u:ije_u,:)*qsat(ijb_u:ije_u,:)*0.01 !hum. rel. en %
977            qgui2(ijb_u:ije_u,:)=qgui2(ijb_u:ije_u,:)*qsat(ijb_u:ije_u,:)*0.01
978        ENDIF
979    ENDIF
980
981    IF (guide_u) THEN
982        ! Calcul des nouvelles valeurs des niveaux de pression du guidage
983        IF (guide_plevs.EQ.1) THEN
984            DO l=1,nlevnc
985                DO j=jjb_u,jje_u
986                    DO i=1,iim
987                        plnc2(i,j,l)=apnc(l)+bpnc(l)*(psnat2(i,j)*aire(i,j)*alpha1p2(i,j) &
988                       &           +psnat2(i+1,j)*aire(i+1,j)*alpha3p4(i+1,j))/aireu(i,j)
989                        plnc1(i,j,l)=apnc(l)+bpnc(l)*(psnat1(i,j)*aire(i,j)*alpha1p2(i,j) &
990                       &           +psnat1(i+1,j)*aire(i+1,j)*alpha3p4(i+1,j))/aireu(i,j)
991                    ENDDO
992                    plnc2(iip1,j,l)=plnc2(1,j,l)
993                    plnc1(iip1,j,l)=plnc1(1,j,l)
994                ENDDO
995            ENDDO
996        ELSE IF (guide_plevs.EQ.2) THEN
997            DO l=1,nlevnc
998                DO j=jjb_u,jje_u
999                    DO i=1,iim
1000                        plnc2(i,j,l)=(pnat2(i,j,l)*aire(i,j)*alpha1p2(i,j) &
1001                       & +pnat2(i+1,j,l)*aire(i,j)*alpha3p4(i+1,j))/aireu(i,j)
1002                        plnc1(i,j,l)=(pnat1(i,j,l)*aire(i,j)*alpha1p2(i,j) &
1003                       & +pnat1(i+1,j,l)*aire(i,j)*alpha3p4(i+1,j))/aireu(i,j)
1004                    ENDDO
1005                    plnc2(iip1,j,l)=plnc2(1,j,l)
1006                    plnc1(iip1,j,l)=plnc1(1,j,l)
1007                ENDDO
1008            ENDDO
1009        ENDIF
1010       
1011        ! Interpolation verticale
1012        CALL pres2lev(unat1(:,jjb_u:jje_u,:),zu1(:,jjb_u:jje_u,:),nlevnc,llm,            &
1013                      plnc1(:,jjb_u:jje_u,:),plunc(:,jjb_u:jje_u,:),iip1,jjn_u,invert_p)
1014        CALL pres2lev(unat2(:,jjb_u:jje_u,:),zu2(:,jjb_u:jje_u,:),nlevnc,llm,            &
1015                      plnc2(:,jjb_u:jje_u,:),plunc(:,jjb_u:jje_u,:),iip1,jjn_u,invert_p)
1016
1017        ! Conversion en variables GCM
1018        do l=1,llm
1019            do j=jjb_u,jje_u
1020                do i=1,iim
1021                    ij=(j-1)*iip1+i
1022                    ugui1(ij,l)=zu1(i,j,l)*cu(i,j)
1023                    ugui2(ij,l)=zu2(i,j,l)*cu(i,j)
1024                enddo
1025                ugui1(j*iip1,l)=ugui1((j-1)*iip1+1,l)   
1026                ugui2(j*iip1,l)=ugui2((j-1)*iip1+1,l)   
1027            enddo
1028            do i=1,iip1
1029                ugui1(i,l)=0.
1030                ugui1(ip1jm+i,l)=0.
1031                ugui2(i,l)=0.
1032                ugui2(ip1jm+i,l)=0.
1033            enddo
1034        enddo
1035    ENDIF
1036   
1037    IF (guide_v) THEN
1038        ! Calcul des nouvelles valeurs des niveaux de pression du guidage
1039        IF (guide_plevs.EQ.1) THEN
1040         CALL Register_SwapFieldHallo(psnat1,psnat1,ip1jmp1,1,jj_Nb_caldyn,1,2,Req)
1041         CALL SendRequest(Req)
1042         CALL WaitRequest(Req)
1043         CALL Register_SwapFieldHallo(psnat2,psnat2,ip1jmp1,1,jj_Nb_caldyn,1,2,Req)
1044         CALL SendRequest(Req)
1045         CALL WaitRequest(Req)
1046            DO l=1,nlevnc
1047                DO j=jjb_v,jje_v
1048                    DO i=1,iip1
1049                        plnc2(i,j,l)=apnc(l)+bpnc(l)*(psnat2(i,j)*aire(i,j)*alpha2p3(i,j) &
1050                       &           +psnat2(i,j+1)*aire(i,j+1)*alpha1p4(i,j+1))/airev(i,j)
1051                        plnc1(i,j,l)=apnc(l)+bpnc(l)*(psnat1(i,j)*aire(i,j)*alpha2p3(i,j) &
1052                       &           +psnat1(i,j+1)*aire(i,j+1)*alpha1p4(i,j+1))/airev(i,j)
1053                    ENDDO
1054                ENDDO
1055            ENDDO
1056        ELSE IF (guide_plevs.EQ.2) THEN
1057         CALL Register_SwapFieldHallo(pnat1,pnat1,ip1jmp1,llm,jj_Nb_caldyn,1,2,Req)
1058         CALL SendRequest(Req)
1059         CALL WaitRequest(Req)
1060         CALL Register_SwapFieldHallo(pnat2,pnat2,ip1jmp1,llm,jj_Nb_caldyn,1,2,Req)
1061         CALL SendRequest(Req)
1062         CALL WaitRequest(Req)
1063            DO l=1,nlevnc
1064                DO j=jjb_v,jje_v
1065                    DO i=1,iip1
1066                        plnc2(i,j,l)=(pnat2(i,j,l)*aire(i,j)*alpha2p3(i,j) &
1067                       & +pnat2(i,j+1,l)*aire(i,j)*alpha1p4(i,j+1))/airev(i,j)
1068                        plnc1(i,j,l)=(pnat1(i,j,l)*aire(i,j)*alpha2p3(i,j) &
1069                       & +pnat1(i,j+1,l)*aire(i,j)*alpha1p4(i,j+1))/airev(i,j)
1070                    ENDDO
1071                ENDDO
1072            ENDDO
1073        ENDIF
1074        ! Interpolation verticale
1075        CALL pres2lev(vnat1(:,jjb_v:jje_v,:),zv1(:,jjb_v:jje_v,:),nlevnc,llm,             &
1076                      plnc1(:,jjb_v:jje_v,:),plvnc(:,jjb_v:jje_v,:),iip1,jjn_v,invert_p)
1077        CALL pres2lev(vnat2(:,jjb_v:jje_v,:),zv2(:,jjb_v:jje_v,:),nlevnc,llm,             &
1078                      plnc2(:,jjb_v:jje_v,:),plvnc(:,jjb_v:jje_v,:),iip1,jjn_v,invert_p)
1079        ! Conversion en variables GCM
1080        do l=1,llm
1081            do j=jjb_v,jje_v
1082                do i=1,iim
1083                    ij=(j-1)*iip1+i
1084                    vgui1(ij,l)=zv1(i,j,l)*cv(i,j)
1085                    vgui2(ij,l)=zv2(i,j,l)*cv(i,j)
1086                enddo
1087                vgui1(j*iip1,l)=vgui1((j-1)*iip1+1,l)   
1088                vgui2(j*iip1,l)=vgui2((j-1)*iip1+1,l)   
1089            enddo
1090        enddo
1091    ENDIF
1092   
1093
1094  END SUBROUTINE guide_interp
1095
1096!=======================================================================
1097  SUBROUTINE tau2alpha(typ,pim,pjm,factt,taumin,taumax,alpha)
1098
1099! Calcul des constantes de rappel alpha (=1/tau)
1100
1101    implicit none
1102
1103    include "dimensions.h"
1104    include "paramet.h"
1105    include "comconst.h"
1106    include "comgeom2.h"
1107    include "serre.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    IMPLICIT NONE
1816
1817    INCLUDE "dimensions.h"
1818    INCLUDE "paramet.h"
1819    INCLUDE "netcdf.inc"
1820    INCLUDE "comgeom2.h"
1821    INCLUDE "comconst.h"
1822    INCLUDE "comvert.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.