source: trunk/libf/dyn3d/startvar.F90 @ 1

Last change on this file since 1 was 1, checked in by emillour, 14 years ago

Import initial LMDZ5
File size: 31.2 KB
RevLine 
[1]1!
2! $Id: startvar.F90 1425 2010-09-02 13:45:23Z lguez $
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, &
646     lon_in2, lat_in2, pls_in, var3d, ibar)
647
648  use pchsp_95_m, only: pchsp_95
649  use pchfe_95_m, only: pchfe_95
650
651! Arguments:
652  CHARACTER(LEN=*),             INTENT(IN)    :: varname
653  INTEGER,                      INTENT(IN)    :: iml, jml, lml
654  REAL, DIMENSION(iml),         INTENT(IN)    :: lon_in
655  REAL, DIMENSION(jml),         INTENT(IN)    :: lat_in
656  INTEGER,                      INTENT(IN)    :: jml2
657  REAL, DIMENSION(iml),         INTENT(IN)    :: lon_in2
658  REAL, DIMENSION(jml2),        INTENT(IN)    :: lat_in2
659  REAL, DIMENSION(iml, jml, lml), INTENT(IN)    :: pls_in
660  REAL, DIMENSION(iml, jml, lml), INTENT(OUT)   :: var3d
661  LOGICAL,                      INTENT(IN)    :: ibar
662!----------------------------------------------------------------------------
663! Local variables:
664#include "iniprint.h"
665  LOGICAL:: check=.TRUE., skip
666  REAL                  chmin, chmax
667  INTEGER ii, ij, il, ierr
668  integer n_extrap ! number of extrapolated points
669  REAL, DIMENSION(iml, jml, llm_dyn):: var_tmp3d
670  REAL,    DIMENSION(:),     ALLOCATABLE :: lon_rad, lat_rad, lon_ini, lat_ini
671  REAL, DIMENSION(llm_dyn):: lev_dyn, ax, ay, yder
672
673!---------------------------------------------------------------------------
674  IF(check) WRITE(lunout, *)'Going into flinget to extract the 3D  field.'
675  IF(check) WRITE(lunout, *) fid_dyn, varname, iml_dyn, jml_dyn, llm_dyn, &
676       ttm_dyn
677  IF(check) WRITE(lunout, *) 'Allocating space for interpolation', iml, jml, &
678       llm_dyn
679
680  IF(.NOT.ALLOCATED(var_ana3d)) ALLOCATE(var_ana3d(iml_dyn, jml_dyn, llm_dyn))
681  CALL flinget(fid_dyn, varname, iml_dyn, jml_dyn, llm_dyn, ttm_dyn, 1, 1, &
682       var_ana3d)
683
684!--- IF ANGLES ARE IN DEGREES, THEY ARE CONVERTED INTO RADIANS
685  ALLOCATE(lon_ini(iml_dyn), lat_ini(jml_dyn))
686  lon_ini(:)=lon_dyn(:, 1); IF(MAXVAL(lon_dyn)>pi) lon_ini=lon_ini*deg2rad
687  lat_ini(:)=lat_dyn(1, :); IF(MAXVAL(lat_dyn)>pi) lat_ini=lat_ini*deg2rad
688
689!--- FIELDS ARE PROCESSED TO BE ON STANDARD ANGULAR DOMAINS
690  ALLOCATE(lon_rad(iml_dyn), lat_rad(jml_dyn))
691  CALL conf_dat3d (varname, iml_dyn, jml_dyn, llm_dyn, lon_ini, lat_ini,      &
692                   levdyn_ini, lon_rad, lat_rad, lev_dyn, var_ana3d, ibar)
693  DEALLOCATE(lon_ini, lat_ini)
694
695!--- COMPUTING THE REQUIRED FIELDS USING ROUTINE grid_noro
696  DO il=1, llm_dyn
697    CALL interp_startvar(varname, ibar, il==1, iml_dyn, jml_dyn, lon_rad, &
698         lat_rad, var_ana3d(:, :, il), iml, jml, jml2, lon_in, lat_in, &
699         lon_in2, lat_in2, var_tmp3d(:, :, il))
700  END DO
701  DEALLOCATE(lon_rad, lat_rad)
702
703!--- VERTICAL INTERPOLATION IS PERFORMED FROM TOP OF ATMOSPHERE TO GROUND
704  ax = lev_dyn(llm_dyn:1:-1)
705  skip = .false.
706  n_extrap = 0
707  DO ij=1, jml
708    DO ii=1, iml-1
709      ay = var_tmp3d(ii, ij, llm_dyn:1:-1)
710      yder = pchsp_95(ax, ay, ibeg=2, iend=2, vc_beg=0., vc_end=0.)
711      CALL pchfe_95(ax, ay, yder, skip, pls_in(ii, ij, lml:1:-1), &
712           var3d(ii, ij, lml:1:-1), ierr)
713      if (ierr < 0) stop 1
714      n_extrap = n_extrap + ierr
715    END DO
716  END DO
717  if (n_extrap /= 0) then
718     print *, "start_inter_3d pchfe_95: n_extrap = ", n_extrap
719  end if
720  var3d(iml, :, :) = var3d(1, :, :)
721
722  DO il=1, lml
723    CALL minmax(iml*jml, var3d(1, 1, il), chmin, chmax)
724    WRITE(lunout, *)' '//TRIM(varname)//'  min max l ', il, chmin, chmax
725  END DO
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.