source: LMDZ6/trunk/libf/dyn3d/guide_mod.F90 @ 5104

Last change on this file since 5104 was 5084, checked in by Laurent Fairhead, 12 months ago

Reverting to r4065. Updating fortran standard broke too much stuff. Will do it by smaller chunks
AB, LF

  • 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: 65.6 KB
RevLine 
[1170]1!
[1186]2! $Id$
[1170]3!
4MODULE guide_mod
5
6!=======================================================================
7!   Auteur:  F.Hourdin
8!            F. Codron 01/09
9!=======================================================================
10
[5084]11  USE getparam, only: ini_getparam, fin_getparam, getpar
[1170]12  USE Write_Field
[5084]13  use netcdf, only: nf90_nowrite, nf90_open, nf90_inq_varid, nf90_close, &
14                    nf90_inq_dimid, nf90_inquire_dimension
15  use pres2lev_mod, only: pres2lev
[1170]16
17  IMPLICIT NONE
18
19! ---------------------------------------------
20! Declarations des cles logiques et parametres
21! ---------------------------------------------
22  INTEGER, PRIVATE, SAVE  :: iguide_read,iguide_int,iguide_sav
[3995]23  INTEGER, PRIVATE, SAVE  :: nlevnc, guide_plevs
[1170]24  LOGICAL, PRIVATE, SAVE  :: guide_u,guide_v,guide_T,guide_Q,guide_P
25  LOGICAL, PRIVATE, SAVE  :: guide_hr,guide_teta 
26  LOGICAL, PRIVATE, SAVE  :: guide_BL,guide_reg,guide_add,gamma4,guide_zon
[3995]27  LOGICAL, PRIVATE, SAVE  :: invert_p,invert_y,ini_anal
28  LOGICAL, PRIVATE, SAVE  :: guide_2D,guide_sav,guide_modele
29!FC
30  LOGICAL, PRIVATE, SAVE  :: convert_Pa
[1170]31 
32  REAL, PRIVATE, SAVE     :: tau_min_u,tau_max_u
33  REAL, PRIVATE, SAVE     :: tau_min_v,tau_max_v
34  REAL, PRIVATE, SAVE     :: tau_min_T,tau_max_T
35  REAL, PRIVATE, SAVE     :: tau_min_Q,tau_max_Q
36  REAL, PRIVATE, SAVE     :: tau_min_P,tau_max_P
37
38  REAL, PRIVATE, SAVE     :: lat_min_g,lat_max_g
39  REAL, PRIVATE, SAVE     :: lon_min_g,lon_max_g
40  REAL, PRIVATE, SAVE     :: tau_lon,tau_lat
41
[4246]42  REAL, PRIVATE, SAVE     :: plim_guide_BL
43
[1170]44  REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE     :: alpha_u,alpha_v
[2740]45  REAL, ALLOCATABLE, DIMENSION(:, :), PRIVATE, SAVE     :: alpha_T,alpha_Q
[1170]46  REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE     :: alpha_P,alpha_pcor
47 
48! ---------------------------------------------
49! Variables de guidage
50! ---------------------------------------------
51! Variables des fichiers de guidage
52  REAL, ALLOCATABLE, DIMENSION(:,:,:), PRIVATE, SAVE   :: unat1,unat2
53  REAL, ALLOCATABLE, DIMENSION(:,:,:), PRIVATE, SAVE   :: vnat1,vnat2
54  REAL, ALLOCATABLE, DIMENSION(:,:,:), PRIVATE, SAVE   :: tnat1,tnat2
55  REAL, ALLOCATABLE, DIMENSION(:,:,:), PRIVATE, SAVE   :: qnat1,qnat2
[3995]56  REAL, ALLOCATABLE, DIMENSION(:,:,:), PRIVATE, SAVE   :: pnat1,pnat2
[1170]57  REAL, ALLOCATABLE, DIMENSION(:,:),   PRIVATE, SAVE   :: psnat1,psnat2
58  REAL, ALLOCATABLE, DIMENSION(:),     PRIVATE, SAVE   :: apnc,bpnc
59! Variables aux dimensions du modele
60  REAL, ALLOCATABLE, DIMENSION(:,:),   PRIVATE, SAVE   :: ugui1,ugui2
61  REAL, ALLOCATABLE, DIMENSION(:,:),   PRIVATE, SAVE   :: vgui1,vgui2
62  REAL, ALLOCATABLE, DIMENSION(:,:),   PRIVATE, SAVE   :: tgui1,tgui2
63  REAL, ALLOCATABLE, DIMENSION(:,:),   PRIVATE, SAVE   :: qgui1,qgui2
64  REAL, ALLOCATABLE, DIMENSION(:),   PRIVATE, SAVE   :: psgui1,psgui2
65
66CONTAINS
67!=======================================================================
68
69  SUBROUTINE guide_init
70
[5084]71    use netcdf, only: nf90_noerr
[2598]72    USE control_mod, ONLY: day_step
73    USE serre_mod, ONLY: grossismx
[1403]74
[1170]75    IMPLICIT NONE
76 
77    INCLUDE "dimensions.h"
78    INCLUDE "paramet.h"
[5084]79    INCLUDE "netcdf.inc"
[1170]80
81    INTEGER                :: error,ncidpl,rid,rcod
82    CHARACTER (len = 80)   :: abort_message
83    CHARACTER (len = 20)   :: modname = 'guide_init'
[3995]84    CHARACTER (len = 20)   :: namedim
[1170]85
86! ---------------------------------------------
87! Lecture des parametres: 
88! ---------------------------------------------
[2263]89    call ini_getparam("nudging_parameters_out.txt")
[1170]90! Variables guidees
91    CALL getpar('guide_u',.true.,guide_u,'guidage de u')
92    CALL getpar('guide_v',.true.,guide_v,'guidage de v')
93    CALL getpar('guide_T',.true.,guide_T,'guidage de T')
94    CALL getpar('guide_P',.true.,guide_P,'guidage de P')
95    CALL getpar('guide_Q',.true.,guide_Q,'guidage de Q')
96    CALL getpar('guide_hr',.true.,guide_hr,'guidage de Q par H.R')
97    CALL getpar('guide_teta',.false.,guide_teta,'guidage de T par Teta')
98
99    CALL getpar('guide_add',.false.,guide_add,'for�age constant?')
100    CALL getpar('guide_zon',.false.,guide_zon,'guidage moy zonale')
[2134]101    if (guide_zon .and. abs(grossismx - 1.) > 0.01) &
102         call abort_gcm("guide_init", &
103         "zonal nudging requires grid regular in longitude", 1)
[1170]104
105!   Constantes de rappel. Unite : fraction de jour
106    CALL getpar('tau_min_u',0.02,tau_min_u,'Cste de rappel min, u')
107    CALL getpar('tau_max_u', 10.,tau_max_u,'Cste de rappel max, u')
108    CALL getpar('tau_min_v',0.02,tau_min_v,'Cste de rappel min, v')
109    CALL getpar('tau_max_v', 10.,tau_max_v,'Cste de rappel max, v')
110    CALL getpar('tau_min_T',0.02,tau_min_T,'Cste de rappel min, T')
111    CALL getpar('tau_max_T', 10.,tau_max_T,'Cste de rappel max, T')
112    CALL getpar('tau_min_Q',0.02,tau_min_Q,'Cste de rappel min, Q')
113    CALL getpar('tau_max_Q', 10.,tau_max_Q,'Cste de rappel max, Q')
114    CALL getpar('tau_min_P',0.02,tau_min_P,'Cste de rappel min, P')
115    CALL getpar('tau_max_P', 10.,tau_max_P,'Cste de rappel max, P')
116    CALL getpar('gamma4',.false.,gamma4,'Zone sans rappel elargie')
117    CALL getpar('guide_BL',.true.,guide_BL,'guidage dans C.Lim')
[4246]118    CALL getpar('plim_guide_BL',85000.,plim_guide_BL,'BL top presnivs value')
[2249]119
[4246]120
[1170]121! Sauvegarde du for�age
122    CALL getpar('guide_sav',.false.,guide_sav,'sauvegarde guidage')
123    CALL getpar('iguide_sav',4,iguide_sav,'freq. sauvegarde guidage')
124    ! frequences f>0: fx/jour; f<0: tous les f jours; f=0: 1 seule fois.
[5084]125    IF (iguide_sav.GT.0) THEN
[2124]126       iguide_sav=day_step/iguide_sav
127    ELSE if (iguide_sav == 0) then
128       iguide_sav = huge(0)
[2134]129    ELSE
[2124]130       iguide_sav=day_step*iguide_sav
[1170]131    ENDIF
132
133! Guidage regional seulement (sinon constant ou suivant le zoom)
134    CALL getpar('guide_reg',.false.,guide_reg,'guidage regional')
135    CALL getpar('lat_min_g',-90.,lat_min_g,'Latitude mini guidage ')
136    CALL getpar('lat_max_g', 90.,lat_max_g,'Latitude maxi guidage ')
137    CALL getpar('lon_min_g',-180.,lon_min_g,'longitude mini guidage ')
138    CALL getpar('lon_max_g', 180.,lon_max_g,'longitude maxi guidage ')
139    CALL getpar('tau_lat', 5.,tau_lat,'raideur lat guide regional ')
140    CALL getpar('tau_lon', 5.,tau_lon,'raideur lon guide regional ')
141
142! Parametres pour lecture des fichiers
143    CALL getpar('iguide_read',4,iguide_read,'freq. lecture guidage')
[2134]144    CALL getpar('iguide_int',4,iguide_int,'freq. interpolation vert')
[5084]145    IF (iguide_int.EQ.0) THEN
[2134]146        iguide_int=1
[5084]147    ELSEIF (iguide_int.GT.0) THEN
[1170]148        iguide_int=day_step/iguide_int
149    ELSE
150        iguide_int=day_step*iguide_int
151    ENDIF
[3995]152    CALL getpar('guide_plevs',0,guide_plevs,'niveaux pression fichiers guidage')
153    ! Pour compatibilite avec ancienne version avec guide_modele
154    CALL getpar('guide_modele',.false.,guide_modele,'niveaux pression ap+bp*psol')
155    IF (guide_modele) THEN
156        guide_plevs=1
157    ENDIF
158!FC
159    CALL getpar('convert_Pa',.true.,convert_Pa,'Convert Pressure levels in Pa')
160    ! Fin raccord
[1170]161    CALL getpar('ini_anal',.false.,ini_anal,'Etat initial = analyse')
162    CALL getpar('guide_invertp',.true.,invert_p,'niveaux p inverses')
163    CALL getpar('guide_inverty',.true.,invert_y,'inversion N-S')
164    CALL getpar('guide_2D',.false.,guide_2D,'fichier guidage lat-P')
165
[2263]166    call fin_getparam
167   
[1170]168! ---------------------------------------------
169! Determination du nombre de niveaux verticaux
170! des fichiers guidage
171! ---------------------------------------------
172    ncidpl=-99
[5084]173    if (guide_plevs.EQ.1) then
174       if (ncidpl.eq.-99) then
[1902]175          rcod=nf90_open('apbp.nc',Nf90_NOWRITe, ncidpl)
[5084]176          if (rcod.NE.NF90_NOERR) THEN
[3995]177             abort_message=' Nudging error -> no file apbp.nc'
178             CALL abort_gcm(modname,abort_message,1)
[1902]179          endif
180       endif
[5084]181    elseif (guide_plevs.EQ.2) then
182       if (ncidpl.EQ.-99) then
[3995]183          rcod=nf90_open('P.nc',Nf90_NOWRITe,ncidpl)
[5084]184          if (rcod.NE.NF90_NOERR) THEN
[3995]185             abort_message=' Nudging error -> no file P.nc'
186             CALL abort_gcm(modname,abort_message,1)
187          endif
188       endif
189
190    elseif (guide_u) then
[5084]191           if (ncidpl.eq.-99) then
[1902]192               rcod=nf90_open('u.nc',Nf90_NOWRITe,ncidpl)
[5084]193               if (rcod.NE.NF90_NOERR) THEN
[2249]194                  CALL abort_gcm(modname, &
[3995]195                       ' Nudging error -> no file u.nc',1)
[1902]196               endif
197           endif
[3995]198
199    elseif (guide_v) then
[5084]200           if (ncidpl.eq.-99) then
[1902]201               rcod=nf90_open('v.nc',nf90_nowrite,ncidpl)
[5084]202               if (rcod.NE.NF90_NOERR) THEN
[2249]203                  CALL abort_gcm(modname, &
[3995]204                       ' Nudging error -> no file v.nc',1)
[1902]205               endif
206           endif
[3995]207    elseif (guide_T) then
[5084]208           if (ncidpl.eq.-99) then
[1902]209               rcod=nf90_open('T.nc',nf90_nowrite,ncidpl)
[5084]210               if (rcod.NE.NF90_NOERR) THEN
[2249]211                  CALL abort_gcm(modname, &
[3995]212                       ' Nudging error -> no file T.nc',1)
[1902]213               endif
214           endif
[3995]215    elseif (guide_Q) then
[5084]216           if (ncidpl.eq.-99) then
[1902]217               rcod=nf90_open('hur.nc',nf90_nowrite, ncidpl)
[5084]218               if (rcod.NE.NF90_NOERR) THEN
[2249]219                  CALL abort_gcm(modname, &
[3995]220                       ' Nudging error -> no file hur.nc',1)
[1902]221               endif
222           endif
[3995]223
224
[1170]225    endif
226    error=NF_INQ_DIMID(ncidpl,'LEVEL',rid)
[5084]227    IF (error.NE.NF90_NOERR) error=NF_INQ_DIMID(ncidpl,'PRESSURE',rid)
228    IF (error.NE.NF90_NOERR) THEN
[3995]229        CALL abort_gcm(modname,'Nudging: error reading pressure levels',1)
[1170]230    ENDIF
231    error=NF_INQ_DIMLEN(ncidpl,rid,nlevnc)
[3995]232    write(*,*)trim(modname)//' : number of vertical levels nlevnc', nlevnc
[1188]233    rcod = nf90_close(ncidpl)
[1170]234
235! ---------------------------------------------
236! Allocation des variables
237! ---------------------------------------------
[3995]238    abort_message='nudging allocation error'
[1170]239
240    ALLOCATE(apnc(nlevnc), stat = error)
241    IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
242    ALLOCATE(bpnc(nlevnc), stat = error)
243    IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
244    apnc=0.;bpnc=0.
245
246    ALLOCATE(alpha_pcor(llm), stat = error)
247    IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
248    ALLOCATE(alpha_u(ip1jmp1), stat = error)
249    IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
250    ALLOCATE(alpha_v(ip1jm), stat = error)
251    IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
[2740]252    ALLOCATE(alpha_T(iip1, jjp1), stat = error)
[1170]253    IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
[2740]254    ALLOCATE(alpha_Q(iip1, jjp1), stat = error)
[1170]255    IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
256    ALLOCATE(alpha_P(ip1jmp1), stat = error)
257    IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
258    alpha_u=0.;alpha_v=0;alpha_T=0;alpha_Q=0;alpha_P=0
259   
260    IF (guide_u) THEN
261        ALLOCATE(unat1(iip1,jjp1,nlevnc), stat = error)
262        IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
263        ALLOCATE(ugui1(ip1jmp1,llm), stat = error)
264        IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
265        ALLOCATE(unat2(iip1,jjp1,nlevnc), stat = error)
266        IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
267        ALLOCATE(ugui2(ip1jmp1,llm), stat = error)
268        IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
269        unat1=0.;unat2=0.;ugui1=0.;ugui2=0.
270    ENDIF
271
272    IF (guide_T) THEN
273        ALLOCATE(tnat1(iip1,jjp1,nlevnc), stat = error)
274        IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
275        ALLOCATE(tgui1(ip1jmp1,llm), stat = error)
276        IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
277        ALLOCATE(tnat2(iip1,jjp1,nlevnc), stat = error)
278        IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
279        ALLOCATE(tgui2(ip1jmp1,llm), stat = error)
280        IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
281        tnat1=0.;tnat2=0.;tgui1=0.;tgui2=0.
282    ENDIF
283     
284    IF (guide_Q) THEN
285        ALLOCATE(qnat1(iip1,jjp1,nlevnc), stat = error)
286        IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
287        ALLOCATE(qgui1(ip1jmp1,llm), stat = error)
288        IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
289        ALLOCATE(qnat2(iip1,jjp1,nlevnc), stat = error)
290        IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
291        ALLOCATE(qgui2(ip1jmp1,llm), stat = error)
292        IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
293        qnat1=0.;qnat2=0.;qgui1=0.;qgui2=0.
294    ENDIF
295
296    IF (guide_v) THEN
297        ALLOCATE(vnat1(iip1,jjm,nlevnc), stat = error)
298        IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
299        ALLOCATE(vgui1(ip1jm,llm), stat = error)
300        IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
301        ALLOCATE(vnat2(iip1,jjm,nlevnc), stat = error)
302        IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
303        ALLOCATE(vgui2(ip1jm,llm), stat = error)
304        IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
305        vnat1=0.;vnat2=0.;vgui1=0.;vgui2=0.
306    ENDIF
307
[5084]308    IF (guide_plevs.EQ.2) THEN
[3995]309        ALLOCATE(pnat1(iip1,jjp1,nlevnc), stat = error)
310        IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
311        ALLOCATE(pnat2(iip1,jjp1,nlevnc), stat = error)
312        IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
313        pnat1=0.;pnat2=0.;
314    ENDIF
315
[5084]316    IF (guide_P.OR.guide_plevs.EQ.1) THEN
[1170]317        ALLOCATE(psnat1(iip1,jjp1), stat = error)
318        IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
319        ALLOCATE(psnat2(iip1,jjp1), stat = error)
320        IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
321        psnat1=0.;psnat2=0.;
322    ENDIF
323    IF (guide_P) THEN
324        ALLOCATE(psgui2(ip1jmp1), stat = error)
325        IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
326        ALLOCATE(psgui1(ip1jmp1), stat = error)
327        IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
328        psgui1=0.;psgui2=0.
329    ENDIF
330
331! ---------------------------------------------
332!   Lecture du premier etat de guidage.
333! ---------------------------------------------
334    IF (guide_2D) THEN
335        CALL guide_read2D(1)
336    ELSE
337        CALL guide_read(1)
338    ENDIF
339    IF (guide_v) vnat1=vnat2
340    IF (guide_u) unat1=unat2
341    IF (guide_T) tnat1=tnat2
342    IF (guide_Q) qnat1=qnat2
[5084]343    IF (guide_plevs.EQ.2) pnat1=pnat2
344    IF (guide_P.OR.guide_plevs.EQ.1) psnat1=psnat2
[1170]345
346  END SUBROUTINE guide_init
347
348!=======================================================================
349  SUBROUTINE guide_main(itau,ucov,vcov,teta,q,masse,ps)
[1403]350
[3995]351    USE exner_hyb_m, ONLY: exner_hyb
352    USE exner_milieu_m, ONLY: exner_milieu
[2597]353    USE control_mod, ONLY: day_step, iperiod
[3995]354    USE comconst_mod, ONLY: cpp, dtvr, daysec,kappa
355    USE comvert_mod, ONLY: ap, bp, preff, presnivs, pressure_exner
[1170]356 
357    IMPLICIT NONE
358 
359    INCLUDE "dimensions.h"
360    INCLUDE "paramet.h"
[4470]361    INCLUDE "iniprint.h"
[1170]362
[4470]363
[1170]364    ! Variables entree
365    INTEGER,                       INTENT(IN)    :: itau !pas de temps
366    REAL, DIMENSION (ip1jmp1,llm), INTENT(INOUT) :: ucov,teta,q,masse
367    REAL, DIMENSION (ip1jm,llm),   INTENT(INOUT) :: vcov
368    REAL, DIMENSION (ip1jmp1),     INTENT(INOUT) :: ps
369
370    ! Variables locales
371    LOGICAL, SAVE :: first=.TRUE.
372    LOGICAL       :: f_out ! sortie guidage
373    REAL, DIMENSION (ip1jmp1,llm) :: f_add ! var aux: champ de guidage
[3995]374    REAL :: pk(ip1jmp1,llm) ! Exner at mid-layers
375    REAL :: pks(ip1jmp1) ! Exner at the surface
376    REAL :: unskap ! 1./kappa
377    REAL, DIMENSION (ip1jmp1,llmp1) :: p ! Pressure at inter-layers
[1170]378    ! Compteurs temps:
379    INTEGER, SAVE :: step_rea,count_no_rea,itau_test ! lecture guidage
380    REAL          :: ditau, dday_step
381    REAL          :: tau,reste ! position entre 2 etats de guidage
382    REAL, SAVE    :: factt ! pas de temps en fraction de jour
383   
384    INTEGER       :: l
[3995]385    CHARACTER(LEN=20) :: modname="guide_main"
[4470]386    CHARACTER (len = 80)   :: abort_message
[1170]387
[4470]388
[1170]389!-----------------------------------------------------------------------
390! Initialisations au premier passage
391!-----------------------------------------------------------------------
392    IF (first) THEN
393        first=.FALSE.
394        CALL guide_init
395        itau_test=1001
396        step_rea=1
397        count_no_rea=0
398! Calcul des constantes de rappel
399        factt=dtvr*iperiod/daysec
400        call tau2alpha(3,iip1,jjm ,factt,tau_min_v,tau_max_v,alpha_v)
401        call tau2alpha(2,iip1,jjp1,factt,tau_min_u,tau_max_u,alpha_u)
402        call tau2alpha(1,iip1,jjp1,factt,tau_min_T,tau_max_T,alpha_T)
403        call tau2alpha(1,iip1,jjp1,factt,tau_min_P,tau_max_P,alpha_P)
404        call tau2alpha(1,iip1,jjp1,factt,tau_min_Q,tau_max_Q,alpha_Q)
405! correction de rappel dans couche limite
406        if (guide_BL) then
407             alpha_pcor(:)=1.
408        else
409            do l=1,llm
[4246]410               alpha_pcor(l)=(1.+tanh(((plim_guide_BL-presnivs(l))/preff)/0.05))/2.
[1170]411            enddo
412        endif
413! ini_anal: etat initial egal au guidage       
414        IF (ini_anal) THEN
415            CALL guide_interp(ps,teta)
416            IF (guide_u) ucov=ugui2
417            IF (guide_v) vcov=ugui2
418            IF (guide_T) teta=tgui2
419            IF (guide_Q) q=qgui2
420            IF (guide_P) THEN
421                ps=psgui2
422                CALL pression(ip1jmp1,ap,bp,ps,p)
423                CALL massdair(p,masse)
424            ENDIF
425            RETURN
426        ENDIF
427! Verification structure guidage
[3995]428!        IF (guide_u) THEN
429!            CALL writefield('unat',unat1)
430!            CALL writefield('ucov',RESHAPE(ucov,(/iip1,jjp1,llm/)))
431!        ENDIF
432!        IF (guide_T) THEN
433!            CALL writefield('tnat',tnat1)
434!            CALL writefield('teta',RESHAPE(teta,(/iip1,jjp1,llm/)))
435!        ENDIF
[1170]436
437    ENDIF !first
438
439!-----------------------------------------------------------------------
440! Lecture des fichiers de guidage ?
441!-----------------------------------------------------------------------
[5084]442    IF (iguide_read.NE.0) THEN
[1170]443      ditau=real(itau)
444      dday_step=real(day_step)
[5084]445      IF (iguide_read.LT.0) THEN
[1403]446          tau=ditau/dday_step/REAL(iguide_read)
[1170]447      ELSE
[1403]448          tau=REAL(iguide_read)*ditau/dday_step
[1170]449      ENDIF
450      reste=tau-AINT(tau)
[5084]451      IF (reste.EQ.0.) THEN
452          IF (itau_test.EQ.itau) THEN
[4470]453            write(lunout,*)trim(modname)//' second pass in advreel at itau=',&
[3995]454            itau
[4470]455              abort_message='stopped'
456              CALL abort_gcm(modname,abort_message,1)
[1170]457          ELSE
458              IF (guide_v) vnat1=vnat2
459              IF (guide_u) unat1=unat2
460              IF (guide_T) tnat1=tnat2
461              IF (guide_Q) qnat1=qnat2
[5084]462              IF (guide_plevs.EQ.2) pnat1=pnat2
463              IF (guide_P.OR.guide_plevs.EQ.1) psnat1=psnat2
[1170]464              step_rea=step_rea+1
465              itau_test=itau
[3995]466              write(*,*)trim(modname)//' Reading nudging files, step ',&
467                     step_rea,'after ',count_no_rea,' skips'
[1170]468              IF (guide_2D) THEN
469                  CALL guide_read2D(step_rea)
470              ELSE
471                  CALL guide_read(step_rea)
472              ENDIF
473              count_no_rea=0
474          ENDIF
475      ELSE
476        count_no_rea=count_no_rea+1
477
478      ENDIF
479    ENDIF !iguide_read=0
480
481!-----------------------------------------------------------------------
482! Interpolation et conversion des champs de guidage
483!-----------------------------------------------------------------------
[5084]484    IF (MOD(itau,iguide_int).EQ.0) THEN
[1170]485        CALL guide_interp(ps,teta)
486    ENDIF
487! Repartition entre 2 etats de guidage
[5084]488    IF (iguide_read.NE.0) THEN
[1170]489        tau=reste
490    ELSE
491        tau=1.
492    ENDIF
493
494!-----------------------------------------------------------------------
495!   Ajout des champs de guidage
496!-----------------------------------------------------------------------
497! Sauvegarde du guidage?
[5084]498    f_out=((MOD(itau,iguide_sav).EQ.0).AND.guide_sav) 
[3995]499    IF (f_out) THEN
500      ! compute pressures at layer interfaces
501      CALL pression(ip1jmp1,ap,bp,ps,p)
502      if (pressure_exner) then
503        call exner_hyb(ip1jmp1,ps,p,pks,pk)
504      else
505        call exner_milieu(ip1jmp1,ps,p,pks,pk)
506      endif
507      unskap=1./kappa
508      ! Now compute pressures at mid-layer
509      do l=1,llm
510        p(:,l)=preff*(pk(:,l)/cpp)**unskap
511      enddo
512      CALL guide_out("SP",jjp1,llm,p(:,1:llm))
513    ENDIF
[1170]514   
515    if (guide_u) then
516        if (guide_add) then
517           f_add=(1.-tau)*ugui1+tau*ugui2
518        else
519           f_add=(1.-tau)*ugui1+tau*ugui2-ucov
520        endif
521        if (guide_zon) CALL guide_zonave(1,jjp1,llm,f_add)
522        CALL guide_addfield(ip1jmp1,llm,f_add,alpha_u)
[2025]523        IF (f_out) CALL guide_out("ua",jjp1,llm,(1.-tau)*ugui1+tau*ugui2)
524        IF (f_out) CALL guide_out("u",jjp1,llm,ucov)
525        IF (f_out) CALL guide_out("ucov",jjp1,llm,f_add/factt)
[1170]526        ucov=ucov+f_add
527    endif
528
529    if (guide_T) then
530        if (guide_add) then
531           f_add=(1.-tau)*tgui1+tau*tgui2
532        else
533           f_add=(1.-tau)*tgui1+tau*tgui2-teta
534        endif
535        if (guide_zon) CALL guide_zonave(2,jjp1,llm,f_add)
536        CALL guide_addfield(ip1jmp1,llm,f_add,alpha_T)
[2025]537        IF (f_out) CALL guide_out("teta",jjp1,llm,f_add/factt)
[1170]538        teta=teta+f_add
539    endif
540
541    if (guide_P) then
542        if (guide_add) then
543           f_add(1:ip1jmp1,1)=(1.-tau)*psgui1+tau*psgui2
544        else
545           f_add(1:ip1jmp1,1)=(1.-tau)*psgui1+tau*psgui2-ps
546        endif
547        if (guide_zon) CALL guide_zonave(2,jjp1,1,f_add(1:ip1jmp1,1))
548        CALL guide_addfield(ip1jmp1,1,f_add(1:ip1jmp1,1),alpha_P)
[3995]549!        IF (f_out) CALL guide_out("ps",jjp1,1,f_add(1:ip1jmp1,1)/factt)
[1170]550        ps=ps+f_add(1:ip1jmp1,1)
551        CALL pression(ip1jmp1,ap,bp,ps,p)
552        CALL massdair(p,masse)
553    endif
554
555    if (guide_Q) then
556        if (guide_add) then
557           f_add=(1.-tau)*qgui1+tau*qgui2
558        else
559           f_add=(1.-tau)*qgui1+tau*qgui2-q
560        endif
561        if (guide_zon) CALL guide_zonave(2,jjp1,llm,f_add)
562        CALL guide_addfield(ip1jmp1,llm,f_add,alpha_Q)
[2025]563        IF (f_out) CALL guide_out("q",jjp1,llm,f_add/factt)
[1170]564        q=q+f_add
565    endif
566
567    if (guide_v) then
568        if (guide_add) then
569           f_add(1:ip1jm,:)=(1.-tau)*vgui1+tau*vgui2
570        else
571           f_add(1:ip1jm,:)=(1.-tau)*vgui1+tau*vgui2-vcov
572        endif
573        if (guide_zon) CALL guide_zonave(2,jjm,llm,f_add(1:ip1jm,:))
574        CALL guide_addfield(ip1jm,llm,f_add(1:ip1jm,:),alpha_v)
[2025]575        IF (f_out) CALL guide_out("v",jjm,llm,vcov)
576        IF (f_out) CALL guide_out("va",jjm,llm,(1.-tau)*vgui1+tau*vgui2)
577        IF (f_out) CALL guide_out("vcov",jjm,llm,f_add(1:ip1jm,:)/factt)
[1170]578        vcov=vcov+f_add(1:ip1jm,:)
579    endif
580  END SUBROUTINE guide_main
581
582!=======================================================================
583  SUBROUTINE guide_addfield(hsize,vsize,field,alpha)
584! field1=a*field1+alpha*field2
585
586    IMPLICIT NONE
587
588    ! input variables
589    INTEGER,                      INTENT(IN)    :: hsize
590    INTEGER,                      INTENT(IN)    :: vsize
591    REAL, DIMENSION(hsize),       INTENT(IN)    :: alpha
592    REAL, DIMENSION(hsize,vsize), INTENT(INOUT) :: field
593
594    ! Local variables
595    INTEGER :: l
596
597    do l=1,vsize
598        field(:,l)=alpha*field(:,l)*alpha_pcor(l)
599    enddo
600
601  END SUBROUTINE guide_addfield
602
603!=======================================================================
604  SUBROUTINE guide_zonave(typ,hsize,vsize,field)
605
[2597]606    USE comconst_mod, ONLY: pi
607   
[1170]608    IMPLICIT NONE
609
610    INCLUDE "dimensions.h"
611    INCLUDE "paramet.h"
612    INCLUDE "comgeom.h"
613   
614    ! input/output variables
615    INTEGER,                           INTENT(IN)    :: typ
616    INTEGER,                           INTENT(IN)    :: vsize
617    INTEGER,                           INTENT(IN)    :: hsize
618    REAL, DIMENSION(hsize*iip1,vsize), INTENT(INOUT) :: field
619
620    ! Local variables
621    LOGICAL, SAVE                :: first=.TRUE.
622    INTEGER, DIMENSION (2), SAVE :: imin, imax ! averaging domain
623    INTEGER                      :: i,j,l,ij
624    REAL, DIMENSION (iip1)       :: lond       ! longitude in Deg.
625    REAL, DIMENSION (hsize,vsize):: fieldm     ! zon-averaged field
626
627    IF (first) THEN
628        first=.FALSE.
629!Compute domain for averaging
630        lond=rlonu*180./pi
631        imin(1)=1;imax(1)=iip1;
632        imin(2)=1;imax(2)=iip1;
633        IF (guide_reg) THEN
634            DO i=1,iim
[5084]635                IF (lond(i).LT.lon_min_g) imin(1)=i
636                IF (lond(i).LE.lon_max_g) imax(1)=i
[1170]637            ENDDO
638            lond=rlonv*180./pi
639            DO i=1,iim
[5084]640                IF (lond(i).LT.lon_min_g) imin(2)=i
641                IF (lond(i).LE.lon_max_g) imax(2)=i
[1170]642            ENDDO
643        ENDIF
644    ENDIF
645
646    fieldm=0.
647    DO l=1,vsize
648    ! Compute zonal average
649        DO j=1,hsize
650            DO i=imin(typ),imax(typ)
651                ij=(j-1)*iip1+i
652                fieldm(j,l)=fieldm(j,l)+field(ij,l)
653            ENDDO
654        ENDDO
[1403]655        fieldm(:,l)=fieldm(:,l)/REAL(imax(typ)-imin(typ)+1)
[1170]656    ! Compute forcing
657        DO j=1,hsize
658            DO i=1,iip1
659                ij=(j-1)*iip1+i
660                field(ij,l)=fieldm(j,l)
661            ENDDO
662        ENDDO
663    ENDDO
664
665  END SUBROUTINE guide_zonave
666
667!=======================================================================
668  SUBROUTINE guide_interp(psi,teta)
669 
[2021]670  use exner_hyb_m, only: exner_hyb
671  use exner_milieu_m, only: exner_milieu
[2597]672  use comconst_mod, only: kappa, cpp
[2600]673  use comvert_mod, only: preff, pressure_exner, bp, ap
[1170]674  IMPLICIT NONE
675
676  include "dimensions.h"
677  include "paramet.h"
678  include "comgeom2.h"
679
680  REAL, DIMENSION (iip1,jjp1),     INTENT(IN) :: psi ! Psol gcm
681  REAL, DIMENSION (iip1,jjp1,llm), INTENT(IN) :: teta ! Temp. Pot. gcm
682
683  LOGICAL, SAVE                      :: first=.TRUE.
684  ! Variables pour niveaux pression:
685  REAL, DIMENSION (iip1,jjp1,nlevnc) :: plnc1,plnc2 !niveaux pression guidage
686  REAL, DIMENSION (iip1,jjp1,llm)    :: plunc,plsnc !niveaux pression modele
687  REAL, DIMENSION (iip1,jjm,llm)     :: plvnc       !niveaux pression modele
688  REAL, DIMENSION (iip1,jjp1,llmp1)  :: p           ! pression intercouches
689  REAL, DIMENSION (iip1,jjp1,llm)    :: pls, pext   ! var intermediaire
690  REAL, DIMENSION (iip1,jjp1,llm)    :: pbarx
691  REAL, DIMENSION (iip1,jjm,llm)     :: pbary
692  ! Variables pour fonction Exner (P milieu couche)
[2021]693  REAL, DIMENSION (iip1,jjp1,llm)    :: pk
[1170]694  REAL, DIMENSION (iip1,jjp1)        :: pks   
695  REAL                               :: prefkap,unskap
696  ! Pression de vapeur saturante
697  REAL, DIMENSION (ip1jmp1,llm)      :: qsat
698  !Variables intermediaires interpolation
699  REAL, DIMENSION (iip1,jjp1,llm)    :: zu1,zu2
700  REAL, DIMENSION (iip1,jjm,llm)     :: zv1,zv2
701 
702  INTEGER                            :: i,j,l,ij
[3995]703  CHARACTER(LEN=20),PARAMETER :: modname="guide_interp"
[1170]704 
[3995]705    write(*,*)trim(modname)//': interpolate nudging variables'
[1170]706! -----------------------------------------------------------------
707! Calcul des niveaux de pression champs guidage
708! -----------------------------------------------------------------
709if (guide_modele) then
710    do i=1,iip1
711        do j=1,jjp1
712            do l=1,nlevnc
713                plnc2(i,j,l)=apnc(l)+bpnc(l)*psnat2(i,j)
714                plnc1(i,j,l)=apnc(l)+bpnc(l)*psnat1(i,j)
715            enddo
716        enddo
717    enddo
718else
719    do i=1,iip1
720        do j=1,jjp1
721            do l=1,nlevnc
722                plnc2(i,j,l)=apnc(l)
723                plnc1(i,j,l)=apnc(l)
724           enddo
725        enddo
726    enddo
727
728endif
729    if (first) then
730        first=.FALSE.
[3995]731        write(*,*)trim(modname)//' : check vertical level order'
732        write(*,*)trim(modname)//' LMDZ :'
[1170]733        do l=1,llm
[3995]734          write(*,*)trim(modname)//' PL(',l,')=',(ap(l)+ap(l+1))/2. &
[1170]735                  +psi(1,jjp1)*(bp(l)+bp(l+1))/2.
736        enddo
[3995]737        write(*,*)trim(modname)//' nudging file :'
[1170]738        do l=1,nlevnc
[3995]739          write(*,*)trim(modname)//' PL(',l,')=',plnc2(1,1,l)
[1170]740        enddo
[3995]741        write(*,*)trim(modname)//' invert ordering: invert_p=',invert_p
[1170]742        if (guide_u) then
743            do l=1,nlevnc
[3995]744              write(*,*)trim(modname)//' U(',l,')=',unat2(1,1,l)
[1170]745            enddo
746        endif
747        if (guide_T) then
748            do l=1,nlevnc
[3995]749              write(*,*)trim(modname)//' T(',l,')=',tnat2(1,1,l)
[1170]750            enddo
751        endif
752    endif
753   
754! -----------------------------------------------------------------
755! Calcul niveaux pression modele
756! -----------------------------------------------------------------
757    CALL pression( ip1jmp1, ap, bp, psi, p )
[1625]758    if (pressure_exner) then
[2021]759      CALL exner_hyb(ip1jmp1,psi,p,pks,pk)
[1625]760    else
[2021]761      CALL exner_milieu(ip1jmp1,psi,p,pks,pk)
[1520]762    endif
[1170]763!    ....  Calcul de pls , pression au milieu des couches ,en Pascals
764    unskap=1./kappa
765    prefkap =  preff  ** kappa
766    DO l = 1, llm
767        DO j=1,jjp1
768            DO i =1, iip1
769                pls(i,j,l) = preff * ( pk(i,j,l)/cpp) ** unskap
770            ENDDO
771        ENDDO
772    ENDDO
773
774!   calcul des pressions pour les grilles u et v
775    do l=1,llm
776        do j=1,jjp1
777            do i=1,iip1
778                pext(i,j,l)=pls(i,j,l)*aire(i,j)
779            enddo
780        enddo
781    enddo
782    call massbar(pext, pbarx, pbary )
783    do l=1,llm
784        do j=1,jjp1
785            do i=1,iip1
786                plunc(i,j,l)=pbarx(i,j,l)/aireu(i,j)
787                plsnc(i,j,l)=pls(i,j,l)
788            enddo
789        enddo
790    enddo
791    do l=1,llm
792        do j=1,jjm
793            do i=1,iip1
794                plvnc(i,j,l)=pbary(i,j,l)/airev(i,j)
795            enddo
796        enddo
797    enddo
798
799! -----------------------------------------------------------------
800! Interpolation champs guidage sur niveaux modele (+inversion N/S)
801! Conversion en variables gcm (ucov, vcov...)
802! -----------------------------------------------------------------
803    if (guide_P) then
804        do j=1,jjp1
805            do i=1,iim
806                ij=(j-1)*iip1+i
807                psgui1(ij)=psnat1(i,j)
808                psgui2(ij)=psnat2(i,j)
809            enddo
810            psgui1(iip1*j)=psnat1(1,j)
811            psgui2(iip1*j)=psnat2(1,j)
812        enddo
813    endif
814
815    IF (guide_u) THEN
816        CALL pres2lev(unat1,zu1,nlevnc,llm,plnc1,plunc,iip1,jjp1,invert_p)
817        CALL pres2lev(unat2,zu2,nlevnc,llm,plnc2,plunc,iip1,jjp1,invert_p)
818        do l=1,llm
819            do j=1,jjp1
820                do i=1,iim
821                    ij=(j-1)*iip1+i
822                    ugui1(ij,l)=zu1(i,j,l)*cu(i,j)
823                    ugui2(ij,l)=zu2(i,j,l)*cu(i,j)
824                enddo
825                ugui1(j*iip1,l)=ugui1((j-1)*iip1+1,l)   
826                ugui2(j*iip1,l)=ugui2((j-1)*iip1+1,l)   
827            enddo
828            do i=1,iip1
829                ugui1(i,l)=0.
830                ugui1(ip1jm+i,l)=0.
831                ugui2(i,l)=0.
832                ugui2(ip1jm+i,l)=0.
833            enddo
834        enddo
835    ENDIF
836   
837    IF (guide_T) THEN
838        CALL pres2lev(tnat1,zu1,nlevnc,llm,plnc1,plsnc,iip1,jjp1,invert_p)
839        CALL pres2lev(tnat2,zu2,nlevnc,llm,plnc2,plsnc,iip1,jjp1,invert_p)
840        do l=1,llm
841            do j=1,jjp1
842                IF (guide_teta) THEN
[2597]843                    do i=1,iim
844                        ij=(j-1)*iip1+i
845                        tgui1(ij,l)=zu1(i,j,l)
846                        tgui2(ij,l)=zu2(i,j,l)
847                    enddo
[1170]848                ELSE
[2597]849                    do i=1,iim
850                        ij=(j-1)*iip1+i
851                        tgui1(ij,l)=zu1(i,j,l)*cpp/pk(i,j,l)
852                        tgui2(ij,l)=zu2(i,j,l)*cpp/pk(i,j,l)
853                    enddo
[1170]854                ENDIF
855                tgui1(j*iip1,l)=tgui1((j-1)*iip1+1,l)   
856                tgui2(j*iip1,l)=tgui2((j-1)*iip1+1,l)   
857            enddo
858            do i=1,iip1
859                tgui1(i,l)=tgui1(1,l)
860                tgui1(ip1jm+i,l)=tgui1(ip1jm+1,l)
861                tgui2(i,l)=tgui2(1,l)
862                tgui2(ip1jm+i,l)=tgui2(ip1jm+1,l)
863            enddo
864        enddo
865    ENDIF
866
867    IF (guide_v) THEN
868
869        CALL pres2lev(vnat1,zv1,nlevnc,llm,plnc1(:,:jjm,:),plvnc,iip1,jjm,invert_p)
870        CALL pres2lev(vnat2,zv2,nlevnc,llm,plnc2(:,:jjm,:),plvnc,iip1,jjm,invert_p)
871
872        do l=1,llm
873            do j=1,jjm
874                do i=1,iim
875                    ij=(j-1)*iip1+i
876                    vgui1(ij,l)=zv1(i,j,l)*cv(i,j)
877                    vgui2(ij,l)=zv2(i,j,l)*cv(i,j)
878                enddo
879                vgui1(j*iip1,l)=vgui1((j-1)*iip1+1,l)   
880                vgui2(j*iip1,l)=vgui2((j-1)*iip1+1,l)   
881            enddo
882        enddo
883    ENDIF
884   
885    IF (guide_Q) THEN
886        ! On suppose qu'on a la bonne variable dans le fichier de guidage:
887        ! Hum.Rel si guide_hr, Hum.Spec. sinon.
888        CALL pres2lev(qnat1,zu1,nlevnc,llm,plnc1,plsnc,iip1,jjp1,invert_p)
889        CALL pres2lev(qnat2,zu2,nlevnc,llm,plnc2,plsnc,iip1,jjp1,invert_p)
890        do l=1,llm
891            do j=1,jjp1
892                do i=1,iim
893                    ij=(j-1)*iip1+i
894                    qgui1(ij,l)=zu1(i,j,l)
895                    qgui2(ij,l)=zu2(i,j,l)
896                enddo
897                qgui1(j*iip1,l)=qgui1((j-1)*iip1+1,l)   
898                qgui2(j*iip1,l)=qgui2((j-1)*iip1+1,l)   
899            enddo
900            do i=1,iip1
901                qgui1(i,l)=qgui1(1,l)
902                qgui1(ip1jm+i,l)=qgui1(ip1jm+1,l)
903                qgui2(i,l)=qgui2(1,l)
904                qgui2(ip1jm+i,l)=qgui2(ip1jm+1,l)
905            enddo
906        enddo
907        IF (guide_hr) THEN
908            CALL q_sat(iip1*jjp1*llm,teta*pk/cpp,plsnc,qsat)
909            qgui1=qgui1*qsat*0.01 !hum. rel. en %
910            qgui2=qgui2*qsat*0.01
911        ENDIF
912    ENDIF
913
914  END SUBROUTINE guide_interp
915
916!=======================================================================
917  SUBROUTINE tau2alpha(typ,pim,pjm,factt,taumin,taumax,alpha)
918
919! Calcul des constantes de rappel alpha (=1/tau)
920
[2597]921    use comconst_mod, only: pi
[2598]922    use serre_mod, only: clon, clat, grossismx, grossismy
[2597]923   
[1170]924    implicit none
925
926    include "dimensions.h"
927    include "paramet.h"
928    include "comgeom2.h"
929
930! input arguments :
931    INTEGER, INTENT(IN) :: typ    ! u(2),v(3), ou scalaire(1)
932    INTEGER, INTENT(IN) :: pim,pjm ! dimensions en lat, lon
933    REAL, INTENT(IN)    :: factt   ! pas de temps en fraction de jour
934    REAL, INTENT(IN)    :: taumin,taumax
935! output arguments:
936    REAL, DIMENSION(pim,pjm), INTENT(OUT) :: alpha
937 
938!  local variables:
939    LOGICAL, SAVE               :: first=.TRUE.
940    REAL, SAVE                  :: gamma,dxdy_min,dxdy_max
941    REAL, DIMENSION (iip1,jjp1) :: zdx,zdy
942    REAL, DIMENSION (iip1,jjp1) :: dxdys,dxdyu
943    REAL, DIMENSION (iip1,jjm)  :: dxdyv
944    real dxdy_
945    real zlat,zlon
946    real alphamin,alphamax,xi
947    integer i,j,ilon,ilat
[3995]948    character(len=20),parameter :: modname="tau2alpha"
[4470]949    CHARACTER (len = 80)   :: abort_message
[1170]950
951
952    alphamin=factt/taumax
953    alphamax=factt/taumin
954    IF (guide_reg.OR.guide_add) THEN
955        alpha=alphamax
956!-----------------------------------------------------------------------
957! guide_reg: alpha=alpha_min dans region, 0. sinon.
958!-----------------------------------------------------------------------
959        IF (guide_reg) THEN
960            do j=1,pjm
961                do i=1,pim
[5084]962                    if (typ.eq.2) then
[1170]963                       zlat=rlatu(j)*180./pi
964                       zlon=rlonu(i)*180./pi
[5084]965                    elseif (typ.eq.1) then
[1170]966                       zlat=rlatu(j)*180./pi
967                       zlon=rlonv(i)*180./pi
[5084]968                    elseif (typ.eq.3) then
[1170]969                       zlat=rlatv(j)*180./pi
970                       zlon=rlonv(i)*180./pi
971                    endif
972                    alpha(i,j)=alphamax/16.* &
973                              (1.+tanh((zlat-lat_min_g)/tau_lat))* &
974                              (1.+tanh((lat_max_g-zlat)/tau_lat))* &
975                              (1.+tanh((zlon-lon_min_g)/tau_lon))* &
976                              (1.+tanh((lon_max_g-zlon)/tau_lon))
977                enddo
978            enddo
979        ENDIF
980    ELSE
981!-----------------------------------------------------------------------
982! Sinon, alpha varie entre alpha_min et alpha_max suivant le zoom.
983!-----------------------------------------------------------------------
984!Calcul de l'aire des mailles
985        do j=2,jjm
986            do i=2,iip1
987               zdx(i,j)=0.5*(cu(i-1,j)+cu(i,j))/cos(rlatu(j))
988            enddo
989            zdx(1,j)=zdx(iip1,j)
990        enddo
991        do j=2,jjm
992            do i=1,iip1
993               zdy(i,j)=0.5*(cv(i,j-1)+cv(i,j))
994            enddo
995        enddo
996        do i=1,iip1
997            zdx(i,1)=zdx(i,2)
998            zdx(i,jjp1)=zdx(i,jjm)
999            zdy(i,1)=zdy(i,2)
1000            zdy(i,jjp1)=zdy(i,jjm)
1001        enddo
1002        do j=1,jjp1
1003            do i=1,iip1
1004               dxdys(i,j)=sqrt(zdx(i,j)*zdx(i,j)+zdy(i,j)*zdy(i,j))
1005            enddo
1006        enddo
[5084]1007        IF (typ.EQ.2) THEN
[1170]1008            do j=1,jjp1
1009                do i=1,iim
1010                   dxdyu(i,j)=0.5*(dxdys(i,j)+dxdys(i+1,j))
1011                enddo
1012                dxdyu(iip1,j)=dxdyu(1,j)
1013            enddo
1014        ENDIF
[5084]1015        IF (typ.EQ.3) THEN
[1170]1016            do j=1,jjm
1017                do i=1,iip1
1018                   dxdyv(i,j)=0.5*(dxdys(i,j)+dxdys(i,j+1))
1019                enddo
1020            enddo
1021        ENDIF
1022! Premier appel: calcul des aires min et max et de gamma.
1023        IF (first) THEN
1024            first=.FALSE.
1025            ! coordonnees du centre du zoom
1026            CALL coordij(clon,clat,ilon,ilat)
1027            ! aire de la maille au centre du zoom
1028            dxdy_min=dxdys(ilon,ilat)
1029            ! dxdy maximale de la maille
1030            dxdy_max=0.
1031            do j=1,jjp1
1032                do i=1,iip1
1033                     dxdy_max=max(dxdy_max,dxdys(i,j))
1034                enddo
1035            enddo
1036            ! Calcul de gamma
[5084]1037            if (abs(grossismx-1.).lt.0.1.or.abs(grossismy-1.).lt.0.1) then
[3995]1038              write(*,*)trim(modname)//' ATTENTION modele peu zoome'
1039              write(*,*)trim(modname)//' ATTENTION on prend une constante de guidage cste'
1040              gamma=0.
[1170]1041            else
[3995]1042              gamma=(dxdy_max-2.*dxdy_min)/(dxdy_max-dxdy_min)
1043              write(*,*)trim(modname)//' gamma=',gamma
[5084]1044              if (gamma.lt.1.e-5) then
[3995]1045                write(*,*)trim(modname)//' gamma =',gamma,'<1e-5'
[4470]1046                abort_message='stopped'
1047                CALL abort_gcm(modname,abort_message,1)
[3995]1048              endif
1049              gamma=log(0.5)/log(gamma)
1050              if (gamma4) then
1051                gamma=min(gamma,4.)
1052              endif
1053              write(*,*)trim(modname)//' gamma=',gamma
[1170]1054            endif
1055        ENDIF !first
1056
1057        do j=1,pjm
1058            do i=1,pim
[5084]1059                if (typ.eq.1) then
[1170]1060                   dxdy_=dxdys(i,j)
1061                   zlat=rlatu(j)*180./pi
[5084]1062                elseif (typ.eq.2) then
[1170]1063                   dxdy_=dxdyu(i,j)
1064                   zlat=rlatu(j)*180./pi
[5084]1065                elseif (typ.eq.3) then
[1170]1066                   dxdy_=dxdyv(i,j)
1067                   zlat=rlatv(j)*180./pi
1068                endif
[5084]1069                if (abs(grossismx-1.).lt.0.1.or.abs(grossismy-1.).lt.0.1) then
[1170]1070                ! pour une grille reguliere, xi=xxx**0=1 -> alpha=alphamin
1071                    alpha(i,j)=alphamin
1072                else
1073                    xi=((dxdy_max-dxdy_)/(dxdy_max-dxdy_min))**gamma
1074                    xi=min(xi,1.)
[5084]1075                    if(lat_min_g.le.zlat .and. zlat.le.lat_max_g) then
[1170]1076                        alpha(i,j)=xi*alphamin+(1.-xi)*alphamax
1077                    else
1078                        alpha(i,j)=0.
1079                    endif
1080                endif
1081            enddo
1082        enddo
1083    ENDIF ! guide_reg
1084
[2134]1085    if (.not. guide_add) alpha = 1. - exp(- alpha)
1086
[1170]1087  END SUBROUTINE tau2alpha
1088
1089!=======================================================================
1090  SUBROUTINE guide_read(timestep)
[5084]1091
1092    use netcdf, only: NF90_GET_VAR, nf90_noerr
1093
[1170]1094    IMPLICIT NONE
1095
[3995]1096    include "dimensions.h"
1097    include "paramet.h"
[1170]1098
1099    INTEGER, INTENT(IN)   :: timestep
1100
1101    LOGICAL, SAVE         :: first=.TRUE.
1102! Identification fichiers et variables NetCDF:
[3995]1103    INTEGER, SAVE         :: ncidu,varidu,ncidv,varidv,ncidp,varidp
1104    INTEGER, SAVE         :: ncidQ,varidQ,ncidt,varidt,ncidps,varidps
1105    INTEGER               :: ncidpl,varidpl,varidap,varidbp,dimid,lendim
[1170]1106! Variables auxiliaires NetCDF:
1107    INTEGER, DIMENSION(4) :: start,count
1108    INTEGER               :: status,rcode
[1902]1109    CHARACTER (len = 80)   :: abort_message
1110    CHARACTER (len = 20)   :: modname = 'guide_read'
[3995]1111    CHARACTER (len = 20)   :: namedim
1112
[1170]1113! -----------------------------------------------------------------
1114! Premier appel: initialisation de la lecture des fichiers
1115! -----------------------------------------------------------------
1116    if (first) then
1117         ncidpl=-99
[4252]1118         write(*,*) trim(modname)//': opening nudging files '
[1170]1119! Niveaux de pression si non constants
[5084]1120         if (guide_plevs.EQ.1) then
[4252]1121             write(*,*) trim(modname)//' Reading nudging on model levels'
[1188]1122             rcode = nf90_open('apbp.nc', nf90_nowrite, ncidpl)
[5084]1123             IF (rcode.NE.NF90_NOERR) THEN
[3995]1124              abort_message='Nudging: error -> no file apbp.nc'
[1902]1125              CALL abort_gcm(modname,abort_message,1)
1126             ENDIF
[1188]1127             rcode = nf90_inq_varid(ncidpl, 'AP', varidap)
[5084]1128             IF (rcode.NE.NF90_NOERR) THEN
[3995]1129              abort_message='Nudging: error -> no AP variable in file apbp.nc'
[1902]1130              CALL abort_gcm(modname,abort_message,1)
1131             ENDIF
[1188]1132             rcode = nf90_inq_varid(ncidpl, 'BP', varidbp)
[5084]1133             IF (rcode.NE.NF90_NOERR) THEN
[3995]1134              abort_message='Nudging: error -> no BP variable in file apbp.nc'
[1902]1135              CALL abort_gcm(modname,abort_message,1)
1136             ENDIF
[4252]1137             write(*,*) trim(modname)//' ncidpl,varidap',ncidpl,varidap
[1170]1138         endif
[3995]1139
1140! Pression si guidage sur niveaux P variables
[5084]1141         if (guide_plevs.EQ.2) then
[3995]1142             rcode = nf90_open('P.nc', nf90_nowrite, ncidp)
[5084]1143             IF (rcode.NE.NF90_NOERR) THEN
[3995]1144              abort_message='Nudging: error -> no file P.nc'
1145              CALL abort_gcm(modname,abort_message,1)
1146             ENDIF
1147             rcode = nf90_inq_varid(ncidp, 'PRES', varidp)
[5084]1148             IF (rcode.NE.NF90_NOERR) THEN
[3995]1149              abort_message='Nudging: error -> no PRES variable in file P.nc'
1150              CALL abort_gcm(modname,abort_message,1)
1151             ENDIF
[4252]1152             write(*,*) trim(modname)//' ncidp,varidp',ncidp,varidp
[5084]1153             if (ncidpl.eq.-99) ncidpl=ncidp
[3995]1154         endif
1155
[1170]1156! Vent zonal
1157         if (guide_u) then
[1188]1158             rcode = nf90_open('u.nc', nf90_nowrite, ncidu)
[5084]1159             IF (rcode.NE.NF90_NOERR) THEN
[3995]1160              abort_message='Nudging: error -> no file u.nc'
[1902]1161              CALL abort_gcm(modname,abort_message,1)
1162             ENDIF
[1188]1163             rcode = nf90_inq_varid(ncidu, 'UWND', varidu)
[5084]1164             IF (rcode.NE.NF90_NOERR) THEN
[3995]1165              abort_message='Nudging: error -> no UWND variable in file u.nc'
[1902]1166              CALL abort_gcm(modname,abort_message,1)
1167             ENDIF
[4252]1168             write(*,*) trim(modname)//' ncidu,varidu',ncidu,varidu
[5084]1169             if (ncidpl.eq.-99) ncidpl=ncidu
[3995]1170
1171             status=NF90_INQ_DIMID(ncidu, "LONU", dimid)
1172             status=NF90_INQUIRE_DIMENSION(ncidu,dimid,namedim,lendim)
[5084]1173             IF (lendim .NE. iip1) THEN
[3995]1174                abort_message='dimension LONU different from iip1 in u.nc'
1175                CALL abort_gcm(modname,abort_message,1)
1176             ENDIF
1177
1178             status=NF90_INQ_DIMID(ncidu, "LATU", dimid)
1179             status=NF90_INQUIRE_DIMENSION(ncidu,dimid,namedim,lendim)
[5084]1180             IF (lendim .NE. jjp1) THEN
[3995]1181                abort_message='dimension LATU different from jjp1 in u.nc'
1182                CALL abort_gcm(modname,abort_message,1)
1183             ENDIF
1184
[1170]1185         endif
[3995]1186
[1170]1187! Vent meridien
1188         if (guide_v) then
[1188]1189             rcode = nf90_open('v.nc', nf90_nowrite, ncidv)
[5084]1190             IF (rcode.NE.NF90_NOERR) THEN
[3995]1191              abort_message='Nudging: error -> no file v.nc'
[1902]1192              CALL abort_gcm(modname,abort_message,1)
1193             ENDIF
[1188]1194             rcode = nf90_inq_varid(ncidv, 'VWND', varidv)
[5084]1195             IF (rcode.NE.NF90_NOERR) THEN
[3995]1196              abort_message='Nudging: error -> no VWND variable in file v.nc'
[1902]1197              CALL abort_gcm(modname,abort_message,1)
1198             ENDIF
[4252]1199             write(*,*) trim(modname)//' ncidv,varidv',ncidv,varidv
[5084]1200             if (ncidpl.eq.-99) ncidpl=ncidv
[3995]1201             
1202             status=NF90_INQ_DIMID(ncidv, "LONV", dimid)
1203             status=NF90_INQUIRE_DIMENSION(ncidv,dimid,namedim,lendim)
1204             
[5084]1205                IF (lendim .NE. iip1) THEN
[3995]1206                abort_message='dimension LONV different from iip1 in v.nc'
1207                CALL abort_gcm(modname,abort_message,1)
1208             ENDIF
1209
1210
1211             status=NF90_INQ_DIMID(ncidv, "LATV", dimid)
1212             status=NF90_INQUIRE_DIMENSION(ncidv,dimid,namedim,lendim)
[5084]1213             IF (lendim .NE. jjm) THEN
[3995]1214                abort_message='dimension LATV different from jjm in v.nc'
1215                CALL abort_gcm(modname,abort_message,1)
1216             ENDIF
1217       
[1170]1218         endif
[3995]1219
[1170]1220! Temperature
1221         if (guide_T) then
[1188]1222             rcode = nf90_open('T.nc', nf90_nowrite, ncidt)
[5084]1223             IF (rcode.NE.NF90_NOERR) THEN
[3995]1224              abort_message='Nudging: error -> no file T.nc'
[1902]1225              CALL abort_gcm(modname,abort_message,1)
1226             ENDIF
[1188]1227             rcode = nf90_inq_varid(ncidt, 'AIR', varidt)
[5084]1228             IF (rcode.NE.NF90_NOERR) THEN
[3995]1229              abort_message='Nudging: error -> no AIR variable in file T.nc'
[1902]1230              CALL abort_gcm(modname,abort_message,1)
1231             ENDIF
[4252]1232             write(*,*) trim(modname)//' ncidT,varidT',ncidt,varidt
[5084]1233             if (ncidpl.eq.-99) ncidpl=ncidt
[3995]1234
1235             status=NF90_INQ_DIMID(ncidt, "LONV", dimid)
1236             status=NF90_INQUIRE_DIMENSION(ncidt,dimid,namedim,lendim)
[5084]1237             IF (lendim .NE. iip1) THEN
[3995]1238                abort_message='dimension LONV different from iip1 in T.nc'
1239                CALL abort_gcm(modname,abort_message,1)
1240             ENDIF
1241
1242             status=NF90_INQ_DIMID(ncidt, "LATU", dimid)
1243             status=NF90_INQUIRE_DIMENSION(ncidt,dimid,namedim,lendim)
[5084]1244             IF (lendim .NE. jjp1) THEN
[3995]1245                abort_message='dimension LATU different from jjp1 in T.nc'
1246                CALL abort_gcm(modname,abort_message,1)
1247             ENDIF
1248
[1170]1249         endif
[3995]1250
[1170]1251! Humidite
1252         if (guide_Q) then
[1188]1253             rcode = nf90_open('hur.nc', nf90_nowrite, ncidQ)
[5084]1254             IF (rcode.NE.NF90_NOERR) THEN
[3995]1255              abort_message='Nudging: error -> no file hur.nc'
[1902]1256              CALL abort_gcm(modname,abort_message,1)
1257             ENDIF
[1188]1258             rcode = nf90_inq_varid(ncidQ, 'RH', varidQ)
[5084]1259             IF (rcode.NE.NF90_NOERR) THEN
[3995]1260              abort_message='Nudging: error -> no RH variable in file hur.nc'
[1902]1261              CALL abort_gcm(modname,abort_message,1)
1262             ENDIF
[4252]1263             write(*,*) trim(modname)//' ncidQ,varidQ',ncidQ,varidQ
[5084]1264             if (ncidpl.eq.-99) ncidpl=ncidQ
[3995]1265
1266             status=NF90_INQ_DIMID(ncidQ, "LONV", dimid)
1267             status=NF90_INQUIRE_DIMENSION(ncidQ,dimid,namedim,lendim)
[5084]1268             IF (lendim .NE. iip1) THEN
[3995]1269                abort_message='dimension LONV different from iip1 in hur.nc'
1270                CALL abort_gcm(modname,abort_message,1)
1271             ENDIF
1272
1273             status=NF90_INQ_DIMID(ncidQ, "LATU", dimid)
1274             status=NF90_INQUIRE_DIMENSION(ncidQ,dimid,namedim,lendim)
[5084]1275             IF (lendim .NE. jjp1) THEN
[3995]1276                abort_message='dimension LATU different from jjp1 in hur.nc'
1277                CALL abort_gcm(modname,abort_message,1)
1278             ENDIF
1279
[1170]1280         endif
[3995]1281
[1170]1282! Pression de surface
1283         if ((guide_P).OR.(guide_modele)) then
[1188]1284             rcode = nf90_open('ps.nc', nf90_nowrite, ncidps)
[5084]1285             IF (rcode.NE.NF90_NOERR) THEN
[3995]1286              abort_message='Nudging: error -> no file ps.nc'
[1902]1287              CALL abort_gcm(modname,abort_message,1)
1288             ENDIF
[1188]1289             rcode = nf90_inq_varid(ncidps, 'SP', varidps)
[5084]1290             IF (rcode.NE.NF90_NOERR) THEN
[3995]1291              abort_message='Nudging: error -> no SP variable in file ps.nc'
[1902]1292              CALL abort_gcm(modname,abort_message,1)
1293             ENDIF
[4252]1294             write(*,*) trim(modname)//' ncidps,varidps',ncidps,varidps
[1170]1295         endif
1296! Coordonnee verticale
[5084]1297         if (guide_plevs.EQ.0) then
[1188]1298              rcode = nf90_inq_varid(ncidpl, 'LEVEL', varidpl)
[5084]1299              IF (rcode.NE.0) rcode = nf90_inq_varid(ncidpl, 'PRESSURE', varidpl)
[4252]1300              write(*,*) trim(modname)//' ncidpl,varidpl',ncidpl,varidpl
[1170]1301         endif
1302! Coefs ap, bp pour calcul de la pression aux differents niveaux
[5084]1303         if (guide_plevs.EQ.1) then
[4254]1304             status=NF90_GET_VAR(ncidpl,varidap,apnc,[1],[nlevnc])
1305             status=NF90_GET_VAR(ncidpl,varidbp,bpnc,[1],[nlevnc])
[5084]1306         ELSEIF (guide_plevs.EQ.0) THEN
[4254]1307             status=NF90_GET_VAR(ncidpl,varidpl,apnc,[1],[nlevnc])
[3995]1308!FC Pour les corrections la pression est deja en Pascals on commente la ligne ci-dessous
1309             IF(convert_Pa) apnc=apnc*100.! conversion en Pascals
[1170]1310             bpnc(:)=0.
1311         endif
1312         first=.FALSE.
1313     endif ! (first)
1314
1315! -----------------------------------------------------------------
1316!   lecture des champs u, v, T, Q, ps
1317! -----------------------------------------------------------------
1318
1319!  dimensions pour les champs scalaires et le vent zonal
1320     start(1)=1
1321     start(2)=1
1322     start(3)=1
1323     start(4)=timestep
1324
1325     count(1)=iip1
1326     count(2)=jjp1
1327     count(3)=nlevnc
1328     count(4)=1
1329
[3995]1330! Pression
[5084]1331     if (guide_plevs.EQ.2) then
[4254]1332         status=NF90_GET_VAR(ncidp,varidp,pnat2,start,count)
[3995]1333         IF (invert_y) THEN
1334!           PRINT*,"Invertion impossible actuellement"
1335!           CALL abort_gcm(modname,abort_message,1)
1336           CALL invert_lat(iip1,jjp1,nlevnc,pnat2)
1337         ENDIF
1338     endif
1339
[1170]1340!  Vent zonal
1341     if (guide_u) then
[4254]1342         status=NF90_GET_VAR(ncidu,varidu,unat2,start,count)
[1170]1343         IF (invert_y) THEN
[1304]1344           CALL invert_lat(iip1,jjp1,nlevnc,unat2)
[1170]1345         ENDIF
1346     endif
1347
1348!  Temperature
1349     if (guide_T) then
[4254]1350         status=NF90_GET_VAR(ncidt,varidt,tnat2,start,count)
[1170]1351         IF (invert_y) THEN
[1304]1352           CALL invert_lat(iip1,jjp1,nlevnc,tnat2)
[1170]1353         ENDIF
1354     endif
1355
1356!  Humidite
1357     if (guide_Q) then
[4254]1358         status=NF90_GET_VAR(ncidQ,varidQ,qnat2,start,count)
[1170]1359         IF (invert_y) THEN
[1304]1360           CALL invert_lat(iip1,jjp1,nlevnc,qnat2)
[1170]1361         ENDIF
1362         
1363     endif
1364
1365!  Vent meridien
1366     if (guide_v) then
1367         count(2)=jjm
[4254]1368         status=NF90_GET_VAR(ncidv,varidv,vnat2,start,count)
[1170]1369         IF (invert_y) THEN
[1304]1370           CALL invert_lat(iip1,jjm,nlevnc,vnat2)
[1170]1371         ENDIF
1372     endif
1373
1374!  Pression de surface
1375     if ((guide_P).OR.(guide_modele))  then
1376         start(3)=timestep
1377         start(4)=0
1378         count(2)=jjp1
1379         count(3)=1
1380         count(4)=0
[4254]1381         status=NF90_GET_VAR(ncidps,varidps,psnat2,start,count)
[1170]1382         IF (invert_y) THEN
1383           CALL invert_lat(iip1,jjp1,1,psnat2)
1384         ENDIF
1385     endif
1386
1387  END SUBROUTINE guide_read
1388
1389!=======================================================================
1390  SUBROUTINE guide_read2D(timestep)
[5084]1391
1392    use netcdf, only: nf90_get_var, nf90_noerr
1393
[1170]1394    IMPLICIT NONE
1395
[3995]1396    include "dimensions.h"
1397    include "paramet.h"
[1170]1398
1399    INTEGER, INTENT(IN)   :: timestep
1400
1401    LOGICAL, SAVE         :: first=.TRUE.
1402! Identification fichiers et variables NetCDF:
[3995]1403    INTEGER, SAVE         :: ncidu,varidu,ncidv,varidv,ncidp,varidp
1404    INTEGER, SAVE         :: ncidQ,varidQ,ncidt,varidt,ncidps,varidps
[1170]1405    INTEGER               :: ncidpl,varidpl,varidap,varidbp
1406! Variables auxiliaires NetCDF:
1407    INTEGER, DIMENSION(4) :: start,count
1408    INTEGER               :: status,rcode
1409! Variables for 3D extension:
1410    REAL, DIMENSION (jjp1,llm) :: zu
1411    REAL, DIMENSION (jjm,llm)  :: zv
1412    INTEGER               :: i
1413
[1902]1414    CHARACTER (len = 80)   :: abort_message
1415    CHARACTER (len = 20)   :: modname = 'guide_read2D'
[1170]1416! -----------------------------------------------------------------
1417! Premier appel: initialisation de la lecture des fichiers
1418! -----------------------------------------------------------------
1419    if (first) then
1420         ncidpl=-99
[3995]1421         write(*,*)trim(modname)//' : opening nudging files '
1422! Ap et Bp si niveaux de pression hybrides
[5084]1423         if (guide_plevs.EQ.1) then
[3995]1424           write(*,*)trim(modname)//' Reading nudging on model levels'
1425           rcode = nf90_open('apbp.nc', nf90_nowrite, ncidpl)
[5084]1426           IF (rcode.NE.NF90_NOERR) THEN
[3995]1427             abort_message='Nudging: error -> no file apbp.nc'
1428           CALL abort_gcm(modname,abort_message,1)
1429           ENDIF
1430           rcode = nf90_inq_varid(ncidpl, 'AP', varidap)
[5084]1431           IF (rcode.NE.NF90_NOERR) THEN
[3995]1432             abort_message='Nudging: error -> no AP variable in file apbp.nc'
1433           CALL abort_gcm(modname,abort_message,1)
1434           ENDIF
1435           rcode = nf90_inq_varid(ncidpl, 'BP', varidbp)
[5084]1436           IF (rcode.NE.NF90_NOERR) THEN
[3995]1437             abort_message='Nudging: error -> no BP variable in file apbp.nc'
1438             CALL abort_gcm(modname,abort_message,1)
1439           ENDIF
1440           write(*,*)trim(modname)//'ncidpl,varidap',ncidpl,varidap
[1170]1441         endif
[3995]1442! Pression
[5084]1443         if (guide_plevs.EQ.2) then
[3995]1444           rcode = nf90_open('P.nc', nf90_nowrite, ncidp)
[5084]1445           IF (rcode.NE.NF90_NOERR) THEN
[3995]1446             abort_message='Nudging: error -> no file P.nc'
1447             CALL abort_gcm(modname,abort_message,1)
1448           ENDIF
1449           rcode = nf90_inq_varid(ncidp, 'PRES', varidp)
[5084]1450           IF (rcode.NE.NF90_NOERR) THEN
[3995]1451             abort_message='Nudging: error -> no PRES variable in file P.nc'
1452             CALL abort_gcm(modname,abort_message,1)
1453           ENDIF
1454           write(*,*)trim(modname)//' ncidp,varidp',ncidp,varidp
[5084]1455           if (ncidpl.eq.-99) ncidpl=ncidp
[3995]1456         endif
[1170]1457! Vent zonal
1458         if (guide_u) then
[3995]1459           rcode = nf90_open('u.nc', nf90_nowrite, ncidu)
[5084]1460           IF (rcode.NE.NF90_NOERR) THEN
[3995]1461             abort_message='Nudging: error -> no file u.nc'
1462             CALL abort_gcm(modname,abort_message,1)
1463           ENDIF
1464           rcode = nf90_inq_varid(ncidu, 'UWND', varidu)
[5084]1465           IF (rcode.NE.NF90_NOERR) THEN
[3995]1466             abort_message='Nudging: error -> no UWND variable in file u.nc'
1467             CALL abort_gcm(modname,abort_message,1)
1468           ENDIF
1469           write(*,*)trim(modname)//' ncidu,varidu',ncidu,varidu
[5084]1470           if (ncidpl.eq.-99) ncidpl=ncidu
[1170]1471         endif
1472! Vent meridien
1473         if (guide_v) then
[3995]1474           rcode = nf90_open('v.nc', nf90_nowrite, ncidv)
[5084]1475           IF (rcode.NE.NF90_NOERR) THEN
[3995]1476             abort_message='Nudging: error -> no file v.nc'
1477             CALL abort_gcm(modname,abort_message,1)
1478           ENDIF
1479           rcode = nf90_inq_varid(ncidv, 'VWND', varidv)
[5084]1480           IF (rcode.NE.NF90_NOERR) THEN
[3995]1481             abort_message='Nudging: error -> no VWND variable in file v.nc'
1482             CALL abort_gcm(modname,abort_message,1)
1483           ENDIF
1484           write(*,*)trim(modname)//' ncidv,varidv',ncidv,varidv
[5084]1485           if (ncidpl.eq.-99) ncidpl=ncidv
[1170]1486         endif
1487! Temperature
1488         if (guide_T) then
[3995]1489           rcode = nf90_open('T.nc', nf90_nowrite, ncidt)
[5084]1490           IF (rcode.NE.NF90_NOERR) THEN
[3995]1491             abort_message='Nudging: error -> no file T.nc'
1492             CALL abort_gcm(modname,abort_message,1)
1493           ENDIF
1494           rcode = nf90_inq_varid(ncidt, 'AIR', varidt)
[5084]1495           IF (rcode.NE.NF90_NOERR) THEN
[3995]1496             abort_message='Nudging: error -> no AIR variable in file T.nc'
1497             CALL abort_gcm(modname,abort_message,1)
1498           ENDIF
1499           write(*,*)trim(modname)//' ncidT,varidT',ncidt,varidt
[5084]1500           if (ncidpl.eq.-99) ncidpl=ncidt
[1170]1501         endif
1502! Humidite
1503         if (guide_Q) then
[3995]1504           rcode = nf90_open('hur.nc', nf90_nowrite, ncidQ)
[5084]1505           IF (rcode.NE.NF90_NOERR) THEN
[3995]1506             abort_message='Nudging: error -> no file hur.nc'
1507             CALL abort_gcm(modname,abort_message,1)
1508           ENDIF
1509           rcode = nf90_inq_varid(ncidQ, 'RH', varidQ)
[5084]1510           IF (rcode.NE.NF90_NOERR) THEN
[3995]1511             abort_message='Nudging: error -> no RH,variable in file hur.nc'
1512             CALL abort_gcm(modname,abort_message,1)
1513           ENDIF
1514           write(*,*)trim(modname)//' ncidQ,varidQ',ncidQ,varidQ
[5084]1515           if (ncidpl.eq.-99) ncidpl=ncidQ
[1170]1516         endif
1517! Pression de surface
1518         if ((guide_P).OR.(guide_modele)) then
[3995]1519           rcode = nf90_open('ps.nc', nf90_nowrite, ncidps)
[5084]1520           IF (rcode.NE.NF90_NOERR) THEN
[3995]1521             abort_message='Nudging: error -> no file ps.nc'
1522             CALL abort_gcm(modname,abort_message,1)
1523           ENDIF
1524           rcode = nf90_inq_varid(ncidps, 'SP', varidps)
[5084]1525           IF (rcode.NE.NF90_NOERR) THEN
[3995]1526             abort_message='Nudging: error -> no SP variable in file ps.nc'
1527             CALL abort_gcm(modname,abort_message,1)
1528           ENDIF
1529           write(*,*)trim(modname)//' ncidps,varidps',ncidps,varidps
[1170]1530         endif
1531! Coordonnee verticale
[5084]1532         if (guide_plevs.EQ.0) then
[3995]1533           rcode = nf90_inq_varid(ncidpl, 'LEVEL', varidpl)
[5084]1534           IF (rcode.NE.0) rcode = nf90_inq_varid(ncidpl, 'PRESSURE', varidpl)
[3995]1535           write(*,*)trim(modname)//' ncidpl,varidpl',ncidpl,varidpl
[1170]1536         endif
1537! Coefs ap, bp pour calcul de la pression aux differents niveaux
[5084]1538         if (guide_plevs.EQ.1) then
[4254]1539             status=NF90_GET_VAR(ncidpl,varidap,apnc,[1],[nlevnc])
1540             status=NF90_GET_VAR(ncidpl,varidbp,bpnc,[1],[nlevnc])
[5084]1541         elseif (guide_plevs.EQ.0) THEN
[4254]1542             status=NF90_GET_VAR(ncidpl,varidpl,apnc,[1],[nlevnc])
[1170]1543             apnc=apnc*100.! conversion en Pascals
1544             bpnc(:)=0.
1545         endif
1546         first=.FALSE.
1547     endif ! (first)
1548
1549! -----------------------------------------------------------------
1550!   lecture des champs u, v, T, Q, ps
1551! -----------------------------------------------------------------
1552
1553!  dimensions pour les champs scalaires et le vent zonal
1554     start(1)=1
1555     start(2)=1
1556     start(3)=1
1557     start(4)=timestep
1558
1559     count(1)=1
1560     count(2)=jjp1
1561     count(3)=nlevnc
1562     count(4)=1
1563
[3995]1564!  Pression
[5084]1565     if (guide_plevs.EQ.2) then
[4254]1566         status=NF90_GET_VAR(ncidp,varidp,zu,start,count)
[3995]1567         DO i=1,iip1
1568             pnat2(i,:,:)=zu(:,:)
1569         ENDDO
1570
1571         IF (invert_y) THEN
1572!           PRINT*,"Invertion impossible actuellement"
1573!           CALL abort_gcm(modname,abort_message,1)
1574           CALL invert_lat(iip1,jjp1,nlevnc,pnat2)
1575         ENDIF
1576     endif
[1170]1577!  Vent zonal
1578     if (guide_u) then
[4254]1579         status=NF90_GET_VAR(ncidu,varidu,zu,start,count)
[1170]1580         DO i=1,iip1
1581             unat2(i,:,:)=zu(:,:)
1582         ENDDO
1583
1584         IF (invert_y) THEN
[1304]1585           CALL invert_lat(iip1,jjp1,nlevnc,unat2)
[1170]1586         ENDIF
1587
1588     endif
1589
1590!  Temperature
1591     if (guide_T) then
[4254]1592         status=NF90_GET_VAR(ncidt,varidt,zu,start,count)
[1170]1593         DO i=1,iip1
1594             tnat2(i,:,:)=zu(:,:)
1595         ENDDO
1596
1597         IF (invert_y) THEN
[1304]1598           CALL invert_lat(iip1,jjp1,nlevnc,tnat2)
[1170]1599         ENDIF
1600
1601     endif
1602
1603!  Humidite
1604     if (guide_Q) then
[4254]1605         status=NF90_GET_VAR(ncidQ,varidQ,zu,start,count)
[1170]1606         DO i=1,iip1
1607             qnat2(i,:,:)=zu(:,:)
1608         ENDDO
1609
1610         IF (invert_y) THEN
[1304]1611           CALL invert_lat(iip1,jjp1,nlevnc,qnat2)
[1170]1612         ENDIF
1613
1614     endif
1615
1616!  Vent meridien
1617     if (guide_v) then
1618         count(2)=jjm
[4254]1619         status=NF90_GET_VAR(ncidv,varidv,zv,start,count)
[1170]1620         DO i=1,iip1
1621             vnat2(i,:,:)=zv(:,:)
1622         ENDDO
1623
1624         IF (invert_y) THEN
[1304]1625           CALL invert_lat(iip1,jjm,nlevnc,vnat2)
[1170]1626         ENDIF
1627
1628     endif
1629
1630!  Pression de surface
[5084]1631     if ((guide_P).OR.(guide_plevs.EQ.1))  then
[1170]1632         start(3)=timestep
1633         start(4)=0
1634         count(2)=jjp1
1635         count(3)=1
1636         count(4)=0
[4254]1637         status=NF90_GET_VAR(ncidps,varidps,zu(:,1),start,count)
[1170]1638         DO i=1,iip1
1639             psnat2(i,:)=zu(:,1)
1640         ENDDO
1641
1642         IF (invert_y) THEN
1643           CALL invert_lat(iip1,jjp1,1,psnat2)
1644         ENDIF
1645
1646     endif
1647
1648  END SUBROUTINE guide_read2D
1649 
1650!=======================================================================
1651  SUBROUTINE guide_out(varname,hsize,vsize,field)
1652
[2597]1653    USE comconst_mod, ONLY: pi
[2600]1654    USE comvert_mod, ONLY: presnivs
[5084]1655    use netcdf95, only: nf95_def_var, nf95_put_var
1656    use netcdf, only: nf90_float, nf90_def_var
1657   
[1170]1658    IMPLICIT NONE
1659
1660    INCLUDE "dimensions.h"
1661    INCLUDE "paramet.h"
[5084]1662    INCLUDE "netcdf.inc"
[1170]1663    INCLUDE "comgeom2.h"
1664   
1665    ! Variables entree
[2025]1666    CHARACTER*(*), INTENT(IN)                          :: varname
[1170]1667    INTEGER,   INTENT (IN)                         :: hsize,vsize
1668    REAL, DIMENSION (iip1,hsize,vsize), INTENT(IN) :: field
1669
1670    ! Variables locales
1671    INTEGER, SAVE :: timestep=0
1672    ! Identites fichier netcdf
1673    INTEGER       :: nid, id_lonu, id_lonv, id_latu, id_latv, id_tim, id_lev
1674    INTEGER       :: vid_lonu,vid_lonv,vid_latu,vid_latv,vid_cu,vid_cv,vid_lev
[2740]1675    INTEGER       :: vid_au,vid_av, varid_alpha_t, varid_alpha_q
[1170]1676    INTEGER, DIMENSION (3) :: dim3
1677    INTEGER, DIMENSION (4) :: dim4,count,start
[2025]1678    INTEGER                :: ierr, varid,l
1679    REAL, DIMENSION (iip1,hsize,vsize) :: field2
[3995]1680    CHARACTER(LEN=20),PARAMETER :: modname="guide_out"
[1170]1681
[3995]1682    write(*,*)trim(modname)//': output timestep',timestep,'var ',varname
[5084]1683    IF (timestep.EQ.0) THEN
[1170]1684! ----------------------------------------------
1685! initialisation fichier de sortie
1686! ----------------------------------------------
1687! Ouverture du fichier
[3803]1688        ierr=NF_CREATE("guide_ins.nc",IOR(NF_CLOBBER,NF_64BIT_OFFSET),nid)
[1170]1689! Definition des dimensions
1690        ierr=NF_DEF_DIM(nid,"LONU",iip1,id_lonu)
1691        ierr=NF_DEF_DIM(nid,"LONV",iip1,id_lonv)
1692        ierr=NF_DEF_DIM(nid,"LATU",jjp1,id_latu)
1693        ierr=NF_DEF_DIM(nid,"LATV",jjm,id_latv)
1694        ierr=NF_DEF_DIM(nid,"LEVEL",llm,id_lev)
1695        ierr=NF_DEF_DIM(nid,"TIME",NF_UNLIMITED,id_tim)
1696
1697! Creation des variables dimensions
[4253]1698        ierr=NF90_DEF_VAR(nid,"LONU",NF90_FLOAT,id_lonu,vid_lonu)
1699        ierr=NF90_DEF_VAR(nid,"LONV",NF90_FLOAT,id_lonv,vid_lonv)
1700        ierr=NF90_DEF_VAR(nid,"LATU",NF90_FLOAT,id_latu,vid_latu)
1701        ierr=NF90_DEF_VAR(nid,"LATV",NF90_FLOAT,id_latv,vid_latv)
1702        ierr=NF90_DEF_VAR(nid,"LEVEL",NF90_FLOAT,id_lev,vid_lev)
1703        ierr=NF90_DEF_VAR(nid,"cu",NF90_FLOAT,(/id_lonu,id_latu/),vid_cu)
1704        ierr=NF90_DEF_VAR(nid,"cv",NF90_FLOAT,(/id_lonv,id_latv/),vid_cv)
1705        ierr=NF90_DEF_VAR(nid,"au",NF90_FLOAT,(/id_lonu,id_latu/),vid_au)
1706        ierr=NF90_DEF_VAR(nid,"av",NF90_FLOAT,(/id_lonv,id_latv/),vid_av)
[2740]1707        call nf95_def_var(nid, "alpha_T", nf90_float, (/id_lonv, id_latu/), &
1708             varid_alpha_t)
1709        call nf95_def_var(nid, "alpha_q", nf90_float, (/id_lonv, id_latu/), &
1710             varid_alpha_q)
[1170]1711       
1712        ierr=NF_ENDDEF(nid)
1713
1714! Enregistrement des variables dimensions
[5084]1715#ifdef NC_DOUBLE
1716        ierr = NF_PUT_VAR_DOUBLE(nid,vid_lonu,rlonu*180./pi)
1717        ierr = NF_PUT_VAR_DOUBLE(nid,vid_lonv,rlonv*180./pi)
1718        ierr = NF_PUT_VAR_DOUBLE(nid,vid_latu,rlatu*180./pi)
1719        ierr = NF_PUT_VAR_DOUBLE(nid,vid_latv,rlatv*180./pi)
1720        ierr = NF_PUT_VAR_DOUBLE(nid,vid_lev,presnivs)
1721        ierr = NF_PUT_VAR_DOUBLE(nid,vid_cu,cu)
1722        ierr = NF_PUT_VAR_DOUBLE(nid,vid_cv,cv)
1723        ierr = NF_PUT_VAR_DOUBLE(nid,vid_au,alpha_u)
1724        ierr = NF_PUT_VAR_DOUBLE(nid,vid_av,alpha_v)
1725#else
1726        ierr = NF_PUT_VAR_REAL(nid,vid_lonu,rlonu*180./pi)
1727        ierr = NF_PUT_VAR_REAL(nid,vid_lonv,rlonv*180./pi)
1728        ierr = NF_PUT_VAR_REAL(nid,vid_latu,rlatu*180./pi)
1729        ierr = NF_PUT_VAR_REAL(nid,vid_latv,rlatv*180./pi)
1730        ierr = NF_PUT_VAR_REAL(nid,vid_lev,presnivs)
1731        ierr = NF_PUT_VAR_REAL(nid,vid_cu,cu)
1732        ierr = NF_PUT_VAR_REAL(nid,vid_cv,cv)
1733        ierr = NF_PUT_VAR_REAL(nid,vid_au,alpha_u)
1734        ierr = NF_PUT_VAR_REAL(nid,vid_av,alpha_v)
1735#endif
[2740]1736        call nf95_put_var(nid, varid_alpha_t, alpha_t)
1737        call nf95_put_var(nid, varid_alpha_q, alpha_q)
[1170]1738! --------------------------------------------------------------------
1739! Cr�ation des variables sauvegard�es
1740! --------------------------------------------------------------------
1741        ierr = NF_REDEF(nid)
[3995]1742! Pressure (GCM)
1743        dim4=(/id_lonv,id_latu,id_lev,id_tim/)
[4253]1744        ierr = NF90_DEF_VAR(nid,"SP",NF90_FLOAT,dim4,varid)
[1170]1745! Surface pressure (guidage)
1746        IF (guide_P) THEN
1747            dim3=(/id_lonv,id_latu,id_tim/)
[4253]1748            ierr = NF90_DEF_VAR(nid,"ps",NF90_FLOAT,dim3,varid)
[1170]1749        ENDIF
1750! Zonal wind
1751        IF (guide_u) THEN
1752            dim4=(/id_lonu,id_latu,id_lev,id_tim/)
[4253]1753            ierr = NF90_DEF_VAR(nid,"u",NF90_FLOAT,dim4,varid)
1754            ierr = NF90_DEF_VAR(nid,"ua",NF90_FLOAT,dim4,varid)
1755            ierr = NF90_DEF_VAR(nid,"ucov",NF90_FLOAT,dim4,varid)
[1170]1756        ENDIF
1757! Merid. wind
1758        IF (guide_v) THEN
1759            dim4=(/id_lonv,id_latv,id_lev,id_tim/)
[4253]1760            ierr = NF90_DEF_VAR(nid,"v",NF90_FLOAT,dim4,varid)
1761            ierr = NF90_DEF_VAR(nid,"va",NF90_FLOAT,dim4,varid)
1762            ierr = NF90_DEF_VAR(nid,"vcov",NF90_FLOAT,dim4,varid)
[1170]1763        ENDIF
1764! Pot. Temperature
1765        IF (guide_T) THEN
1766            dim4=(/id_lonv,id_latu,id_lev,id_tim/)
[4253]1767            ierr = NF90_DEF_VAR(nid,"teta",NF90_FLOAT,dim4,varid)
[1170]1768        ENDIF
1769! Specific Humidity
1770        IF (guide_Q) THEN
1771            dim4=(/id_lonv,id_latu,id_lev,id_tim/)
[4253]1772            ierr = NF90_DEF_VAR(nid,"q",NF90_FLOAT,dim4,varid)
[1170]1773        ENDIF
1774       
1775        ierr = NF_ENDDEF(nid)
1776        ierr = NF_CLOSE(nid)
1777    ENDIF ! timestep=0
1778
1779! --------------------------------------------------------------------
1780! Enregistrement du champ
1781! --------------------------------------------------------------------
1782    ierr=NF_OPEN("guide_ins.nc",NF_WRITE,nid)
1783
[2025]1784    IF (varname=="SP") timestep=timestep+1
1785
1786    ierr = NF_INQ_VARID(nid,varname,varid)
[1170]1787    SELECT CASE (varname)
[2025]1788    CASE ("SP","ps")
[3995]1789        start=(/1,1,1,timestep/)
1790        count=(/iip1,jjp1,llm,1/)
[2025]1791    CASE ("v","va","vcov")
[1170]1792        start=(/1,1,1,timestep/)
1793        count=(/iip1,jjm,llm,1/)
[2025]1794    CASE DEFAULT
[1170]1795        start=(/1,1,1,timestep/)
1796        count=(/iip1,jjp1,llm,1/)
[2025]1797    END SELECT
1798
1799    SELECT CASE (varname)
1800    CASE("u","ua")
1801        DO l=1,llm ; field2(:,2:jjm,l)=field(:,2:jjm,l)/cu(:,2:jjm) ; ENDDO
1802        field2(:,1,:)=0. ; field2(:,jjp1,:)=0.
1803    CASE("v","va")
1804        DO l=1,llm ; field2(:,:,l)=field(:,:,l)/cv(:,:) ; ENDDO
1805    CASE DEFAULT
1806        field2=field
1807    END SELECT
1808
[5084]1809
1810#ifdef NC_DOUBLE
1811    ierr = NF_PUT_VARA_DOUBLE(nid,varid,start,count,field2)
1812#else
1813    ierr = NF_PUT_VARA_REAL(nid,varid,start,count,field2)
1814#endif
1815
[1170]1816    ierr = NF_CLOSE(nid)
1817
1818  END SUBROUTINE guide_out
1819   
1820 
1821!===========================================================================
1822  subroutine correctbid(iim,nl,x)
1823    integer iim,nl
1824    real x(iim+1,nl)
1825    integer i,l
1826    real zz
1827
1828    do l=1,nl
1829        do i=2,iim-1
[5084]1830            if(abs(x(i,l)).gt.1.e10) then
[1170]1831               zz=0.5*(x(i-1,l)+x(i+1,l))
1832              print*,'correction ',i,l,x(i,l),zz
1833               x(i,l)=zz
1834            endif
1835         enddo
1836     enddo
1837     return
1838  end subroutine correctbid
1839
1840!===========================================================================
1841END MODULE guide_mod
Note: See TracBrowser for help on using the repository browser.