source: LMDZ4/branches/LMDZ4-dev-20091210/libf/dyn3dpar/startvar.F @ 1533

Last change on this file since 1533 was 1222, checked in by Ehouarn Millour, 15 years ago

Changes and cleanups to enable compiling without physics
and without ioipsl.

IOIPSL related cleanups:

  • bibio/writehist.F encapsulate the routine (which needs IOIPSL to function)

with #ifdef IOIPSL flag.

  • dyn3d/abort_gcm.F, dyn3dpar/abort_gcm.F and dyn3dpar/getparam.F90: use ioipsl_getincom module when not compiling with IOIPSL library, in order to always be able to use getin() routine.
  • removed unused "use IOIPSL" in dyn3dpar/guide_p_mod.F90
  • calendar related issue: Initialize day_ref and annee_ref in iniacademic.F (i.e. when they are not read from start.nc file)

Earth-specific programs/routines/modules:
create_etat0.F, fluxstokenc.F, limit_netcdf.F, startvar.F
(versions in dyn3d and dyn3dpar)
These routines and modules, which by design and porpose are made to function with
Earth physics are encapsulated with #CPP_EARTH cpp flag.

Earth-specific instructions:

  • calls to qminimum (specific treatment of first 2 tracers, i.e. water) in dyn3d/caladvtrac.F, dyn3d/integrd.F, dyn3dpar/caladvtrac_p.F, dyn3dpar/integrd_p.F only if (planet_type == 'earth')

Interaction with parallel physics:

  • routine dyn3dpar/parallel.F90 uses "surface_data" module (which is in the physics ...) to know value of "type_ocean" . Encapsulated that with #ifdef CPP_EARTH and set to a default type_ocean="dummy" otherwise.
  • So far, only Earth physics are parallelized, so all the interaction between parallel dynamics and parallel physics are encapsulated with #ifdef CCP_EARTH (this way we can run parallel without any physics). The (dyn3dpar) routines which contains such interaction are: bands.F90, gr_dyn_fi_p.F, gr_fi_dyn_p.F, mod_interface_dyn_phys.F90 This should later (when improving dyn/phys interface) be encapsulated with a more general and appropriate #ifdef CPP_PHYS cpp flag.

I checked that these changes do not alter results (on the simple
32x24x11 bench) on Ciclad (seq & mpi), Brodie (seq, mpi & omp) and
Vargas (seq, mpi & omp).

EM

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 39.9 KB
Line 
1!
2! $Id: startvar.F 1222 2009-08-07 11:48:33Z musat $
3!
4      MODULE startvar
5#ifdef CPP_EARTH
6! This module is designed to work for Earth (and with ioipsl)
7    !
8    !
9    !      There are three ways to access data from the database of atmospheric data which
10    !       can be used to initialize the model. This depends on the type of field which needs
11    !       to be extracted. In any case the call should come after a restget and should be of the type :
12    !                CALL startget(...)
13    !
14    !       We will details the possible arguments to startget here :
15    !
16    !        - A 2D variable on the dynamical grid :
17    !           CALL startget(varname, iml, jml, lon_in, lat_in, champ, val_ex, jml2, lon_in2, lat_in2, interbar )             
18    !
19    !        - A 1D variable on the physical grid :
20    !            CALL startget(varname, iml, jml, lon_in, lat_in, nbindex, champ, val_exp, jml2, lon_in2, lat_in2, interbar )
21    !
22    !
23    !         - A 3D variable on the dynamical grid :
24    !            CALL startget(varname, iml, jml, lon_in, lat_in, lml, pls, workvar, champ, val_exp, jml2, lon_in2, lat_in2, interbar )
25    !
26    !
27    !         There is special constraint on the atmospheric data base except that the
28    !         the data needs to be in netCDF and the variables should have the the following
29    !        names in the file :
30    !
31    !      'RELIEF'  : High resolution orography
32    !       'ST'            : Surface temperature
33    !       'CDSW'     : Soil moisture
34    !       'Z'               : Surface geopotential
35    !       'SP'            : Surface pressure
36    !        'U'              : East ward wind
37    !        'V'              : Northward wind
38    !        'TEMP'             : Temperature
39    !        'R'             : Relative humidity
40    !     
41      USE ioipsl
42    !
43    !
44      IMPLICIT NONE
45    !
46    !
47      PRIVATE
48      PUBLIC startget
49    !
50    !
51      INTERFACE startget
52        MODULE PROCEDURE startget_phys2d, startget_phys1d, startget_dyn
53      END INTERFACE
54    !
55      INTEGER, SAVE :: fid_phys, fid_dyn
56      INTEGER, SAVE  :: iml_phys, iml_rel, iml_dyn
57      INTEGER, SAVE :: jml_phys,  jml_rel, jml_dyn
58      INTEGER, SAVE ::  llm_dyn, ttm_dyn
59      REAL, ALLOCATABLE, SAVE, DIMENSION (:,:)  :: lon_phys, lon_rug,
60     . lon_alb, lon_rel, lon_dyn
61      REAL, ALLOCATABLE, SAVE, DIMENSION (:,:)  :: lat_phys, lat_rug,
62     . lat_alb, lat_rel, lat_dyn
63      REAL, ALLOCATABLE, SAVE, DIMENSION (:)  :: levdyn_ini
64      REAL, ALLOCATABLE, SAVE, DIMENSION (:,:)  :: relief, zstd, zsig,
65     . zgam, zthe, zpic, zval
66      REAL, ALLOCATABLE, SAVE, DIMENSION (:,:)  :: rugo, masque, phis
67    !
68      REAL, ALLOCATABLE, SAVE, DIMENSION (:,:)  :: tsol, qsol, psol_dyn
69      REAL, ALLOCATABLE, SAVE, DIMENSION (:,:,:)  ::   var_ana3d
70    !
71      CONTAINS
72    !
73    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
74    !
75      SUBROUTINE startget_phys2d(varname, iml, jml, lon_in, lat_in,
76     . champ, val_exp, jml2, lon_in2, lat_in2 , interbar, masque_lu )
77    !
78    !    There is a big mess with the size in logitude, should it be iml or iml+1.
79    !    I have chosen to use the iml+1 as an argument to this routine and we declare
80    !   internaly smaler fields when needed. This needs to be cleared once and for all in LMDZ.
81    !  A convention is required.
82    !
83    !
84      CHARACTER*(*), INTENT(in) :: varname
85      INTEGER, INTENT(in) :: iml, jml ,jml2
86      REAL, INTENT(in) :: lon_in(iml), lat_in(jml)
87      REAL, INTENT(in) :: lon_in2(iml), lat_in2(jml2)
88      REAL, INTENT(inout) :: champ(iml,jml)
89      REAL, INTENT(in) :: val_exp
90      REAL, INTENT(in), optional :: masque_lu(iml,jml)
91      LOGICAL interbar
92    !
93    !   This routine only works if the variable does not exist or is constant
94    !
95      IF ( MINVAL(champ(:,:)).EQ.MAXVAL(champ(:,:)) .AND.
96     .MINVAL(champ(:,:)).EQ.val_exp ) THEN
97          !
98          SELECTCASE(varname)
99              !
100              CASE ('relief')
101                  !
102                  !  If we do not have the orography we need to get it
103                  !
104                  IF ( .NOT.ALLOCATED(relief)) THEN
105                      !
106                    if (present(masque_lu)) then
107                      CALL start_init_orog( iml, jml, lon_in, lat_in,
108     .                    jml2,lon_in2,lat_in2, interbar, masque_lu )
109                    else
110                      CALL start_init_orog( iml, jml, lon_in, lat_in,
111     .                    jml2,lon_in2,lat_in2, interbar)
112                    endif
113                      !
114                  ENDIF
115                  !
116                  IF (SIZE(relief) .NE. SIZE(champ)) THEN
117                      !
118                      WRITE(*,*) 'STARTVAR module has been',
119     .' initialized to the wrong size'
120                      STOP
121                      !
122                  ENDIF
123                  !
124                  champ(:,:) = relief(:,:)
125                  !
126              CASE ('rugosite')
127                  !
128                  !  If we do not have the orography we need to get it
129                  !
130                  IF ( .NOT.ALLOCATED(rugo)) THEN
131                      !
132                      CALL start_init_orog( iml, jml, lon_in, lat_in,
133     .                    jml2,lon_in2,lat_in2 , interbar )
134                      !
135                  ENDIF
136                  !
137                  IF (SIZE(rugo) .NE. SIZE(champ)) THEN
138                      !
139                      WRITE(*,*)
140     .  'STARTVAR module has been initialized to the wrong size'
141                      STOP
142                      !
143                  ENDIF
144                  !
145                  champ(:,:) = rugo(:,:)
146                  !
147              CASE ('masque')
148                  !
149                  !  If we do not have the orography we need to get it
150                  !
151                  IF ( .NOT.ALLOCATED(masque)) THEN
152                      !
153                      CALL start_init_orog( iml, jml, lon_in, lat_in,
154     .                     jml2,lon_in2,lat_in2 , interbar )
155                      !
156                  ENDIF
157                  !
158                  IF (SIZE(masque) .NE. SIZE(champ)) THEN
159                      !
160                      WRITE(*,*)
161     .   'STARTVAR module has been initialized to the wrong size'
162                      STOP
163                      !
164                  ENDIF
165                  !
166                  champ(:,:) = masque(:,:)
167                  !
168              CASE ('surfgeo')
169                  !
170                  !  If we do not have the orography we need to get it
171                  !
172                  IF ( .NOT.ALLOCATED(phis)) THEN
173                      !
174                      CALL start_init_orog( iml, jml, lon_in, lat_in,
175     .                   jml2,lon_in2, lat_in2 , interbar )
176                      !
177                  ENDIF
178                  !
179                  IF (SIZE(phis) .NE. SIZE(champ)) THEN
180                      !
181                      WRITE(*,*)
182     . 'STARTVAR module has been initialized to the wrong size'
183                      STOP
184                      !
185                  ENDIF
186                  !
187                  champ(:,:) = phis(:,:)
188                  !
189              CASE ('psol')
190                  !
191                  !  If we do not have the orography we need to get it
192                  !
193                  IF ( .NOT.ALLOCATED(psol_dyn)) THEN
194                      !
195                      CALL start_init_dyn( iml, jml, lon_in, lat_in,
196     .                   jml2,lon_in2, lat_in2 , interbar )
197                      !
198                  ENDIF
199                  !
200                  IF (SIZE(psol_dyn) .NE. SIZE(champ)) THEN
201                      !
202                      WRITE(*,*)
203     . 'STARTVAR module has been initialized to the wrong size'
204                      STOP
205                      !
206                  ENDIF
207                  !
208                  champ(:,:) = psol_dyn(:,:)
209                  !
210              CASE DEFAULT
211                  !
212                  WRITE(*,*) 'startget_phys2d'
213                  WRITE(*,*) 'No rule is present to extract variable',
214     .                 varname(:LEN_TRIM(varname)),' from any data set'
215                  STOP
216                  !
217          END SELECT
218          !
219      ELSE
220          !
221          ! There are a few fields we might need if we need to interpolate 3D filed. Thus if they come through here we
222          ! will catch them
223          !
224          SELECTCASE(varname)
225              !
226              CASE ('surfgeo')
227                  !
228                  IF ( .NOT.ALLOCATED(phis)) THEN
229                      ALLOCATE(phis(iml,jml))
230                  ENDIF
231                  !
232                  IF (SIZE(phis) .NE. SIZE(champ)) THEN
233                      !
234                      WRITE(*,*)
235     .  'STARTVAR module has been initialized to the wrong size'
236                      STOP
237                      !
238                  ENDIF
239                  !
240                  phis(:,:) = champ(:,:)
241                  !
242          END SELECT
243          !
244      ENDIF
245    !
246      END SUBROUTINE startget_phys2d
247    !
248    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
249    !
250      SUBROUTINE start_init_orog ( iml,jml,lon_in, lat_in,jml2,lon_in2 ,
251     ,   lat_in2 , interbar, masque_lu )
252    !
253      INTEGER, INTENT(in) :: iml, jml, jml2
254      REAL, INTENT(in) :: lon_in(iml), lat_in(jml)
255      REAL, INTENT(in) :: lon_in2(iml), lat_in2(jml2)
256      REAL, intent(in), optional :: masque_lu(iml,jml)
257      LOGICAL interbar
258    !
259    !  LOCAL
260    !
261      LOGICAL interbar2
262      REAL :: lev(1), date, dt,chmin,chmax
263      INTEGER :: itau(1), fid
264      INTEGER ::  llm_tmp, ttm_tmp
265      INTEGER :: i, j
266      INTEGER :: iret
267      CHARACTER*25 title
268      REAL, ALLOCATABLE :: relief_hi(:,:)
269      REAL, ALLOCATABLE :: lon_rad(:), lat_rad(:)
270      REAL, ALLOCATABLE :: lon_ini(:), lat_ini(:)
271      REAL, ALLOCATABLE :: tmp_var(:,:)
272      INTEGER, ALLOCATABLE :: tmp_int(:,:)
273    !
274      CHARACTER*120 :: orogfname
275      LOGICAL :: check=.TRUE.
276    !
277    !
278      orogfname = 'Relief.nc'
279    !
280      IF ( check ) WRITE(*,*) 'Reading the high resolution orography'
281    !
282      CALL flininfo(orogfname,iml_rel, jml_rel, llm_tmp, ttm_tmp, fid)
283    !
284      ALLOCATE (lat_rel(iml_rel,jml_rel), stat=iret)
285      ALLOCATE (lon_rel(iml_rel,jml_rel), stat=iret)
286      ALLOCATE (relief_hi(iml_rel,jml_rel), stat=iret)
287    !
288      CALL flinopen(orogfname, .FALSE., iml_rel, jml_rel,
289     .llm_tmp, lon_rel, lat_rel, lev, ttm_tmp,
290     .      itau, date, dt, fid)
291    !
292      CALL flinget(fid, 'RELIEF', iml_rel, jml_rel, llm_tmp,
293     . ttm_tmp, 1, 1, relief_hi)
294    !
295      CALL flinclo(fid)
296    !
297    !   In case we have a file which is in degrees we do the transformation
298    !
299      ALLOCATE(lon_rad(iml_rel))
300      ALLOCATE(lon_ini(iml_rel))
301
302      IF ( MAXVAL(lon_rel(:,:)) .GT. 2.0 * ASIN(1.0) ) THEN
303          lon_ini(:) = lon_rel(:,1) * 2.0 * ASIN(1.0) / 180.0
304      ELSE
305          lon_ini(:) = lon_rel(:,1)
306      ENDIF
307
308      ALLOCATE(lat_rad(jml_rel))
309      ALLOCATE(lat_ini(jml_rel))
310
311      IF ( MAXVAL(lat_rel(:,:)) .GT. 2.0 * ASIN(1.0) ) THEN
312          lat_ini(:) = lat_rel(1,:) * 2.0 * ASIN(1.0) / 180.0
313      ELSE
314          lat_ini(:) = lat_rel(1,:)
315      ENDIF
316    !
317    !
318
319      title='RELIEF'
320
321      interbar2 = .FALSE.
322      CALL conf_dat2d(title,iml_rel, jml_rel, lon_ini, lat_ini,
323     . lon_rad, lat_rad, relief_hi , interbar2  )
324
325      IF ( check ) WRITE(*,*) 'Computes all the parameters needed',
326     .' for the gravity wave drag code'
327    !
328    !    Allocate the data we need to put in the interpolated fields
329    !
330    !            RELIEF:  orographie moyenne
331      ALLOCATE(relief(iml,jml))
332    !            zphi :  orographie moyenne
333      ALLOCATE(phis(iml,jml))
334    !             zstd:  deviation standard de l'orographie sous-maille
335      ALLOCATE(zstd(iml,jml))
336    !             zsig:  pente de l'orographie sous-maille
337      ALLOCATE(zsig(iml,jml))
338    !             zgam:  anisotropy de l'orographie sous maille
339      ALLOCATE(zgam(iml,jml))
340    !             zthe:  orientation de l'axe oriente dans la direction
341    !                    de plus grande pente de l'orographie sous maille
342      ALLOCATE(zthe(iml,jml))
343    !             zpic:  hauteur pics de la SSO
344      ALLOCATE(zpic(iml,jml))
345    !             zval:  hauteur vallees de la SSO
346      ALLOCATE(zval(iml,jml))
347    !             masque : Masque terre ocean
348      ALLOCATE(tmp_int(iml,jml))
349      ALLOCATE(masque(iml,jml))
350
351      masque = -99999.
352      if (present(masque_lu)) then
353        masque = masque_lu
354      endif
355    !
356      CALL grid_noro(iml_rel, jml_rel, lon_rad, lat_rad, relief_hi,
357     . iml-1, jml, lon_in, lat_in,
358     . phis, relief, zstd, zsig, zgam, zthe, zpic, zval, masque)
359      phis = phis * 9.81
360    !
361!      masque(:,:) = FLOAT(tmp_int(:,:))
362    !
363    !  Compute surface roughness
364    !
365      IF ( check ) WRITE(*,*)
366     .'Compute surface roughness induced by the orography'
367    !
368      ALLOCATE(rugo(iml,jml))
369      ALLOCATE(tmp_var(iml-1,jml))
370    !
371      CALL rugsoro(iml_rel, jml_rel, lon_rad, lat_rad, relief_hi,
372     . iml-1, jml, lon_in, lat_in, tmp_var)
373    !
374      DO j = 1, jml
375        DO i = 1, iml-1
376          rugo(i,j) = tmp_var(i,j)
377        ENDDO
378        rugo(iml,j) = tmp_var(1,j)
379      ENDDO
380c
381cc   ***   rugo  n'est pas utilise pour l'instant  ******
382    !
383    !   Build land-sea mask
384    !
385    !
386      RETURN
387    !
388      END SUBROUTINE start_init_orog
389    !
390    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
391    !
392      SUBROUTINE startget_phys1d(varname, iml, jml, lon_in,
393     .lat_in, nbindex, champ, val_exp ,jml2, lon_in2, lat_in2,interbar)
394    !
395      CHARACTER*(*), INTENT(in) :: varname
396      INTEGER, INTENT(in) :: iml, jml, nbindex, jml2
397      REAL, INTENT(in) :: lon_in(iml), lat_in(jml)
398      REAL, INTENT(in) :: lon_in2(iml), lat_in2(jml2)
399      REAL, INTENT(inout) :: champ(nbindex)
400      REAL, INTENT(in) :: val_exp
401      LOGICAL interbar
402    !
403    !
404    !   This routine only works if the variable does not exist or is constant
405    !
406      IF ( MINVAL(champ(:)).EQ.MAXVAL(champ(:)) .AND.
407     .MINVAL(champ(:)).EQ.val_exp ) THEN
408          SELECTCASE(varname)
409            CASE ('tsol')
410              IF ( .NOT.ALLOCATED(tsol)) THEN
411                CALL start_init_phys( iml, jml, lon_in, lat_in,
412     .              jml2, lon_in2, lat_in2, interbar )
413              ENDIF
414              IF ( SIZE(tsol) .NE. SIZE(lon_in)*SIZE(lat_in) ) THEN
415                WRITE(*,*)
416     . 'STARTVAR module has been initialized to the wrong size'
417                 STOP
418              ENDIF
419              CALL gr_dyn_fi(1, iml, jml, nbindex, tsol, champ)
420            CASE ('qsol')
421              IF ( .NOT.ALLOCATED(qsol)) THEN
422                CALL start_init_phys( iml, jml, lon_in, lat_in,
423     .              jml2, lon_in2,lat_in2 , interbar )
424              ENDIF
425              IF ( SIZE(qsol) .NE. SIZE(lon_in)*SIZE(lat_in) ) THEN
426                WRITE(*,*)
427     . 'STARTVAR module has been initialized to the wrong size'
428                STOP
429              ENDIF
430              CALL gr_dyn_fi(1, iml, jml, nbindex, qsol, champ)
431            CASE ('psol')
432              IF ( .NOT.ALLOCATED(psol_dyn)) THEN
433                CALL start_init_dyn( iml, jml, lon_in, lat_in,
434     .              jml2, lon_in2,lat_in2 , interbar )
435              ENDIF
436              IF (SIZE(psol_dyn) .NE. SIZE(lon_in)*SIZE(lat_in)) THEN
437                WRITE(*,*)
438     . 'STARTVAR module has been initialized to the wrong size'
439                STOP
440              ENDIF
441              CALL gr_dyn_fi(1, iml, jml, nbindex, psol_dyn, champ)
442            CASE ('zmea')
443              IF ( .NOT.ALLOCATED(relief)) THEN
444                CALL start_init_orog( iml, jml, lon_in, lat_in,
445     .            jml2, lon_in2,lat_in2 , interbar )
446              ENDIF
447              IF ( SIZE(relief) .NE. SIZE(lon_in)*SIZE(lat_in) ) THEN
448                WRITE(*,*)
449     . 'STARTVAR module has been initialized to the wrong size'
450                 STOP
451              ENDIF
452              CALL gr_dyn_fi(1, iml, jml, nbindex, relief, champ)
453            CASE ('zstd')
454              IF ( .NOT.ALLOCATED(zstd)) THEN
455                CALL start_init_orog( iml, jml, lon_in, lat_in,
456     .              jml2, lon_in2,lat_in2 , interbar )
457              ENDIF
458              IF ( SIZE(zstd) .NE. SIZE(lon_in)*SIZE(lat_in) ) THEN
459                WRITE(*,*)
460     . 'STARTVAR module has been initialized to the wrong size'
461                 STOP
462              ENDIF
463              CALL gr_dyn_fi(1, iml, jml, nbindex,zstd, champ)
464            CASE ('zsig')
465              IF ( .NOT.ALLOCATED(zsig)) THEN
466                CALL start_init_orog( iml, jml, lon_in, lat_in,
467     .               jml2, lon_in2,lat_in2 , interbar )
468              ENDIF
469              IF ( SIZE(zsig) .NE. SIZE(lon_in)*SIZE(lat_in) ) THEN
470                WRITE(*,*)
471     . 'STARTVAR module has been initialized to the wrong size'
472                 STOP
473              ENDIF
474              CALL gr_dyn_fi(1, iml, jml, nbindex,zsig, champ)
475            CASE ('zgam')
476              IF ( .NOT.ALLOCATED(zgam)) THEN
477                CALL start_init_orog( iml, jml, lon_in, lat_in,
478     .            jml2, lon_in2,lat_in2 , interbar )
479              ENDIF
480              IF ( SIZE(zgam) .NE. SIZE(lon_in)*SIZE(lat_in) ) THEN
481                WRITE(*,*)
482     . 'STARTVAR module has been initialized to the wrong size'
483                 STOP
484              ENDIF
485              CALL gr_dyn_fi(1, iml, jml, nbindex,zgam, champ)
486            CASE ('zthe')
487              IF ( .NOT.ALLOCATED(zthe)) THEN
488                CALL start_init_orog( iml, jml, lon_in, lat_in,
489     .            jml2, lon_in2,lat_in2 , interbar )
490              ENDIF
491              IF ( SIZE(zthe) .NE. SIZE(lon_in)*SIZE(lat_in) ) THEN
492                WRITE(*,*)
493     . 'STARTVAR module has been initialized to the wrong size'
494                 STOP
495              ENDIF
496              CALL gr_dyn_fi(1, iml, jml, nbindex,zthe, champ)
497            CASE ('zpic')
498              IF ( .NOT.ALLOCATED(zpic)) THEN
499                CALL start_init_orog( iml, jml, lon_in, lat_in,
500     .            jml2, lon_in2,lat_in2 , interbar )
501              ENDIF
502              IF ( SIZE(zpic) .NE. SIZE(lon_in)*SIZE(lat_in) ) THEN
503                WRITE(*,*)
504     . 'STARTVAR module has been initialized to the wrong size'
505                 STOP
506              ENDIF
507              CALL gr_dyn_fi(1, iml, jml, nbindex,zpic, champ)
508            CASE ('zval')
509              IF ( .NOT.ALLOCATED(zval)) THEN
510                CALL start_init_orog( iml, jml, lon_in, lat_in,
511     .            jml2, lon_in2,lat_in2 , interbar )
512              ENDIF
513              IF ( SIZE(zval) .NE. SIZE(lon_in)*SIZE(lat_in) ) THEN
514                WRITE(*,*)
515     . 'STARTVAR module has been initialized to the wrong size'
516                 STOP
517              ENDIF
518              CALL gr_dyn_fi(1, iml, jml, nbindex,zval, champ)
519            CASE ('rads')
520                  champ(:) = 0.0
521            CASE ('snow')
522                  champ(:) = 0.0
523cIM "slab" ocean
524            CASE ('tslab')
525                   champ(:) = 0.0
526            CASE ('seaice')
527                  champ(:) = 0.0
528            CASE ('rugmer')
529                  champ(:) = 0.001
530            CASE ('agsno')
531                  champ(:) = 50.0
532            CASE DEFAULT
533              WRITE(*,*) 'startget_phys1d'
534              WRITE(*,*) 'No rule is present to extract variable  ',
535     . varname(:LEN_TRIM(varname)),' from any data set'
536              STOP
537          END SELECT
538      ELSE
539        !
540        ! If we see tsol we catch it as we may need it for a 3D interpolation
541        !
542        SELECTCASE(varname)
543          CASE ('tsol')
544            IF ( .NOT.ALLOCATED(tsol)) THEN
545              ALLOCATE(tsol(SIZE(lon_in),SIZE(lat_in) ))
546            ENDIF
547            CALL gr_fi_dyn(1, iml, jml, nbindex, champ, tsol)
548        END SELECT
549      ENDIF
550      END SUBROUTINE startget_phys1d
551    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
552    !
553      SUBROUTINE start_init_phys( iml, jml, lon_in, lat_in, jml2,
554     .                 lon_in2, lat_in2 , interbar )
555    !
556      INTEGER, INTENT(in) :: iml, jml ,jml2
557      REAL, INTENT(in) :: lon_in(iml), lat_in(jml)
558      REAL, INTENT(in) :: lon_in2(iml), lat_in2(jml2)
559      LOGICAL interbar
560    !
561    !  LOCAL
562    !
563!ac     REAL :: lev(1), date, dt
564      REAL :: date, dt
565      REAL, DIMENSION(:), ALLOCATABLE :: levphys_ini
566!ac
567      INTEGER :: itau(1)
568      INTEGER ::  llm_tmp, ttm_tmp
569      INTEGER :: i, j
570    !
571      CHARACTER*25 title
572      CHARACTER*120 :: physfname
573      LOGICAL :: check=.TRUE.
574    !
575      REAL, ALLOCATABLE :: lon_rad(:), lat_rad(:)
576      REAL, ALLOCATABLE :: lon_ini(:), lat_ini(:)
577      REAL, ALLOCATABLE :: var_ana(:,:), tmp_var(:,:)
578    !
579      physfname = 'ECPHY.nc'
580    !
581      IF ( check ) WRITE(*,*) 'Opening the surface analysis'
582    !
583      CALL flininfo(physfname, iml_phys, jml_phys, llm_tmp,
584     . ttm_tmp, fid_phys)
585    !
586      ALLOCATE (lat_phys(iml_phys,jml_phys))
587      ALLOCATE (lon_phys(iml_phys,jml_phys))
588!ac
589      ALLOCATE (levphys_ini(llm_tmp))
590    !
591!      CALL flinopen(physfname, .FALSE., iml_phys, jml_phys,
592!     . llm_tmp, lon_phys, lat_phys, lev, ttm_tmp,
593!     . itau, date, dt, fid_phys)
594    !
595      CALL flinopen(physfname, .FALSE., iml_phys, jml_phys,
596     . llm_tmp, lon_phys, lat_phys, levphys_ini, ttm_tmp,
597     . itau, date, dt, fid_phys)
598    !
599      DEALLOCATE (levphys_ini)
600!ac
601    !
602    ! Allocate the space we will need to get the data out of this file
603    !
604      ALLOCATE(var_ana(iml_phys, jml_phys))
605    !
606    !   In case we have a file which is in degrees we do the transformation
607    !
608      ALLOCATE(lon_rad(iml_phys))
609      ALLOCATE(lon_ini(iml_phys))
610
611      IF ( MAXVAL(lon_phys(:,:)) .GT. 2.0 * ASIN(1.0) ) THEN
612          lon_ini(:) = lon_phys(:,1) * 2.0 * ASIN(1.0) / 180.0
613      ELSE
614          lon_ini(:) = lon_phys(:,1)
615      ENDIF
616
617      ALLOCATE(lat_rad(jml_phys))
618      ALLOCATE(lat_ini(jml_phys))
619
620      IF ( MAXVAL(lat_phys(:,:)) .GT. 2.0 * ASIN(1.0) ) THEN
621          lat_ini(:) = lat_phys(1,:) * 2.0 * ASIN(1.0) / 180.0
622      ELSE
623          lat_ini(:) = lat_phys(1,:)
624      ENDIF
625
626
627    !
628    !   We get the two standard varibales
629    !   Surface temperature
630    !
631      ALLOCATE(tsol(iml,jml))
632      ALLOCATE(tmp_var(iml-1,jml))
633    !
634    !
635
636      CALL flinget(fid_phys, 'ST', iml_phys, jml_phys,
637     .llm_tmp, ttm_tmp, 1, 1, var_ana)
638
639      title='ST'
640      CALL conf_dat2d(title,iml_phys, jml_phys, lon_ini, lat_ini,
641     . lon_rad, lat_rad, var_ana , interbar  )
642
643      IF ( interbar )   THEN
644        WRITE(6,*) '-------------------------------------------------',
645     ,'--------------'
646        WRITE(6,*) '$$$ Utilisation de l interpolation barycentrique ',
647     , ' pour  ST $$$ '
648        WRITE(6,*) '-------------------------------------------------',
649     ,'--------------'
650        CALL inter_barxy ( iml_phys,jml_phys -1,lon_rad,lat_rad ,
651     ,   var_ana, iml-1, jml-1, lon_in2, lat_in2, jml, tmp_var   )
652      ELSE
653        CALL grille_m(iml_phys, jml_phys, lon_rad, lat_rad,
654     .    var_ana, iml-1, jml, lon_in, lat_in, tmp_var     )
655      ENDIF
656
657      CALL gr_int_dyn(tmp_var, tsol, iml-1, jml)
658    !
659    ! Soil moisture
660    !
661      ALLOCATE(qsol(iml,jml))
662      CALL flinget(fid_phys, 'CDSW', iml_phys, jml_phys,
663     . llm_tmp, ttm_tmp, 1, 1, var_ana)
664
665      title='CDSW'
666      CALL conf_dat2d(title,iml_phys, jml_phys, lon_ini, lat_ini,
667     . lon_rad, lat_rad, var_ana, interbar  )
668
669      IF ( interbar )   THEN
670        WRITE(6,*) '-------------------------------------------------',
671     ,'--------------'
672        WRITE(6,*) '$$$ Utilisation de l interpolation barycentrique ',
673     , ' pour  CDSW $$$ '
674        WRITE(6,*) '-------------------------------------------------',
675     ,'--------------'
676        CALL inter_barxy ( iml_phys,jml_phys -1,lon_rad,lat_rad ,
677     ,   var_ana, iml-1, jml-1, lon_in2, lat_in2, jml, tmp_var  )
678      ELSE
679        CALL grille_m(iml_phys, jml_phys, lon_rad, lat_rad,
680     .    var_ana, iml-1, jml, lon_in, lat_in, tmp_var     )
681      ENDIF
682c
683        CALL gr_int_dyn(tmp_var, qsol, iml-1, jml)
684    !
685       CALL flinclo(fid_phys)
686    !
687      END SUBROUTINE start_init_phys
688    !
689    !
690    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
691    !
692    !
693      SUBROUTINE startget_dyn(varname, iml, jml, lon_in, lat_in,
694     . lml, pls, workvar, champ, val_exp,jml2, lon_in2, lat_in2 ,
695     ,  interbar )
696    !
697    !   ARGUMENTS
698    !
699      CHARACTER*(*), INTENT(in) :: varname
700      INTEGER, INTENT(in) :: iml, jml, lml, jml2
701      REAL, INTENT(in) :: lon_in(iml), lat_in(jml)
702      REAL, INTENT(in) :: lon_in2(iml), lat_in2(jml2)
703      REAL, INTENT(in) :: pls(iml, jml, lml)
704      REAL, INTENT(in) :: workvar(iml, jml, lml)
705      REAL, INTENT(inout) :: champ(iml, jml, lml)
706      REAL, INTENT(in) :: val_exp
707      LOGICAL interbar
708    !
709    !    LOCAL
710    !
711      INTEGER :: il, ij, ii
712      REAL :: xppn, xpps
713    !
714#include "dimensions.h"
715#include "paramet.h"
716#include "comgeom2.h"
717#include "comconst.h"
718    !
719    !   This routine only works if the variable does not exist or is constant
720    !
721      IF ( MINVAL(champ(:,:,:)).EQ.MAXVAL(champ(:,:,:)) .AND.
722     . MINVAL(champ(:,:,:)).EQ.val_exp ) THEN
723        !
724        SELECTCASE(varname)
725          CASE ('u')
726            IF ( .NOT.ALLOCATED(psol_dyn)) THEN
727              CALL start_init_dyn( iml, jml, lon_in, lat_in, jml2 ,
728     .          lon_in2,lat_in2 , interbar )
729            ENDIF
730            CALL start_inter_3d('U', iml, jml, lml, lon_in,
731     .       lat_in, jml2, lon_in2, lat_in2,  pls, champ,interbar )
732            DO il=1,lml
733              DO ij=1,jml
734                DO ii=1,iml-1
735                  champ(ii,ij,il) = champ(ii,ij,il) * cu(ii,ij)
736                ENDDO
737                champ(iml,ij, il) = champ(1,ij, il)
738              ENDDO
739            ENDDO
740          CASE ('v')
741            IF ( .NOT.ALLOCATED(psol_dyn)) THEN
742              CALL start_init_dyn( iml, jml, lon_in, lat_in , jml2,
743     .           lon_in2, lat_in2 , interbar )
744            ENDIF
745            CALL start_inter_3d('V', iml, jml, lml, lon_in,
746     .       lat_in, jml2, lon_in2, lat_in2,  pls, champ, interbar )
747            DO il=1,lml
748              DO ij=1,jml
749                DO ii=1,iml-1
750                  champ(ii,ij,il) = champ(ii,ij,il) * cv(ii,ij)
751                ENDDO
752                champ(iml,ij, il) = champ(1,ij, il)
753              ENDDO
754            ENDDO
755          CASE ('t')
756            IF ( .NOT.ALLOCATED(psol_dyn)) THEN
757              CALL start_init_dyn( iml, jml, lon_in, lat_in, jml2 ,
758     .           lon_in2, lat_in2 ,interbar )
759            ENDIF
760            CALL start_inter_3d('TEMP', iml, jml, lml, lon_in,
761     .       lat_in, jml2, lon_in2, lat_in2,  pls, champ, interbar )
762 
763          CASE ('tpot')
764            IF ( .NOT.ALLOCATED(psol_dyn)) THEN
765              CALL start_init_dyn( iml, jml, lon_in, lat_in , jml2 ,
766     .            lon_in2, lat_in2 , interbar )
767            ENDIF
768            CALL start_inter_3d('TEMP', iml, jml, lml, lon_in,
769     .       lat_in, jml2, lon_in2, lat_in2,  pls, champ, interbar )
770            IF ( MINVAL(workvar(:,:,:)) .NE. MAXVAL(workvar(:,:,:)) )
771     .                                    THEN
772              DO il=1,lml
773                DO ij=1,jml
774                  DO ii=1,iml-1
775                    champ(ii,ij,il) = champ(ii,ij,il) * cpp
776     .                                 / workvar(ii,ij,il)
777                  ENDDO
778                  champ(iml,ij,il) = champ(1,ij,il)
779                ENDDO
780              ENDDO
781              DO il=1,lml
782                xppn = SUM(aire(:,1)*champ(:,1,il))/apoln
783                xpps = SUM(aire(:,jml)*champ(:,jml,il))/apols
784                champ(:,1,il) = xppn
785                champ(:,jml,il) = xpps
786              ENDDO
787            ELSE
788              WRITE(*,*)'Could not compute potential temperature as the'
789              WRITE(*,*)'Exner function is missing or constant.'
790              STOP
791            ENDIF
792          CASE ('q')
793            IF ( .NOT.ALLOCATED(psol_dyn)) THEN
794              CALL start_init_dyn( iml, jml, lon_in, lat_in, jml2 ,
795     .           lon_in2, lat_in2 , interbar )
796            ENDIF
797            CALL start_inter_3d('R', iml, jml, lml, lon_in, lat_in,
798     .        jml2, lon_in2, lat_in2,  pls, champ, interbar )
799            IF ( MINVAL(workvar(:,:,:)) .NE. MAXVAL(workvar(:,:,:)) )
800     .                                     THEN
801              DO il=1,lml
802                DO ij=1,jml
803                  DO ii=1,iml-1
804                    champ(ii,ij,il) = 0.01 * champ(ii,ij,il) *
805     .                                       workvar(ii,ij,il)
806                  ENDDO
807                  champ(iml,ij,il) = champ(1,ij,il)
808                ENDDO
809              ENDDO
810              WHERE ( champ .LT. 0.) champ = 1.0E-10
811              DO il=1,lml
812                xppn = SUM(aire(:,1)*champ(:,1,il))/apoln
813                xpps = SUM(aire(:,jml)*champ(:,jml,il))/apols
814                champ(:,1,il) = xppn
815                champ(:,jml,il) = xpps
816              ENDDO
817            ELSE
818              WRITE(*,*)'Could not compute specific humidity as the'
819              WRITE(*,*)'saturated humidity is missing or constant.'
820              STOP
821            ENDIF
822          CASE DEFAULT
823            WRITE(*,*) 'startget_dyn'
824            WRITE(*,*) 'No rule is present to extract variable  ',
825     . varname(:LEN_TRIM(varname)),' from any data set'
826            STOP
827          END SELECT
828      ENDIF
829      END SUBROUTINE startget_dyn
830    !
831    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
832    !
833      SUBROUTINE start_init_dyn( iml, jml, lon_in, lat_in,jml2,lon_in2 ,
834     ,             lat_in2 , interbar )
835    !
836      INTEGER, INTENT(in) :: iml, jml, jml2
837      REAL, INTENT(in) :: lon_in(iml), lat_in(jml)
838      REAL, INTENT(in) :: lon_in2(iml), lat_in2(jml2)
839      LOGICAL interbar
840    !
841    !  LOCAL
842    !
843      REAL :: lev(1), date, dt
844      INTEGER :: itau(1)
845      INTEGER :: i, j
846      integer :: iret
847    !
848      CHARACTER*120 :: physfname
849      LOGICAL :: check=.TRUE.
850    !
851      REAL, ALLOCATABLE :: lon_rad(:), lat_rad(:)
852      REAL, ALLOCATABLE :: lon_ini(:), lat_ini(:)
853      REAL, ALLOCATABLE :: var_ana(:,:), tmp_var(:,:), z(:,:)
854      REAL, ALLOCATABLE :: xppn(:), xpps(:)
855      LOGICAL :: allo
856    !
857    !
858#include "dimensions.h"
859#include "paramet.h"
860#include "comgeom2.h"
861
862      CHARACTER*25 title
863
864    !
865      physfname = 'ECDYN.nc'
866    !
867      IF ( check ) WRITE(*,*) 'Opening the surface analysis'
868    !
869      CALL flininfo(physfname, iml_dyn, jml_dyn, llm_dyn,
870     .                            ttm_dyn, fid_dyn)
871      IF ( check ) WRITE(*,*) 'Values read: ', iml_dyn, jml_dyn,
872     .                                         llm_dyn, ttm_dyn
873    !
874      ALLOCATE (lat_dyn(iml_dyn,jml_dyn), stat=iret)
875      ALLOCATE (lon_dyn(iml_dyn,jml_dyn), stat=iret)
876      ALLOCATE (levdyn_ini(llm_dyn), stat=iret)
877    !
878      CALL flinopen(physfname, .FALSE., iml_dyn, jml_dyn, llm_dyn,
879     . lon_dyn, lat_dyn, levdyn_ini, ttm_dyn,
880     . itau, date, dt, fid_dyn)
881    !
882
883      allo = allocated (var_ana)
884      if (allo) then
885        DEALLOCATE(var_ana, stat=iret)
886      endif
887      ALLOCATE(var_ana(iml_dyn, jml_dyn), stat=iret)
888
889      allo = allocated (lon_rad)
890      if (allo) then
891        DEALLOCATE(lon_rad, stat=iret)
892      endif
893
894      ALLOCATE(lon_rad(iml_dyn), stat=iret)
895      ALLOCATE(lon_ini(iml_dyn))
896       
897      IF ( MAXVAL(lon_dyn(:,:)) .GT. 2.0 * ASIN(1.0) ) THEN
898          lon_ini(:) = lon_dyn(:,1) * 2.0 * ASIN(1.0) / 180.0
899      ELSE
900          lon_ini(:) = lon_dyn(:,1)
901      ENDIF
902
903      ALLOCATE(lat_rad(jml_dyn))
904      ALLOCATE(lat_ini(jml_dyn))
905
906      IF ( MAXVAL(lat_dyn(:,:)) .GT. 2.0 * ASIN(1.0) ) THEN
907          lat_ini(:) = lat_dyn(1,:) * 2.0 * ASIN(1.0) / 180.0
908      ELSE
909          lat_ini(:) = lat_dyn(1,:)
910      ENDIF
911    !
912
913
914      ALLOCATE(z(iml, jml))
915      ALLOCATE(tmp_var(iml-1,jml))
916    !
917      CALL flinget(fid_dyn, 'Z', iml_dyn, jml_dyn, 0, ttm_dyn,
918     .              1, 1, var_ana)
919c
920      title='Z'
921      CALL conf_dat2d( title,iml_dyn, jml_dyn,lon_ini, lat_ini,
922     . lon_rad, lat_rad, var_ana, interbar  )
923c
924      IF ( interbar )   THEN
925        WRITE(6,*) '-------------------------------------------------',
926     ,'--------------'
927        WRITE(6,*) '$$$ Utilisation de l interpolation barycentrique ',
928     , ' pour  Z  $$$ '
929        WRITE(6,*) '-------------------------------------------------',
930     ,'--------------'
931        CALL inter_barxy ( iml_dyn,jml_dyn -1,lon_rad,lat_rad ,
932     ,    var_ana, iml-1, jml-1, lon_in2, lat_in2, jml, tmp_var)
933      ELSE
934        CALL grille_m(iml_dyn, jml_dyn , lon_rad, lat_rad, var_ana,
935     .               iml-1, jml, lon_in, lat_in, tmp_var)
936      ENDIF
937
938      CALL gr_int_dyn(tmp_var, z, iml-1, jml)
939    !
940      ALLOCATE(psol_dyn(iml, jml))
941    !
942      CALL flinget(fid_dyn, 'SP', iml_dyn, jml_dyn, 0, ttm_dyn,
943     .              1, 1, var_ana)
944
945       title='SP'
946      CALL conf_dat2d( title,iml_dyn, jml_dyn,lon_ini, lat_ini,
947     . lon_rad, lat_rad, var_ana, interbar  )
948
949      IF ( interbar )   THEN
950        WRITE(6,*) '-------------------------------------------------',
951     ,'--------------'
952        WRITE(6,*) '$$$ Utilisation de l interpolation barycentrique ',
953     , ' pour  SP  $$$ '
954        WRITE(6,*) '-------------------------------------------------',
955     ,'--------------'
956        CALL inter_barxy ( iml_dyn,jml_dyn -1,lon_rad,lat_rad ,
957     ,    var_ana, iml-1, jml-1, lon_in2, lat_in2, jml, tmp_var)
958      ELSE
959        CALL grille_m(iml_dyn, jml_dyn , lon_rad, lat_rad, var_ana,
960     .             iml-1, jml, lon_in, lat_in, tmp_var  )
961      ENDIF
962
963      CALL gr_int_dyn(tmp_var, psol_dyn, iml-1, jml)
964    !
965      IF ( .NOT.ALLOCATED(tsol)) THEN
966    !   These variables may have been allocated by the need to
967    !   create a start field for them or by the varibale
968    !   coming out of the restart file. In case we dor have it we will initialize it.
969    !
970        CALL start_init_phys( iml, jml, lon_in, lat_in,jml2,lon_in2,
971     .                 lat_in2 , interbar )
972      ELSE
973        IF ( SIZE(tsol) .NE. SIZE(psol_dyn) ) THEN
974        WRITE(*,*) 'start_init_dyn :'
975        WRITE(*,*) 'The temperature field we have does not ',
976     .             'have the right size'
977        STOP
978      ENDIF
979      ENDIF
980      IF ( .NOT.ALLOCATED(phis)) THEN
981            !
982            !    These variables may have been allocated by the need to create a start field for them or by the varibale
983            !     coming out of the restart file. In case we dor have it we will initialize it.
984            !
985        CALL start_init_orog( iml, jml, lon_in, lat_in, jml2, lon_in2 ,
986     .      lat_in2 , interbar )
987            !
988      ELSE
989            !
990          IF (SIZE(phis) .NE. SIZE(psol_dyn)) THEN
991                !
992              WRITE(*,*) 'start_init_dyn :'
993              WRITE(*,*) 'The orography field we have does not ',
994     .                   ' have the right size'
995              STOP
996          ENDIF
997            !
998      ENDIF
999    !
1000    !     PSOL is computed in Pascals
1001    !
1002    !
1003      DO j = 1, jml
1004        DO i = 1, iml-1
1005          psol_dyn(i,j) = psol_dyn(i,j)*(1.0+(z(i,j)-phis(i,j))
1006     .                    /287.0/tsol(i,j))
1007        ENDDO
1008        psol_dyn(iml,j) = psol_dyn(1,j)
1009      ENDDO
1010    !
1011    !
1012      ALLOCATE(xppn(iml-1))
1013      ALLOCATE(xpps(iml-1))
1014    !
1015      DO  i   = 1, iml-1
1016        xppn(i) = aire( i,1) * psol_dyn( i,1)
1017        xpps(i) = aire( i,jml) * psol_dyn( i,jml)
1018      ENDDO
1019    !
1020      DO i   = 1, iml
1021        psol_dyn(i,1    )  = SUM(xppn)/apoln
1022        psol_dyn(i,jml)  = SUM(xpps)/apols
1023      ENDDO
1024    !
1025      RETURN
1026    !
1027      END SUBROUTINE start_init_dyn
1028    !
1029    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1030    !
1031      SUBROUTINE start_inter_3d(varname, iml, jml, lml, lon_in,
1032     .      lat_in, jml2, lon_in2, lat_in2, pls_in, var3d, interbar )
1033    !
1034    !    This subroutine gets a variables from a 3D file and does the interpolations needed
1035    !
1036    !
1037    !    ARGUMENTS
1038    !
1039      CHARACTER*(*) :: varname
1040      INTEGER :: iml, jml, lml, jml2
1041      REAL :: lon_in(iml), lat_in(jml), pls_in(iml, jml, lml)
1042      REAL :: lon_in2(iml) , lat_in2(jml2)
1043      REAL :: var3d(iml, jml, lml)
1044      LOGICAL interbar
1045      real chmin,chmax
1046    !
1047    !  LOCAL
1048    !
1049      CHARACTER*25 title
1050      INTEGER :: ii, ij, il, jsort,i,j,l
1051      REAL :: bx, by
1052      REAL, ALLOCATABLE :: lon_rad(:), lat_rad(:)
1053      REAL, ALLOCATABLE :: lon_ini(:), lat_ini(:) , lev_dyn(:)
1054      REAL, ALLOCATABLE :: var_tmp2d(:,:), var_tmp3d(:,:,:)
1055      REAL, ALLOCATABLE :: ax(:), ay(:), yder(:)
1056!       REAL, ALLOCATABLE :: varrr(:,:,:)
1057      INTEGER, ALLOCATABLE :: lind(:)
1058    !
1059      LOGICAL :: check = .TRUE.
1060    !
1061      IF ( .NOT. ALLOCATED(var_ana3d)) THEN
1062          ALLOCATE(var_ana3d(iml_dyn, jml_dyn, llm_dyn))
1063      ENDIF
1064!          ALLOCATE(varrr(iml_dyn, jml_dyn, llm_dyn))
1065    !
1066    !
1067      IF ( check) WRITE(*,*) 'Going into flinget to extract the 3D ',
1068     .  ' field.', fid_dyn
1069      IF ( check) WRITE(*,*) fid_dyn, varname, iml_dyn, jml_dyn,
1070     .                        llm_dyn,ttm_dyn
1071    !
1072      CALL flinget(fid_dyn, varname, iml_dyn, jml_dyn, llm_dyn,
1073     . ttm_dyn, 1, 1, var_ana3d)
1074    !
1075      IF ( check) WRITE(*,*) 'Allocating space for the interpolation',
1076     . iml, jml, llm_dyn
1077    !
1078      ALLOCATE(lon_rad(iml_dyn))
1079      ALLOCATE(lon_ini(iml_dyn))
1080
1081      IF ( MAXVAL(lon_dyn(:,:)) .GT. 2.0 * ASIN(1.0) ) THEN
1082          lon_ini(:) = lon_dyn(:,1) * 2.0 * ASIN(1.0) / 180.0
1083      ELSE
1084          lon_ini(:) = lon_dyn(:,1)
1085      ENDIF
1086
1087      ALLOCATE(lat_rad(jml_dyn))
1088      ALLOCATE(lat_ini(jml_dyn))
1089
1090      ALLOCATE(lev_dyn(llm_dyn))
1091
1092      IF ( MAXVAL(lat_dyn(:,:)) .GT. 2.0 * ASIN(1.0) ) THEN
1093          lat_ini(:) = lat_dyn(1,:) * 2.0 * ASIN(1.0) / 180.0
1094      ELSE
1095          lat_ini(:) = lat_dyn(1,:)
1096      ENDIF
1097    !
1098
1099      CALL conf_dat3d ( varname,iml_dyn, jml_dyn, llm_dyn, lon_ini,
1100     . lat_ini, levdyn_ini, lon_rad, lat_rad, lev_dyn, var_ana3d  ,
1101     ,  interbar                                                   )
1102
1103      ALLOCATE(var_tmp2d(iml-1, jml))
1104      ALLOCATE(var_tmp3d(iml, jml, llm_dyn))
1105      ALLOCATE(ax(llm_dyn))
1106      ALLOCATE(ay(llm_dyn))
1107      ALLOCATE(yder(llm_dyn))
1108      ALLOCATE(lind(llm_dyn))
1109    !
1110 
1111      DO il=1,llm_dyn
1112        !
1113      IF( interbar )  THEN
1114       IF( il.EQ.1 )  THEN
1115        WRITE(6,*) '-------------------------------------------------',
1116     ,'--------------'
1117        WRITE(6,*) '$$$ Utilisation de l interpolation barycentrique ',
1118     , ' pour ', varname
1119        WRITE(6,*) '-------------------------------------------------',
1120     ,'--------------'
1121       ENDIF
1122       CALL inter_barxy ( iml_dyn, jml_dyn -1,lon_rad, lat_rad,
1123     , var_ana3d(:,:,il),iml-1, jml2, lon_in2, lat_in2,jml,var_tmp2d )
1124      ELSE
1125       CALL grille_m(iml_dyn, jml_dyn, lon_rad, lat_rad,
1126     .  var_ana3d(:,:,il), iml-1, jml, lon_in, lat_in, var_tmp2d )
1127      ENDIF
1128        !
1129        CALL gr_int_dyn(var_tmp2d, var_tmp3d(:,:,il), iml-1, jml)
1130        !
1131       ENDDO
1132       !
1133          DO il=1,llm_dyn
1134            lind(il) = llm_dyn-il+1
1135          ENDDO
1136    !
1137c
1138c  ... Pour l'interpolation verticale ,on interpole du haut de l'atmosphere
1139c                    vers  le  sol  ...
1140c
1141      DO ij=1,jml
1142        DO ii=1,iml-1
1143          !
1144          ax(:) = lev_dyn(lind(:))
1145          ay(:) = var_tmp3d(ii, ij, lind(:))
1146          !
1147         
1148          CALL SPLINE(ax, ay, llm_dyn, 1.e30, 1.e30, yder)
1149          !
1150          DO il=1,lml
1151            bx = pls_in(ii, ij, il)
1152            CALL SPLINT(ax, ay, yder, llm_dyn, bx, by)
1153            var3d(ii, ij, il) = by
1154          ENDDO
1155          !
1156        ENDDO
1157        var3d(iml, ij, :) = var3d(1, ij, :)
1158      ENDDO
1159
1160      do il=1,lml
1161        call minmax(iml*jml,var3d(1,1,il),chmin,chmax)
1162      SELECTCASE(varname)
1163       CASE('U')
1164          WRITE(*,*) ' U  min max l ',il,chmin,chmax
1165       CASE('V')
1166          WRITE(*,*) ' V  min max l ',il,chmin,chmax
1167       CASE('TEMP')
1168          WRITE(*,*) ' TEMP  min max l ',il,chmin,chmax
1169       CASE('R')
1170          WRITE(*,*) ' R  min max l ',il,chmin,chmax
1171      END SELECT
1172           enddo
1173
1174      DEALLOCATE(lon_rad)
1175      DEALLOCATE(lon_ini)
1176      DEALLOCATE(lat_rad)
1177      DEALLOCATE(lat_ini)
1178      DEALLOCATE(lev_dyn)
1179      DEALLOCATE(var_tmp2d)
1180      DEALLOCATE(var_tmp3d)
1181      DEALLOCATE(ax)
1182      DEALLOCATE(ay)
1183      DEALLOCATE(yder)
1184      DEALLOCATE(lind)
1185
1186    !
1187      RETURN
1188    !
1189      END SUBROUTINE start_inter_3d
1190    !
1191#endif
1192! of #ifdef CPP_EARTH
1193      END MODULE startvar
Note: See TracBrowser for help on using the repository browser.