source: LMDZ5/trunk/libf/dyn3dpar/startvar.F90 @ 1495

Last change on this file since 1495 was 1425, checked in by lguez, 14 years ago

Replaced Numerical Recipes procedures for spline interpolation (not in
the public domain) by procedures from the Pchip package in the Slatec
library. This only affects the program "ce0l", not the program
"gcm". Tested on Brodie SX8 with "-debug" and "-prod", "-parallel
none" and "-parallel mpi". "start.nc" and "limit.nc" are
changed. "startphy.nc" is not changed. The relative change is of order
1e-7 or less. The revision makes the program faster (tested on Brodie
with "-prod -d 144x142x39", CPU time is 38 s, instead of 54
s). Procedures from Slatec are untouched, except for
"i1mach.F". Created procedures "pchfe_95" and "pchsp_95" which are
wrappers for "pchfe" and "pchsp" from Slatec. "pchfe_95" and
"pchsp_95" have a safer and simpler interface.

Replaced "make" by "sxgmake" in "arch-SX8_BRODIE.fcm". Added files for
compilation by FCM with "g95".

In "arch-linux-32bit.fcm", replaced "pgf90" by "pgf95". There was no
difference between "dev" and "debug" so added "-O1" to "dev". Added
debugging options. Removed "-Wl,-Bstatic
-L/usr/lib/gcc-lib/i386-linux/2.95.2", which usually produces an error
at link-time.

Bash is now ubiquitous while KornShell? is not so use Bash instead of
KornShell? in FCM.

Replaced some statements "write(6,*)" by "write(lunout,*)". Replaced
"stop" by "stop 1" in the case where "abort_gcm" is called with "ierr
/= 0". Removed "stop" statements at the end of procedures
"limit_netcdf" and main program "ce0l" (why not let the program end
normally?).

Made some arrays automatic instead of allocatable in "start_inter_3d".

Zeroed "wake_pe", "fm_therm", "entr_therm" and "detr_therm" in
"dyn3dpar/etat0_netcdf.F90". The parallel and sequential results of
"ce0l" are thus identical.

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 31.2 KB
Line 
1!
2! $Id: startvar.F90 1425 2010-09-02 13:45:23Z 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, &
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.