source: LMDZ5/branches/LF-private/libf/dyn3dmem/guide_loc_mod.F90 @ 5017

Last change on this file since 5017 was 1848, checked in by yann meurdesoif, 11 years ago

Solve performance problem comming from declarations of a derived type.

YM

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