source: LMDZ6/branches/DYNAMICO-conv/libf/dyn3dmem/guide_loc_mod.F90 @ 3301

Last change on this file since 3301 was 2740, checked in by lguez, 8 years ago

Added output of nudging coefficients for temperature and humidity,
along those for the wind.

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