source: LMDZ6/branches/LMDZ-INCA-Dyn/libf/dyn3dmem/guide_loc_mod.F90 @ 5428

Last change on this file since 5428 was 3803, checked in by Ehouarn Millour, 4 years ago

Make dynamics restart.nc file be in "64-bit offset" mode to enable generating large files. While at it also enforce this for limit.nc, guide_ins.nc and grilles_gcm.nc generated by the model.
EM

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