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

Last change on this file since 112 was 99, checked in by lmdzadmin, 25 years ago

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