source: LMDZ6/trunk/libf/dyn3dmem/guide_loc_mod.f90 @ 5408

Last change on this file since 5408 was 5285, checked in by abarral, 7 weeks ago

As discussed internally, remove generic ONLY: ... for new _mod_h modules

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