source: LMDZ5/trunk/libf/dyn3dmem/guide_loc_mod.F90 @ 2397

Last change on this file since 2397 was 2263, checked in by lguez, 10 years ago

In the sequential version, ini_getparam and fin_getparam were not
called so guide_init created a file "fort.99". Added call to
ini_getparam and fin_getparam in guide_init. The created file is now
"nudging_parameters_out.txt".

In the parallel version, ini_getparam was called from gcm, with the
argument "out.def". So "out.def" was created even without nudging, it
remained empty. Moved the call to ini_getparam from gcm to
guide_init. Also changed the name of the created file to
"nudging_parameters_out.txt", as in the sequential version. "out.def"
was too vague and confusing a name.

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