source: LMDZ.3.3/branches/rel-LF/libf/dyn3d/startvar.F @ 278

Last change on this file since 278 was 177, checked in by lmdzadmin, 25 years ago

Lots of stuff, plus particulierement:

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