source: LMDZ4/branches/LMDZ4V5.0-dev/libf/dyn3dpar/startvar.F @ 1356

Last change on this file since 1356 was 1299, checked in by Laurent Fairhead, 15 years ago

Nettoyage general pour se rapprocher des normes et éviter des erreurs a la
compilation:

  • tous les FLOAT() sont remplacés par des REAL()
  • tous les STOP dans phylmd sont remplacés par des appels à abort_gcm
  • le common control défini dans le fichier control.h est remplacé par le module control_mod pour éviter des messages sur l'alignement des variables dans les déclarations
  • des $Header$ remplacés par des $Id$ pour svn

Quelques remplacements à faire ont pu m'échapper

General cleanup of the code to try and adhere to norms and to prevent some
compilation errors:

  • all FLOAT() instructions have been replaced by REAL() instructions
  • all STOP instructions in phylmd have been replaced by calls to abort_gcm
  • the common block control defined in the control.h file has been replaced by the control_mod to prevent compilation warnings on the alignement of declared variables
  • $Header$ replaced by $Id$ for svn

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