source: LMDZ5/trunk/libf/dynphy_lonlat/phylmd/etat0dyn_netcdf.F90 @ 2588

Last change on this file since 2588 was 2361, checked in by dcugnet, 9 years ago

In etat0dyn: removed few useless lines: "masque" is always known because etat0dyn is called after etat0phys.
In grid_atob: shape error in routine fine2coarse now fixed: "msk" argument and local variable mask must have the dimensions of the output grid. Working unit of dist_sphe is no longer degree, but radian.

File size: 21.7 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  IMPLICIT NONE
39
40  PRIVATE
41  PUBLIC :: etat0dyn_netcdf
42
43  include "iniprint.h"
44  include "dimensions.h"
45  include "paramet.h"
46  include "comgeom2.h"
47  include "comvert.h"
48  include "comconst.h"
49  include "temps.h"
50  include "comdissnew.h"
51  include "serre.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
103! Compute psol AND tsol, knowing phis.
104!*******************************************************************************
105  CALL start_init_dyn(rlonv, rlatu, rlonu, rlatv, phis, psol)
106
107! Mid-levels pressure computation
108!*******************************************************************************
109  CALL pression(ip1jmp1, ap, bp, psol, p3d)             !--- Update p3d
110  IF(pressure_exner) THEN                               !--- Update pk, pks
111    CALL exner_hyb   (ip1jmp1,psol,p3d,pks,pk)
112  ELSE
113    CALL exner_milieu(ip1jmp1,psol,p3d,pks,pk)
114  END IF
115  pls(:,:,:)=preff*(pk(:,:,:)/cpp)**(1./kappa)          !--- Update pls
116
117! Update uvent, vvent, t3d and tpot
118!*******************************************************************************
119  uvent(:,:,:) = 0.0 ; vvent(:,:,:) = 0.0 ; t3d (:,:,:) = 0.0
120  CALL startget_dyn3d('u'   ,rlonu,rlatu,pls,y ,uvent,rlonv,rlatv)
121  CALL startget_dyn3d('v'   ,rlonv,rlatv,pls(:,:jjm,:),y(:,:jjm,:),vvent,      &
122 &                           rlonu,rlatu(:jjm))
123  CALL startget_dyn3d('t'   ,rlonv,rlatu,pls,y ,t3d ,rlonu,rlatv)
124  tpot(:,:,:)=t3d(:,:,:)
125  CALL startget_dyn3d('tpot',rlonv,rlatu,pls,pk,tpot,rlonu,rlatv)
126
127  WRITE(lunout,*) 'T3D min,max:',MINVAL(t3d(:,:,:)),MAXVAL(t3d(:,:,:))
128  WRITE(lunout,*) 'PLS min,max:',MINVAL(pls(:,:,:)),MAXVAL(pls(:,:,:))
129
130! Humidity at saturation computation
131!*******************************************************************************
132  WRITE(lunout,*) 'avant q_sat'
133  CALL q_sat(llm*jjp1*iip1, t3d, pls, qsat)
134  WRITE(lunout,*) 'apres q_sat'
135  WRITE(lunout,*) 'QSAT min,max:',MINVAL(qsat(:,:,:)),MAXVAL(qsat(:,:,:))
136!  WRITE(lunout,*) 'QSAT :',qsat(10,20,:)
137  qd (:,:,:) = 0.0
138  CALL startget_dyn3d('q',rlonv,rlatu,pls,qsat,qd,rlonu,rlatv)
139  ALLOCATE(q3d(iip1,jjp1,llm,nqtot)); q3d(:,:,:,:)=0.0 ; q3d(:,:,:,1)=qd(:,:,:)
140  CALL flinclo(fid_dyn)
141
142! Parameterization of ozone chemistry:
143!*******************************************************************************
144! Look for ozone tracer:
145  DO i=1,nqtot; IF(ANY(["O3","o3"]==tname(i))) EXIT; END DO
146  IF(i/=nqtot+1) THEN
147    CALL regr_lat_time_coefoz
148    CALL press_coefoz
149    CALL regr_pr_o3(p3d, q3d(:,:,:,i))
150    q3d(:,:,:,i)=q3d(:,:,:,i)*48./ 29.                  !--- Mole->mass fraction         
151  END IF
152  q3d(iip1,:,:,:)=q3d(1,:,:,:)
153
154! Writing
155!*******************************************************************************
156  CALL inidissip(lstardis, nitergdiv, nitergrot, niterh, tetagdiv, tetagrot,   &
157       tetatemp, vert_prof_dissip)
158  WRITE(lunout,*)'sortie inidissip'
159  itau=0
160  itau_dyn=0
161  itau_phy=0
162  iday=dayref+itau/day_step
163  time=FLOAT(itau-(iday-dayref)*day_step)/day_step
164  IF(time>1.) THEN
165   time=time-1
166   iday=iday+1
167  END IF
168  day_ref=dayref
169  annee_ref=anneeref
170  CALL geopot( ip1jmp1, tpot, pk, pks, phis, phi )
171  WRITE(lunout,*)'sortie geopot'
172  CALL caldyn0( itau, uvent, vvent, tpot, psol, masse, pk, phis,               &
173                phi,  w, pbaru, pbarv, time+iday-dayref)
174  WRITE(lunout,*)'sortie caldyn0'     
175#ifdef CPP_PARA
176  CALL dynredem0_loc( "start.nc", dayref, phis)
177#else
178  CALL dynredem0( "start.nc", dayref, phis)
179#endif
180  WRITE(lunout,*)'sortie dynredem0'
181#ifdef CPP_PARA
182  CALL dynredem1_loc( "start.nc", 0.0, vvent, uvent, tpot, q3d, masse, psol)
183#else
184  CALL dynredem1( "start.nc", 0.0, vvent, uvent, tpot, q3d, masse, psol)
185#endif
186  WRITE(lunout,*)'sortie dynredem1'
187  CALL histclo()
188
189END SUBROUTINE etat0dyn_netcdf
190!
191!-------------------------------------------------------------------------------
192
193
194!-------------------------------------------------------------------------------
195!
196SUBROUTINE startget_dyn3d(var, lon_in,  lat_in,  pls,  workvar,&
197                        champ, lon_in2, lat_in2)
198!-------------------------------------------------------------------------------
199  IMPLICIT NONE
200!===============================================================================
201! Purpose: Compute some quantities (u,v,t,q,tpot) using variables U,V,TEMP and R
202!     (3D fields) of file dynfname.
203!-------------------------------------------------------------------------------
204! Note: An input auxilliary field "workvar" has to be specified in two cases:
205!     * for "q":    the saturated humidity.
206!     * for "tpot": the Exner function.
207!===============================================================================
208! Arguments:
209  CHARACTER(LEN=*), INTENT(IN)    :: var
210  REAL,             INTENT(IN)    :: lon_in(:)        ! dim (iml)
211  REAL,             INTENT(IN)    :: lat_in(:)        ! dim (jml)
212  REAL,             INTENT(IN)    :: pls    (:, :, :) ! dim (iml, jml, lml)
213  REAL,             INTENT(IN)    :: workvar(:, :, :) ! dim (iml, jml, lml)
214  REAL,             INTENT(INOUT) :: champ  (:, :, :) ! dim (iml, jml, lml)
215  REAL,             INTENT(IN)    :: lon_in2(:)       ! dim (iml)
216  REAL,             INTENT(IN)    :: lat_in2(:)       ! dim (jml2)
217!-------------------------------------------------------------------------------
218! Local variables:
219  CHARACTER(LEN=10)  :: vname
220  CHARACTER(LEN=256) :: msg, modname="startget_dyn3d"
221  INTEGER            :: iml, jml, jml2, lml, il
222  REAL               :: xppn, xpps
223!-------------------------------------------------------------------------------
224  iml=assert_eq([SIZE(lon_in),SIZE(pls,1),SIZE(workvar,1),SIZE(champ,1),       &
225     &                                    SIZE(lon_in2)], TRIM(modname)//" iml")
226  jml=assert_eq( SIZE(lat_in),SIZE(pls,2),SIZE(workvar,2),SIZE(champ,2),       &
227     &                                                    TRIM(modname)//" jml")
228  lml=assert_eq(              SIZE(pls,3),SIZE(workvar,3),SIZE(champ,3),       &
229     &                                                    TRIM(modname)//" lml")
230  jml2=SIZE(lat_in2)
231
232!--- CHECK IF THE FIELD IS KNOWN
233   SELECT CASE(var)
234    CASE('u');    vname='U'
235    CASE('v');    vname='V'
236    CASE('t');    vname='TEMP'
237    CASE('q');    vname='R';    msg='humidity as the saturated humidity'
238    CASE('tpot'); msg='potential temperature as the Exner function'
239    CASE DEFAULT; msg='No rule to extract variable '//TRIM(var)
240      CALL abort_gcm(modname,TRIM(msg)//' from any data set',1)
241  END SELECT
242
243!--- CHECK IF SOMETHING IS MISSING
244  IF((var=='tpot'.OR.var=='q').AND.MINVAL(workvar)==MAXVAL(workvar)) THEN
245    msg='Could not compute '//TRIM(msg)//' is missing or constant.'
246    CALL abort_gcm(modname,TRIM(msg),1)
247  END IF
248
249!--- INTERPOLATE 3D FIELD IF NEEDED
250  IF(var/='tpot') CALL start_inter_3d(TRIM(vname),lon_in,lat_in,lon_in2,      &
251                                                  lat_in2,pls,champ)
252
253!--- COMPUTE THE REQUIRED FILED
254  SELECT CASE(var)
255    CASE('u'); DO il=1,lml; champ(:,:,il)=champ(:,:,il)*cu(:,1:jml); END DO
256      champ(iml,:,:)=champ(1,:,:)                   !--- Eastward wind
257
258    CASE('v'); DO il=1,lml; champ(:,:,il)=champ(:,:,il)*cv(:,1:jml); END DO
259      champ(iml,:,:)=champ(1,:,:)                   !--- Northward wind
260
261    CASE('tpot','q')
262      IF(var=='tpot') THEN; champ=champ*cpp/workvar !--- Potential temperature
263      ELSE;                 champ=champ*.01*workvar !--- Relative humidity
264        WHERE(champ<0.) champ=1.0E-10
265      END IF
266      DO il=1,lml
267        xppn = SUM(aire(:,1  )*champ(:,1  ,il))/apoln
268        xpps = SUM(aire(:,jml)*champ(:,jml,il))/apols
269        champ(:,1  ,il) = xppn
270        champ(:,jml,il) = xpps
271      END DO
272  END SELECT
273
274END SUBROUTINE startget_dyn3d
275!
276!-------------------------------------------------------------------------------
277
278
279!-------------------------------------------------------------------------------
280!
281SUBROUTINE start_init_dyn(lon_in,lat_in,lon_in2,lat_in2,zs,psol)
282!
283!-------------------------------------------------------------------------------
284  IMPLICIT NONE
285!===============================================================================
286! Purpose:   Compute psol, knowing phis.
287!===============================================================================
288! Arguments:
289  REAL,    INTENT(IN)  :: lon_in (:),  lat_in (:)    ! dim (iml) (jml)
290  REAL,    INTENT(IN)  :: lon_in2(:),  lat_in2(:)    ! dim (iml) (jml2)
291  REAL,    INTENT(IN)  :: zs  (:,:)                  ! dim (iml,jml)
292  REAL,    INTENT(OUT) :: psol(:,:)                  ! dim (iml,jml)
293!-------------------------------------------------------------------------------
294! Local variables:
295  CHARACTER(LEN=256) :: modname='start_init_dyn'
296  REAL               :: date, dt
297  INTEGER            :: iml, jml, jml2, itau(1)
298  REAL, ALLOCATABLE  :: lon_rad(:), lon_ini(:), var_ana(:,:)
299  REAL, ALLOCATABLE  :: lat_rad(:), lat_ini(:)
300  REAL, ALLOCATABLE  :: z(:,:), ps(:,:), ts(:,:)
301!-------------------------------------------------------------------------------
302  iml=assert_eq(SIZE(lon_in),SIZE(zs,1),SIZE(psol,1),SIZE(lon_in2),            &
303      &                                              TRIM(modname)//" iml")
304  jml=assert_eq(SIZE(lat_in),SIZE(zs,2),SIZE(psol,2),TRIM(modname)//" jml")
305  jml2=SIZE(lat_in2)
306
307  WRITE(lunout,*) 'Opening the surface analysis'
308  CALL flininfo(dynfname, iml_dyn, jml_dyn, llm_dyn, ttm_dyn, fid_dyn)
309  WRITE(lunout,*) 'Values read: ', iml_dyn, jml_dyn, llm_dyn, ttm_dyn
310
311  ALLOCATE(lon_dyn(iml_dyn,jml_dyn), lat_dyn(iml_dyn,jml_dyn))
312  ALLOCATE(levdyn_ini(llm_dyn))
313  CALL flinopen(dynfname, .FALSE., iml_dyn, jml_dyn, llm_dyn,                  &
314                lon_dyn,lat_dyn,levdyn_ini,ttm_dyn,itau,date,dt,fid_dyn)
315
316!--- IF ANGLES ARE IN DEGREES, THEY ARE CONVERTED INTO RADIANS
317  ALLOCATE(lon_ini(iml_dyn),lat_ini(jml_dyn))
318  lon_ini(:)=lon_dyn(:,1); IF(MAXVAL(lon_dyn)>pi) lon_ini=lon_ini*deg2rad
319  lat_ini(:)=lat_dyn(1,:); IF(MAXVAL(lat_dyn)>pi) lat_ini=lat_ini*deg2rad
320
321  ALLOCATE(var_ana(iml_dyn,jml_dyn),lon_rad(iml_dyn),lat_rad(jml_dyn))
322  CALL get_var_dyn('Z',z)                        !--- SURFACE GEOPOTENTIAL
323  CALL get_var_dyn('SP',ps)                      !--- SURFACE PRESSURE
324  CALL get_var_dyn('ST',ts)                      !--- SURFACE TEMPERATURE
325!  CALL flinclo(fid_dyn)
326  DEALLOCATE(var_ana,lon_rad,lat_rad,lon_ini,lat_ini)
327
328!--- PSOL IS COMPUTED IN PASCALS
329  psol(:iml-1,:) = ps(:iml-1,:)*(1.0+(z(:iml-1,:)-zs(:iml-1,:))/287.0          &
330     &            /ts(:iml-1,:))
331  psol(iml,:)=psol(1,:)
332  DEALLOCATE(z,ps,ts)
333  psol(:,1  )=SUM(aire(1:iml-1,1  )*psol(1:iml-1,1  ))/apoln  !--- NORTH POLE
334  psol(:,jml)=SUM(aire(1:iml-1,jml)*psol(1:iml-1,jml))/apols  !--- SOUTH POLE
335
336CONTAINS
337
338!-------------------------------------------------------------------------------
339!
340SUBROUTINE get_var_dyn(title,field)
341!
342!-------------------------------------------------------------------------------
343  USE conf_dat_m, ONLY: conf_dat2d
344  IMPLICIT NONE
345!-------------------------------------------------------------------------------
346! Arguments:
347  CHARACTER(LEN=*),  INTENT(IN)    :: title
348  REAL, ALLOCATABLE, INTENT(INOUT) :: field(:,:)
349!-------------------------------------------------------------------------------
350! Local variables:
351  CHARACTER(LEN=256) :: msg
352  INTEGER :: tllm
353!-------------------------------------------------------------------------------
354  SELECT CASE(title)
355    CASE('Z');     tllm=0;       msg='geopotential'
356    CASE('SP');    tllm=0;       msg='surface pressure'
357    CASE('ST');    tllm=llm_dyn; msg='temperature'
358  END SELECT
359  IF(.NOT.ALLOCATED(field)) THEN
360    ALLOCATE(field(iml,jml))
361    CALL flinget(fid_dyn, title, iml_dyn,jml_dyn, tllm, ttm_dyn, 1, 1, var_ana)
362    CALL conf_dat2d(title, lon_ini, lat_ini, lon_rad, lat_rad, var_ana, .TRUE.)
363    CALL interp_startvar(title, .TRUE., lon_rad,lat_rad, var_ana,              &
364                                        lon_in, lat_in, lon_in2, lat_in2, field)
365  ELSE IF(SIZE(field)/=SIZE(z)) THEN
366    msg='The '//TRIM(msg)//' field we have does not have the right size'
367    CALL abort_gcm(TRIM(modname),msg,1)
368  END IF
369
370END SUBROUTINE get_var_dyn
371!
372!-------------------------------------------------------------------------------
373
374END SUBROUTINE start_init_dyn
375!
376!-------------------------------------------------------------------------------
377
378
379!-------------------------------------------------------------------------------
380!
381SUBROUTINE start_inter_3d(var,lon_in,lat_in,lon_in2,lat_in2,pls_in,var3d)
382!
383!-------------------------------------------------------------------------------
384  USE conf_dat_m, ONLY: conf_dat3d
385  USE pchsp_95_m, ONLY: pchsp_95
386  USE pchfe_95_m, ONLY: pchfe_95
387  IMPLICIT NONE
388!-------------------------------------------------------------------------------
389! Arguments:
390  CHARACTER(LEN=*), INTENT(IN) :: var
391  REAL,    INTENT(IN)  :: lon_in(:),  lat_in(:)   ! dim (iml) (jml)
392  REAL,    INTENT(IN)  :: lon_in2(:), lat_in2(:)  ! dim (iml) (jml2)
393  REAL,    INTENT(IN)  :: pls_in(:,:,:)           ! dim (iml,jml,lml)
394  REAL,    INTENT(OUT) :: var3d (:,:,:)           ! dim (iml,jml,lml)
395!-------------------------------------------------------------------------------
396! Local variables:
397  CHARACTER(LEN=256) :: modname='start_inter_3d'
398  LOGICAL :: skip
399  REAL    :: chmin, chmax
400  INTEGER :: iml, jml, lml, jml2, ii, ij, il, ierr
401  INTEGER :: n_extrap                             ! Extrapolated points number
402  REAL, ALLOCATABLE :: ax(:), lon_rad(:), lon_ini(:), lev_dyn(:), yder(:)
403  REAL, ALLOCATABLE :: ay(:), lat_rad(:), lat_ini(:), var_tmp3d(:,:,:)
404  REAL, ALLOCATABLE, SAVE :: var_ana3d(:,:,:)
405!-------------------------------------------------------------------------------
406  iml=assert_eq(SIZE(lon_in),SIZE(lon_in2),SIZE(pls_in,1),SIZE(var3d,1),TRIM(modname)//" iml")
407  jml=assert_eq(SIZE(lat_in),              SIZE(pls_in,2),SIZE(var3d,2),TRIM(modname)//" jml")
408  lml=assert_eq(SIZE(pls_in,3),SIZE(var3d,3),TRIM(modname)//" lml"); jml2=SIZE(lat_in2)
409
410  WRITE(lunout, *)'Going into flinget to extract the 3D field.'
411  IF(.NOT.ALLOCATED(var_ana3d)) ALLOCATE(var_ana3d(iml_dyn, jml_dyn, llm_dyn))
412  CALL flinget(fid_dyn,var,iml_dyn,jml_dyn,llm_dyn,ttm_dyn,1,1,var_ana3d)
413
414!--- ANGLES IN DEGREES ARE CONVERTED INTO RADIANS
415  ALLOCATE(lon_ini(iml_dyn), lat_ini(jml_dyn))
416  lon_ini(:)=lon_dyn(:,1); IF(MAXVAL(lon_dyn)>pi) lon_ini=lon_ini*deg2rad
417  lat_ini(:)=lat_dyn(1,:); IF(MAXVAL(lat_dyn)>pi) lat_ini=lat_ini*deg2rad
418
419!--- FIELDS ARE PROCESSED TO BE ON STANDARD ANGULAR DOMAINS
420  ALLOCATE(lon_rad(iml_dyn), lat_rad(jml_dyn), lev_dyn(llm_dyn))
421  CALL conf_dat3d(var, lon_ini, lat_ini, levdyn_ini,                           &
422                       lon_rad, lat_rad, lev_dyn, var_ana3d, .TRUE.)
423  DEALLOCATE(lon_ini, lat_ini)
424
425!--- COMPUTE THE REQUIRED FIELDS USING ROUTINE grid_noro
426  ALLOCATE(var_tmp3d(iml,jml,llm_dyn))
427  DO il = 1,llm_dyn
428    CALL interp_startvar(var,il==1,lon_rad,lat_rad,var_ana3d(:,:,il),          &
429                          lon_in,lat_in,lon_in2,lat_in2,var_tmp3d(:,:,il))
430  END DO
431  DEALLOCATE(lon_rad, lat_rad)
432
433!--- VERTICAL INTERPOLATION FROM TOP OF ATMOSPHERE TO GROUND
434  ALLOCATE(ax(llm_dyn),ay(llm_dyn),yder(llm_dyn))
435  ax = lev_dyn(llm_dyn:1:-1)
436  skip = .FALSE.
437  n_extrap = 0
438  DO ij=1, jml
439    DO ii=1, iml-1
440      ay = var_tmp3d(ii, ij, llm_dyn:1:-1)
441      yder = pchsp_95(ax, ay, ibeg=2, iend=2, vc_beg=0., vc_end=0.)
442      CALL pchfe_95(ax, ay, yder, skip, pls_in(ii, ij, lml:1:-1),              &
443           var3d(ii, ij, lml:1:-1), ierr)
444      IF(ierr<0) CALL abort_gcm(TRIM(modname),'error in pchfe_95',1)
445      n_extrap = n_extrap + ierr
446    END DO
447  END DO
448  IF(n_extrap/=0) WRITE(lunout,*)TRIM(modname)//" pchfe_95: n_extrap=", n_extrap
449  var3d(iml, :, :) = var3d(1, :, :)
450
451  DO il=1, lml
452    CALL minmax(iml*jml, var3d(1, 1, il), chmin, chmax)
453    WRITE(lunout, *)' '//TRIM(var)//'  min max l ', il, chmin, chmax
454  END DO
455
456END SUBROUTINE start_inter_3d
457!
458!-------------------------------------------------------------------------------
459
460
461!-------------------------------------------------------------------------------
462!
463SUBROUTINE interp_startvar(nam,ibeg,lon,lat,vari,lon1,lat1,lon2,lat2,varo)
464!
465!-------------------------------------------------------------------------------
466  USE inter_barxy_m, ONLY: inter_barxy
467  IMPLICIT NONE
468!-------------------------------------------------------------------------------
469! Arguments:
470  CHARACTER(LEN=*), INTENT(IN)  :: nam
471  LOGICAL,          INTENT(IN)  :: ibeg
472  REAL,             INTENT(IN)  :: lon(:), lat(:)   ! dim (ii) (jj)
473  REAL,             INTENT(IN)  :: vari(:,:)        ! dim (ii,jj)
474  REAL,             INTENT(IN)  :: lon1(:), lat1(:) ! dim (i1) (j1)
475  REAL,             INTENT(IN)  :: lon2(:), lat2(:) ! dim (i1) (j2)
476  REAL,             INTENT(OUT) :: varo(:,:)        ! dim (i1) (j1)
477!-------------------------------------------------------------------------------
478! Local variables:
479  CHARACTER(LEN=256) :: modname="interp_startvar"
480  INTEGER            :: ii, jj, i1, j1, j2
481  REAL, ALLOCATABLE  :: vtmp(:,:)
482!-------------------------------------------------------------------------------
483  ii=assert_eq(SIZE(lon),            SIZE(vari,1),TRIM(modname)//" ii")
484  jj=assert_eq(SIZE(lat),            SIZE(vari,2),TRIM(modname)//" jj")
485  i1=assert_eq(SIZE(lon1),SIZE(lon2),SIZE(varo,1),TRIM(modname)//" i1")
486  j1=assert_eq(SIZE(lat1),           SIZE(varo,2),TRIM(modname)//" j1")
487  j2=SIZE(lat2)
488  ALLOCATE(vtmp(i1-1,j1))
489  IF(ibeg.AND.prt_level>1) THEN
490    WRITE(lunout,*)"---------------------------------------------------------"
491    WRITE(lunout,*)"$$$ Interpolation barycentrique pour "//TRIM(nam)//" $$$"
492    WRITE(lunout,*)"---------------------------------------------------------"
493  END IF
494  CALL inter_barxy(lon, lat(:jj-1), vari, lon2(:i1-1), lat2, vtmp)
495  CALL gr_int_dyn(vtmp, varo, i1-1, j1)
496
497END SUBROUTINE interp_startvar
498!
499!-------------------------------------------------------------------------------
500
501END MODULE etat0dyn
502!
503!*******************************************************************************
Note: See TracBrowser for help on using the repository browser.