source: LMDZ4/branches/LMDZ4_AR5/libf/dyn3d/startvar.F90 @ 3400

Last change on this file since 3400 was 1323, checked in by Laurent Fairhead, 15 years ago

Changes made in r1293 are integrated into the trunk
Start files are identical between r1293 and this version


Les modifications de la r1293 sont intégrées à la trunk
Les fichiers start et startphy sont identiques entre la version 1293 et celle-ci

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 31.3 KB
Line 
1!
2! $Id: startvar.F90 1323 2010-03-12 16:19:12Z fairhead $
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_phys1d((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_phys2d(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_dyn((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_phys2d, startget_phys1d, startget_dyn
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,  lon_in,  lat_in, pls,workvar,&
257                     champ, val_exp, lon_in2, lat_in2, ibar)
258
259      use assert_eq_m, only: assert_eq
260
261
262!-------------------------------------------------------------------------------
263! Comment:
264!   This routine only works if the variable does not exist or is constant.
265!-------------------------------------------------------------------------------
266! Arguments:
267  CHARACTER(LEN=*), INTENT(IN)    :: varname
268  REAL, INTENT(IN)    :: lon_in(:) ! dim(iml)
269  REAL, INTENT(IN)    :: lat_in(:) ! dim(jml)
270  REAL, INTENT(IN)    :: pls(:, :, :) ! dim(iml, jml, lml)
271  REAL, INTENT(IN)    :: workvar(:, :, :) ! dim(iml, jml, lml)
272  REAL, INTENT(INOUT) :: champ(:, :, :) ! dim(iml, jml, lml)
273  REAL, INTENT(IN)    :: val_exp
274  REAL, INTENT(IN)    :: lon_in2(:) ! dim(iml)
275  REAL, INTENT(IN)    :: lat_in2(:) ! dim(jml2)
276  LOGICAL,                      INTENT(IN)    :: ibar
277!-------------------------------------------------------------------------------
278! Local variables:
279#include "iniprint.h"
280#include "dimensions.h"
281#include "comconst.h"
282#include "paramet.h"
283#include "comgeom2.h"
284  INTEGER    :: iml, jml
285  INTEGER    :: lml
286  INTEGER    :: jml2
287  REAL, DIMENSION(:,:,:), POINTER :: v3d=>NULL()
288  CHARACTER(LEN=10) :: vname
289  INTEGER :: il
290  REAL    :: xppn, xpps
291!-------------------------------------------------------------------------------
292  NULLIFY(v3d)
293  IF(MINVAL(champ)==MAXVAL(champ).AND.MINVAL(champ)==val_exp) THEN
294
295      iml = assert_eq((/size(lon_in), size(pls, 1), size(workvar, 1), &
296     &     size(champ, 1), size(lon_in2)/), "startget_dyn iml")
297      jml = assert_eq(size(lat_in), size(pls, 2), size(workvar, 2),   &
298     &     size(champ, 2), "startget_dyn jml")
299      lml = assert_eq(size(pls, 3), size(workvar, 3), size(champ, 3), &
300     &     "startget_dyn lml")
301      jml2 = size(lat_in2)
302
303!--- READING UNALLOCATED FILES
304    IF(.NOT.ALLOCATED(psol_dyn))                                               &
305      CALL start_init_dyn(iml,jml,lon_in,lat_in,jml2,lon_in2,lat_in2,ibar)
306
307!--- CHECKING IF THE FIELD IS KNOWN AND INTERPOLATING 3D FIELDS
308    SELECT CASE(varname)
309      CASE('u');        vname='U'
310      CASE('v');        vname='V'
311      CASE('t','tpot'); vname='TEMP'
312      CASE('q');        vname='R'
313      CASE DEFAULT
314        WRITE(lunout,*)'startget_dyn'
315        WRITE(lunout,*)'No rule is present to extract variable '//TRIM(varname)&
316                //' from any data set'; STOP
317    END SELECT
318    CALL start_inter_3d(TRIM(vname), iml, jml, lml, lon_in, lat_in, jml2,      &
319                        lon_in2, lat_in2,  pls, champ,ibar )
320
321!--- COMPUTING THE REQUIRED FILED
322    SELECT CASE(varname)
323      CASE('u')                                            !--- Eastward wind
324        DO il=1,lml; champ(:,:,il)=champ(:,:,il)*cu(:,1:jml); END DO
325        champ(iml,:,:)=champ(1,:,:)
326
327      CASE('v')                                            !--- Northward wind
328        DO il=1,lml; champ(:,:,il)=champ(:,:,il)*cv(:,1:jml); END DO
329        champ(iml,:,:)=champ(1,:,:)
330
331      CASE('tpot')                                         !--- Temperature
332        IF(MINVAL(workvar)/=MAXVAL(workvar)) THEN
333          champ=champ*cpp/workvar
334          DO il=1,lml
335            xppn = SUM(aire(:,1  )*champ(:,1  ,il))/apoln
336            xpps = SUM(aire(:,jml)*champ(:,jml,il))/apols
337            champ(:,1  ,il) = xppn
338            champ(:,jml,il) = xpps
339          END DO
340        ELSE
341          WRITE(lunout,*)'Could not compute potential temperature as the'
342          WRITE(lunout,*)'Exner function is missing or constant.'; STOP
343        END IF
344
345      CASE('q')                                            !--- Relat. humidity
346        IF(MINVAL(workvar)/=MAXVAL(workvar)) THEN
347          champ=0.01*champ*workvar
348          WHERE(champ<0.) champ=1.0E-10
349          DO il=1,lml
350            xppn = SUM(aire(:,1  )*champ(:,1  ,il))/apoln
351            xpps = SUM(aire(:,jml)*champ(:,jml,il))/apols
352            champ(:,1  ,il) = xppn
353            champ(:,jml,il) = xpps
354          END DO
355        ELSE
356          WRITE(lunout,*)'Could not compute specific humidity as the'
357          WRITE(lunout,*)'saturated humidity is missing or constant.'; STOP
358        END IF
359
360    END SELECT
361
362  END IF
363
364END SUBROUTINE startget_dyn
365!
366!-------------------------------------------------------------------------------
367
368
369!-------------------------------------------------------------------------------
370!
371SUBROUTINE start_init_orog(iml,jml,lon_in,lat_in,masque_lu)
372!
373!-------------------------------------------------------------------------------
374! Arguments:
375  INTEGER,                  INTENT(IN)           :: iml, jml
376  REAL, DIMENSION(iml),     INTENT(IN)           :: lon_in
377  REAL, DIMENSION(jml),     INTENT(IN)           :: lat_in
378  REAL, DIMENSION(iml,jml), INTENT(IN), OPTIONAL :: masque_lu
379!-------------------------------------------------------------------------------
380! Local variables:
381#include "iniprint.h"
382  CHARACTER(LEN=25)     :: title
383  CHARACTER(LEN=120)    :: orofname
384  LOGICAL               :: check=.TRUE.
385  REAL,    DIMENSION(1) :: lev
386  REAL                  :: date, dt
387  INTEGER, DIMENSION(1) :: itau
388  INTEGER               :: fid, llm_tmp, ttm_tmp
389  REAL,    DIMENSION(:,:), ALLOCATABLE :: relief_hi, tmp_var
390  REAL,    DIMENSION(:),   ALLOCATABLE :: lon_rad, lat_rad, lon_ini, lat_ini
391!-------------------------------------------------------------------------------
392  pi=2.0*ASIN(1.0); deg2rad=pi/180.0
393
394  orofname = 'Relief.nc'; title='RELIEF'
395  IF(check) WRITE(lunout,*)'Reading the high resolution orography'
396  CALL flininfo(orofname, iml_rel, jml_rel, llm_tmp, ttm_tmp, fid)
397
398  ALLOCATE(lat_rel(iml_rel,jml_rel),lon_rel(iml_rel,jml_rel))
399  CALL flinopen(orofname, .FALSE., iml_rel, jml_rel, llm_tmp, lon_rel, lat_rel,&
400                lev, ttm_tmp, itau, date, dt, fid)
401  ALLOCATE(relief_hi(iml_rel,jml_rel))
402  CALL flinget(fid, title, iml_rel, jml_rel, llm_tmp, ttm_tmp, 1, 1, relief_hi)
403  CALL flinclo(fid)
404
405!--- IF ANGLES ARE IN DEGREES, THEY ARE CONVERTED INTO RADIANS
406  ALLOCATE(lon_ini(iml_rel),lat_ini(jml_rel))
407  lon_ini(:)=lon_rel(:,1); IF(MAXVAL(lon_rel)>pi) lon_ini=lon_ini*deg2rad
408  lat_ini(:)=lat_rel(1,:); IF(MAXVAL(lat_rel)>pi) lat_ini=lat_ini*deg2rad
409
410!--- FIELDS ARE PROCESSED TO BE ON STANDARD ANGULAR DOMAINS
411  ALLOCATE(lon_rad(iml_rel),lat_rad(jml_rel))
412  CALL conf_dat2d(title, iml_rel, jml_rel, lon_ini, lat_ini, lon_rad, lat_rad, &
413                  relief_hi, .FALSE.)
414  DEALLOCATE(lon_ini,lat_ini)
415
416!--- COMPUTING THE REQUIRED FIELDS USING ROUTINE grid_noro
417  IF(check) WRITE(lunout,*)'Computes all parameters needed for gravity wave dra&
418     &g code'
419
420  ALLOCATE(phis(iml,jml))      ! Geopotentiel au sol
421  ALLOCATE(zstd(iml,jml))      ! Deviation standard de l'orographie sous-maille
422  ALLOCATE(zsig(iml,jml))      ! Pente de l'orographie sous-maille
423  ALLOCATE(zgam(iml,jml))      ! Anisotropie de l'orographie sous maille
424  ALLOCATE(zthe(iml,jml))      ! Orientation axe +grande pente d'oro sous maille
425  ALLOCATE(zpic(iml,jml))      ! Hauteur pics de la SSO
426  ALLOCATE(zval(iml,jml))      ! Hauteur vallees de la SSO
427  ALLOCATE(relief(iml,jml))    ! Orographie moyenne
428  ALLOCATE(masque(iml,jml))    ! Masque terre ocean
429  masque = -99999.
430  IF(PRESENT(masque_lu)) masque=masque_lu
431
432  CALL grid_noro(iml_rel, jml_rel, lon_rad, lat_rad, relief_hi, iml-1, jml,    &
433       lon_in, lat_in, phis, relief, zstd, zsig, zgam, zthe, zpic, zval, masque)
434  phis = phis * 9.81
435
436!--- SURFACE ROUGHNESS COMPUTATION (UNUSED FOR THE MOMENT !!! )
437  IF(check) WRITE(lunout,*)'Compute surface roughness induced by the orography'
438  ALLOCATE(rugo   (iml  ,jml))
439  ALLOCATE(tmp_var(iml-1,jml))
440  CALL rugsoro(iml_rel, jml_rel, lon_rad, lat_rad, relief_hi, iml-1, jml,      &
441       lon_in, lat_in, tmp_var)
442  rugo(1:iml-1,:)=tmp_var; rugo(iml,:)=tmp_var(1,:)
443  DEALLOCATE(relief_hi,tmp_var,lon_rad,lat_rad)
444  RETURN
445
446END SUBROUTINE start_init_orog
447!
448!-------------------------------------------------------------------------------
449
450
451!-------------------------------------------------------------------------------
452!
453SUBROUTINE start_init_phys(iml,jml,lon_in,lat_in,jml2,lon_in2,lat_in2,ibar)
454!
455!-------------------------------------------------------------------------------
456! Arguments:
457  INTEGER,               INTENT(IN) :: iml, jml
458  REAL, DIMENSION(iml),  INTENT(IN) :: lon_in
459  REAL, DIMENSION(jml),  INTENT(IN) :: lat_in
460  INTEGER,               INTENT(IN) :: jml2
461  REAL, DIMENSION(iml),  INTENT(IN) :: lon_in2
462  REAL, DIMENSION(jml2), INTENT(IN) :: lat_in2
463  LOGICAL,               INTENT(IN) :: ibar
464!-------------------------------------------------------------------------------
465! Local variables:
466#include "iniprint.h"
467  CHARACTER(LEN=25)     :: title
468  CHARACTER(LEN=120)    :: physfname
469  LOGICAL               :: check=.TRUE.
470  REAL                  :: date, dt
471  INTEGER, DIMENSION(1) :: itau
472  INTEGER               :: llm_tmp, ttm_tmp
473  REAL,    DIMENSION(:,:), ALLOCATABLE :: var_ana
474  REAL,    DIMENSION(:),   ALLOCATABLE :: lon_rad, lat_rad, lon_ini, lat_ini
475  REAL,    DIMENSION(:),   ALLOCATABLE :: levphys_ini
476!-------------------------------------------------------------------------------
477  physfname = 'ECPHY.nc'; pi=2.0*ASIN(1.0); deg2rad=pi/180.0
478  IF(check) WRITE(lunout,*)'Opening the surface analysis'
479  CALL flininfo(physfname, iml_phys, jml_phys, llm_tmp, ttm_tmp, fid_phys)
480
481  ALLOCATE(lat_phys(iml_phys,jml_phys))
482  ALLOCATE(lon_phys(iml_phys,jml_phys))
483  ALLOCATE(levphys_ini(llm_tmp))
484  CALL flinopen(physfname, .FALSE., iml_phys, jml_phys, llm_tmp, lon_phys,     &
485                lat_phys, levphys_ini, ttm_tmp, itau, date, dt, fid_phys)
486  DEALLOCATE(levphys_ini)
487
488!--- IF ANGLES ARE IN DEGREES, THEY ARE CONVERTED INTO RADIANS
489  ALLOCATE(lon_ini(iml_phys),lat_ini(jml_phys))
490  lon_ini(:)=lon_phys(:,1); IF(MAXVAL(lon_phys)>pi) lon_ini=lon_ini*deg2rad
491  lat_ini(:)=lat_phys(1,:); IF(MAXVAL(lat_phys)>pi) lat_ini=lat_ini*deg2rad
492
493  ALLOCATE(var_ana(iml_phys,jml_phys),lon_rad(iml_phys),lat_rad(jml_phys))
494
495!--- SURFACE TEMPERATURE
496  title='ST'
497  ALLOCATE(tsol(iml,jml))
498  CALL flinget(fid_phys,title,iml_phys,jml_phys,llm_tmp,ttm_tmp,1,1,var_ana)
499  CALL conf_dat2d(title,iml_phys, jml_phys, lon_ini, lat_ini, lon_rad, lat_rad,&
500                  var_ana , ibar  )
501  CALL interp_startvar(title, ibar, .TRUE.,                                    &
502      iml_phys, jml_phys, lon_rad, lat_rad, var_ana, iml, jml, jml-1,          &
503      lon_in,   lat_in,   lon_in2, lat_in2, tsol)
504
505!--- SOIL MOISTURE
506  title='CDSW'
507  ALLOCATE(qsol(iml,jml))
508  CALL flinget(fid_phys,title,iml_phys,jml_phys,llm_tmp,ttm_tmp,1,1,var_ana)
509  CALL conf_dat2d(title,iml_phys, jml_phys, lon_ini, lat_ini, lon_rad, lat_rad,&
510                  var_ana, ibar  )
511  CALL interp_startvar(title, ibar, .TRUE.,                                    &
512      iml_phys, jml_phys, lon_rad, lat_rad, var_ana, iml, jml, jml-1,          &
513      lon_in,   lat_in,   lon_in2, lat_in2, qsol)
514
515  CALL flinclo(fid_phys)
516
517  DEALLOCATE(var_ana,lon_rad,lat_rad,lon_ini,lat_ini)
518
519END SUBROUTINE start_init_phys
520!
521!-------------------------------------------------------------------------------
522
523
524!-------------------------------------------------------------------------------
525!
526SUBROUTINE start_init_dyn(iml,jml,lon_in,lat_in,jml2,lon_in2,lat_in2,ibar)
527!
528!-------------------------------------------------------------------------------
529! Arguments:
530  INTEGER,               INTENT(IN) :: iml, jml
531  REAL, DIMENSION(iml),  INTENT(IN) :: lon_in
532  REAL, DIMENSION(jml),  INTENT(IN) :: lat_in
533  INTEGER,               INTENT(IN) :: jml2
534  REAL, DIMENSION(iml),  INTENT(IN) :: lon_in2
535  REAL, DIMENSION(jml2), INTENT(IN) :: lat_in2
536  LOGICAL,               INTENT(IN) :: ibar
537!-------------------------------------------------------------------------------
538! Local variables:
539#include "iniprint.h"
540#include "dimensions.h"
541#include "paramet.h"
542#include "comgeom2.h"
543  CHARACTER(LEN=25)     :: title
544  CHARACTER(LEN=120)    :: physfname
545  LOGICAL               :: check=.TRUE.
546  REAL                  :: date, dt
547  INTEGER, DIMENSION(1) :: itau
548  INTEGER               :: i, j
549  REAL,    DIMENSION(:,:), ALLOCATABLE :: var_ana, z
550  REAL,    DIMENSION(:),   ALLOCATABLE :: lon_rad, lat_rad, lon_ini, lat_ini
551  REAL,    DIMENSION(:),   ALLOCATABLE :: xppn, xpps
552!-------------------------------------------------------------------------------
553
554!--- KINETIC ENERGY
555  physfname = 'ECDYN.nc'; pi=2.0*ASIN(1.0); deg2rad=pi/180.0
556  IF(check) WRITE(lunout,*) 'Opening the surface analysis'
557  CALL flininfo(physfname, iml_dyn, jml_dyn, llm_dyn, ttm_dyn, fid_dyn)
558  IF(check) WRITE(lunout,*) 'Values read: ', iml_dyn, jml_dyn, llm_dyn, ttm_dyn
559
560  ALLOCATE(lat_dyn(iml_dyn,jml_dyn))
561  ALLOCATE(lon_dyn(iml_dyn,jml_dyn))
562  ALLOCATE(levdyn_ini(llm_dyn))
563  CALL flinopen(physfname, .FALSE., iml_dyn, jml_dyn, llm_dyn, lon_dyn,lat_dyn,&
564                levdyn_ini, ttm_dyn, itau, date, dt, fid_dyn)
565
566!--- IF ANGLES ARE IN DEGREES, THEY ARE CONVERTED INTO RADIANS
567  ALLOCATE(lon_ini(iml_dyn),lat_ini(jml_dyn))
568  lon_ini(:)=lon_dyn(:,1); IF(MAXVAL(lon_dyn)>pi) lon_ini=lon_ini*deg2rad
569  lat_ini(:)=lat_dyn(1,:); IF(MAXVAL(lat_dyn)>pi) lat_ini=lat_ini*deg2rad
570
571  ALLOCATE(var_ana(iml_dyn,jml_dyn),lon_rad(iml_dyn),lat_rad(jml_dyn))
572
573!--- SURFACE GEOPOTENTIAL
574  title='Z'
575  ALLOCATE(z(iml,jml))
576  CALL flinget(fid_dyn, title, iml_dyn, jml_dyn, 0, ttm_dyn, 1, 1, var_ana)
577  CALL conf_dat2d(title, iml_dyn, jml_dyn, lon_ini, lat_ini, lon_rad, lat_rad, &
578                  var_ana, ibar)
579  CALL interp_startvar(title, ibar, .TRUE.,                                    &
580      iml_dyn, jml_dyn, lon_rad, lat_rad, var_ana, iml, jml, jml-1,            &
581      lon_in,  lat_in,  lon_in2, lat_in2, z)
582
583!--- SURFACE PRESSURE
584  title='SP'
585  ALLOCATE(psol_dyn(iml,jml))
586  CALL flinget(fid_dyn, title, iml_dyn, jml_dyn, 0, ttm_dyn, 1, 1, var_ana)
587  CALL conf_dat2d(title, iml_dyn, jml_dyn, lon_ini, lat_ini, lon_rad, lat_rad, &
588                  var_ana, ibar)
589  CALL interp_startvar(title, ibar, .TRUE.,                                    &
590      iml_dyn, jml_dyn, lon_rad, lat_rad, var_ana, iml, jml, jml-1,            &
591      lon_in,  lat_in,  lon_in2, lat_in2, psol_dyn)
592
593  DEALLOCATE(var_ana,lon_rad,lat_rad,lon_ini,lat_ini)
594
595!--- ALLOCATION OF VARIABLES CREATED IN OR COMING FROM RESTART FILE
596  IF(.NOT.ALLOCATED(tsol)) THEN
597    CALL start_init_phys(iml,jml,lon_in,lat_in,jml2,lon_in2,lat_in2,ibar)
598  ELSE
599    IF(SIZE(tsol)/=SIZE(psol_dyn)) THEN
600      WRITE(lunout,*)'start_init_dyn :'
601      WRITE(lunout,*)'The temperature field we have does not have the right size'
602      STOP
603    END IF
604  END IF
605
606  IF(.NOT.ALLOCATED(phis)) THEN
607    CALL start_init_orog(iml,jml,lon_in,lat_in)
608  ELSE
609    IF(SIZE(phis)/=SIZE(psol_dyn)) THEN
610      WRITE(lunout,*)'start_init_dyn :'
611      WRITE(lunout,*)'The orography field we have does not have the right size'
612      STOP
613    END IF
614  END IF
615
616!--- PSOL IS COMPUTED IN PASCALS
617  DO j = 1, jml
618    DO i = 1, iml-1
619      psol_dyn(i,j) = psol_dyn(i,j)*(1.0+(z(i,j)-phis(i,j))/287.0/tsol(i,j))
620    END DO
621    psol_dyn(iml,j) = psol_dyn(1,j)
622  END DO
623  DEALLOCATE(z)
624
625  ALLOCATE(xppn(iml-1),xpps(iml-1))
626  DO i = 1, iml-1
627    xppn(i) = aire( i,1) * psol_dyn( i,1)
628    xpps(i) = aire( i,jml) * psol_dyn( i,jml)
629  END DO
630  DO i = 1, iml
631    psol_dyn(i,1  ) = SUM(xppn)/apoln
632    psol_dyn(i,jml) = SUM(xpps)/apols
633  END DO
634  DEALLOCATE(xppn,xpps)
635
636  RETURN
637
638END SUBROUTINE start_init_dyn
639!
640!-------------------------------------------------------------------------------
641
642
643!-------------------------------------------------------------------------------
644!
645SUBROUTINE start_inter_3d(varname, iml, jml, lml, lon_in, lat_in, jml2,lon_in2,&
646                          lat_in2, pls_in, var3d, ibar)
647!
648!-------------------------------------------------------------------------------
649! Arguments:
650  CHARACTER(LEN=*),             INTENT(IN)    :: varname
651  INTEGER,                      INTENT(IN)    :: iml, jml, lml
652  REAL, DIMENSION(iml),         INTENT(IN)    :: lon_in
653  REAL, DIMENSION(jml),         INTENT(IN)    :: lat_in
654  INTEGER,                      INTENT(IN)    :: jml2
655  REAL, DIMENSION(iml),         INTENT(IN)    :: lon_in2
656  REAL, DIMENSION(jml2),        INTENT(IN)    :: lat_in2
657  REAL, DIMENSION(iml,jml,lml), INTENT(IN)    :: pls_in
658  REAL, DIMENSION(iml,jml,lml), INTENT(OUT)   :: var3d
659  LOGICAL,                      INTENT(IN)    :: ibar
660!-------------------------------------------------------------------------------
661! Local variables:
662#include "iniprint.h"
663  LOGICAL               :: check=.TRUE.
664  REAL                  :: bx, by, chmin, chmax
665  INTEGER               :: ii, ij, il
666  REAL,    DIMENSION(:,:,:), ALLOCATABLE :: var_tmp3d
667  REAL,    DIMENSION(:),     ALLOCATABLE :: lon_rad, lat_rad, lon_ini, lat_ini
668  REAL,    DIMENSION(:),     ALLOCATABLE :: lev_dyn, ax, ay, yder
669  INTEGER, DIMENSION(:),     ALLOCATABLE :: lind
670!-------------------------------------------------------------------------------
671  IF(check) WRITE(lunout,*)'Going into flinget to extract the 3D  field.'
672  IF(check) WRITE(lunout,*)fid_dyn, varname, iml_dyn, jml_dyn, llm_dyn, ttm_dyn
673  IF(check) WRITE(lunout,*)'Allocating space for interpolation',iml,jml,llm_dyn
674
675  IF(.NOT.ALLOCATED(var_ana3d)) ALLOCATE(var_ana3d(iml_dyn,jml_dyn,llm_dyn))
676  CALL flinget(fid_dyn,varname,iml_dyn,jml_dyn,llm_dyn,ttm_dyn,1,1,var_ana3d)
677
678!--- IF ANGLES ARE IN DEGREES, THEY ARE CONVERTED INTO RADIANS
679  ALLOCATE(lon_ini(iml_dyn),lat_ini(jml_dyn))
680  lon_ini(:)=lon_dyn(:,1); IF(MAXVAL(lon_dyn)>pi) lon_ini=lon_ini*deg2rad
681  lat_ini(:)=lat_dyn(1,:); IF(MAXVAL(lat_dyn)>pi) lat_ini=lat_ini*deg2rad
682
683!--- FIELDS ARE PROCESSED TO BE ON STANDARD ANGULAR DOMAINS
684  ALLOCATE(lon_rad(iml_dyn),lat_rad(jml_dyn),lev_dyn(llm_dyn))
685  CALL conf_dat3d (varname, iml_dyn, jml_dyn, llm_dyn, lon_ini, lat_ini,       &
686                   levdyn_ini, lon_rad, lat_rad, lev_dyn, var_ana3d, ibar)
687  DEALLOCATE(lon_ini,lat_ini)
688
689!--- COMPUTING THE REQUIRED FIELDS USING ROUTINE grid_noro
690  ALLOCATE(var_tmp3d(iml,jml,llm_dyn))
691  DO il=1,llm_dyn
692    CALL interp_startvar(varname, ibar, il==1,                                 &
693      iml_dyn, jml_dyn, lon_rad, lat_rad, var_ana3d(:,:,il), iml, jml, jml2,   &
694      lon_in,  lat_in,  lon_in2, lat_in2, var_tmp3d(:,:,il))
695  END DO
696  DEALLOCATE(lon_rad,lat_rad)
697
698  ALLOCATE(lind(llm_dyn))
699  DO il=1,llm_dyn
700    lind(il) = llm_dyn-il+1
701  END DO
702
703!--- VERTICAL INTERPOLATION IS PERFORMED FROM TOP OF ATMOSPHERE TO GROUND
704  ALLOCATE(ax(llm_dyn),ay(llm_dyn),yder(llm_dyn))
705  DO ij=1,jml
706    DO ii=1,iml-1
707      ax(:)=lev_dyn(lind(:))
708      ay(:)=var_tmp3d(ii,ij,lind(:))
709      CALL SPLINE(ax, ay, llm_dyn, 1.e30, 1.e30, yder)
710      DO il=1,lml
711        bx=pls_in(ii,ij,il)
712        CALL SPLINT(ax, ay, yder, llm_dyn, bx, by)
713        var3d(ii,ij,il) = by
714      END DO
715    END DO
716    var3d(iml,ij,:) = var3d(1,ij,:)
717  END DO
718  DEALLOCATE(var_tmp3d,lev_dyn,ax,ay,yder,lind)
719
720  DO il=1,lml
721    CALL minmax(iml*jml,var3d(1,1,il),chmin,chmax)
722    WRITE(lunout,*)' '//TRIM(varname)//'  min max l ',il,chmin,chmax
723  END DO
724
725  RETURN
726
727END SUBROUTINE start_inter_3d
728!
729!-------------------------------------------------------------------------------
730
731
732!-------------------------------------------------------------------------------
733!
734SUBROUTINE interp_startvar(vname, ibar, ibeg, ii, jj,    lon,  lat,  vari,     &
735                                 i1, j1, j2, lon1, lat1, lon2, lat2, varo)
736!
737!-------------------------------------------------------------------------------
738
739  USE inter_barxy_m, only: inter_barxy
740
741! Arguments:
742  CHARACTER(LEN=*),       INTENT(IN)  :: vname
743  LOGICAL,                INTENT(IN)  :: ibar, ibeg
744  INTEGER,                INTENT(IN)  :: ii, jj
745  REAL, DIMENSION(ii),    INTENT(IN)  :: lon
746  REAL, DIMENSION(jj),    INTENT(IN)  :: lat
747  REAL, DIMENSION(ii,jj), INTENT(IN)  :: vari
748  INTEGER,                INTENT(IN)  :: i1, j1, j2
749  REAL, DIMENSION(i1),    INTENT(IN)  :: lon1
750  REAL, DIMENSION(j1),    INTENT(IN)  :: lat1
751  REAL, DIMENSION(i1),    INTENT(IN)  :: lon2
752  REAL, DIMENSION(j2),    INTENT(IN)  :: lat2
753  REAL, DIMENSION(i1,j1), INTENT(OUT) :: varo
754!-------------------------------------------------------------------------------
755! Local variables:
756#include "iniprint.h"
757  REAL, DIMENSION(i1-1,j1) :: vtmp
758!-------------------------------------------------------------------------------
759  IF(ibar) THEN
760    IF(ibeg) THEN
761      WRITE(lunout,*)                                                          &
762               '---------------------------------------------------------------'
763      WRITE(lunout,*)                                                          &
764 '$$$ Utilisation de l interpolation barycentrique  pour  '//TRIM(vname)//' $$$'
765      WRITE(lunout,*)                                                          &
766               '---------------------------------------------------------------'
767    END IF
768    CALL inter_barxy(lon, lat(:jj-1), vari, lon2(:i1-1), lat2(:j2), vtmp)
769  ELSE
770    CALL grille_m   (ii, jj,   lon, lat, vari, i1-1, j1, lon1, lat1,     vtmp)
771  END IF
772  CALL gr_int_dyn(vtmp, varo, i1-1, j1)
773
774END SUBROUTINE interp_startvar
775!
776!-------------------------------------------------------------------------------
777
778#endif
779! of #ifdef CPP_EARTH
780
781END MODULE startvar
782!
783!*******************************************************************************
Note: See TracBrowser for help on using the repository browser.