source: LMDZ.3.3/trunk/libf/dyn3d/startvar.F @ 6

Last change on this file since 6 was 2, checked in by lmdz, 25 years ago

Initial revision
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 32.6 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    !
319      masque(:,:) = FLOAT(tmp_int(:,:))
320    !
321    !  Compute surface roughness
322    !
323      IF ( check ) WRITE(*,*)
324     .'Compute surface roughness induced by the orography'
325    !
326      ALLOCATE(rugo(iml,jml))
327      ALLOCATE(tmp_var(iml-1,jml))
328    !
329      CALL rugsoro(iml_rel, jml_rel, lon_rad, lat_rad, relief_hi,
330     . iml-1, jml, lon_in, lat_in, tmp_var)
331    !
332      DO j = 1, jml
333        DO i = 1, iml-1
334          rugo(i,j) = tmp_var(i,j)
335        ENDDO
336        rugo(iml,j) = tmp_var(1,j)
337      ENDDO
338    !
339    !   Build land-sea mask
340    !
341    !
342      RETURN
343    !
344      END SUBROUTINE start_init_orog
345    !
346    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
347    !
348      SUBROUTINE startget_phys1d(varname, iml, jml, lon_in,
349     .lat_in, nbindex, champ, val_exp)
350    !
351      CHARACTER*(*), INTENT(in) :: varname
352      INTEGER, INTENT(in) :: iml, jml, nbindex
353      REAL, INTENT(in) :: lon_in(iml), lat_in(jml)
354      REAL, INTENT(inout) :: champ(nbindex)
355      REAL, INTENT(in) :: val_exp
356    !
357    !
358    !   This routine only works if the variable does not exist or is constant
359    !
360      IF ( MINVAL(champ(:)).EQ.MAXVAL(champ(:)) .AND.
361     .MINVAL(champ(:)).EQ.val_exp ) THEN
362          SELECTCASE(varname)
363            CASE ('tsol')
364              IF ( .NOT.ALLOCATED(tsol)) THEN
365                CALL start_init_phys( iml, jml, lon_in, lat_in)
366              ENDIF
367              IF ( SIZE(tsol) .NE. SIZE(lon_in)*SIZE(lat_in) ) THEN
368                WRITE(*,*)
369     . 'STARTVAR module has been initialized to the wrong size'
370                 STOP
371              ENDIF
372              CALL gr_dyn_fi(1, iml, jml, nbindex, tsol, champ)
373            CASE ('qsol')
374              IF ( .NOT.ALLOCATED(qsol)) THEN
375                CALL start_init_phys( iml, jml, lon_in, lat_in)
376              ENDIF
377              IF ( SIZE(qsol) .NE. SIZE(lon_in)*SIZE(lat_in) ) THEN
378                WRITE(*,*)
379     . 'STARTVAR module has been initialized to the wrong size'
380                STOP
381              ENDIF
382              CALL gr_dyn_fi(1, iml, jml, nbindex, qsol, champ)
383            CASE ('psol')
384              IF ( .NOT.ALLOCATED(psol_dyn)) THEN
385                CALL start_init_dyn( iml, jml, lon_in, lat_in)
386              ENDIF
387              IF (SIZE(psol_dyn) .NE. SIZE(lon_in)*SIZE(lat_in)) THEN
388                WRITE(*,*)
389     . 'STARTVAR module has been initialized to the wrong size'
390                STOP
391              ENDIF
392              CALL gr_dyn_fi(1, iml, jml, nbindex, psol_dyn, champ)
393            CASE ('zmea')
394              IF ( .NOT.ALLOCATED(relief)) THEN
395                CALL start_init_orog( iml, jml, lon_in, lat_in)
396              ENDIF
397              IF ( SIZE(relief) .NE. SIZE(lon_in)*SIZE(lat_in) ) THEN
398                WRITE(*,*)
399     . 'STARTVAR module has been initialized to the wrong size'
400                 STOP
401              ENDIF
402              CALL gr_dyn_fi(1, iml, jml, nbindex, relief, champ)
403            CASE ('zstd')
404              IF ( .NOT.ALLOCATED(zstd)) THEN
405                CALL start_init_orog( iml, jml, lon_in, lat_in)
406              ENDIF
407              IF ( SIZE(zstd) .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,zstd, champ)
413            CASE ('zsig')
414              IF ( .NOT.ALLOCATED(zsig)) THEN
415                CALL start_init_orog( iml, jml, lon_in, lat_in)
416              ENDIF
417              IF ( SIZE(zsig) .NE. SIZE(lon_in)*SIZE(lat_in) ) THEN
418                WRITE(*,*)
419     . 'STARTVAR module has been initialized to the wrong size'
420                 STOP
421              ENDIF
422              CALL gr_dyn_fi(1, iml, jml, nbindex,zsig, champ)
423            CASE ('zgam')
424              IF ( .NOT.ALLOCATED(zgam)) THEN
425                CALL start_init_orog( iml, jml, lon_in, lat_in)
426              ENDIF
427              IF ( SIZE(zgam) .NE. SIZE(lon_in)*SIZE(lat_in) ) THEN
428                WRITE(*,*)
429     . 'STARTVAR module has been initialized to the wrong size'
430                 STOP
431              ENDIF
432              CALL gr_dyn_fi(1, iml, jml, nbindex,zgam, champ)
433            CASE ('zthe')
434              IF ( .NOT.ALLOCATED(zthe)) THEN
435                CALL start_init_orog( iml, jml, lon_in, lat_in)
436              ENDIF
437              IF ( SIZE(zthe) .NE. SIZE(lon_in)*SIZE(lat_in) ) THEN
438                WRITE(*,*)
439     . 'STARTVAR module has been initialized to the wrong size'
440                 STOP
441              ENDIF
442              CALL gr_dyn_fi(1, iml, jml, nbindex,zthe, champ)
443            CASE ('zpic')
444              IF ( .NOT.ALLOCATED(zpic)) THEN
445                CALL start_init_orog( iml, jml, lon_in, lat_in)
446              ENDIF
447              IF ( SIZE(zpic) .NE. SIZE(lon_in)*SIZE(lat_in) ) THEN
448                WRITE(*,*)
449     . 'STARTVAR module has been initialized to the wrong size'
450                 STOP
451              ENDIF
452              CALL gr_dyn_fi(1, iml, jml, nbindex,zpic, champ)
453            CASE ('zval')
454              IF ( .NOT.ALLOCATED(zval)) THEN
455                CALL start_init_orog( iml, jml, lon_in, lat_in)
456              ENDIF
457              IF ( SIZE(zval) .NE. SIZE(lon_in)*SIZE(lat_in) ) THEN
458                WRITE(*,*)
459     . 'STARTVAR module has been initialized to the wrong size'
460                 STOP
461              ENDIF
462              CALL gr_dyn_fi(1, iml, jml, nbindex,zval, champ)
463             CASE ('rads')
464                  champ(:) = 0.0
465            CASE ('snow')
466                  champ(:) = 0.0
467            CASE ('deltat')
468                  champ(:) = 0.0
469            CASE ('rugmer')
470                  champ(:) = 0.001
471            CASE ('agsno')
472                  champ(:) = 50.0
473            CASE DEFAULT
474              WRITE(*,*) 'startget_phys1d'
475              WRITE(*,*) 'No rule is present to extract variable  ',
476     . varname(:LEN_TRIM(varname)),' from any data set'
477              STOP
478          END SELECT
479      ELSE
480        !
481        ! If we see tsol we catch it as we may need it for a 3D interpolation
482        !
483        SELECTCASE(varname)
484          CASE ('tsol')
485            IF ( .NOT.ALLOCATED(tsol)) THEN
486              ALLOCATE(tsol(SIZE(lon_in),SIZE(lat_in) ))
487            ENDIF
488            CALL gr_fi_dyn(1, iml, jml, nbindex, champ, tsol)
489        END SELECT
490      ENDIF
491      END SUBROUTINE startget_phys1d
492    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
493    !
494      SUBROUTINE start_init_phys( iml, jml, lon_in, lat_in)
495    !
496      INTEGER, INTENT(in) :: iml, jml
497      REAL, INTENT(in) :: lon_in(iml), lat_in(jml)
498    !
499    !  LOCAL
500    !
501      REAL :: lev(1), date, dt
502      INTEGER :: itau(1)
503      INTEGER ::  llm_tmp, ttm_tmp
504      INTEGER :: i, j
505    !
506      CHARACTER*120 :: physfname
507      LOGICAL :: check=.TRUE.
508    !
509      REAL, ALLOCATABLE :: lon_rad(:), lat_rad(:)
510      REAL, ALLOCATABLE :: var_ana(:,:), tmp_var(:,:)
511    !
512      physfname = 'ECPHY.nc'
513    !
514      IF ( check ) WRITE(*,*) 'Opening the surface analysis'
515    !
516      CALL flininfo(physfname, iml_phys, jml_phys, llm_tmp,
517     . ttm_tmp, fid_phys)
518    !
519    !
520      ALLOCATE (lat_phys(iml_phys,jml_phys))
521      ALLOCATE (lon_phys(iml_phys,jml_phys))
522    !
523      CALL flinopen(physfname, .FALSE., iml_phys, jml_phys,
524     . llm_tmp, lon_phys, lat_phys, lev, ttm_tmp,
525     . itau, date, dt, fid_phys)
526    !
527    ! Allocate the space we will need to get the data out of this file
528    !
529      ALLOCATE(var_ana(iml_phys, jml_phys))
530    !
531    !   In case we have a file which is in degrees we do the transformation
532    !
533      ALLOCATE(lon_rad(iml_phys))
534      IF ( MAXVAL(lon_phys(:,:)) .GT. 2.0 * ASIN(1.0) ) THEN
535          lon_rad(:) = lon_phys(:,1) * 2.0 * ASIN(1.0) / 180.0
536      ELSE
537          lon_rad(:) = lon_phys(:,1)
538      ENDIF
539      ALLOCATE(lat_rad(jml_phys))
540      IF ( MAXVAL(lat_phys(:,:)) .GT. 2.0 * ASIN(1.0) ) THEN
541          lat_rad(:) = lat_phys(1,:) * 2.0 * ASIN(1.0) / 180.0
542      ELSE
543          lat_rad(:) = lat_phys(1,:)
544      ENDIF
545    !
546    !   We get the two standard varibales
547    !   Surface temperature
548    !
549      ALLOCATE(tsol(iml,jml))
550      ALLOCATE(tmp_var(iml-1,jml))
551    !
552    !
553      CALL flinget(fid_phys, 'ST', iml_phys, jml_phys,
554     .llm_tmp, ttm_tmp, 1, 1, var_ana)
555      CALL grille_m(iml_phys, jml_phys, lon_rad, lat_rad,
556     . var_ana, iml-1, jml, lon_in, lat_in, tmp_var)
557      CALL gr_int_dyn(tmp_var, tsol, iml-1, jml)
558    !
559    ! Soil moisture
560    !
561      ALLOCATE(qsol(iml,jml))
562      CALL flinget(fid_phys, 'CDSW', iml_phys, jml_phys,
563     . llm_tmp, ttm_tmp, 1, 1, var_ana)
564      CALL grille_m(iml_phys, jml_phys, lon_rad, lat_rad,
565     . var_ana, iml-1, jml, lon_in, lat_in, tmp_var)
566      CALL gr_int_dyn(tmp_var, qsol, iml-1, jml)
567    !
568      CALL flinclo(fid_phys)
569    !
570      END SUBROUTINE start_init_phys
571    !
572    !
573    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
574    !
575    !
576      SUBROUTINE startget_dyn(varname, iml, jml, lon_in, lat_in,
577     . lml, pls, workvar, champ, val_exp)
578    !
579    !   ARGUMENTS
580    !
581      CHARACTER*(*), INTENT(in) :: varname
582      INTEGER, INTENT(in) :: iml, jml, lml
583      REAL, INTENT(in) :: lon_in(iml), lat_in(jml)
584      REAL, INTENT(in) :: pls(iml, jml, lml)
585      REAL, INTENT(in) :: workvar(iml, jml, lml)
586      REAL, INTENT(inout) :: champ(iml, jml, lml)
587      REAL, INTENT(in) :: val_exp
588    !
589    !    LOCAL
590    !
591      INTEGER :: il, ij, ii
592      REAL :: xppn, xpps
593    !
594    !   C'est vraiment une galere de devoir rajouter tant de commons just pour avoir les aires.
595    !   Il faudrait mettre une structure plus flexible et moins dangereuse.
596    !
597#include "dimensions.h"
598#include "paramet.h"
599#include "comgeom2.h"
600#include "comconst.h"
601    !
602    !   This routine only works if the variable does not exist or is constant
603    !
604      IF ( MINVAL(champ(:,:,:)).EQ.MAXVAL(champ(:,:,:)) .AND.
605     . MINVAL(champ(:,:,:)).EQ.val_exp ) THEN
606        !
607        SELECTCASE(varname)
608          CASE ('u')
609            IF ( .NOT.ALLOCATED(psol_dyn)) THEN
610              CALL start_init_dyn( iml, jml, lon_in, lat_in)
611            ENDIF
612            CALL start_inter_3d('U', iml, jml, lml, lon_in,
613     .                           lat_in, pls, champ)
614            DO il=1,lml
615              DO ij=1,jml
616                DO ii=1,iml-1
617                  champ(ii,ij,il) = champ(ii,ij,il) * cu(ii,ij)
618                ENDDO
619                champ(iml,ij, il) = champ(1,ij, il)
620              ENDDO
621            ENDDO
622          CASE ('v')
623            IF ( .NOT.ALLOCATED(psol_dyn)) THEN
624              CALL start_init_dyn( iml, jml, lon_in, lat_in)
625            ENDIF
626            CALL start_inter_3d('V', iml, jml, lml, lon_in,
627     .                          lat_in, pls, champ)
628            DO il=1,lml
629              DO ij=1,jml
630                DO ii=1,iml-1
631                  champ(ii,ij,il) = champ(ii,ij,il) * cv(ii,ij)
632                ENDDO
633                champ(iml,ij, il) = champ(1,ij, il)
634              ENDDO
635            ENDDO
636          CASE ('t')
637            IF ( .NOT.ALLOCATED(psol_dyn)) THEN
638              CALL start_init_dyn( iml, jml, lon_in, lat_in)
639            ENDIF
640            CALL start_inter_3d('TEMP', iml, jml, lml, lon_in,
641     .                           lat_in, pls, champ)
642 
643          CASE ('tpot')
644            IF ( .NOT.ALLOCATED(psol_dyn)) THEN
645              CALL start_init_dyn( iml, jml, lon_in, lat_in)
646            ENDIF
647            CALL start_inter_3d('TEMP', iml, jml, lml, lon_in,
648     .                           lat_in, pls, champ)
649            IF ( MINVAL(workvar(:,:,:)) .NE. MAXVAL(workvar(:,:,:)) )
650     .                                    THEN
651              DO il=1,lml
652                DO ij=1,jml
653                  DO ii=1,iml-1
654                    champ(ii,ij,il) = champ(ii,ij,il) * cpp
655     .                                 / workvar(ii,ij,il)
656                  ENDDO
657                  champ(iml,ij,il) = champ(1,ij,il)
658                ENDDO
659              ENDDO
660              DO il=1,lml
661                xppn = SUM(aire(:,1)*champ(:,1,il))/apoln
662                xpps = SUM(aire(:,jml)*champ(:,jml,il))/apoln
663                champ(:,1,il) = xppn
664                champ(:,jml,il) = xpps
665              ENDDO
666            ELSE
667              WRITE(*,*)'Could not compute potential temperature as the'
668              WRITE(*,*)'Exner function is missing or constant.'
669              STOP
670            ENDIF
671          CASE ('q')
672            IF ( .NOT.ALLOCATED(psol_dyn)) THEN
673              CALL start_init_dyn( iml, jml, lon_in, lat_in)
674            ENDIF
675            CALL start_inter_3d('R', iml, jml, lml, lon_in, lat_in,
676     .                           pls, champ)
677            IF ( MINVAL(workvar(:,:,:)) .NE. MAXVAL(workvar(:,:,:)) )
678     .                                     THEN
679              DO il=1,lml
680                DO ij=1,jml
681                  DO ii=1,iml-1
682                    champ(ii,ij,il) = 0.01 * champ(ii,ij,il) *
683     .                                       workvar(ii,ij,il)
684                  ENDDO
685                  champ(iml,ij,il) = champ(1,ij,il)
686                ENDDO
687              ENDDO
688              WHERE ( champ .LT. 0.) champ = 1.0E-10
689              DO il=1,lml
690                xppn = SUM(aire(:,1)*champ(:,1,il))/apoln
691                xpps = SUM(aire(:,jml)*champ(:,jml,il))/apoln
692                champ(:,1,il) = xppn
693                champ(:,jml,il) = xpps
694              ENDDO
695            ELSE
696              WRITE(*,*)'Could not compute specific humidity as the'
697              WRITE(*,*)'saturated humidity is missing or constant.'
698              STOP
699            ENDIF
700          CASE DEFAULT
701            WRITE(*,*) 'startget_dyn'
702            WRITE(*,*) 'No rule is present to extract variable  ',
703     . varname(:LEN_TRIM(varname)),' from any data set'
704            STOP
705          END SELECT
706      ENDIF
707      END SUBROUTINE startget_dyn
708    !
709    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
710    !
711      SUBROUTINE start_init_dyn( iml, jml, lon_in, lat_in)
712    !
713      INTEGER, INTENT(in) :: iml, jml
714      REAL, INTENT(in) :: lon_in(iml), lat_in(jml)
715    !
716    !  LOCAL
717    !
718      REAL :: lev(1), date, dt
719      INTEGER :: itau(1)
720      INTEGER :: i, j
721      integer :: iret
722    !
723      CHARACTER*120 :: physfname
724      LOGICAL :: check=.TRUE.
725    !
726      REAL, ALLOCATABLE :: lon_rad(:), lat_rad(:)
727      REAL, ALLOCATABLE :: var_ana(:,:), tmp_var(:,:), z(:,:)
728      REAL, ALLOCATABLE :: xppn(:), xpps(:)
729      LOGICAL :: allo
730    !
731    !   Ce n'est pas tres pratique d'avoir a charger 3 include pour avoir la grille du modele
732    !
733#include "dimensions.h"
734#include "paramet.h"
735#include "comgeom2.h"
736    !
737      physfname = 'ECDYN.nc'
738    !
739      IF ( check ) WRITE(*,*) 'Opening the surface analysis'
740    !
741      CALL flininfo(physfname, iml_dyn, jml_dyn, llm_dyn,
742     .                            ttm_dyn, fid_dyn)
743      IF ( check ) WRITE(*,*) 'Values read: ', iml_dyn, jml_dyn,
744     .                                         llm_dyn, ttm_dyn
745    !
746      ALLOCATE (lat_dyn(iml_dyn,jml_dyn), stat=iret)
747      ALLOCATE (lon_dyn(iml_dyn,jml_dyn), stat=iret)
748      ALLOCATE (lev_dyn(llm_dyn), stat=iret)
749    !
750      CALL flinopen(physfname, .FALSE., iml_dyn, jml_dyn, llm_dyn,
751     . lon_dyn, lat_dyn, lev_dyn, ttm_dyn,
752     . itau, date, dt, fid_dyn)
753    !
754
755      allo = allocated (var_ana)
756      if (allo) then
757        DEALLOCATE(var_ana, stat=iret)
758      endif
759      ALLOCATE(var_ana(iml_dyn, jml_dyn), stat=iret)
760
761      allo = allocated (lon_rad)
762      if (allo) then
763        DEALLOCATE(lon_rad, stat=iret)
764      endif
765        ALLOCATE(lon_rad(iml_dyn), stat=iret)
766       
767      IF ( MAXVAL(lon_dyn(:,:)) .GT. 2.0 * ASIN(1.0) ) THEN
768          lon_rad(:) = lon_dyn(:,1) * 2.0 * ASIN(1.0) / 180.0
769      ELSE
770          lon_rad(:) = lon_dyn(:,1)
771      ENDIF
772      ALLOCATE(lat_rad(jml_dyn))
773      IF ( MAXVAL(lat_dyn(:,:)) .GT. 2.0 * ASIN(1.0) ) THEN
774          lat_rad(:) = lat_dyn(1,:) * 2.0 * ASIN(1.0) / 180.0
775      ELSE
776          lat_rad(:) = lat_dyn(1,:)
777      ENDIF
778    !
779      ALLOCATE(z(iml, jml))
780      ALLOCATE(tmp_var(iml-1,jml))
781    !
782      CALL flinget(fid_dyn, 'Z', iml_dyn, jml_dyn, 0, ttm_dyn,
783     .              1, 1, var_ana)
784      CALL grille_m(iml_dyn, jml_dyn , lon_rad, lat_rad, var_ana,
785     .               iml-1, jml, lon_in, lat_in, tmp_var)
786      CALL gr_int_dyn(tmp_var, z, iml-1, jml)
787    !
788      ALLOCATE(psol_dyn(iml, jml))
789    !
790      CALL flinget(fid_dyn, 'SP', iml_dyn, jml_dyn, 0, ttm_dyn,
791     .              1, 1, var_ana)
792      CALL grille_m(iml_dyn, jml_dyn , lon_rad, lat_rad, var_ana,
793     .               iml-1, jml, lon_in, lat_in, tmp_var)
794      CALL gr_int_dyn(tmp_var, psol_dyn, iml-1, jml)
795    !
796      IF ( .NOT.ALLOCATED(tsol)) THEN
797    !   These variables may have been allocated by the need to
798    !   create a start field for them or by the varibale
799    !   coming out of the restart file. In case we dor have it we will initialize it.
800    !
801        CALL start_init_phys( iml, jml, lon_in, lat_in)
802      ELSE
803        IF ( SIZE(tsol) .NE. SIZE(psol_dyn) ) THEN
804        WRITE(*,*) 'start_init_dyn :'
805        WRITE(*,*) 'The temperature field we have does not ',
806     .             'have the right size'
807        STOP
808      ENDIF
809      ENDIF
810      IF ( .NOT.ALLOCATED(phis)) THEN
811            !
812            !    These variables may have been allocated by the need to create a start field for them or by the varibale
813            !     coming out of the restart file. In case we dor have it we will initialize it.
814            !
815          CALL start_init_orog( iml, jml, lon_in, lat_in)
816            !
817      ELSE
818            !
819          IF (SIZE(phis) .NE. SIZE(psol_dyn)) THEN
820                !
821              WRITE(*,*) 'start_init_dyn :'
822              WRITE(*,*) 'The orography field we have does not ',
823     .                   ' have the right size'
824              STOP
825          ENDIF
826            !
827      ENDIF
828    !
829    !     PSOL is computed in Pascals
830    !
831    !
832      DO j = 1, jml
833        DO i = 1, iml-1
834          psol_dyn(i,j) = psol_dyn(i,j)*(1.0+(z(i,j)-phis(i,j))
835     .                    /287.0/tsol(i,j))
836        ENDDO
837        psol_dyn(iml,j) = psol_dyn(1,j)
838      ENDDO
839    !
840    !
841      ALLOCATE(xppn(iml-1))
842      ALLOCATE(xpps(iml-1))
843    !
844      DO  i   = 1, iml-1
845        xppn(i) = aire( i,1) * psol_dyn( i,1)
846        xpps(i) = aire( i,jml) * psol_dyn( i,jml)
847      ENDDO
848    !
849      DO i   = 1, iml
850        psol_dyn(i,1    )  = SUM(xppn)/apoln
851        psol_dyn(i,jml)  = SUM(xpps)/apols
852      ENDDO
853    !
854      RETURN
855    !
856      END SUBROUTINE start_init_dyn
857    !
858    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
859    !
860      SUBROUTINE start_inter_3d(varname, iml, jml, lml, lon_in,
861     .                           lat_in, pls_in, var3d)
862    !
863    !    This subroutine gets a variables from a 3D file and does the interpolations needed
864    !
865    !
866    !    ARGUMENTS
867    !
868      CHARACTER*(*) :: varname
869      INTEGER :: iml, jml, lml
870      REAL :: lon_in(iml), lat_in(jml), pls_in(iml, jml, lml)
871      REAL :: var3d(iml, jml, lml)
872    !
873    !  LOCAL
874    !
875      INTEGER :: ii, ij, il
876      REAL :: bx, by
877      REAL, ALLOCATABLE :: lon_rad(:), lat_rad(:)
878      REAL, ALLOCATABLE :: var_tmp2d(:,:), var_tmp3d(:,:,:)
879      REAL, ALLOCATABLE :: ax(:), ay(:), yder(:)
880      INTEGER, ALLOCATABLE :: lind(:)
881    !
882      LOGICAL :: check = .TRUE.
883    !
884      IF ( .NOT. ALLOCATED(var_ana3d)) THEN
885          ALLOCATE(var_ana3d(iml_dyn, jml_dyn, llm_dyn))
886      ENDIF
887    !
888    !
889      IF ( check) WRITE(*,*) 'Going into flinget to extract the 3D ',
890     .  ' field.', fid_dyn
891      IF ( check) WRITE(*,*) fid_dyn, varname, iml_dyn, jml_dyn,
892     .                        llm_dyn,ttm_dyn
893    !
894      CALL flinget(fid_dyn, varname, iml_dyn, jml_dyn, llm_dyn,
895     . ttm_dyn, 1, 1, var_ana3d)
896    !
897      IF ( check) WRITE(*,*) 'Allocating space for the interpolation',
898     . iml, jml, llm_dyn
899    !
900      ALLOCATE(lon_rad(iml_dyn))
901      IF ( MAXVAL(lon_dyn(:,:)) .GT. 2.0 * ASIN(1.0) ) THEN
902          lon_rad(:) = lon_dyn(:,1) * 2.0 * ASIN(1.0) / 180.0
903      ELSE
904          lon_rad(:) = lon_dyn(:,1)
905      ENDIF
906      ALLOCATE(lat_rad(jml_dyn))
907      IF ( MAXVAL(lat_dyn(:,:)) .GT. 2.0 * ASIN(1.0) ) THEN
908          lat_rad(:) = lat_dyn(1,:) * 2.0 * ASIN(1.0) / 180.0
909      ELSE
910          lat_rad(:) = lat_dyn(1,:)
911      ENDIF
912    !
913      ALLOCATE(var_tmp2d(iml-1, jml))
914      ALLOCATE(var_tmp3d(iml, jml, llm_dyn))
915      ALLOCATE(ax(llm_dyn))
916      ALLOCATE(ay(llm_dyn))
917      ALLOCATE(yder(llm_dyn))
918      ALLOCATE(lind(llm_dyn))
919    !
920      DO il=1,llm_dyn
921        !
922        CALL grille_m(iml_dyn, jml_dyn, lon_rad, lat_rad,
923     .var_ana3d(:,:,il), iml-1, jml, lon_in, lat_in, var_tmp2d)
924        !
925        CALL gr_int_dyn(var_tmp2d, var_tmp3d(:,:,il), iml-1, jml)
926        !
927       ENDDO
928       !
929       !  IF needed we return the vertical axis. The spline interpolation
930       !  Requires the coordinate to be in increasing order.
931    !
932      IF ( lev_dyn(1) .LT. lev_dyn(llm_dyn)) THEN
933          DO il=1,llm_dyn
934            lind(il) = il
935          ENDDO       
936      ELSE
937          DO il=1,llm_dyn
938            lind(il) = llm_dyn-il+1
939          ENDDO
940      ENDIF
941    !
942      DO ij=1,jml
943        DO ii=1,iml-1
944          !
945          ax(:) = lev_dyn(lind(:)) * 100
946          ay(:) = var_tmp3d(ii, ij, lind(:))
947          !
948          CALL SPLINE(ax, ay, llm_dyn, 1.e30, 1.e30, yder)
949          !
950          DO il=1,lml
951            bx = pls_in(ii, ij, il)
952            CALL SPLINT(ax, ay, yder, llm_dyn, bx, by)
953            var3d(ii, ij, il) = by
954          ENDDO
955          !
956        ENDDO
957        var3d(iml, ij, :) = var3d(1, ij, :)
958      ENDDO
959
960      DEALLOCATE(lon_rad)
961      DEALLOCATE(lat_rad)
962      DEALLOCATE(var_tmp2d)
963      DEALLOCATE(var_tmp3d)
964      DEALLOCATE(ax)
965      DEALLOCATE(ay)
966      DEALLOCATE(yder)
967      DEALLOCATE(lind)
968
969    !
970      RETURN
971    !
972      END SUBROUTINE start_inter_3d
973    !
974      END MODULE startvar
Note: See TracBrowser for help on using the repository browser.