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

Last change on this file since 5450 was 5285, checked in by abarral, 2 months 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
RevLine 
[5084]1!
2! $Id$
3!
[1632]4MODULE guide_loc_mod
5
6!=======================================================================
7!   Auteur:  F.Hourdin
8!            F. Codron 01/09
9!=======================================================================
10
[3995]11  USE getparam, only: ini_getparam, fin_getparam, getpar
[1632]12  USE Write_Field_loc
[5270]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
[1823]18  USE parallel_lmdz
[3995]19  USE pres2lev_mod, only: pres2lev
[1632]20
[5285]21  USE paramet_mod_h
[5272]22IMPLICIT NONE
[1632]23
24! ---------------------------------------------
[5249]25! Declarations des cles logiques et parametres
[1632]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
[5249]30  LOGICAL, PRIVATE, SAVE  :: guide_hr,guide_teta
31  LOGICAL, PRIVATE, SAVE  :: guide_BL,guide_reg,guide_add,gamma4,guide_zon
[1632]32  LOGICAL, PRIVATE, SAVE  :: invert_p,invert_y,ini_anal
33  LOGICAL, PRIVATE, SAVE  :: guide_2D,guide_sav,guide_modele
[3795]34!FC
35  LOGICAL, PRIVATE, SAVE  :: convert_Pa
[5249]36
[1632]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
[4246]47  REAL, PRIVATE, SAVE     :: plim_guide_BL
48
[5249]49  REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE     :: alpha_u,alpha_v
50  REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE     :: alpha_T,alpha_Q
[1632]51  REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE     :: alpha_P,alpha_pcor
[5249]52
[1632]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
[5249]70
[3995]71  INTEGER,SAVE,PRIVATE :: ijbu,ijbv,ijeu,ijev !,ijnu,ijnv
[1632]72  INTEGER,SAVE,PRIVATE :: jjbu,jjbv,jjeu,jjev,jjnu,jjnv
73
74
75CONTAINS
76!=======================================================================
77
78  SUBROUTINE guide_init
[2134]79
[2598]80    USE control_mod, ONLY: day_step
81    USE serre_mod, ONLY: grossismx
[2134]82
[5271]83    USE dimensions_mod, ONLY: iim, jjm, llm, ndm
84IMPLICIT NONE
[5249]85
[5271]86
[1632]87
[5272]88
[1632]89    INTEGER                :: error,ncidpl,rid,rcod
90    CHARACTER (len = 80)   :: abort_message
91    CHARACTER (len = 20)   :: modname = 'guide_init'
[3984]92    CHARACTER (len = 20)   :: namedim
[1632]93
94! ---------------------------------------------
[5249]95! Lecture des parametres:
[1632]96! ---------------------------------------------
[2263]97    call ini_getparam("nudging_parameters_out.txt")
[1632]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')
[2134]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)
[1632]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')
[4246]126    CALL getpar('plim_guide_BL',85000.,plim_guide_BL,'BL top presnivs value')
127
[1632]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.
[5084]132    IF (iguide_sav.GT.0) THEN
[2124]133       iguide_sav=day_step/iguide_sav
134    ELSE if (iguide_sav == 0) then
135       iguide_sav = huge(0)
[1632]136    ELSE
[2124]137       iguide_sav=day_step*iguide_sav
[1632]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')
[5084]152    IF (iguide_int.EQ.0) THEN
[1632]153        iguide_int=1
[5084]154    ELSEIF (iguide_int.GT.0) THEN
[1632]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
[3795]165!FC
166    CALL getpar('convert_Pa',.true.,convert_Pa,'Convert Pressure levels in Pa')
[1632]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
[2263]173    call fin_getparam
[5249]174
[1632]175! ---------------------------------------------
176! Determination du nombre de niveaux verticaux
177! des fichiers guidage
178! ---------------------------------------------
179    ncidpl=-99
[5084]180    if (guide_plevs.EQ.1) then
181       if (ncidpl.eq.-99) then
[1902]182          rcod=nf90_open('apbp.nc',Nf90_NOWRITe, ncidpl)
[5270]183          if (rcod.NE.nf90_noerr) THEN
[3995]184             abort_message=' Nudging error -> no file apbp.nc'
[1902]185             CALL abort_gcm(modname,abort_message,1)
186          endif
187       endif
[5084]188    elseif (guide_plevs.EQ.2) then
189       if (ncidpl.EQ.-99) then
[1902]190          rcod=nf90_open('P.nc',Nf90_NOWRITe,ncidpl)
[5270]191          if (rcod.NE.nf90_noerr) THEN
[3995]192             abort_message=' Nudging error -> no file P.nc'
[1902]193             CALL abort_gcm(modname,abort_message,1)
194          endif
195       endif
[3984]196
[1632]197    elseif (guide_u) then
[5084]198       if (ncidpl.eq.-99) then
[1902]199          rcod=nf90_open('u.nc',Nf90_NOWRITe,ncidpl)
[5270]200          if (rcod.NE.nf90_noerr) THEN
[3995]201             abort_message=' Nudging error -> no file u.nc'
[1902]202             CALL abort_gcm(modname,abort_message,1)
203          endif
[5249]204
[1902]205       endif
[3984]206
207
[1632]208    elseif (guide_v) then
[5084]209       if (ncidpl.eq.-99) then
[1902]210          rcod=nf90_open('v.nc',nf90_nowrite,ncidpl)
[5270]211          if (rcod.NE.nf90_noerr) THEN
[3995]212             abort_message=' Nudging error -> no file v.nc'
[1902]213             CALL abort_gcm(modname,abort_message,1)
214          endif
215       endif
[3984]216
[5249]217
[1632]218    elseif (guide_T) then
[5084]219       if (ncidpl.eq.-99) then
[1902]220          rcod=nf90_open('T.nc',nf90_nowrite,ncidpl)
[5270]221          if (rcod.NE.nf90_noerr) THEN
[3995]222             abort_message=' Nudging error -> no file T.nc'
[1902]223             CALL abort_gcm(modname,abort_message,1)
224          endif
225       endif
[3984]226
227
228
[1632]229    elseif (guide_Q) then
[5084]230       if (ncidpl.eq.-99) then
[1902]231          rcod=nf90_open('hur.nc',nf90_nowrite, ncidpl)
[5270]232          if (rcod.NE.nf90_noerr) THEN
[3995]233             abort_message=' Nudging error -> no file hur.nc'
[1902]234             CALL abort_gcm(modname,abort_message,1)
235          endif
236       endif
[3984]237
238
[5249]239    endif
[5270]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
[3995]243        abort_message='Nudging: error reading pressure levels'
[1632]244        CALL abort_gcm(modname,abort_message,1)
245    ENDIF
[5270]246    error=nf90_inquire_dimension(ncidpl,rid,len=nlevnc)
[3995]247    write(*,*)trim(modname)//' : number of vertical levels nlevnc', nlevnc
[1632]248    rcod = nf90_close(ncidpl)
249
250! ---------------------------------------------
251! Allocation des variables
252! ---------------------------------------------
[3995]253    abort_message='nudging allocation error'
[1632]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
[5249]274
[1632]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
[5249]298
[1632]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
[5084]323    IF (guide_plevs.EQ.2) THEN
[1632]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
[5084]331    IF (guide_P.OR.guide_plevs.EQ.1) THEN
[1632]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
[5084]358    IF (guide_plevs.EQ.2) pnat1=pnat2
359    IF (guide_P.OR.guide_plevs.EQ.1) psnat1=psnat2
[1632]360
361  END SUBROUTINE guide_init
362
363!=======================================================================
364  SUBROUTINE guide_main(itau,ucov,vcov,teta,q,masse,ps)
[2021]365    use exner_hyb_loc_m, only: exner_hyb_loc
366    use exner_milieu_loc_m, only: exner_milieu_loc
[1823]367    USE parallel_lmdz
[1632]368    USE control_mod
[1806]369    USE write_field_loc
[2597]370    USE comconst_mod, ONLY: cpp, daysec, dtvr, kappa
[2600]371    USE comvert_mod, ONLY: ap, bp, preff, presnivs, pressure_exner
[5249]372
[5271]373    USE dimensions_mod, ONLY: iim, jjm, llm, ndm
[5285]374USE paramet_mod_h
[5271]375IMPLICIT NONE
[5249]376
[5271]377
[1632]378
[5272]379
[1632]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
[1806]390    REAL, ALLOCATABLE, SAVE, DIMENSION (:,:) :: f_addu ! var aux: champ de guidage
391    REAL, ALLOCATABLE, SAVE, DIMENSION (:,:) :: f_addv ! var aux: champ de guidage
[1632]392    ! Variables pour fonction Exner (P milieu couche)
[2021]393    REAL, ALLOCATABLE, SAVE, DIMENSION (:,:,:)    :: pk
[5249]394    REAL, ALLOCATABLE, SAVE, DIMENSION (:,:)        :: pks
[1632]395    REAL                               :: unskap
[1806]396    REAL, ALLOCATABLE, SAVE, DIMENSION (:,:)    :: p ! besoin si guide_P
[1632]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)
[5249]404
[1632]405    INTEGER       :: i,j,l
[3995]406    CHARACTER(LEN=20) :: modname="guide_main"
[5249]407
408!$OMP MASTER
[3995]409    ijbu=ij_begin ; ijeu=ij_end
[5249]410    jjbu=jj_begin ; jjeu=jj_end ; jjnu=jjeu-jjbu+1
[3995]411    ijbv=ij_begin ; ijev=ij_end
[5249]412    jjbv=jj_begin ; jjev=jj_end ; jjnv=jjev-jjbv+1
[1632]413    IF (pole_sud) THEN
[3995]414      ijeu=ij_end-iip1
[1632]415      ijev=ij_end-iip1
416      jjev=jj_end-1
[5249]417      jjnv=jjev-jjbv+1
[1632]418    ENDIF
[3995]419    IF (pole_nord) THEN
420      ijbu=ij_begin+iip1
421      ijbv=ij_begin
422    ENDIF
[1632]423!$OMP END MASTER
424!$OMP BARRIER
[5249]425
[2036]426!    PRINT *,'---> on rentre dans guide_main'
[1632]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)
[5249]432
[1632]433!-----------------------------------------------------------------------
434! Initialisations au premier passage
435!-----------------------------------------------------------------------
436
437    IF (first) THEN
438        first=.FALSE.
439!$OMP MASTER
[5249]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
[1632]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
[5249]452        factt=dtvr*iperiod/daysec
[1632]453!$OMP MASTER
[1806]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)
[1632]459! correction de rappel dans couche limite
460        if (guide_BL) then
461             alpha_pcor(:)=1.
462        else
463            do l=1,llm
[4246]464                alpha_pcor(l)=(1.+tanh(((plim_guide_BL-presnivs(l))/preff)/0.05))/2.
[1632]465            enddo
466        endif
467!$OMP END MASTER
[1806]468!$OMP BARRIER
[5249]469! ini_anal: etat initial egal au guidage
[1632]470        IF (ini_anal) THEN
471            CALL guide_interp(ps,teta)
[5249]472!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[1806]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
[5249]479
[1632]480            IF (guide_P) THEN
[1806]481!$OMP MASTER
[1632]482                ps(ijbu:ijeu)=psgui2(ijbu:ijeu)
[1806]483!$OMP END MASTER
484!$OMP BARRIER
[1632]485                CALL pression_loc(ijnb_u,ap,bp,ps,p)
486                CALL massdair_loc(p,masse)
[1806]487!$OMP BARRIER
[1632]488            ENDIF
489            RETURN
490        ENDIF
491
492    ENDIF !first
493
494!-----------------------------------------------------------------------
495! Lecture des fichiers de guidage ?
496!-----------------------------------------------------------------------
[5084]497    IF (iguide_read.NE.0) THEN
[1632]498      ditau=real(itau)
499      dday_step=real(day_step)
[5084]500      IF (iguide_read.LT.0) THEN
[1632]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)
[5084]506      IF (reste.EQ.0.) THEN
507          IF (itau_test.EQ.itau) THEN
[3995]508            write(*,*)trim(modname)//' second pass in advreel at itau=',&
509            itau
[4469]510            CALL abort_gcm("guide_loc_lod","stopped",1)
[1632]511          ELSE
[1806]512!$OMP MASTER
[1632]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,:)
[5084]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)
[1806]519!$OMP END MASTER
520!$OMP BARRIER
[1632]521              step_rea=step_rea+1
522              itau_test=itau
[3995]523              if (is_master) then
524                write(*,*)trim(modname)//' Reading nudging files, step ',&
525                    step_rea,'after ',count_no_rea,' skips'
526              endif
[1632]527              IF (guide_2D) THEN
[1806]528!$OMP MASTER
[1632]529                  CALL guide_read2D(step_rea)
[1806]530!$OMP END MASTER
531!$OMP BARRIER
[1632]532              ELSE
[1806]533!$OMP MASTER
[1632]534                  CALL guide_read(step_rea)
[1806]535!$OMP END MASTER
536!$OMP BARRIER
[1632]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!-----------------------------------------------------------------------
[5084]549    IF (MOD(itau,iguide_int).EQ.0) THEN
[1632]550        CALL guide_interp(ps,teta)
551    ENDIF
552! Repartition entre 2 etats de guidage
[5084]553    IF (iguide_read.NE.0) THEN
[1632]554        tau=reste
555    ELSE
556        tau=1.
557    ENDIF
558
[1806]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)
[5249]563
564
[3995]565!-----------------------------------------------------------------------
[5249]566!   Ajout des champs de guidage
[1632]567!-----------------------------------------------------------------------
568! Sauvegarde du guidage?
[5249]569    f_out=((MOD(itau,iguide_sav).EQ.0).AND.guide_sav)
[1632]570    IF (f_out) THEN
[2036]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
[1632]584        unskap=1./kappa
[2036]585!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[2028]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
[2036]593
594        CALL guide_out("SP",jjp1,llm,p(ijb_u:ije_u,1:llm),1.)
[1632]595    ENDIF
[5249]596
[1632]597    if (guide_u) then
598        if (guide_add) then
[2036]599!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[1806]600          DO l=1,llm
601           f_addu(ijbu:ijeu,l)=(1.-tau)*ugui1(ijbu:ijeu,l)+tau*ugui2(ijbu:ijeu,l)
602          ENDDO
[1632]603        else
[2036]604!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[1806]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
[5249]608        endif
609
[1806]610!        CALL WriteField_u('f_addu',f_addu)
[1632]611
[1806]612        if (guide_zon) CALL guide_zonave_u(1,llm,f_addu)
613        CALL guide_addfield_u(llm,f_addu,alpha_u)
[2028]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)
[3995]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
[2036]622!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[1806]623        DO l=1,llm
624          ucov(ijbu:ijeu,l)=ucov(ijbu:ijeu,l)+f_addu(ijbu:ijeu,l)
625        ENDDO
626
[1632]627    endif
628
629    if (guide_T) then
630        if (guide_add) then
[2036]631!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[1806]632          DO l=1,llm
633            f_addu(ijbu:ijeu,l)=(1.-tau)*tgui1(ijbu:ijeu,l)+tau*tgui2(ijbu:ijeu,l)
634          ENDDO
[1632]635        else
[2036]636!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[1806]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
[5249]640        endif
[1806]641        if (guide_zon) CALL guide_zonave_u(2,llm,f_addu)
642        CALL guide_addfield_u(llm,f_addu,alpha_T)
[2028]643        IF (f_out) CALL guide_out("teta",jjp1,llm,f_addu(:,:)/factt,factt)
[2036]644!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[1806]645        DO l=1,llm
646          teta(ijbu:ijeu,l)=teta(ijbu:ijeu,l)+f_addu(ijbu:ijeu,l)
647        ENDDO
[1632]648    endif
649
650    if (guide_P) then
651        if (guide_add) then
[1806]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
[1632]656        else
[1806]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
[5249]661        endif
[1806]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)
[2036]664!       IF (f_out) CALL guide_out("ps",jjp1,1,f_addu(ijb_u:ije_u,1)/factt,factt)
[1806]665!$OMP MASTER
666        ps(ijbu:ijeu)=ps(ijbu:ijeu)+f_addu(ijbu:ijeu,1)
667!$OMP END MASTER
668!$OMP BARRIER
[1632]669        CALL pression_loc(ijnb_u,ap,bp,ps,p)
670        CALL massdair_loc(p,masse)
[1806]671!$OMP BARRIER
[1632]672    endif
673
674    if (guide_Q) then
675        if (guide_add) then
[2036]676!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[1806]677          DO l=1,llm
678            f_addu(ijbu:ijeu,l)=(1.-tau)*qgui1(ijbu:ijeu,l)+tau*qgui2(ijbu:ijeu,l)
679          ENDDO
[1632]680        else
[2036]681!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[1806]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
[5249]685        endif
[1806]686        if (guide_zon) CALL guide_zonave_u(2,llm,f_addu)
687        CALL guide_addfield_u(llm,f_addu,alpha_Q)
[2028]688        IF (f_out) CALL guide_out("q",jjp1,llm,f_addu(:,:)/factt,factt)
[1806]689
[2036]690!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[1806]691        DO l=1,llm
692          q(ijbu:ijeu,l)=q(ijbu:ijeu,l)+f_addu(ijbu:ijeu,l)
693        ENDDO
[1632]694    endif
695
696    if (guide_v) then
697        if (guide_add) then
[2036]698!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[1806]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
[1632]703        else
[2036]704!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[1806]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
[5249]709        endif
710
[1806]711        if (guide_zon) CALL guide_zonave_v(2,jjm,llm,f_addv(ijb_v:ije_v,:))
[5249]712
[1806]713        CALL guide_addfield_v(llm,f_addv(ijb_v:ije_v,:),alpha_v)
[2028]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)
[3995]716        IF (f_out) THEN
717          ! Ehouarn: Fill in the gaps adequately
[5249]718          IF (ijbv>ijb_v) f_addv(ijb_v:ijbv-1,:)=0
[3995]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
[1806]722
[2036]723!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[1806]724        DO l=1,llm
725          vcov(ijbv:ijev,l)=vcov(ijbv:ijev,l)+f_addv(ijbv:ijev,l)
726        ENDDO
[1632]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
[5271]735    USE dimensions_mod, ONLY: iim, jjm, llm, ndm
[5285]736USE paramet_mod_h
[5271]737IMPLICIT NONE
738
[1632]739
[5272]740
[1632]741    ! input variables
742    INTEGER,                      INTENT(IN)    :: vsize
[5249]743    REAL, DIMENSION(ijb_u:ije_u),       INTENT(IN)    :: alpha
[1632]744    REAL, DIMENSION(ijb_u:ije_u,vsize), INTENT(INOUT) :: field
745
746    ! Local variables
747    INTEGER :: l
748
[2036]749!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[1632]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
[5271]760    USE dimensions_mod, ONLY: iim, jjm, llm, ndm
[5285]761USE paramet_mod_h
[5271]762IMPLICIT NONE
763
[1632]764
[5272]765
[1632]766    ! input variables
767    INTEGER,                      INTENT(IN)    :: vsize
[5249]768    REAL, DIMENSION(ijb_v:ije_v),       INTENT(IN)    :: alpha
[1632]769    REAL, DIMENSION(ijb_v:ije_v,vsize), INTENT(INOUT) :: field
770
771    ! Local variables
772    INTEGER :: l
773
[2036]774!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[1632]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
[5249]780
[1632]781!=======================================================================
782
783  SUBROUTINE guide_zonave_u(typ,vsize,field)
784
[5281]785    USE comgeom_mod_h
[2597]786    USE comconst_mod, ONLY: pi
[5249]787
[5271]788    USE dimensions_mod, ONLY: iim, jjm, llm, ndm
[5285]789    USE paramet_mod_h
[5271]790IMPLICIT NONE
[1632]791
[5271]792
[5272]793
[5249]794
[1632]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.
[1806]802!$OMP THREADPRIVATE(first)
803
[1632]804    INTEGER, DIMENSION (2), SAVE :: imin, imax ! averaging domain
[5249]805!$OMP THREADPRIVATE(imin,imax)
[1632]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
[5084]818                IF (lond(i).LT.lon_min_g) imin(1)=i
819                IF (lond(i).LE.lon_max_g) imax(1)=i
[1632]820            ENDDO
821            lond=rlonv*180./pi
822            DO i=1,iim
[5084]823                IF (lond(i).LT.lon_min_g) imin(2)=i
824                IF (lond(i).LE.lon_max_g) imax(2)=i
[1632]825            ENDDO
826        ENDIF
827    ENDIF
828
[5249]829
[2036]830!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[1632]831      DO l=1,vsize
[1806]832        fieldm(:,l)=0.
[1632]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
[5249]843          ENDDO
[1632]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
[5281]859    USE comgeom_mod_h
[2597]860    USE comconst_mod, ONLY: pi
[5249]861
[5271]862    USE dimensions_mod, ONLY: iim, jjm, llm, ndm
[5285]863USE paramet_mod_h
[5271]864IMPLICIT NONE
[1632]865
[5271]866
[5272]867
[5249]868
[1632]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.
[1806]877!$OMP THREADPRIVATE(first)
[1632]878    INTEGER, DIMENSION (2), SAVE :: imin, imax ! averaging domain
[1806]879!$OMP THREADPRIVATE(imin, imax)
[1632]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
[5084]892                IF (lond(i).LT.lon_min_g) imin(1)=i
893                IF (lond(i).LE.lon_max_g) imax(1)=i
[1632]894            ENDDO
895            lond=rlonv*180./pi
896            DO i=1,iim
[5084]897                IF (lond(i).LT.lon_min_g) imin(2)=i
898                IF (lond(i).LE.lon_max_g) imax(2)=i
[1632]899            ENDDO
900        ENDIF
901    ENDIF
902
[2036]903!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[1632]904      DO l=1,vsize
905      ! Compute zonal average
[1806]906          fieldm(:,l)=0.
[1632]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
[5249]912          ENDDO
[1632]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
[5249]925
[1632]926!=======================================================================
927  SUBROUTINE guide_interp(psi,teta)
[5281]928    USE comgeom2_mod_h
929  use exner_hyb_loc_m, only: exner_hyb_loc
[2021]930    use exner_milieu_loc_m, only: exner_milieu_loc
[1823]931  USE parallel_lmdz
[1632]932  USE mod_hallo
933  USE Bands
[2597]934  USE comconst_mod, ONLY: cpp, kappa
[2600]935  USE comvert_mod, ONLY: preff, pressure_exner, bp, ap, disvert_type
[5271]936  USE dimensions_mod, ONLY: iim, jjm, llm, ndm
[5285]937USE paramet_mod_h
[5271]938IMPLICIT NONE
[1632]939
[5271]940
[5272]941
[1632]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.
[1806]947!$OMP THREADPRIVATE(first)
[1632]948  ! Variables pour niveaux pression:
[1806]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
[5249]952  REAL, ALLOCATABLE, SAVE,DIMENSION (:,:,:)  :: p           ! pression intercouches
[1806]953  REAL, ALLOCATABLE, SAVE,DIMENSION (:,:,:)    :: pls, pext   ! var intermediaire
[5249]954  REAL, ALLOCATABLE, SAVE,DIMENSION (:,:,:)    :: pbarx
955  REAL, ALLOCATABLE, SAVE,DIMENSION (:,:,:)     :: pbary
[1632]956  ! Variables pour fonction Exner (P milieu couche)
[2021]957  REAL, ALLOCATABLE, SAVE,DIMENSION (:,:,:)    :: pk
[5249]958  REAL ,ALLOCATABLE, SAVE, DIMENSION (:,:)        :: pks
[1632]959  REAL                               :: unskap
960  ! Pression de vapeur saturante
[1806]961  REAL, ALLOCATABLE, SAVE,DIMENSION (:,:)      :: qsat
[1632]962  !Variables intermediaires interpolation
[5249]963  REAL, ALLOCATABLE, SAVE,DIMENSION (:,:,:)    :: zu1,zu2
[1806]964  REAL, ALLOCATABLE, SAVE,DIMENSION (:,:,:)     :: zv1,zv2
[5249]965
[1632]966  INTEGER                            :: i,j,l,ij
[3995]967  CHARACTER(LEN=20),PARAMETER :: modname="guide_interp"
[5249]968  TYPE(Request),SAVE :: Req
[1848]969!$OMP THREADPRIVATE(Req)
[5249]970
[3995]971    if (is_master) write(*,*)trim(modname)//': interpolate nudging variables'
[1632]972! -----------------------------------------------------------------
973! Calcul des niveaux de pression champs guidage (pour T et Q)
974! -----------------------------------------------------------------
[1806]975    IF (first) THEN
976!$OMP MASTER
[5249]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) )
[1806]993      ALLOCATE(zv2(iip1,jjb_v:jje_v,llm) )
994!$OMP END MASTER
995!$OMP BARRIER
[5249]996    ENDIF
[1806]997
[5249]998
999
1000
[5084]1001    IF (guide_plevs.EQ.0) THEN
[1806]1002!$OMP DO
[1632]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
[5249]1011    ENDIF
[1632]1012
1013    if (first) then
1014        first=.FALSE.
[1806]1015!$OMP MASTER
[3995]1016        write(*,*)trim(modname)//' : check vertical level order'
1017        write(*,*)trim(modname)//' LMDZ :'
[1632]1018        do l=1,llm
[3995]1019          write(*,*)trim(modname)//' PL(',l,')=',(ap(l)+ap(l+1))/2. &
[1632]1020                  +psi(1,jjeu)*(bp(l)+bp(l+1))/2.
1021        enddo
[3995]1022        write(*,*)trim(modname)//' nudging file :'
[1632]1023        SELECT CASE (guide_plevs)
[5249]1024        CASE (0)
[1632]1025            do l=1,nlevnc
[3995]1026              write(*,*)trim(modname)//' PL(',l,')=',plnc2(1,jjbu,l)
[1632]1027            enddo
1028        CASE (1)
1029            DO l=1,nlevnc
[3995]1030              write(*,*)trim(modname)//' PL(',l,')=',&
[4041]1031                        apnc(l)+bpnc(l)*psnat2(1,jjbu)
[3995]1032            ENDDO
[1632]1033        CASE (2)
1034            do l=1,nlevnc
[3995]1035              write(*,*)trim(modname)//' PL(',l,')=',pnat2(1,jjbu,l)
[1632]1036            enddo
1037        END SELECT
[3995]1038        write(*,*)trim(modname)//' invert ordering: invert_p=',invert_p
[1632]1039        if (guide_u) then
1040            do l=1,nlevnc
[3995]1041              write(*,*)trim(modname)//' U(',l,')=',unat2(1,jjbu,l)
[1632]1042            enddo
1043        endif
1044        if (guide_T) then
1045            do l=1,nlevnc
[3995]1046              write(*,*)trim(modname)//' T(',l,')=',tnat2(1,jjbu,l)
[1632]1047            enddo
1048        endif
[1806]1049!$OMP END MASTER
[3995]1050    endif ! of if (first)
[4042]1051
[4054]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
[1632]1062! -----------------------------------------------------------------
[5249]1063! Calcul niveaux pression modele
[1632]1064! -----------------------------------------------------------------
1065
1066!    ....  Calcul de pls , pression au milieu des couches ,en Pascals
[5084]1067    IF (guide_plevs.EQ.1) THEN
[2036]1068!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[1632]1069        DO l=1,llm
[2597]1070            DO j=jjbu,jjeu
1071                DO i =1, iip1
[1632]1072                    pls(i,j,l)=(ap(l)+ap(l+1))/2.+psi(i,j)*(bp(l)+bp(l+1))/2.
[2028]1073                ENDDO
1074            ENDDO
[1632]1075        ENDDO
1076    ELSE
[2028]1077        unskap=1./kappa
[1806]1078!$OMP BARRIER
[2036]1079!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[2028]1080   DO l = 1, llm
1081       DO j=jjbu,jjeu
[2597]1082           DO i =1, iip1
1083               pls(i,j,l) = preff * ( pk(i,j,l)/cpp) ** unskap
1084           ENDDO
[2028]1085       ENDDO
1086   ENDDO
[1632]1087    ENDIF
1088
1089!   calcul des pressions pour les grilles u et v
[2036]1090!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[1632]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)
[1806]1101!$OMP BARRIER
[1632]1102     CALL WaitRequest(Req)
[1806]1103!$OMP BARRIER
[1632]1104
1105    call massbar_loc(pext, pbarx, pbary )
[1806]1106!$OMP BARRIER
[2036]1107!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[1632]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
[2036]1116!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[1632]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
[1806]1130!$OMP MASTER
[1632]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
[1806]1140!$OMP END MASTER
1141!$OMP BARRIER
[1632]1142    endif
1143
1144    IF (guide_T) THEN
1145        ! Calcul des nouvelles valeurs des niveaux de pression du guidage
[5084]1146        IF (guide_plevs.EQ.1) THEN
[1806]1147!$OMP DO
[1632]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
[5084]1156        ELSE IF (guide_plevs.EQ.2) THEN
[1806]1157!$OMP DO
[1632]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
[1806]1169!$OMP MASTER
[1632]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)
[1806]1174!$OMP END MASTER
1175!$OMP BARRIER
[1632]1176        ! Conversion en variables GCM
[2036]1177!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[1632]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
[5249]1193                tgui1(j*iip1,l)=tgui1((j-1)*iip1+1,l)
1194                tgui2(j*iip1,l)=tgui2((j-1)*iip1+1,l)
[1632]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
[5249]1204                tgui1(ip1jm+i,l)=tgui1(ip1jm+1,l)
1205                tgui2(ip1jm+i,l)=tgui2(ip1jm+1,l)
[1632]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
[5084]1213        IF (guide_plevs.EQ.1) THEN
[1806]1214!$OMP DO
[1632]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
[5084]1223        ELSE IF (guide_plevs.EQ.2) THEN
[1806]1224!$OMP DO
[1632]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
[1806]1236!$OMP MASTER
[1632]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)
[1806]1241!$OMP END MASTER
1242!$OMP BARRIER
[1632]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.
[2036]1247!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[1632]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
[5249]1255                qgui1(j*iip1,l)=qgui1((j-1)*iip1+1,l)
1256                qgui2(j*iip1,l)=qgui2((j-1)*iip1+1,l)
[1632]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
[3547]1264            if (pole_sud) then
[1632]1265              do i=1,iip1
[5249]1266                qgui1(ip1jm+i,l)=qgui1(ip1jm+1,l)
1267                qgui2(ip1jm+i,l)=qgui2(ip1jm+1,l)
[1632]1268              enddo
1269            endif
1270        enddo
1271        IF (guide_hr) THEN
[2036]1272!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[1806]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 %
[5249]1277            qgui2(ijbu:ijeu,l)=qgui2(ijbu:ijeu,l)*qsat(ijbu:ijeu,l)*0.01
[1806]1278          enddo
1279
[1632]1280        ENDIF
1281    ENDIF
1282
1283    IF (guide_u) THEN
1284        ! Calcul des nouvelles valeurs des niveaux de pression du guidage
[5084]1285        IF (guide_plevs.EQ.1) THEN
[1806]1286!$OMP DO
[1632]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
[5084]1299        ELSE IF (guide_plevs.EQ.2) THEN
[1806]1300!$OMP DO
[1632]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
[5249]1314
[1632]1315        ! Interpolation verticale
[1806]1316!$OMP MASTER
[1632]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)
[1806]1321!$OMP END MASTER
1322!$OMP BARRIER
[1632]1323
1324        ! Conversion en variables GCM
[2036]1325!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[1632]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
[5249]1333                ugui1(j*iip1,l)=ugui1((j-1)*iip1+1,l)
1334                ugui2(j*iip1,l)=ugui2((j-1)*iip1+1,l)
[1632]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
[5249]1350
[1632]1351    IF (guide_v) THEN
1352        ! Calcul des nouvelles valeurs des niveaux de pression du guidage
[5084]1353        IF (guide_plevs.EQ.1) THEN
[2597]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)
[1806]1357!$OMP BARRIER
[2597]1358         CALL WaitRequest(Req)
[1806]1359!$OMP BARRIER
1360!$OMP DO
[1632]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
[5084]1371        ELSE IF (guide_plevs.EQ.2) THEN
[2597]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)
[1806]1375!$OMP BARRIER
[2597]1376         CALL WaitRequest(Req)
[1806]1377!$OMP BARRIER
1378!$OMP DO
[1632]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
[1806]1391
1392!$OMP MASTER
[1632]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)
[1806]1397!$OMP END MASTER
1398!$OMP BARRIER
[1632]1399        ! Conversion en variables GCM
[2036]1400!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[1632]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
[5249]1408                vgui1(j*iip1,l)=vgui1((j-1)*iip1+1,l)
1409                vgui2(j*iip1,l)=vgui2((j-1)*iip1+1,l)
[1632]1410            enddo
1411        enddo
1412    ENDIF
1413
[5249]1414
[1632]1415  END SUBROUTINE guide_interp
1416
1417!=======================================================================
[1806]1418  SUBROUTINE tau2alpha(typ,pim,jjb,jje,factt,taumin,taumax,alpha)
[1632]1419
1420! Calcul des constantes de rappel alpha (=1/tau)
1421
[2597]1422    use comconst_mod, only: pi
[2598]1423    use serre_mod, only: clat, clon, grossismx, grossismy
[5249]1424
[5271]1425    USE dimensions_mod, ONLY: iim, jjm, llm, ndm
[5285]1426    USE paramet_mod_h
[5281]1427    USE comgeom2_mod_h
1428    implicit none
[1632]1429
[5271]1430
[1632]1431! input arguments :
1432    INTEGER, INTENT(IN) :: typ    ! u(2),v(3), ou scalaire(1)
[1806]1433    INTEGER, INTENT(IN) :: pim ! dimensions en lon
1434    INTEGER, INTENT(IN) :: jjb,jje ! dimensions en lat
[1632]1435    REAL, INTENT(IN)    :: factt   ! pas de temps en fraction de jour
1436    REAL, INTENT(IN)    :: taumin,taumax
1437! output arguments:
[5249]1438    REAL, DIMENSION(pim,jjb:jje), INTENT(OUT) :: alpha
1439
[1632]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
[3995]1450    character(len=20),parameter :: modname="tau2alpha"
[1632]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
[1806]1461            do j=jjb,jje
[1632]1462                do i=1,pim
[5084]1463                    if (typ.eq.2) then
[1632]1464                       zlat=rlatu(j)*180./pi
1465                       zlon=rlonu(i)*180./pi
[5084]1466                    elseif (typ.eq.1) then
[1632]1467                       zlat=rlatu(j)*180./pi
1468                       zlon=rlonv(i)*180./pi
[5084]1469                    elseif (typ.eq.3) then
[1632]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
[5084]1508        IF (typ.EQ.2) THEN
[1632]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
[5084]1516        IF (typ.EQ.3) THEN
[1632]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.
[5249]1524        IF (first) THEN
[1632]1525            first=.FALSE.
1526            ! coordonnees du centre du zoom
[5249]1527            CALL coordij(clon,clat,ilon,ilat)
[1632]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
[5084]1538            if (abs(grossismx-1.).lt.0.1.or.abs(grossismy-1.).lt.0.1) then
[3995]1539              write(*,*)trim(modname)//' ATTENTION modele peu zoome'
1540              write(*,*)trim(modname)//' ATTENTION on prend une constante de guidage cste'
1541              gamma=0.
[1632]1542            else
[3995]1543              gamma=(dxdy_max-2.*dxdy_min)/(dxdy_max-dxdy_min)
1544              write(*,*)trim(modname)//' gamma=',gamma
[5084]1545              if (gamma.lt.1.e-5) then
[3995]1546                write(*,*)trim(modname)//' gamma =',gamma,'<1e-5'
[4469]1547                CALL abort_gcm("guide_loc_mod","stopped",1)
[3995]1548              endif
1549              gamma=log(0.5)/log(gamma)
[5249]1550              if (gamma4) then
[3995]1551                gamma=min(gamma,4.)
1552              endif
1553              write(*,*)trim(modname)//' gamma=',gamma
[1632]1554            endif
1555        ENDIF !first
1556
[1806]1557        do j=jjb,jje
[1632]1558            do i=1,pim
[5084]1559                if (typ.eq.1) then
[1632]1560                   dxdy_=dxdys(i,j)
1561                   zlat=rlatu(j)*180./pi
[5084]1562                elseif (typ.eq.2) then
[1632]1563                   dxdy_=dxdyu(i,j)
1564                   zlat=rlatu(j)*180./pi
[5084]1565                elseif (typ.eq.3) then
[1632]1566                   dxdy_=dxdyv(i,j)
1567                   zlat=rlatv(j)*180./pi
1568                endif
[5084]1569                if (abs(grossismx-1.).lt.0.1.or.abs(grossismy-1.).lt.0.1) then
[1632]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.)
[5084]1575                    if(lat_min_g.le.zlat .and. zlat.le.lat_max_g) then
[1632]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
[2134]1585    if (.not. guide_add) alpha = 1. - exp(- alpha)
1586
[1632]1587  END SUBROUTINE tau2alpha
1588
1589!=======================================================================
1590  SUBROUTINE guide_read(timestep)
[5249]1591    USE netcdf, ONLY: nf90_put_var
[5271]1592    USE dimensions_mod, ONLY: iim, jjm, llm, ndm
[5285]1593USE paramet_mod_h
[5271]1594IMPLICIT NONE
1595
[1632]1596
[5272]1597
[1632]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
[3984]1604    INTEGER               :: ncidpl,varidpl,varidap,varidbp,dimid,lendim
[1632]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'
[3984]1610    CHARACTER (len = 20)   :: namedim
[1632]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
[5279]1618         write(*,*) trim(modname)//': opening nudging files '
[1632]1619! Ap et Bp si Niveaux de pression hybrides
[5084]1620         if (guide_plevs.EQ.1) then
[5279]1621             write(*,*) trim(modname)//' Reading nudging on model levels'
[1632]1622             rcode = nf90_open('apbp.nc', nf90_nowrite, ncidpl)
[5270]1623             IF (rcode.NE.nf90_noerr) THEN
[3995]1624              abort_message='Nudging: error -> no file apbp.nc'
[1902]1625              CALL abort_gcm(modname,abort_message,1)
1626             ENDIF
[1632]1627             rcode = nf90_inq_varid(ncidpl, 'AP', varidap)
[5270]1628             IF (rcode.NE.nf90_noerr) THEN
[3995]1629              abort_message='Nudging: error -> no AP variable in file apbp.nc'
[1902]1630              CALL abort_gcm(modname,abort_message,1)
1631             ENDIF
[1632]1632             rcode = nf90_inq_varid(ncidpl, 'BP', varidbp)
[5270]1633             IF (rcode.NE.nf90_noerr) THEN
[3995]1634              abort_message='Nudging: error -> no BP variable in file apbp.nc'
[1902]1635              CALL abort_gcm(modname,abort_message,1)
1636             ENDIF
[5279]1637             write(*,*) trim(modname)//' ncidpl,varidap',ncidpl,varidap
[1632]1638         endif
[5249]1639
[1632]1640! Pression si guidage sur niveaux P variables
[5084]1641         if (guide_plevs.EQ.2) then
[1632]1642             rcode = nf90_open('P.nc', nf90_nowrite, ncidp)
[5270]1643             IF (rcode.NE.nf90_noerr) THEN
[3995]1644              abort_message='Nudging: error -> no file P.nc'
[1902]1645              CALL abort_gcm(modname,abort_message,1)
1646             ENDIF
[1632]1647             rcode = nf90_inq_varid(ncidp, 'PRES', varidp)
[5270]1648             IF (rcode.NE.nf90_noerr) THEN
[3995]1649              abort_message='Nudging: error -> no PRES variable in file P.nc'
[1902]1650              CALL abort_gcm(modname,abort_message,1)
1651             ENDIF
[5279]1652             write(*,*) trim(modname)//' ncidp,varidp',ncidp,varidp
[5084]1653             if (ncidpl.eq.-99) ncidpl=ncidp
[1632]1654         endif
[3984]1655
[1632]1656! Vent zonal
1657         if (guide_u) then
1658             rcode = nf90_open('u.nc', nf90_nowrite, ncidu)
[5270]1659             IF (rcode.NE.nf90_noerr) THEN
[3995]1660              abort_message='Nudging: error -> no file u.nc'
[1902]1661              CALL abort_gcm(modname,abort_message,1)
1662             ENDIF
[1632]1663             rcode = nf90_inq_varid(ncidu, 'UWND', varidu)
[5270]1664             IF (rcode.NE.nf90_noerr) THEN
[3995]1665              abort_message='Nudging: error -> no UWND variable in file u.nc'
[1902]1666              CALL abort_gcm(modname,abort_message,1)
1667             ENDIF
[5279]1668             write(*,*) trim(modname)//' ncidu,varidu',ncidu,varidu
[5084]1669             if (ncidpl.eq.-99) ncidpl=ncidu
[3984]1670
[5249]1671
[3984]1672             status=NF90_INQ_DIMID(ncidu, "LONU", dimid)
1673             status=NF90_INQUIRE_DIMENSION(ncidu,dimid,namedim,lendim)
[5084]1674             IF (lendim .NE. iip1) THEN
[3995]1675                abort_message='dimension LONU different from iip1 in u.nc'
[3984]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)
[5084]1681             IF (lendim .NE. jjp1) THEN
[3995]1682                abort_message='dimension LATU different from jjp1 in u.nc'
[3984]1683                CALL abort_gcm(modname,abort_message,1)
1684             ENDIF
[5249]1685
[1632]1686         endif
[3984]1687
[1632]1688! Vent meridien
1689         if (guide_v) then
1690             rcode = nf90_open('v.nc', nf90_nowrite, ncidv)
[5270]1691             IF (rcode.NE.nf90_noerr) THEN
[3995]1692              abort_message='Nudging: error -> no file v.nc'
[1902]1693              CALL abort_gcm(modname,abort_message,1)
1694             ENDIF
[1632]1695             rcode = nf90_inq_varid(ncidv, 'VWND', varidv)
[5270]1696             IF (rcode.NE.nf90_noerr) THEN
[3995]1697              abort_message='Nudging: error -> no VWND variable in file v.nc'
[1902]1698              CALL abort_gcm(modname,abort_message,1)
1699             ENDIF
[5279]1700             write(*,*) trim(modname)//' ncidv,varidv',ncidv,varidv
[5084]1701             if (ncidpl.eq.-99) ncidpl=ncidv
[5249]1702
[3984]1703             status=NF90_INQ_DIMID(ncidv, "LONV", dimid)
1704             status=NF90_INQUIRE_DIMENSION(ncidv,dimid,namedim,lendim)
[5249]1705
[5084]1706                IF (lendim .NE. iip1) THEN
[3995]1707                abort_message='dimension LONV different from iip1 in v.nc'
[3984]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)
[5084]1714             IF (lendim .NE. jjm) THEN
[3995]1715                abort_message='dimension LATV different from jjm in v.nc'
[3984]1716                CALL abort_gcm(modname,abort_message,1)
1717             ENDIF
[5249]1718
[3984]1719        endif
1720
[1632]1721! Temperature
1722         if (guide_T) then
1723             rcode = nf90_open('T.nc', nf90_nowrite, ncidt)
[5270]1724             IF (rcode.NE.nf90_noerr) THEN
[3995]1725              abort_message='Nudging: error -> no file T.nc'
[1902]1726              CALL abort_gcm(modname,abort_message,1)
1727             ENDIF
[1632]1728             rcode = nf90_inq_varid(ncidt, 'AIR', varidt)
[5270]1729             IF (rcode.NE.nf90_noerr) THEN
[3995]1730              abort_message='Nudging: error -> no AIR variable in file T.nc'
[1902]1731              CALL abort_gcm(modname,abort_message,1)
1732             ENDIF
[5279]1733             write(*,*) trim(modname)//' ncidT,varidT',ncidt,varidt
[5084]1734             if (ncidpl.eq.-99) ncidpl=ncidt
[3984]1735
1736             status=NF90_INQ_DIMID(ncidt, "LONV", dimid)
1737             status=NF90_INQUIRE_DIMENSION(ncidt,dimid,namedim,lendim)
[5084]1738             IF (lendim .NE. iip1) THEN
[3995]1739                abort_message='dimension LONV different from iip1 in T.nc'
[3984]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)
[5084]1745             IF (lendim .NE. jjp1) THEN
[3995]1746                abort_message='dimension LATU different from jjp1 in T.nc'
[3984]1747                CALL abort_gcm(modname,abort_message,1)
1748             ENDIF
1749
[1632]1750         endif
[3984]1751
[1632]1752! Humidite
1753         if (guide_Q) then
1754             rcode = nf90_open('hur.nc', nf90_nowrite, ncidQ)
[5270]1755             IF (rcode.NE.nf90_noerr) THEN
[3995]1756              abort_message='Nudging: error -> no file hur.nc'
[1902]1757              CALL abort_gcm(modname,abort_message,1)
1758             ENDIF
[1632]1759             rcode = nf90_inq_varid(ncidQ, 'RH', varidQ)
[5270]1760             IF (rcode.NE.nf90_noerr) THEN
[3995]1761              abort_message='Nudging: error -> no RH variable in file hur.nc'
[1902]1762              CALL abort_gcm(modname,abort_message,1)
1763             ENDIF
[5279]1764             write(*,*) trim(modname)//' ncidQ,varidQ',ncidQ,varidQ
[5084]1765             if (ncidpl.eq.-99) ncidpl=ncidQ
[3984]1766
1767
1768             status=NF90_INQ_DIMID(ncidQ, "LONV", dimid)
1769             status=NF90_INQUIRE_DIMENSION(ncidQ,dimid,namedim,lendim)
[5084]1770             IF (lendim .NE. iip1) THEN
[3995]1771                abort_message='dimension LONV different from iip1 in hur.nc'
[3984]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)
[5084]1777             IF (lendim .NE. jjp1) THEN
[3995]1778                abort_message='dimension LATU different from jjp1 in hur.nc'
[3984]1779                CALL abort_gcm(modname,abort_message,1)
1780             ENDIF
1781
1782
[1632]1783         endif
1784! Pression de surface
[5084]1785         if ((guide_P).OR.(guide_plevs.EQ.1)) then
[1632]1786             rcode = nf90_open('ps.nc', nf90_nowrite, ncidps)
[5270]1787             IF (rcode.NE.nf90_noerr) THEN
[3995]1788              abort_message='Nudging: error -> no file ps.nc'
[1902]1789              CALL abort_gcm(modname,abort_message,1)
1790             ENDIF
[1632]1791             rcode = nf90_inq_varid(ncidps, 'SP', varidps)
[5270]1792             IF (rcode.NE.nf90_noerr) THEN
[3995]1793              abort_message='Nudging: error -> no SP variable in file ps.nc'
[1902]1794              CALL abort_gcm(modname,abort_message,1)
1795             ENDIF
[5279]1796             write(*,*) trim(modname)//' ncidps,varidps',ncidps,varidps
[1632]1797         endif
1798! Coordonnee verticale
[5084]1799         if (guide_plevs.EQ.0) then
[1632]1800              rcode = nf90_inq_varid(ncidpl, 'LEVEL', varidpl)
[5084]1801              IF (rcode.NE.0) rcode = nf90_inq_varid(ncidpl, 'PRESSURE', varidpl)
[5279]1802              write(*,*) trim(modname)//' ncidpl,varidpl',ncidpl,varidpl
[1632]1803         endif
1804! Coefs ap, bp pour calcul de la pression aux differents niveaux
[5084]1805         IF (guide_plevs.EQ.1) THEN
[5249]1806             status = nf90_put_var(ncidpl, varidap, apnc, [1], [nlevnc])
1807             status = nf90_put_var(ncidpl, varidbp, bpnc, [1], [nlevnc])
[5084]1808         ELSEIF (guide_plevs.EQ.0) THEN
[5249]1809             status = nf90_put_var(ncidpl, varidpl, apnc, [1], [nlevnc])
[3795]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
[1632]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
[5249]1833! Pression
[5084]1834     if (guide_plevs.EQ.2) then
[5249]1835         status = nf90_put_var(ncidp, varidp, pnat2, start, count)
[1632]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
[5249]1845         status = nf90_put_var(ncidu, varidu, unat2, start, count)
[1632]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
[2028]1854
[1632]1855!  Temperature
1856     if (guide_T) then
[5249]1857         status = nf90_put_var(ncidt, varidt, tnat2, start, count)
[1632]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
[5249]1867         status = nf90_put_var(ncidQ, varidQ, qnat2, start, count)
[1632]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
[5249]1881         status = nf90_put_var(ncidv, varidv, vnat2, start, count)
[1632]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
[5084]1890     if ((guide_P).OR.(guide_plevs.EQ.1))  then
[1632]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
[5249]1898         status = nf90_put_var(ncidps, varidps, psnat2, start, count)
[1632]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)
[5249]1910    USE netcdf, ONLY: nf90_put_var
[5271]1911    USE dimensions_mod, ONLY: iim, jjm, llm, ndm
[5285]1912USE paramet_mod_h
[5271]1913IMPLICIT NONE
1914
[1632]1915
[5272]1916
[1632]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
[1906]1932    CHARACTER (len = 20)   :: modname = 'guide_read2D'
[1632]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
[3995]1940         write(*,*)trim(modname)//' : opening nudging files '
[1632]1941! Ap et Bp si niveaux de pression hybrides
[5084]1942         if (guide_plevs.EQ.1) then
[3995]1943           write(*,*)trim(modname)//' Reading nudging on model levels'
1944           rcode = nf90_open('apbp.nc', nf90_nowrite, ncidpl)
[5270]1945           IF (rcode.NE.nf90_noerr) THEN
[3995]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)
[5270]1950           IF (rcode.NE.nf90_noerr) THEN
[3995]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)
[5270]1955           IF (rcode.NE.nf90_noerr) THEN
[3995]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
[1632]1960         endif
1961! Pression
[5084]1962         if (guide_plevs.EQ.2) then
[3995]1963           rcode = nf90_open('P.nc', nf90_nowrite, ncidp)
[5270]1964           IF (rcode.NE.nf90_noerr) THEN
[3995]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)
[5270]1969           IF (rcode.NE.nf90_noerr) THEN
[3995]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
[5084]1974           if (ncidpl.eq.-99) ncidpl=ncidp
[1632]1975         endif
1976! Vent zonal
1977         if (guide_u) then
[3995]1978           rcode = nf90_open('u.nc', nf90_nowrite, ncidu)
[5270]1979           IF (rcode.NE.nf90_noerr) THEN
[3995]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)
[5270]1984           IF (rcode.NE.nf90_noerr) THEN
[3995]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
[5084]1989           if (ncidpl.eq.-99) ncidpl=ncidu
[1632]1990         endif
[2028]1991
[1632]1992! Vent meridien
1993         if (guide_v) then
[3995]1994           rcode = nf90_open('v.nc', nf90_nowrite, ncidv)
[5270]1995           IF (rcode.NE.nf90_noerr) THEN
[3995]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)
[5270]2000           IF (rcode.NE.nf90_noerr) THEN
[3995]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
[5084]2005           if (ncidpl.eq.-99) ncidpl=ncidv
[3995]2006        endif
[1632]2007! Temperature
2008         if (guide_T) then
[3995]2009           rcode = nf90_open('T.nc', nf90_nowrite, ncidt)
[5270]2010           IF (rcode.NE.nf90_noerr) THEN
[3995]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)
[5270]2015           IF (rcode.NE.nf90_noerr) THEN
[3995]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
[5084]2020           if (ncidpl.eq.-99) ncidpl=ncidt
[1632]2021         endif
2022! Humidite
2023         if (guide_Q) then
[3995]2024           rcode = nf90_open('hur.nc', nf90_nowrite, ncidQ)
[5270]2025           IF (rcode.NE.nf90_noerr) THEN
[3995]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)
[5270]2030           IF (rcode.NE.nf90_noerr) THEN
[3995]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
[5084]2035           if (ncidpl.eq.-99) ncidpl=ncidQ
[1632]2036         endif
2037! Pression de surface
[5084]2038         if ((guide_P).OR.(guide_plevs.EQ.1)) then
[3995]2039           rcode = nf90_open('ps.nc', nf90_nowrite, ncidps)
[5270]2040           IF (rcode.NE.nf90_noerr) THEN
[3995]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)
[5270]2045           IF (rcode.NE.nf90_noerr) THEN
[3995]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
[1632]2050         endif
2051! Coordonnee verticale
[5084]2052         if (guide_plevs.EQ.0) then
[3995]2053           rcode = nf90_inq_varid(ncidpl, 'LEVEL', varidpl)
[5084]2054           IF (rcode.NE.0) rcode = nf90_inq_varid(ncidpl, 'PRESSURE', varidpl)
[3995]2055           write(*,*)trim(modname)//' ncidpl,varidpl',ncidpl,varidpl
[1632]2056         endif
2057! Coefs ap, bp pour calcul de la pression aux differents niveaux
[5084]2058         if (guide_plevs.EQ.1) then
[5249]2059             status = nf90_put_var(ncidpl, varidap, apnc, [1], [nlevnc])
2060             status = nf90_put_var(ncidpl, varidbp, bpnc, [1], [nlevnc])
[5084]2061         elseif (guide_plevs.EQ.0) THEN
[5249]2062             status = nf90_put_var(ncidpl, varidpl, apnc, [1], [nlevnc])
[1632]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
[5084]2086     if (guide_plevs.EQ.2) then
[5249]2087         status = nf90_put_var(ncidp, varidp, zu, start, count)
[1632]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
[5249]2100         status = nf90_put_var(ncidu, varidu, zu, start, count)
[1632]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
[2028]2112
[1632]2113!  Temperature
2114     if (guide_T) then
[5249]2115         status = nf90_put_var(ncidt, varidt, zu, start, count)
[1632]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
[5249]2129         status = nf90_put_var(ncidQ, varidQ, zu, start, count)
[1632]2130         DO i=1,iip1
2131             qnat2(i,:,:)=zu(:,:)
2132         ENDDO
[5249]2133
[1632]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
[5249]2146         status = nf90_put_var(ncidv, varidv, zv, start, count)
[1632]2147         DO i=1,iip1
2148             vnat2(i,:,:)=zv(:,:)
2149         ENDDO
2150
2151         IF (invert_y) THEN
[5249]2152
[1632]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
[5084]2160     if ((guide_P).OR.(guide_plevs.EQ.1))  then
[1632]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
[5249]2168         status = nf90_put_var(ncidps, varidps, zu(:, 1), start, count)
[1632]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
[5249]2181
[1632]2182!=======================================================================
[2028]2183  SUBROUTINE guide_out(varname,hsize,vsize,field_loc,factt)
[1823]2184    USE parallel_lmdz
[2028]2185    USE mod_hallo, ONLY : gather_field_u, gather_field_v
[2597]2186    USE comconst_mod, ONLY: pi
[2600]2187    USE comvert_mod, ONLY: presnivs
[2740]2188    use netcdf95, only: nf95_def_var, nf95_put_var
[5249]2189    use netcdf, only: nf90_float, nf90_put_var
[2740]2190
[5271]2191    USE dimensions_mod, ONLY: iim, jjm, llm, ndm
[5285]2192    USE paramet_mod_h
[5281]2193    USE comgeom2_mod_h
[5271]2194IMPLICIT NONE
[1632]2195
[5271]2196
[5272]2197
[5249]2198
[1632]2199    ! Variables entree
[2028]2200    CHARACTER*(*), INTENT(IN)                      :: varname
[1632]2201    INTEGER,   INTENT (IN)                         :: hsize,vsize
[2036]2202!   REAL, DIMENSION (iip1,hsize,vsize), INTENT(IN) :: field_loc
2203    REAL, DIMENSION (:,:), INTENT(IN) :: field_loc
2204    REAL factt
[1632]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
[2740]2211    INTEGER       :: vid_au,vid_av, varid_alpha_t, varid_alpha_q
[1632]2212    INTEGER, DIMENSION (3) :: dim3
2213    INTEGER, DIMENSION (4) :: dim4,count,start
[2028]2214    INTEGER                :: ierr, varid,l
[2740]2215    REAL zu(ip1jmp1),zv(ip1jm), zt(iip1, jjp1), zq(iip1, jjp1)
[2036]2216    REAL, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: field_glo
[3995]2217    CHARACTER(LEN=20),PARAMETER :: modname="guide_out"
[5249]2218
[2028]2219!$OMP MASTER
[2036]2220    ALLOCATE(field_glo(iip1,hsize,vsize))
2221!$OMP END MASTER
2222!$OMP BARRIER
[2028]2223
[3995]2224!    write(*,*)trim(modname)//' after allocation ',hsize,vsize
[2036]2225
[2028]2226    IF (hsize==jjp1) THEN
[2036]2227        CALL gather_field_u(field_loc,field_glo,vsize)
[2028]2228    ELSE IF (hsize==jjm) THEN
2229       CALL gather_field_v(field_loc,field_glo, vsize)
2230    ENDIF
2231
[3995]2232!    write(*,*)trim(modname)//' after gather '
[2028]2233    CALL Gather_field_u(alpha_u,zu,1)
[2740]2234    CALL Gather_field_u(alpha_t,zt,1)
2235    CALL Gather_field_u(alpha_q,zq,1)
[2028]2236    CALL Gather_field_v(alpha_v,zv,1)
[2036]2237
2238    IF (mpi_rank >  0) THEN
2239!$OMP MASTER
2240       DEALLOCATE(field_glo)
[2034]2241!$OMP END MASTER
2242!$OMP BARRIER
[2028]2243
[2036]2244       RETURN
2245    ENDIF
[5249]2246
[2034]2247!$OMP MASTER
[5249]2248    IF (timestep.EQ.0) THEN
[1632]2249! ----------------------------------------------
2250! initialisation fichier de sortie
2251! ----------------------------------------------
2252! Ouverture du fichier
[5270]2253        ierr=nf90_create("guide_ins.nc",IOR(nf90_clobber,nf90_64bit_offset),nid)
[1632]2254! Definition des dimensions
[5270]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)
[1632]2261
2262! Creation des variables dimensions
[5270]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)
[2740]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)
[5249]2276
[5270]2277        ierr=nf90_enddef(nid)
[1632]2278
2279! Enregistrement des variables dimensions
[5249]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)
[2740]2289        call nf95_put_var(nid, varid_alpha_t, zt)
2290        call nf95_put_var(nid, varid_alpha_q, zq)
[1632]2291! --------------------------------------------------------------------
2292! Cr�ation des variables sauvegard�es
2293! --------------------------------------------------------------------
[5270]2294        ierr = nf90_redef(nid)
[1632]2295! Pressure (GCM)
2296        dim4=(/id_lonv,id_latu,id_lev,id_tim/)
[5270]2297        ierr = nf90_def_var(nid,"SP",nf90_float,dim4,varid)
[1632]2298! Surface pressure (guidage)
2299        IF (guide_P) THEN
2300            dim3=(/id_lonv,id_latu,id_tim/)
[5270]2301            ierr = nf90_def_var(nid,"ps",nf90_float,dim3,varid)
[1632]2302        ENDIF
2303! Zonal wind
2304        IF (guide_u) THEN
2305            dim4=(/id_lonu,id_latu,id_lev,id_tim/)
[5270]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)
[1632]2309        ENDIF
2310! Merid. wind
2311        IF (guide_v) THEN
2312            dim4=(/id_lonv,id_latv,id_lev,id_tim/)
[5270]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)
[1632]2316        ENDIF
2317! Pot. Temperature
2318        IF (guide_T) THEN
2319            dim4=(/id_lonv,id_latu,id_lev,id_tim/)
[5270]2320            ierr = nf90_def_var(nid,"teta",nf90_float,dim4,varid)
[1632]2321        ENDIF
2322! Specific Humidity
2323        IF (guide_Q) THEN
2324            dim4=(/id_lonv,id_latu,id_lev,id_tim/)
[5270]2325            ierr = nf90_def_var(nid,"q",nf90_float,dim4,varid)
[1632]2326        ENDIF
[5249]2327
[5270]2328        ierr = nf90_enddef(nid)
2329        ierr = nf90_close(nid)
[1632]2330    ENDIF ! timestep=0
2331
2332! --------------------------------------------------------------------
2333! Enregistrement du champ
2334! --------------------------------------------------------------------
[5249]2335
[5270]2336    ierr=nf90_open("guide_ins.nc",nf90_write,nid)
[1632]2337
[2028]2338    IF (varname=="SP") timestep=timestep+1
2339
[5270]2340    ierr = nf90_inq_varid(nid,varname,varid)
[1632]2341    SELECT CASE (varname)
[2028]2342    CASE ("SP","ps")
[1632]2343        start=(/1,1,1,timestep/)
2344        count=(/iip1,jjp1,llm,1/)
[2028]2345    CASE ("v","va","vcov")
[1632]2346        start=(/1,1,1,timestep/)
2347        count=(/iip1,jjm,llm,1/)
[2028]2348    CASE DEFAULT
[1632]2349        start=(/1,1,1,timestep/)
2350        count=(/iip1,jjp1,llm,1/)
[2028]2351    END SELECT
2352
[2036]2353!$OMP END MASTER
2354!$OMP BARRIER
2355
[2028]2356    SELECT CASE (varname)
2357
2358    CASE("u","ua")
[2036]2359!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[5249]2360        DO l=1,llm
[2039]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
[2028]2364    CASE("v","va")
[2036]2365!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[2039]2366        DO l=1,llm
2367           field_glo(:,:,l)=field_glo(:,:,l)/cv(:,:)
2368        ENDDO
[2028]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
[2036]2375
2376!$OMP MASTER
2377
[5249]2378    ierr = nf90_put_var(nid, varid, field_glo, start, count)
[5270]2379    ierr = nf90_close(nid)
[1632]2380
[2036]2381       DEALLOCATE(field_glo)
[2028]2382!$OMP END MASTER
2383!$OMP BARRIER
2384
[1632]2385  END SUBROUTINE guide_out
[5249]2386
2387
[1632]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
[5084]2397            if(abs(x(i,l)).gt.1.e10) then
[1632]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
[2028]2407
2408!====================================================================
2409! Ascii debug output. Could be reactivated
2410!====================================================================
2411
2412subroutine dump2du(var,varname)
[5281]2413     use parallel_lmdz
[2028]2414use mod_hallo
[5271]2415USE dimensions_mod, ONLY: iim, jjm, llm, ndm
[5285]2416USE paramet_mod_h
[2028]2417implicit none
[5271]2418
[2028]2419
[5272]2420
[2028]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
[5271]2447     USE dimensions_mod, ONLY: iim, jjm, llm, ndm
[5285]2448     USE paramet_mod_h
[5281]2449     USE comgeom_mod_h
[5271]2450implicit none
2451
[5272]2452
[2028]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
[1632]2460!===========================================================================
2461END MODULE guide_loc_mod
Note: See TracBrowser for help on using the repository browser.