source: LMDZ6/branches/Amaury_dev/libf/dyn3d/guide_mod.F90 @ 5736

Last change on this file since 5736 was 5159, checked in by abarral, 16 months ago

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