My Project
 All Classes Files Functions Variables Macros
startvar.F90
Go to the documentation of this file.
1 !
2 ! $Id: startvar.F90 1425 2010-09-02 13:45:23Z lguez $
3 !
4 !*******************************************************************************
5 !
6 MODULE 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 !
80 SUBROUTINE 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 
164 END SUBROUTINE startget_phys1d
165 !
166 !-------------------------------------------------------------------------------
167 
168 
169 !-------------------------------------------------------------------------------
170 !
171 SUBROUTINE 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 
249 END SUBROUTINE startget_phys2d
250 !
251 !-------------------------------------------------------------------------------
252 
253 
254 !-------------------------------------------------------------------------------
255 !
256 SUBROUTINE 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 
364 END SUBROUTINE startget_dyn
365 !
366 !-------------------------------------------------------------------------------
367 
368 
369 !-------------------------------------------------------------------------------
370 !
371 SUBROUTINE 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,*)
418 'Computes all parameters needed for gravity wave dra& &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 
446 END SUBROUTINE start_init_orog
447 !
448 !-------------------------------------------------------------------------------
449 
450 
451 !-------------------------------------------------------------------------------
452 !
453 SUBROUTINE 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 
519 END SUBROUTINE start_init_phys
520 !
521 !-------------------------------------------------------------------------------
522 
523 
524 !-------------------------------------------------------------------------------
525 !
526 SUBROUTINE 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 
638 END SUBROUTINE start_init_dyn
639 !
640 !-------------------------------------------------------------------------------
641 
642 
643 !-------------------------------------------------------------------------------
644 !
645 SUBROUTINE 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 
727 END SUBROUTINE start_inter_3d
728 !
729 !-------------------------------------------------------------------------------
730 
731 
732 !-------------------------------------------------------------------------------
733 !
734 SUBROUTINE 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 
774 END SUBROUTINE interp_startvar
775 !
776 !-------------------------------------------------------------------------------
777 
778 #endif
779 ! of #ifdef CPP_EARTH
780 
781 END MODULE startvar
782 !
783 !*******************************************************************************