source: LMDZ6/branches/LMDZ_DECOUPLE/libf/dynphy_lonlat/phylmd/etat0dyn_netcdf.F90

Last change on this file was 3435, checked in by Laurent Fairhead, 6 years ago

"Historic" :-) commit merging the physics branch used for DYNAMICO with the LMDZ trunk.
The same physics branch can now be used seamlessly with the traditional lon-lat LMDZ
dynamical core and DYNAMICO.
Testing consisted in running a lon-lat LMDZ bucket simulation with the NPv6.1 physics package
with the original trunk sources and the merged sources. Tests were succesful in the sense that
numeric continuity was preserved in the restart files from both simulation. Further tests
included running both versions of the physics codes for one year in a LMDZOR setting in which
the restart files also came out identical.

Caution:

  • as the physics package now manages unstructured grids, grid information needs to be transmitted

to the surface scheme ORCHIDEE. This means that the interface defined in surf_land_orchidee_mod.F90
is only compatible with ORCHIDEE version orchidee2.1 and later versions. If previous versions of
ORCHIDEE need to be used, the CPP key ORCHIDEE_NOUNSTRUCT needs to be set at compilation time.
This is done automatically if makelmdz/makelmdz_fcm are called with the veget orchidee2.0 switch

  • due to a limitation in XIOS, the time at which limit conditions will be read in by DYNAMICO will be

delayed by one physic timestep with respect to the time it is read in by the lon-lat model. This is caused
by the line

IF (MOD(itime-1, lmt_pas) == 0 .OR. (jour_lu /= jour .AND. grid_type /= unstructured)) THEN ! time to read

in limit_read_mod.F90

Work still needed on COSP integration and XML files for DYNAMICO

EM, YM, LF

File size: 21.8 KB
Line 
1MODULE etat0dyn
2!
3!*******************************************************************************
4! Purpose: Create dynamical initial state using atmospheric fields from a
5!          database of atmospheric to initialize the model.
6!-------------------------------------------------------------------------------
7! Comments:
8!
9!    *  This module is designed to work for Earth (and with ioipsl)
10!
11!    *  etat0dyn_netcdf routine can access to NetCDF data through the following
12!  routine (to be called after restget):
13!    CALL startget_dyn3d(varname, lon_in,  lat_in, pls, workvar,&
14!                          champ, lon_in2, lat_in2)
15!
16!    *  Variables should have the following names in the NetCDF files:
17!            'U'      : East ward wind              (in "ECDYN.nc")
18!            'V'      : Northward wind              (in "ECDYN.nc")
19!            'TEMP'   : Temperature                 (in "ECDYN.nc")
20!            'R'      : Relative humidity           (in "ECDYN.nc")
21!            'RELIEF' : High resolution orography   (in "Relief.nc")
22!
23!    * The land mask and corresponding weights can be:
24!      1) already known (in particular if etat0dyn has been called before) ;
25!         in this case, ANY(masque(:,:)/=-99999.) = .TRUE.
26!      2) computed using the ocean mask from the ocean model (to ensure ocean
27!         fractions are the same for atmosphere and ocean) for coupled runs.
28!         File name: "o2a.nc"  ;  variable name: "OceMask"
29!      3) computed from topography file "Relief.nc" for forced runs.
30!
31!   *   There is a big mess with the longitude size. Should it be iml or iml+1 ?
32!  I have chosen to use the iml+1 as an argument to this routine and we declare
33!  internaly smaller fields when needed. This needs to be cleared once and for
34!  all in LMDZ. A convention is required.
35!-------------------------------------------------------------------------------
36  USE ioipsl,         ONLY: flininfo, flinopen, flinget, flinclo, histclo
37  USE assert_eq_m,    ONLY: assert_eq
38  USE comconst_mod, ONLY: pi, cpp, kappa
39  USE comvert_mod, ONLY: ap, bp, preff, pressure_exner
40  USE temps_mod, ONLY: annee_ref, day_ref, itau_dyn, itau_phy, start_time
41 
42  IMPLICIT NONE
43
44  PRIVATE
45  PUBLIC :: etat0dyn_netcdf
46
47  include "iniprint.h"
48  include "dimensions.h"
49  include "paramet.h"
50  include "comgeom2.h"
51  include "comdissnew.h"
52  REAL, SAVE :: deg2rad
53  INTEGER,            SAVE      :: iml_dyn, jml_dyn, llm_dyn, ttm_dyn, fid_dyn
54  REAL, ALLOCATABLE,  SAVE      :: lon_dyn(:,:), lat_dyn(:,:), levdyn_ini(:)
55  CHARACTER(LEN=120), PARAMETER :: dynfname='ECDYN.nc'
56
57CONTAINS
58
59!-------------------------------------------------------------------------------
60!
61SUBROUTINE etat0dyn_netcdf(masque, phis)
62!
63!-------------------------------------------------------------------------------
64! Purpose: Create dynamical initial states.
65!-------------------------------------------------------------------------------
66! Notes:  1) This routine is designed to work for Earth
67!         2) If masque(:,:)/=-99999., masque and phis are already known.
68!         Otherwise: compute it.
69!-------------------------------------------------------------------------------
70  USE control_mod
71  USE regr_lat_time_coefoz_m, ONLY: regr_lat_time_coefoz
72  USE regr_pr_o3_m,   ONLY: regr_pr_o3
73  USE press_coefoz_m, ONLY: press_coefoz
74  USE exner_hyb_m,    ONLY: exner_hyb
75  USE exner_milieu_m, ONLY: exner_milieu
76  USE infotrac,       ONLY: nqtot, tname
77  USE filtreg_mod
78  IMPLICIT NONE
79!-------------------------------------------------------------------------------
80! Arguments:
81  REAL,    INTENT(INOUT) :: masque(iip1,jjp1)   !--- Land-ocean mask
82  REAL,    INTENT(INOUT) :: phis  (iip1,jjp1)   !--- Ground geopotential
83!-------------------------------------------------------------------------------
84! Local variables:
85  CHARACTER(LEN=256) :: modname, fmt
86  INTEGER            :: i, j, l, ji, itau, iday
87  REAL               :: xpn, xps, time, phystep
88  REAL, DIMENSION(iip1,jjp1)       :: psol
89  REAL, DIMENSION(iip1,jjp1,llm+1) :: p3d
90  REAL, DIMENSION(iip1,jjp1,llm)   :: uvent, t3d, tpot, qsat, qd
91  REAL, DIMENSION(iip1,jjp1,llm)   :: pk, pls, y, masse
92  REAL, DIMENSION(iip1,jjm ,llm)   :: vvent
93  REAL, DIMENSION(ip1jm    ,llm)   :: pbarv
94  REAL, DIMENSION(ip1jmp1  ,llm)   :: pbaru, phi, w
95  REAL, DIMENSION(ip1jmp1)         :: pks
96  REAL, DIMENSION(iim)             :: xppn, xpps
97  REAL, ALLOCATABLE                :: q3d(:,:,:,:)
98!-------------------------------------------------------------------------------
99  modname='etat0dyn_netcdf'
100
101  deg2rad = pi/180.0
102  y(:,:,:)=0  !ym warning unitialized variable
103 
104! Compute psol AND tsol, knowing phis.
105!*******************************************************************************
106  CALL start_init_dyn(rlonv, rlatu, rlonu, rlatv, phis, psol)
107
108! Mid-levels pressure computation
109!*******************************************************************************
110  CALL pression(ip1jmp1, ap, bp, psol, p3d)             !--- Update p3d
111  IF(pressure_exner) THEN                               !--- Update pk, pks
112    CALL exner_hyb   (ip1jmp1,psol,p3d,pks,pk)
113  ELSE
114    CALL exner_milieu(ip1jmp1,psol,p3d,pks,pk)
115  END IF
116  pls(:,:,:)=preff*(pk(:,:,:)/cpp)**(1./kappa)          !--- Update pls
117
118! Update uvent, vvent, t3d and tpot
119!*******************************************************************************
120  uvent(:,:,:) = 0.0 ; vvent(:,:,:) = 0.0 ; t3d (:,:,:) = 0.0
121  CALL startget_dyn3d('u'   ,rlonu,rlatu,pls,y ,uvent,rlonv,rlatv)
122  CALL startget_dyn3d('v'   ,rlonv,rlatv,pls(:,:jjm,:),y(:,:jjm,:),vvent,      &
123 &                           rlonu,rlatu(:jjm))
124  CALL startget_dyn3d('t'   ,rlonv,rlatu,pls,y ,t3d ,rlonu,rlatv)
125  tpot(:,:,:)=t3d(:,:,:)
126  CALL startget_dyn3d('tpot',rlonv,rlatu,pls,pk,tpot,rlonu,rlatv)
127
128  WRITE(lunout,*) 'T3D min,max:',MINVAL(t3d(:,:,:)),MAXVAL(t3d(:,:,:))
129  WRITE(lunout,*) 'PLS min,max:',MINVAL(pls(:,:,:)),MAXVAL(pls(:,:,:))
130
131! Humidity at saturation computation
132!*******************************************************************************
133  WRITE(lunout,*) 'avant q_sat'
134  CALL q_sat(llm*jjp1*iip1, t3d, pls, qsat)
135  WRITE(lunout,*) 'apres q_sat'
136  WRITE(lunout,*) 'QSAT min,max:',MINVAL(qsat(:,:,:)),MAXVAL(qsat(:,:,:))
137!  WRITE(lunout,*) 'QSAT :',qsat(10,20,:)
138  qd (:,:,:) = 0.0
139  CALL startget_dyn3d('q',rlonv,rlatu,pls,qsat,qd,rlonu,rlatv)
140  ALLOCATE(q3d(iip1,jjp1,llm,nqtot)); q3d(:,:,:,:)=0.0 ; q3d(:,:,:,1)=qd(:,:,:)
141  CALL flinclo(fid_dyn)
142
143! Parameterization of ozone chemistry:
144!*******************************************************************************
145! Look for ozone tracer:
146#ifndef INCA
147  DO i=1,nqtot; IF(ANY(["O3","o3"]==tname(i))) EXIT; END DO
148  IF(i/=nqtot+1) THEN
149    CALL regr_lat_time_coefoz
150    CALL press_coefoz
151    CALL regr_pr_o3(p3d, q3d(:,:,:,i))
152    q3d(:,:,:,i)=q3d(:,:,:,i)*48./ 29.                  !--- Mole->mass fraction         
153  END IF
154#endif
155  q3d(iip1,:,:,:)=q3d(1,:,:,:)
156
157! Writing
158!*******************************************************************************
159  CALL inidissip(lstardis, nitergdiv, nitergrot, niterh, tetagdiv, tetagrot,   &
160       tetatemp, vert_prof_dissip)
161  WRITE(lunout,*)'sortie inidissip'
162  itau=0
163  itau_dyn=0
164  itau_phy=0
165  iday=dayref+itau/day_step
166  time=FLOAT(itau-(iday-dayref)*day_step)/day_step
167  IF(time>1.) THEN
168   time=time-1
169   iday=iday+1
170  END IF
171  day_ref=dayref
172  annee_ref=anneeref
173  CALL geopot( ip1jmp1, tpot, pk, pks, phis, phi )
174  WRITE(lunout,*)'sortie geopot'
175  CALL caldyn0( itau, uvent, vvent, tpot, psol, masse, pk, phis,               &
176                phi,  w, pbaru, pbarv, time+iday-dayref)
177  WRITE(lunout,*)'sortie caldyn0'
178  start_time = 0.
179#ifdef CPP_PARA
180  CALL dynredem0_loc( "start.nc", dayref, phis)
181#else
182  CALL dynredem0( "start.nc", dayref, phis)
183#endif
184  WRITE(lunout,*)'sortie dynredem0'
185#ifdef CPP_PARA
186  CALL dynredem1_loc( "start.nc", 0.0, vvent, uvent, tpot, q3d, masse, psol)
187#else
188  CALL dynredem1( "start.nc", 0.0, vvent, uvent, tpot, q3d, masse, psol)
189#endif
190  WRITE(lunout,*)'sortie dynredem1'
191  CALL histclo()
192
193END SUBROUTINE etat0dyn_netcdf
194!
195!-------------------------------------------------------------------------------
196
197
198!-------------------------------------------------------------------------------
199!
200SUBROUTINE startget_dyn3d(var, lon_in,  lat_in,  pls,  workvar,&
201                        champ, lon_in2, lat_in2)
202!-------------------------------------------------------------------------------
203  IMPLICIT NONE
204!===============================================================================
205! Purpose: Compute some quantities (u,v,t,q,tpot) using variables U,V,TEMP and R
206!     (3D fields) of file dynfname.
207!-------------------------------------------------------------------------------
208! Note: An input auxilliary field "workvar" has to be specified in two cases:
209!     * for "q":    the saturated humidity.
210!     * for "tpot": the Exner function.
211!===============================================================================
212! Arguments:
213  CHARACTER(LEN=*), INTENT(IN)    :: var
214  REAL,             INTENT(IN)    :: lon_in(:)        ! dim (iml)
215  REAL,             INTENT(IN)    :: lat_in(:)        ! dim (jml)
216  REAL,             INTENT(IN)    :: pls    (:, :, :) ! dim (iml, jml, lml)
217  REAL,             INTENT(IN)    :: workvar(:, :, :) ! dim (iml, jml, lml)
218  REAL,             INTENT(INOUT) :: champ  (:, :, :) ! dim (iml, jml, lml)
219  REAL,             INTENT(IN)    :: lon_in2(:)       ! dim (iml)
220  REAL,             INTENT(IN)    :: lat_in2(:)       ! dim (jml2)
221!-------------------------------------------------------------------------------
222! Local variables:
223  CHARACTER(LEN=10)  :: vname
224  CHARACTER(LEN=256) :: msg, modname="startget_dyn3d"
225  INTEGER            :: iml, jml, jml2, lml, il
226  REAL               :: xppn, xpps
227!-------------------------------------------------------------------------------
228  iml=assert_eq([SIZE(lon_in),SIZE(pls,1),SIZE(workvar,1),SIZE(champ,1),       &
229     &                                    SIZE(lon_in2)], TRIM(modname)//" iml")
230  jml=assert_eq( SIZE(lat_in),SIZE(pls,2),SIZE(workvar,2),SIZE(champ,2),       &
231     &                                                    TRIM(modname)//" jml")
232  lml=assert_eq(              SIZE(pls,3),SIZE(workvar,3),SIZE(champ,3),       &
233     &                                                    TRIM(modname)//" lml")
234  jml2=SIZE(lat_in2)
235
236!--- CHECK IF THE FIELD IS KNOWN
237   SELECT CASE(var)
238    CASE('u');    vname='U'
239    CASE('v');    vname='V'
240    CASE('t');    vname='TEMP'
241    CASE('q');    vname='R';    msg='humidity as the saturated humidity'
242    CASE('tpot'); msg='potential temperature as the Exner function'
243    CASE DEFAULT; msg='No rule to extract variable '//TRIM(var)
244      CALL abort_gcm(modname,TRIM(msg)//' from any data set',1)
245  END SELECT
246
247!--- CHECK IF SOMETHING IS MISSING
248  IF((var=='tpot'.OR.var=='q').AND.MINVAL(workvar)==MAXVAL(workvar)) THEN
249    msg='Could not compute '//TRIM(msg)//' is missing or constant.'
250    CALL abort_gcm(modname,TRIM(msg),1)
251  END IF
252
253!--- INTERPOLATE 3D FIELD IF NEEDED
254  IF(var/='tpot') CALL start_inter_3d(TRIM(vname),lon_in,lat_in,lon_in2,      &
255                                                  lat_in2,pls,champ)
256
257!--- COMPUTE THE REQUIRED FILED
258  SELECT CASE(var)
259    CASE('u'); DO il=1,lml; champ(:,:,il)=champ(:,:,il)*cu(:,1:jml); END DO
260      champ(iml,:,:)=champ(1,:,:)                   !--- Eastward wind
261
262    CASE('v'); DO il=1,lml; champ(:,:,il)=champ(:,:,il)*cv(:,1:jml); END DO
263      champ(iml,:,:)=champ(1,:,:)                   !--- Northward wind
264
265    CASE('tpot','q')
266      IF(var=='tpot') THEN; champ=champ*cpp/workvar !--- Potential temperature
267      ELSE;                 champ=champ*.01*workvar !--- Relative humidity
268        WHERE(champ<0.) champ=1.0E-10
269      END IF
270      DO il=1,lml
271        xppn = SUM(aire(:,1  )*champ(:,1  ,il))/apoln
272        xpps = SUM(aire(:,jml)*champ(:,jml,il))/apols
273        champ(:,1  ,il) = xppn
274        champ(:,jml,il) = xpps
275      END DO
276  END SELECT
277
278END SUBROUTINE startget_dyn3d
279!
280!-------------------------------------------------------------------------------
281
282
283!-------------------------------------------------------------------------------
284!
285SUBROUTINE start_init_dyn(lon_in,lat_in,lon_in2,lat_in2,zs,psol)
286!
287!-------------------------------------------------------------------------------
288  IMPLICIT NONE
289!===============================================================================
290! Purpose:   Compute psol, knowing phis.
291!===============================================================================
292! Arguments:
293  REAL,    INTENT(IN)  :: lon_in (:),  lat_in (:)    ! dim (iml) (jml)
294  REAL,    INTENT(IN)  :: lon_in2(:),  lat_in2(:)    ! dim (iml) (jml2)
295  REAL,    INTENT(IN)  :: zs  (:,:)                  ! dim (iml,jml)
296  REAL,    INTENT(OUT) :: psol(:,:)                  ! dim (iml,jml)
297!-------------------------------------------------------------------------------
298! Local variables:
299  CHARACTER(LEN=256) :: modname='start_init_dyn'
300  REAL               :: date, dt
301  INTEGER            :: iml, jml, jml2, itau(1)
302  REAL, ALLOCATABLE  :: lon_rad(:), lon_ini(:), var_ana(:,:)
303  REAL, ALLOCATABLE  :: lat_rad(:), lat_ini(:)
304  REAL, ALLOCATABLE  :: z(:,:), ps(:,:), ts(:,:)
305!-------------------------------------------------------------------------------
306  iml=assert_eq(SIZE(lon_in),SIZE(zs,1),SIZE(psol,1),SIZE(lon_in2),            &
307      &                                              TRIM(modname)//" iml")
308  jml=assert_eq(SIZE(lat_in),SIZE(zs,2),SIZE(psol,2),TRIM(modname)//" jml")
309  jml2=SIZE(lat_in2)
310
311  WRITE(lunout,*) 'Opening the surface analysis'
312  CALL flininfo(dynfname, iml_dyn, jml_dyn, llm_dyn, ttm_dyn, fid_dyn)
313  WRITE(lunout,*) 'Values read: ', iml_dyn, jml_dyn, llm_dyn, ttm_dyn
314
315  ALLOCATE(lon_dyn(iml_dyn,jml_dyn), lat_dyn(iml_dyn,jml_dyn))
316  ALLOCATE(levdyn_ini(llm_dyn))
317  CALL flinopen(dynfname, .FALSE., iml_dyn, jml_dyn, llm_dyn,                  &
318                lon_dyn,lat_dyn,levdyn_ini,ttm_dyn,itau,date,dt,fid_dyn)
319
320!--- IF ANGLES ARE IN DEGREES, THEY ARE CONVERTED INTO RADIANS
321  ALLOCATE(lon_ini(iml_dyn),lat_ini(jml_dyn))
322  lon_ini(:)=lon_dyn(:,1); IF(MAXVAL(lon_dyn)>pi) lon_ini=lon_ini*deg2rad
323  lat_ini(:)=lat_dyn(1,:); IF(MAXVAL(lat_dyn)>pi) lat_ini=lat_ini*deg2rad
324
325  ALLOCATE(var_ana(iml_dyn,jml_dyn),lon_rad(iml_dyn),lat_rad(jml_dyn))
326  CALL get_var_dyn('Z',z)                        !--- SURFACE GEOPOTENTIAL
327  CALL get_var_dyn('SP',ps)                      !--- SURFACE PRESSURE
328  CALL get_var_dyn('ST',ts)                      !--- SURFACE TEMPERATURE
329!  CALL flinclo(fid_dyn)
330  DEALLOCATE(var_ana,lon_rad,lat_rad,lon_ini,lat_ini)
331
332!--- PSOL IS COMPUTED IN PASCALS
333  psol(:iml-1,:) = ps(:iml-1,:)*(1.0+(z(:iml-1,:)-zs(:iml-1,:))/287.0          &
334     &            /ts(:iml-1,:))
335  psol(iml,:)=psol(1,:)
336  DEALLOCATE(z,ps,ts)
337  psol(:,1  )=SUM(aire(1:iml-1,1  )*psol(1:iml-1,1  ))/apoln  !--- NORTH POLE
338  psol(:,jml)=SUM(aire(1:iml-1,jml)*psol(1:iml-1,jml))/apols  !--- SOUTH POLE
339
340CONTAINS
341
342!-------------------------------------------------------------------------------
343!
344SUBROUTINE get_var_dyn(title,field)
345!
346!-------------------------------------------------------------------------------
347  USE conf_dat_m, ONLY: conf_dat2d
348  IMPLICIT NONE
349!-------------------------------------------------------------------------------
350! Arguments:
351  CHARACTER(LEN=*),  INTENT(IN)    :: title
352  REAL, ALLOCATABLE, INTENT(INOUT) :: field(:,:)
353!-------------------------------------------------------------------------------
354! Local variables:
355  CHARACTER(LEN=256) :: msg
356  INTEGER :: tllm
357!-------------------------------------------------------------------------------
358  SELECT CASE(title)
359    CASE('Z');     tllm=0;       msg='geopotential'
360    CASE('SP');    tllm=0;       msg='surface pressure'
361    CASE('ST');    tllm=llm_dyn; msg='temperature'
362  END SELECT
363  IF(.NOT.ALLOCATED(field)) THEN
364    ALLOCATE(field(iml,jml))
365    CALL flinget(fid_dyn, title, iml_dyn,jml_dyn, tllm, ttm_dyn, 1, 1, var_ana)
366    CALL conf_dat2d(title, lon_ini, lat_ini, lon_rad, lat_rad, var_ana, .TRUE.)
367    CALL interp_startvar(title, .TRUE., lon_rad,lat_rad, var_ana,              &
368                                        lon_in, lat_in, lon_in2, lat_in2, field)
369  ELSE IF(SIZE(field)/=SIZE(z)) THEN
370    msg='The '//TRIM(msg)//' field we have does not have the right size'
371    CALL abort_gcm(TRIM(modname),msg,1)
372  END IF
373
374END SUBROUTINE get_var_dyn
375!
376!-------------------------------------------------------------------------------
377
378END SUBROUTINE start_init_dyn
379!
380!-------------------------------------------------------------------------------
381
382
383!-------------------------------------------------------------------------------
384!
385SUBROUTINE start_inter_3d(var,lon_in,lat_in,lon_in2,lat_in2,pls_in,var3d)
386!
387!-------------------------------------------------------------------------------
388  USE conf_dat_m, ONLY: conf_dat3d
389  USE pchsp_95_m, ONLY: pchsp_95
390  USE pchfe_95_m, ONLY: pchfe_95
391  IMPLICIT NONE
392!-------------------------------------------------------------------------------
393! Arguments:
394  CHARACTER(LEN=*), INTENT(IN) :: var
395  REAL,    INTENT(IN)  :: lon_in(:),  lat_in(:)   ! dim (iml) (jml)
396  REAL,    INTENT(IN)  :: lon_in2(:), lat_in2(:)  ! dim (iml) (jml2)
397  REAL,    INTENT(IN)  :: pls_in(:,:,:)           ! dim (iml,jml,lml)
398  REAL,    INTENT(OUT) :: var3d (:,:,:)           ! dim (iml,jml,lml)
399!-------------------------------------------------------------------------------
400! Local variables:
401  CHARACTER(LEN=256) :: modname='start_inter_3d'
402  LOGICAL :: skip
403  REAL    :: chmin, chmax
404  INTEGER :: iml, jml, lml, jml2, ii, ij, il, ierr
405  INTEGER :: n_extrap                             ! Extrapolated points number
406  REAL, ALLOCATABLE :: ax(:), lon_rad(:), lon_ini(:), lev_dyn(:), yder(:)
407  REAL, ALLOCATABLE :: ay(:), lat_rad(:), lat_ini(:), var_tmp3d(:,:,:)
408  REAL, ALLOCATABLE, SAVE :: var_ana3d(:,:,:)
409!-------------------------------------------------------------------------------
410  iml=assert_eq(SIZE(lon_in),SIZE(lon_in2),SIZE(pls_in,1),SIZE(var3d,1),TRIM(modname)//" iml")
411  jml=assert_eq(SIZE(lat_in),              SIZE(pls_in,2),SIZE(var3d,2),TRIM(modname)//" jml")
412  lml=assert_eq(SIZE(pls_in,3),SIZE(var3d,3),TRIM(modname)//" lml"); jml2=SIZE(lat_in2)
413
414  WRITE(lunout, *)'Going into flinget to extract the 3D field.'
415  IF(.NOT.ALLOCATED(var_ana3d)) ALLOCATE(var_ana3d(iml_dyn, jml_dyn, llm_dyn))
416  CALL flinget(fid_dyn,var,iml_dyn,jml_dyn,llm_dyn,ttm_dyn,1,1,var_ana3d)
417
418!--- ANGLES IN DEGREES ARE CONVERTED INTO RADIANS
419  ALLOCATE(lon_ini(iml_dyn), lat_ini(jml_dyn))
420  lon_ini(:)=lon_dyn(:,1); IF(MAXVAL(lon_dyn)>pi) lon_ini=lon_ini*deg2rad
421  lat_ini(:)=lat_dyn(1,:); IF(MAXVAL(lat_dyn)>pi) lat_ini=lat_ini*deg2rad
422
423!--- FIELDS ARE PROCESSED TO BE ON STANDARD ANGULAR DOMAINS
424  ALLOCATE(lon_rad(iml_dyn), lat_rad(jml_dyn), lev_dyn(llm_dyn))
425  CALL conf_dat3d(var, lon_ini, lat_ini, levdyn_ini,                           &
426                       lon_rad, lat_rad, lev_dyn, var_ana3d, .TRUE.)
427  DEALLOCATE(lon_ini, lat_ini)
428
429!--- COMPUTE THE REQUIRED FIELDS USING ROUTINE grid_noro
430  ALLOCATE(var_tmp3d(iml,jml,llm_dyn))
431  DO il = 1,llm_dyn
432    CALL interp_startvar(var,il==1,lon_rad,lat_rad,var_ana3d(:,:,il),          &
433                          lon_in,lat_in,lon_in2,lat_in2,var_tmp3d(:,:,il))
434  END DO
435  DEALLOCATE(lon_rad, lat_rad)
436
437!--- VERTICAL INTERPOLATION FROM TOP OF ATMOSPHERE TO GROUND
438  ALLOCATE(ax(llm_dyn),ay(llm_dyn),yder(llm_dyn))
439  ax = lev_dyn(llm_dyn:1:-1)
440  skip = .FALSE.
441  n_extrap = 0
442  DO ij=1, jml
443    DO ii=1, iml-1
444      ay = var_tmp3d(ii, ij, llm_dyn:1:-1)
445      yder = pchsp_95(ax, ay, ibeg=2, iend=2, vc_beg=0., vc_end=0.)
446      CALL pchfe_95(ax, ay, yder, skip, pls_in(ii, ij, lml:1:-1),              &
447           var3d(ii, ij, lml:1:-1), ierr)
448      IF(ierr<0) CALL abort_gcm(TRIM(modname),'error in pchfe_95',1)
449      n_extrap = n_extrap + ierr
450    END DO
451  END DO
452  IF(n_extrap/=0) WRITE(lunout,*)TRIM(modname)//" pchfe_95: n_extrap=", n_extrap
453  var3d(iml, :, :) = var3d(1, :, :)
454
455  DO il=1, lml
456    CALL minmax(iml*jml, var3d(1, 1, il), chmin, chmax)
457    WRITE(lunout, *)' '//TRIM(var)//'  min max l ', il, chmin, chmax
458  END DO
459
460END SUBROUTINE start_inter_3d
461!
462!-------------------------------------------------------------------------------
463
464
465!-------------------------------------------------------------------------------
466!
467SUBROUTINE interp_startvar(nam,ibeg,lon,lat,vari,lon1,lat1,lon2,lat2,varo)
468!
469!-------------------------------------------------------------------------------
470  USE inter_barxy_m, ONLY: inter_barxy
471  IMPLICIT NONE
472!-------------------------------------------------------------------------------
473! Arguments:
474  CHARACTER(LEN=*), INTENT(IN)  :: nam
475  LOGICAL,          INTENT(IN)  :: ibeg
476  REAL,             INTENT(IN)  :: lon(:), lat(:)   ! dim (ii) (jj)
477  REAL,             INTENT(IN)  :: vari(:,:)        ! dim (ii,jj)
478  REAL,             INTENT(IN)  :: lon1(:), lat1(:) ! dim (i1) (j1)
479  REAL,             INTENT(IN)  :: lon2(:), lat2(:) ! dim (i1) (j2)
480  REAL,             INTENT(OUT) :: varo(:,:)        ! dim (i1) (j1)
481!-------------------------------------------------------------------------------
482! Local variables:
483  CHARACTER(LEN=256) :: modname="interp_startvar"
484  INTEGER            :: ii, jj, i1, j1, j2
485  REAL, ALLOCATABLE  :: vtmp(:,:)
486!-------------------------------------------------------------------------------
487  ii=assert_eq(SIZE(lon),            SIZE(vari,1),TRIM(modname)//" ii")
488  jj=assert_eq(SIZE(lat),            SIZE(vari,2),TRIM(modname)//" jj")
489  i1=assert_eq(SIZE(lon1),SIZE(lon2),SIZE(varo,1),TRIM(modname)//" i1")
490  j1=assert_eq(SIZE(lat1),           SIZE(varo,2),TRIM(modname)//" j1")
491  j2=SIZE(lat2)
492  ALLOCATE(vtmp(i1-1,j1))
493  IF(ibeg.AND.prt_level>1) THEN
494    WRITE(lunout,*)"---------------------------------------------------------"
495    WRITE(lunout,*)"$$$ Interpolation barycentrique pour "//TRIM(nam)//" $$$"
496    WRITE(lunout,*)"---------------------------------------------------------"
497  END IF
498  CALL inter_barxy(lon, lat(:jj-1), vari, lon2(:i1-1), lat2, vtmp)
499  CALL gr_int_dyn(vtmp, varo, i1-1, j1)
500
501END SUBROUTINE interp_startvar
502!
503!-------------------------------------------------------------------------------
504
505END MODULE etat0dyn
506!
507!*******************************************************************************
Note: See TracBrowser for help on using the repository browser.