source: LMDZ4/trunk/libf/dyn3d/startvar.F90 @ 1320

Last change on this file since 1320 was 1319, checked in by Laurent Fairhead, 15 years ago
  • Modifications to the start and limit creation routines to account for different

calendars

  • Modification to phyetat0 to force the mask read in the start files to match the

surface fractions read in the limit file

  • Force readaerosol.F90 to read in aerosols file with 12 timesteps

  • Modifications aux routines de créations des fichiers start et limit pour prendre

en compte différents calendriers

  • Modification à phyetat0 pour forcer le masque lu dans le fichier start à être

compatible avec les fractions de surface lu dans le fichier limit

  • Forcer readaerosol à ne lire que des fichiers à 12 pas de temps
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 30.9 KB
Line 
1!
2! $Id: startvar.F90 1319 2010-02-23 21:29:54Z idelkadi $
3!
4!*******************************************************************************
5!
6MODULE startvar
7!
8!*******************************************************************************
9!
10!-------------------------------------------------------------------------------
11! Purpose: Access data from the database of atmospheric to initialize the model.
12!-------------------------------------------------------------------------------
13! Comments:
14!
15!    *  This module is designed to work for Earth (and with ioipsl)
16!
17!    *  There are three ways to acces data, depending on the type of field
18!  which needs to be extracted. In any case the call should come after a restget
19!  and should be of the type :                     CALL startget(...)
20!
21!  - A 1D variable on the physical grid :
22!    CALL startget(varname, iml, jml,  lon_in,  lat_in,  nbindex,              &
23!           champ, val_exp,      jml2, lon_in2, lat_in2, ibar )
24!
25!  - A 2D variable on the dynamical grid :
26!    CALL startget(varname, iml, jml,  lon_in,  lat_in,                        &
27!           champ, val_exp,      jml2, lon_in2, lat_in2, ibar )             
28!
29!  - A 3D variable on the dynamical grid :
30!    CALL startget(varname, iml, jml,  lon_in,  lat_in,  lml, pls, workvar,    &
31!           champ, val_exp,      jml2, lon_in2, lat_in2, ibar )
32!
33!    *  Data needs to be in NetCDF format
34!
35!    *  Variables should have the following names in the files:
36!            'RELIEF' : High resolution orography
37!            'ST'     : Surface temperature
38!            'CDSW'   : Soil moisture
39!            'Z'      : Surface geopotential
40!            'SP'     : Surface pressure
41!            'U'      : East ward wind
42!            'V'      : Northward wind
43!            'TEMP'   : Temperature
44!            'R'      : Relative humidity
45!
46!   *   There is a big mess with the longitude size. Should it be iml or iml+1 ?
47!  I have chosen to use the iml+1 as an argument to this routine and we declare
48!  internaly smaller fields when needed. This needs to be cleared once and for
49!  all in LMDZ. A convention is required.
50!-------------------------------------------------------------------------------
51#ifdef CPP_EARTH
52  USE ioipsl
53  IMPLICIT NONE
54
55  PRIVATE
56  PUBLIC startget
57  INTERFACE startget
58    MODULE PROCEDURE startget_phys1d, startget_phys2d, startget_dyn
59  END INTERFACE
60
61  REAL,    SAVE :: deg2rad,  pi
62  INTEGER, SAVE ::           iml_rel,  jml_rel
63  INTEGER, SAVE :: fid_phys, iml_phys, jml_phys
64  INTEGER, SAVE :: fid_dyn,  iml_dyn,  jml_dyn,  llm_dyn,  ttm_dyn
65  REAL, DIMENSION(:,:),   ALLOCATABLE, TARGET, SAVE :: lon_phys, lon_dyn
66  REAL, DIMENSION(:,:),   ALLOCATABLE, TARGET, SAVE :: lat_phys, lat_dyn
67  REAL, DIMENSION(:,:),   ALLOCATABLE, TARGET, SAVE :: lon_rug, lon_alb, lon_rel
68  REAL, DIMENSION(:,:),   ALLOCATABLE, TARGET, SAVE :: lat_rug, lat_alb, lat_rel
69  REAL, DIMENSION(:),     ALLOCATABLE, TARGET, SAVE :: levdyn_ini
70  REAL, DIMENSION(:,:),   ALLOCATABLE, TARGET, SAVE :: relief, zstd, zsig, zgam
71  REAL, DIMENSION(:,:),   ALLOCATABLE, TARGET, SAVE :: masque, zthe, zpic, zval
72  REAL, DIMENSION(:,:),   ALLOCATABLE, TARGET, SAVE :: rugo, phis, tsol, qsol
73  REAL, DIMENSION(:,:),   ALLOCATABLE, TARGET, SAVE :: psol_dyn
74  REAL, DIMENSION(:,:,:), ALLOCATABLE, TARGET, SAVE :: var_ana3d
75
76   CONTAINS
77
78!-------------------------------------------------------------------------------
79!
80SUBROUTINE startget_phys1d(varname, iml, jml, lon_in, lat_in, nbindex, champ,  &
81                           val_exp ,jml2, lon_in2, lat_in2, ibar)
82!
83!-------------------------------------------------------------------------------
84! Comment:
85!   This routine only works if the variable does not exist or is constant.
86!-------------------------------------------------------------------------------
87! Arguments:
88  CHARACTER(LEN=*),         INTENT(IN)    :: varname
89  INTEGER,                  INTENT(IN)    :: iml, jml
90  REAL, DIMENSION(iml),     INTENT(IN)    :: lon_in
91  REAL, DIMENSION(jml),     INTENT(IN)    :: lat_in
92  INTEGER,                  INTENT(IN)    :: nbindex
93  REAL, DIMENSION(nbindex), INTENT(INOUT) :: champ
94  REAL,                     INTENT(IN)    :: val_exp
95  INTEGER,                  INTENT(IN)    :: jml2
96  REAL, DIMENSION(iml),     INTENT(IN)    :: lon_in2
97  REAL, DIMENSION(jml2),    INTENT(IN)    :: lat_in2
98  LOGICAL,                  INTENT(IN)    :: ibar
99!-------------------------------------------------------------------------------
100! Local variables:
101#include "iniprint.h"
102  REAL, DIMENSION(:,:), POINTER :: v2d
103!-------------------------------------------------------------------------------
104  v2d=>NULL()
105  IF(MINVAL(champ)==MAXVAL(champ).AND.MINVAL(champ)==val_exp) THEN
106
107!--- CHECKING IF THE FIELD IS KNOWN ; READING UNALLOCATED FILES
108    SELECT CASE(varname)
109      CASE('tsol')
110        IF(.NOT.ALLOCATED(tsol))                                               &
111         CALL start_init_phys(iml,jml,lon_in,lat_in,jml2,lon_in2,lat_in2,ibar)
112      CASE('qsol')
113        IF(.NOT.ALLOCATED(qsol))                                               &
114         CALL start_init_phys(iml,jml,lon_in,lat_in,jml2,lon_in2,lat_in2,ibar)
115      CASE('psol')
116        IF(.NOT.ALLOCATED(psol_dyn))                                           &
117         CALL start_init_dyn (iml,jml,lon_in,lat_in,jml2,lon_in2,lat_in2,ibar)
118      CASE('zmea','zstd','zsig','zgam','zthe','zpic','zval')
119        IF(.NOT.ALLOCATED(relief))                                             &
120         CALL start_init_orog(iml,jml,lon_in,lat_in)
121      CASE('rads','snow','tslab','seaice','rugmer','agsno')
122      CASE DEFAULT
123        WRITE(lunout,*)'startget_phys1d'
124        WRITE(lunout,*)'No rule is present to extract variable '//TRIM(varname)&
125                     //' from any data set'; STOP
126    END SELECT
127
128!--- SELECTING v2d FOR WANTED VARIABLE AND CHEKING ITS SIZE
129    SELECT CASE(varname)
130      CASE('rads','snow','tslab','seaice');  champ=0.0
131      CASE('rugmer');                        champ(:)=0.001
132      CASE('agsno');                         champ(:)=50.0
133      CASE DEFAULT
134        SELECT CASE(varname)
135          CASE('tsol'); v2d=>tsol
136          CASE('qsol'); v2d=>qsol
137          CASE('psol'); v2d=>psol_dyn
138          CASE('zmea'); v2d=>relief
139          CASE('zstd'); v2d=>zstd
140          CASE('zsig'); v2d=>zsig
141          CASE('zgam'); v2d=>zgam
142          CASE('zthe'); v2d=>zthe
143          CASE('zpic'); v2d=>zpic
144          CASE('zval'); v2d=>zval
145        END SELECT
146        IF(SIZE(v2d)/=SIZE(lon_in)*SIZE(lat_in)) THEN
147         WRITE(lunout,*)'STARTVAR module has been initialized to the wrong size'
148         STOP
149        END IF
150        CALL gr_dyn_fi(1,iml,jml,nbindex,v2d,champ)
151    END SELECT
152
153  ELSE
154
155!--- SOME FIELDS ARE CAUGHT: MAY BE NEEDED FOR A 3D INTEPROLATION
156    SELECT CASE(varname)
157      CASE('tsol')
158        IF(.NOT.ALLOCATED(tsol)) ALLOCATE(tsol(iml,jml))
159        CALL gr_fi_dyn(1,iml,jml,nbindex,champ,tsol)
160    END SELECT
161
162  END IF
163
164END SUBROUTINE  startget_phys1d
165!
166!-------------------------------------------------------------------------------
167
168
169!-------------------------------------------------------------------------------
170!
171SUBROUTINE startget_phys2d(varname, iml, jml, lon_in, lat_in, champ, val_exp,  &
172                           jml2, lon_in2, lat_in2 , ibar, msk)
173!
174!-------------------------------------------------------------------------------
175! Comment:
176!   This routine only works if the variable does not exist or is constant.
177!-------------------------------------------------------------------------------
178! Arguments:
179  CHARACTER(LEN=*),         INTENT(IN)           :: varname
180  INTEGER,                  INTENT(IN)           :: iml, jml
181  REAL, DIMENSION(iml),     INTENT(IN)           :: lon_in
182  REAL, DIMENSION(jml),     INTENT(IN)           :: lat_in
183  REAL, DIMENSION(iml,jml), INTENT(INOUT)        :: champ
184  REAL,                     INTENT(IN)           :: val_exp
185  INTEGER,                  INTENT(IN)           :: jml2
186  REAL, DIMENSION(iml),     INTENT(IN)           :: lon_in2
187  REAL, DIMENSION(jml2),    INTENT(IN)           :: lat_in2
188  LOGICAL,                  INTENT(IN)           :: ibar
189  REAL, DIMENSION(iml,jml), INTENT(IN), OPTIONAL :: msk
190!-------------------------------------------------------------------------------
191! Local variables:
192#include "iniprint.h"
193  REAL, DIMENSION(:,:), POINTER :: v2d=>NULL()
194  LOGICAL                       :: lrelief1, lrelief2
195!-------------------------------------------------------------------------------
196  v2d=>NULL()
197  lrelief1=(.NOT.ALLOCATED(relief).AND.     PRESENT(msk))
198  lrelief2=(.NOT.ALLOCATED(relief).AND..NOT.PRESENT(msk))
199  IF(MINVAL(champ)==MAXVAL(champ).AND.MINVAL(champ)==val_exp) THEN
200
201!--- CHECKING IF THE FIELD IS KNOWN ; READING UNALLOCATED FILES
202    SELECT CASE(varname)
203      CASE('psol')
204        IF(.NOT.ALLOCATED(psol_dyn))                                           &
205          CALL start_init_dyn (iml,jml,lon_in,lat_in,jml2,lon_in2,lat_in2,ibar)
206      CASE('relief')
207        IF(lrelief1)             CALL start_init_orog(iml,jml,lon_in,lat_in,msk)
208        IF(lrelief2)             CALL start_init_orog(iml,jml,lon_in,lat_in)
209      CASE('surfgeo')
210        IF(.NOT.ALLOCATED(phis)) CALL start_init_orog(iml,jml,lon_in,lat_in)
211      CASE('rugosite','masque')
212        IF(.NOT.ALLOCATED(rugo)) CALL start_init_orog(iml,jml,lon_in,lat_in)
213      CASE DEFAULT
214        WRITE(lunout,*)'startget_phys2d'
215        WRITE(lunout,*)'No rule is present to extract variable '//TRIM(varname)&
216                     //' from any data set'; STOP
217    END SELECT
218
219!--- SELECTING v2d FOR WANTED VARIABLE AND CHEKING ITS SIZE
220    SELECT CASE(varname)
221      CASE('psol');     v2d=>psol_dyn
222      CASE('relief');   v2d=>relief
223      CASE('rugosite'); v2d=>rugo
224      CASE('masque');   v2d=>masque
225      CASE('surfgeo');  v2d=>phis
226    END SELECT
227    IF(SIZE(champ)/=SIZE(v2d)) THEN
228      WRITE(lunout,*) 'STARTVAR module has been initialized to the wrong size'
229      STOP
230    END IF
231
232    champ(:,:)=v2d(:,:)
233
234  ELSE
235
236!--- SOME FIELDS ARE CAUGHT: MAY BE NEEDED FOR A 3D INTEPROLATION
237    SELECT CASE(varname)
238      CASE ('surfgeo')
239        IF(.NOT.ALLOCATED(phis)) ALLOCATE(phis(iml,jml))
240        IF(SIZE(phis)/=SIZE(champ)) THEN
241         WRITE(lunout,*)'STARTVAR module has been initialized to the wrong size'
242         STOP
243        END IF
244        phis(:,:)=champ(:,:)
245    END SELECT
246
247  END IF
248
249END SUBROUTINE startget_phys2d
250!
251!-------------------------------------------------------------------------------
252
253
254!-------------------------------------------------------------------------------
255!
256SUBROUTINE startget_dyn(varname, iml, jml,  lon_in,  lat_in,  lml, pls,workvar,&
257                     champ, val_exp,  jml2, lon_in2, lat_in2, ibar)
258!-------------------------------------------------------------------------------
259! Comment:
260!   This routine only works if the variable does not exist or is constant.
261!-------------------------------------------------------------------------------
262! Arguments:
263  CHARACTER(LEN=*),             INTENT(IN)    :: varname
264  INTEGER,                      INTENT(IN)    :: iml, jml
265  REAL, DIMENSION(iml),         INTENT(IN)    :: lon_in
266  REAL, DIMENSION(jml),         INTENT(IN)    :: lat_in
267  INTEGER,                      INTENT(IN)    :: lml
268  REAL, DIMENSION(iml,jml,lml), INTENT(IN)    :: pls, workvar
269  REAL, DIMENSION(iml,jml,lml), INTENT(INOUT) :: champ
270  REAL,                         INTENT(IN)    :: val_exp
271  INTEGER,                      INTENT(IN)    :: jml2
272  REAL, DIMENSION(iml),         INTENT(IN)    :: lon_in2
273  REAL, DIMENSION(jml2),        INTENT(IN)    :: lat_in2
274  LOGICAL,                      INTENT(IN)    :: ibar
275!-------------------------------------------------------------------------------
276! Local variables:
277#include "iniprint.h"
278#include "dimensions.h"
279#include "comconst.h"
280#include "paramet.h"
281#include "comgeom2.h"
282  REAL, DIMENSION(:,:,:), POINTER :: v3d=>NULL()
283  CHARACTER(LEN=10) :: vname
284  INTEGER :: il
285  REAL    :: xppn, xpps
286!-------------------------------------------------------------------------------
287  NULLIFY(v3d)
288  IF(MINVAL(champ)==MAXVAL(champ).AND.MINVAL(champ)==val_exp) THEN
289
290!--- READING UNALLOCATED FILES
291    IF(.NOT.ALLOCATED(psol_dyn))                                               &
292      CALL start_init_dyn(iml,jml,lon_in,lat_in,jml2,lon_in2,lat_in2,ibar)
293
294!--- CHECKING IF THE FIELD IS KNOWN AND INTERPOLATING 3D FIELDS
295    SELECT CASE(varname)
296      CASE('u');        vname='U'
297      CASE('v');        vname='V'
298      CASE('t','tpot'); vname='TEMP'
299      CASE('q');        vname='R'
300      CASE DEFAULT
301        WRITE(lunout,*)'startget_dyn'
302        WRITE(lunout,*)'No rule is present to extract variable '//TRIM(varname)&
303                //' from any data set'; STOP
304    END SELECT
305    CALL start_inter_3d(TRIM(vname), iml, jml, lml, lon_in, lat_in, jml2,      &
306                        lon_in2, lat_in2,  pls, champ,ibar )
307
308!--- COMPUTING THE REQUIRED FILED
309    SELECT CASE(varname)
310      CASE('u')                                            !--- Eastward wind
311        DO il=1,lml; champ(:,:,il)=champ(:,:,il)*cu(:,1:jml); END DO
312        champ(iml,:,:)=champ(1,:,:)
313
314      CASE('v')                                            !--- Northward wind
315        DO il=1,lml; champ(:,:,il)=champ(:,:,il)*cv(:,1:jml); END DO
316        champ(iml,:,:)=champ(1,:,:)
317
318      CASE('tpot')                                         !--- Temperature
319        IF(MINVAL(workvar)/=MAXVAL(workvar)) THEN
320          champ=champ*cpp/workvar
321          DO il=1,lml
322            xppn = SUM(aire(:,1  )*champ(:,1  ,il))/apoln
323            xpps = SUM(aire(:,jml)*champ(:,jml,il))/apols
324            champ(:,1  ,il) = xppn
325            champ(:,jml,il) = xpps
326          END DO
327        ELSE
328          WRITE(lunout,*)'Could not compute potential temperature as the'
329          WRITE(lunout,*)'Exner function is missing or constant.'; STOP
330        END IF
331
332      CASE('q')                                            !--- Relat. humidity
333        IF(MINVAL(workvar)/=MAXVAL(workvar)) THEN
334          champ=0.01*champ*workvar
335          WHERE(champ<0.) champ=1.0E-10
336          DO il=1,lml
337            xppn = SUM(aire(:,1  )*champ(:,1  ,il))/apoln
338            xpps = SUM(aire(:,jml)*champ(:,jml,il))/apols
339            champ(:,1  ,il) = xppn
340            champ(:,jml,il) = xpps
341          END DO
342        ELSE
343          WRITE(lunout,*)'Could not compute specific humidity as the'
344          WRITE(lunout,*)'saturated humidity is missing or constant.'; STOP
345        END IF
346
347    END SELECT
348
349  END IF
350
351END SUBROUTINE startget_dyn
352!
353!-------------------------------------------------------------------------------
354
355
356!-------------------------------------------------------------------------------
357!
358SUBROUTINE start_init_orog(iml,jml,lon_in,lat_in,masque_lu)
359!
360!-------------------------------------------------------------------------------
361! Arguments:
362  INTEGER,                  INTENT(IN)           :: iml, jml
363  REAL, DIMENSION(iml),     INTENT(IN)           :: lon_in
364  REAL, DIMENSION(jml),     INTENT(IN)           :: lat_in
365  REAL, DIMENSION(iml,jml), INTENT(IN), OPTIONAL :: masque_lu
366!-------------------------------------------------------------------------------
367! Local variables:
368#include "iniprint.h"
369  CHARACTER(LEN=25)     :: title
370  CHARACTER(LEN=120)    :: orofname
371  LOGICAL               :: check=.TRUE.
372  REAL,    DIMENSION(1) :: lev
373  REAL                  :: date, dt
374  INTEGER, DIMENSION(1) :: itau
375  INTEGER               :: fid, llm_tmp, ttm_tmp
376  REAL,    DIMENSION(:,:), ALLOCATABLE :: relief_hi, tmp_var
377  REAL,    DIMENSION(:),   ALLOCATABLE :: lon_rad, lat_rad, lon_ini, lat_ini
378!-------------------------------------------------------------------------------
379  pi=2.0*ASIN(1.0); deg2rad=pi/180.0
380
381  orofname = 'Relief.nc'; title='RELIEF'
382  IF(check) WRITE(lunout,*)'Reading the high resolution orography'
383  CALL flininfo(orofname, iml_rel, jml_rel, llm_tmp, ttm_tmp, fid)
384
385  ALLOCATE(lat_rel(iml_rel,jml_rel),lon_rel(iml_rel,jml_rel))
386  CALL flinopen(orofname, .FALSE., iml_rel, jml_rel, llm_tmp, lon_rel, lat_rel,&
387                lev, ttm_tmp, itau, date, dt, fid)
388  ALLOCATE(relief_hi(iml_rel,jml_rel))
389  CALL flinget(fid, title, iml_rel, jml_rel, llm_tmp, ttm_tmp, 1, 1, relief_hi)
390  CALL flinclo(fid)
391
392!--- IF ANGLES ARE IN DEGREES, THEY ARE CONVERTED INTO RADIANS
393  ALLOCATE(lon_ini(iml_rel),lat_ini(jml_rel))
394  lon_ini(:)=lon_rel(:,1); IF(MAXVAL(lon_rel)>pi) lon_ini=lon_ini*deg2rad
395  lat_ini(:)=lat_rel(1,:); IF(MAXVAL(lat_rel)>pi) lat_ini=lat_ini*deg2rad
396
397!--- FIELDS ARE PROCESSED TO BE ON STANDARD ANGULAR DOMAINS
398  ALLOCATE(lon_rad(iml_rel),lat_rad(jml_rel))
399  CALL conf_dat2d(title, iml_rel, jml_rel, lon_ini, lat_ini, lon_rad, lat_rad, &
400                  relief_hi, .FALSE.)
401  DEALLOCATE(lon_ini,lat_ini)
402
403!--- COMPUTING THE REQUIRED FIELDS USING ROUTINE grid_noro
404  IF(check) WRITE(lunout,*)'Computes all parameters needed for gravity wave dra&
405     &g code'
406
407  ALLOCATE(phis(iml,jml))      ! Geopotentiel au sol
408  ALLOCATE(zstd(iml,jml))      ! Deviation standard de l'orographie sous-maille
409  ALLOCATE(zsig(iml,jml))      ! Pente de l'orographie sous-maille
410  ALLOCATE(zgam(iml,jml))      ! Anisotropie de l'orographie sous maille
411  ALLOCATE(zthe(iml,jml))      ! Orientation axe +grande pente d'oro sous maille
412  ALLOCATE(zpic(iml,jml))      ! Hauteur pics de la SSO
413  ALLOCATE(zval(iml,jml))      ! Hauteur vallees de la SSO
414  ALLOCATE(relief(iml,jml))    ! Orographie moyenne
415  ALLOCATE(masque(iml,jml))    ! Masque terre ocean
416  masque = -99999.
417  IF(PRESENT(masque_lu)) masque=masque_lu
418
419  CALL grid_noro(iml_rel, jml_rel, lon_rad, lat_rad, relief_hi, iml-1, jml,    &
420       lon_in, lat_in, phis, relief, zstd, zsig, zgam, zthe, zpic, zval, masque)
421  phis = phis * 9.81
422
423!--- SURFACE ROUGHNESS COMPUTATION (UNUSED FOR THE MOMENT !!! )
424  IF(check) WRITE(lunout,*)'Compute surface roughness induced by the orography'
425  ALLOCATE(rugo   (iml  ,jml))
426  ALLOCATE(tmp_var(iml-1,jml))
427  CALL rugsoro(iml_rel, jml_rel, lon_rad, lat_rad, relief_hi, iml-1, jml,      &
428       lon_in, lat_in, tmp_var)
429  rugo(1:iml-1,:)=tmp_var; rugo(iml,:)=tmp_var(1,:)
430  DEALLOCATE(relief_hi,tmp_var,lon_rad,lat_rad)
431  RETURN
432
433END SUBROUTINE start_init_orog
434!
435!-------------------------------------------------------------------------------
436
437
438!-------------------------------------------------------------------------------
439!
440SUBROUTINE start_init_phys(iml,jml,lon_in,lat_in,jml2,lon_in2,lat_in2,ibar)
441!
442!-------------------------------------------------------------------------------
443! Arguments:
444  INTEGER,               INTENT(IN) :: iml, jml
445  REAL, DIMENSION(iml),  INTENT(IN) :: lon_in
446  REAL, DIMENSION(jml),  INTENT(IN) :: lat_in
447  INTEGER,               INTENT(IN) :: jml2
448  REAL, DIMENSION(iml),  INTENT(IN) :: lon_in2
449  REAL, DIMENSION(jml2), INTENT(IN) :: lat_in2
450  LOGICAL,               INTENT(IN) :: ibar
451!-------------------------------------------------------------------------------
452! Local variables:
453#include "iniprint.h"
454  CHARACTER(LEN=25)     :: title
455  CHARACTER(LEN=120)    :: physfname
456  LOGICAL               :: check=.TRUE.
457  REAL                  :: date, dt
458  INTEGER, DIMENSION(1) :: itau
459  INTEGER               :: llm_tmp, ttm_tmp
460  REAL,    DIMENSION(:,:), ALLOCATABLE :: var_ana
461  REAL,    DIMENSION(:),   ALLOCATABLE :: lon_rad, lat_rad, lon_ini, lat_ini
462  REAL,    DIMENSION(:),   ALLOCATABLE :: levphys_ini
463!-------------------------------------------------------------------------------
464  physfname = 'ECPHY.nc'; pi=2.0*ASIN(1.0); deg2rad=pi/180.0
465  IF(check) WRITE(lunout,*)'Opening the surface analysis'
466  CALL flininfo(physfname, iml_phys, jml_phys, llm_tmp, ttm_tmp, fid_phys)
467
468  ALLOCATE(lat_phys(iml_phys,jml_phys))
469  ALLOCATE(lon_phys(iml_phys,jml_phys))
470  ALLOCATE(levphys_ini(llm_tmp))
471  CALL flinopen(physfname, .FALSE., iml_phys, jml_phys, llm_tmp, lon_phys,     &
472                lat_phys, levphys_ini, ttm_tmp, itau, date, dt, fid_phys)
473  DEALLOCATE(levphys_ini)
474
475!--- IF ANGLES ARE IN DEGREES, THEY ARE CONVERTED INTO RADIANS
476  ALLOCATE(lon_ini(iml_phys),lat_ini(jml_phys))
477  lon_ini(:)=lon_phys(:,1); IF(MAXVAL(lon_phys)>pi) lon_ini=lon_ini*deg2rad
478  lat_ini(:)=lat_phys(1,:); IF(MAXVAL(lat_phys)>pi) lat_ini=lat_ini*deg2rad
479
480  ALLOCATE(var_ana(iml_phys,jml_phys),lon_rad(iml_phys),lat_rad(jml_phys))
481
482!--- SURFACE TEMPERATURE
483  title='ST'
484  ALLOCATE(tsol(iml,jml))
485  CALL flinget(fid_phys,title,iml_phys,jml_phys,llm_tmp,ttm_tmp,1,1,var_ana)
486  CALL conf_dat2d(title,iml_phys, jml_phys, lon_ini, lat_ini, lon_rad, lat_rad,&
487                  var_ana , ibar  )
488  CALL interp_startvar(title, ibar, .TRUE.,                                    &
489      iml_phys, jml_phys, lon_rad, lat_rad, var_ana, iml, jml, jml-1,          &
490      lon_in,   lat_in,   lon_in2, lat_in2, tsol)
491
492!--- SOIL MOISTURE
493  title='CDSW'
494  ALLOCATE(qsol(iml,jml))
495  CALL flinget(fid_phys,title,iml_phys,jml_phys,llm_tmp,ttm_tmp,1,1,var_ana)
496  CALL conf_dat2d(title,iml_phys, jml_phys, lon_ini, lat_ini, lon_rad, lat_rad,&
497                  var_ana, ibar  )
498  CALL interp_startvar(title, ibar, .TRUE.,                                    &
499      iml_phys, jml_phys, lon_rad, lat_rad, var_ana, iml, jml, jml-1,          &
500      lon_in,   lat_in,   lon_in2, lat_in2, qsol)
501
502  CALL flinclo(fid_phys)
503
504  DEALLOCATE(var_ana,lon_rad,lat_rad,lon_ini,lat_ini)
505
506END SUBROUTINE start_init_phys
507!
508!-------------------------------------------------------------------------------
509
510
511!-------------------------------------------------------------------------------
512!
513SUBROUTINE start_init_dyn(iml,jml,lon_in,lat_in,jml2,lon_in2,lat_in2,ibar)
514!
515!-------------------------------------------------------------------------------
516! Arguments:
517  INTEGER,               INTENT(IN) :: iml, jml
518  REAL, DIMENSION(iml),  INTENT(IN) :: lon_in
519  REAL, DIMENSION(jml),  INTENT(IN) :: lat_in
520  INTEGER,               INTENT(IN) :: jml2
521  REAL, DIMENSION(iml),  INTENT(IN) :: lon_in2
522  REAL, DIMENSION(jml2), INTENT(IN) :: lat_in2
523  LOGICAL,               INTENT(IN) :: ibar
524!-------------------------------------------------------------------------------
525! Local variables:
526#include "iniprint.h"
527#include "dimensions.h"
528#include "paramet.h"
529#include "comgeom2.h"
530  CHARACTER(LEN=25)     :: title
531  CHARACTER(LEN=120)    :: physfname
532  LOGICAL               :: check=.TRUE.
533  REAL                  :: date, dt
534  INTEGER, DIMENSION(1) :: itau
535  INTEGER               :: i, j
536  REAL,    DIMENSION(:,:), ALLOCATABLE :: var_ana, z
537  REAL,    DIMENSION(:),   ALLOCATABLE :: lon_rad, lat_rad, lon_ini, lat_ini
538  REAL,    DIMENSION(:),   ALLOCATABLE :: xppn, xpps
539!-------------------------------------------------------------------------------
540
541!--- KINETIC ENERGY
542  physfname = 'ECDYN.nc'; pi=2.0*ASIN(1.0); deg2rad=pi/180.0
543  IF(check) WRITE(lunout,*) 'Opening the surface analysis'
544  CALL flininfo(physfname, iml_dyn, jml_dyn, llm_dyn, ttm_dyn, fid_dyn)
545  IF(check) WRITE(lunout,*) 'Values read: ', iml_dyn, jml_dyn, llm_dyn, ttm_dyn
546
547  ALLOCATE(lat_dyn(iml_dyn,jml_dyn))
548  ALLOCATE(lon_dyn(iml_dyn,jml_dyn))
549  ALLOCATE(levdyn_ini(llm_dyn))
550  CALL flinopen(physfname, .FALSE., iml_dyn, jml_dyn, llm_dyn, lon_dyn,lat_dyn,&
551                levdyn_ini, ttm_dyn, itau, date, dt, fid_dyn)
552
553!--- IF ANGLES ARE IN DEGREES, THEY ARE CONVERTED INTO RADIANS
554  ALLOCATE(lon_ini(iml_dyn),lat_ini(jml_dyn))
555  lon_ini(:)=lon_dyn(:,1); IF(MAXVAL(lon_dyn)>pi) lon_ini=lon_ini*deg2rad
556  lat_ini(:)=lat_dyn(1,:); IF(MAXVAL(lat_dyn)>pi) lat_ini=lat_ini*deg2rad
557
558  ALLOCATE(var_ana(iml_dyn,jml_dyn),lon_rad(iml_dyn),lat_rad(jml_dyn))
559
560!--- SURFACE GEOPOTENTIAL
561  title='Z'
562  ALLOCATE(z(iml,jml))
563  CALL flinget(fid_dyn, title, iml_dyn, jml_dyn, 0, ttm_dyn, 1, 1, var_ana)
564  CALL conf_dat2d(title, iml_dyn, jml_dyn, lon_ini, lat_ini, lon_rad, lat_rad, &
565                  var_ana, ibar)
566  CALL interp_startvar(title, ibar, .TRUE.,                                    &
567      iml_dyn, jml_dyn, lon_rad, lat_rad, var_ana, iml, jml, jml-1,            &
568      lon_in,  lat_in,  lon_in2, lat_in2, z)
569
570!--- SURFACE PRESSURE
571  title='SP'
572  ALLOCATE(psol_dyn(iml,jml))
573  CALL flinget(fid_dyn, title, iml_dyn, jml_dyn, 0, ttm_dyn, 1, 1, var_ana)
574  CALL conf_dat2d(title, iml_dyn, jml_dyn, lon_ini, lat_ini, lon_rad, lat_rad, &
575                  var_ana, ibar)
576  CALL interp_startvar(title, ibar, .TRUE.,                                    &
577      iml_dyn, jml_dyn, lon_rad, lat_rad, var_ana, iml, jml, jml-1,            &
578      lon_in,  lat_in,  lon_in2, lat_in2, psol_dyn)
579
580  DEALLOCATE(var_ana,lon_rad,lat_rad,lon_ini,lat_ini)
581
582!--- ALLOCATION OF VARIABLES CREATED IN OR COMING FROM RESTART FILE
583  IF(.NOT.ALLOCATED(tsol)) THEN
584    CALL start_init_phys(iml,jml,lon_in,lat_in,jml2,lon_in2,lat_in2,ibar)
585  ELSE
586    IF(SIZE(tsol)/=SIZE(psol_dyn)) THEN
587      WRITE(lunout,*)'start_init_dyn :'
588      WRITE(lunout,*)'The temperature field we have does not have the right size'
589      STOP
590    END IF
591  END IF
592
593  IF(.NOT.ALLOCATED(phis)) THEN
594    CALL start_init_orog(iml,jml,lon_in,lat_in)
595  ELSE
596    IF(SIZE(phis)/=SIZE(psol_dyn)) THEN
597      WRITE(lunout,*)'start_init_dyn :'
598      WRITE(lunout,*)'The orography field we have does not have the right size'
599      STOP
600    END IF
601  END IF
602
603!--- PSOL IS COMPUTED IN PASCALS
604  DO j = 1, jml
605    DO i = 1, iml-1
606      psol_dyn(i,j) = psol_dyn(i,j)*(1.0+(z(i,j)-phis(i,j))/287.0/tsol(i,j))
607    END DO
608    psol_dyn(iml,j) = psol_dyn(1,j)
609  END DO
610  DEALLOCATE(z)
611
612  ALLOCATE(xppn(iml-1),xpps(iml-1))
613  DO i = 1, iml-1
614    xppn(i) = aire( i,1) * psol_dyn( i,1)
615    xpps(i) = aire( i,jml) * psol_dyn( i,jml)
616  END DO
617  DO i = 1, iml
618    psol_dyn(i,1  ) = SUM(xppn)/apoln
619    psol_dyn(i,jml) = SUM(xpps)/apols
620  END DO
621  DEALLOCATE(xppn,xpps)
622
623  RETURN
624
625END SUBROUTINE start_init_dyn
626!
627!-------------------------------------------------------------------------------
628
629
630!-------------------------------------------------------------------------------
631!
632SUBROUTINE start_inter_3d(varname, iml, jml, lml, lon_in, lat_in, jml2,lon_in2,&
633                          lat_in2, pls_in, var3d, ibar)
634!
635!-------------------------------------------------------------------------------
636! Arguments:
637  CHARACTER(LEN=*),             INTENT(IN)    :: varname
638  INTEGER,                      INTENT(IN)    :: iml, jml, lml
639  REAL, DIMENSION(iml),         INTENT(IN)    :: lon_in
640  REAL, DIMENSION(jml),         INTENT(IN)    :: lat_in
641  INTEGER,                      INTENT(IN)    :: jml2
642  REAL, DIMENSION(iml),         INTENT(IN)    :: lon_in2
643  REAL, DIMENSION(jml2),        INTENT(IN)    :: lat_in2
644  REAL, DIMENSION(iml,jml,lml), INTENT(IN)    :: pls_in
645  REAL, DIMENSION(iml,jml,lml), INTENT(OUT)   :: var3d
646  LOGICAL,                      INTENT(IN)    :: ibar
647!-------------------------------------------------------------------------------
648! Local variables:
649#include "iniprint.h"
650  LOGICAL               :: check=.TRUE.
651  REAL                  :: bx, by, chmin, chmax
652  INTEGER               :: ii, ij, il
653  REAL,    DIMENSION(:,:,:), ALLOCATABLE :: var_tmp3d
654  REAL,    DIMENSION(:),     ALLOCATABLE :: lon_rad, lat_rad, lon_ini, lat_ini
655  REAL,    DIMENSION(:),     ALLOCATABLE :: lev_dyn, ax, ay, yder
656  INTEGER, DIMENSION(:),     ALLOCATABLE :: lind
657!-------------------------------------------------------------------------------
658  IF(check) WRITE(lunout,*)'Going into flinget to extract the 3D  field.'
659  IF(check) WRITE(lunout,*)fid_dyn, varname, iml_dyn, jml_dyn, llm_dyn, ttm_dyn
660  IF(check) WRITE(lunout,*)'Allocating space for interpolation',iml,jml,llm_dyn
661
662  IF(.NOT.ALLOCATED(var_ana3d)) ALLOCATE(var_ana3d(iml_dyn,jml_dyn,llm_dyn))
663  CALL flinget(fid_dyn,varname,iml_dyn,jml_dyn,llm_dyn,ttm_dyn,1,1,var_ana3d)
664
665!--- IF ANGLES ARE IN DEGREES, THEY ARE CONVERTED INTO RADIANS
666  ALLOCATE(lon_ini(iml_dyn),lat_ini(jml_dyn))
667  lon_ini(:)=lon_dyn(:,1); IF(MAXVAL(lon_dyn)>pi) lon_ini=lon_ini*deg2rad
668  lat_ini(:)=lat_dyn(1,:); IF(MAXVAL(lat_dyn)>pi) lat_ini=lat_ini*deg2rad
669
670!--- FIELDS ARE PROCESSED TO BE ON STANDARD ANGULAR DOMAINS
671  ALLOCATE(lon_rad(iml_dyn),lat_rad(jml_dyn),lev_dyn(llm_dyn))
672  CALL conf_dat3d (varname, iml_dyn, jml_dyn, llm_dyn, lon_ini, lat_ini,       &
673                   levdyn_ini, lon_rad, lat_rad, lev_dyn, var_ana3d, ibar)
674  DEALLOCATE(lon_ini,lat_ini)
675
676!--- COMPUTING THE REQUIRED FIELDS USING ROUTINE grid_noro
677  ALLOCATE(var_tmp3d(iml,jml,llm_dyn))
678  DO il=1,llm_dyn
679    CALL interp_startvar(varname, ibar, il==1,                                 &
680      iml_dyn, jml_dyn, lon_rad, lat_rad, var_ana3d(:,:,il), iml, jml, jml2,   &
681      lon_in,  lat_in,  lon_in2, lat_in2, var_tmp3d(:,:,il))
682  END DO
683  DEALLOCATE(lon_rad,lat_rad)
684
685  ALLOCATE(lind(llm_dyn))
686  DO il=1,llm_dyn
687    lind(il) = llm_dyn-il+1
688  END DO
689
690!--- VERTICAL INTERPOLATION IS PERFORMED FROM TOP OF ATMOSPHERE TO GROUND
691  ALLOCATE(ax(llm_dyn),ay(llm_dyn),yder(llm_dyn))
692  DO ij=1,jml
693    DO ii=1,iml-1
694      ax(:)=lev_dyn(lind(:))
695      ay(:)=var_tmp3d(ii,ij,lind(:))
696      CALL SPLINE(ax, ay, llm_dyn, 1.e30, 1.e30, yder)
697      DO il=1,lml
698        bx=pls_in(ii,ij,il)
699        CALL SPLINT(ax, ay, yder, llm_dyn, bx, by)
700        var3d(ii,ij,il) = by
701      END DO
702    END DO
703    var3d(iml,ij,:) = var3d(1,ij,:)
704  END DO
705  DEALLOCATE(var_tmp3d,lev_dyn,ax,ay,yder,lind)
706
707  DO il=1,lml
708    CALL minmax(iml*jml,var3d(1,1,il),chmin,chmax)
709    WRITE(lunout,*)' '//TRIM(varname)//'  min max l ',il,chmin,chmax
710  END DO
711
712  RETURN
713
714END SUBROUTINE start_inter_3d
715!
716!-------------------------------------------------------------------------------
717
718
719!-------------------------------------------------------------------------------
720!
721SUBROUTINE interp_startvar(vname, ibar, ibeg, ii, jj,    lon,  lat,  vari,     &
722                                 i1, j1, j2, lon1, lat1, lon2, lat2, varo)
723!
724!-------------------------------------------------------------------------------
725! Arguments:
726  CHARACTER(LEN=*),       INTENT(IN)  :: vname
727  LOGICAL,                INTENT(IN)  :: ibar, ibeg
728  INTEGER,                INTENT(IN)  :: ii, jj
729  REAL, DIMENSION(ii),    INTENT(IN)  :: lon
730  REAL, DIMENSION(jj),    INTENT(IN)  :: lat
731  REAL, DIMENSION(ii,jj), INTENT(IN)  :: vari
732  INTEGER,                INTENT(IN)  :: i1, j1, j2
733  REAL, DIMENSION(i1),    INTENT(IN)  :: lon1
734  REAL, DIMENSION(j1),    INTENT(IN)  :: lat1
735  REAL, DIMENSION(i1),    INTENT(IN)  :: lon2
736  REAL, DIMENSION(j2),    INTENT(IN)  :: lat2
737  REAL, DIMENSION(i1,j1), INTENT(OUT) :: varo
738!-------------------------------------------------------------------------------
739! Local variables:
740#include "iniprint.h"
741  REAL, DIMENSION(i1-1,j1) :: vtmp
742!-------------------------------------------------------------------------------
743  IF(ibar) THEN
744    IF(ibeg) THEN
745      WRITE(lunout,*)                                                          &
746               '---------------------------------------------------------------'
747      WRITE(lunout,*)                                                          &
748 '$$$ Utilisation de l interpolation barycentrique  pour  '//TRIM(vname)//' $$$'
749      WRITE(lunout,*)                                                          &
750               '---------------------------------------------------------------'
751    END IF
752    CALL inter_barxy(ii, jj-1, lon, lat, vari, i1-1, j2, lon2, lat2, j1, vtmp)
753  ELSE
754    CALL grille_m   (ii, jj,   lon, lat, vari, i1-1, j1, lon1, lat1,     vtmp)
755  END IF
756  CALL gr_int_dyn(vtmp, varo, i1-1, j1)
757
758END SUBROUTINE interp_startvar
759!
760!-------------------------------------------------------------------------------
761
762#endif
763! of #ifdef CPP_EARTH
764
765END MODULE startvar
766!
767!*******************************************************************************
Note: See TracBrowser for help on using the repository browser.