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

Last change on this file since 248 was 24, checked in by lmdz, 25 years ago

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