source: LMDZ6/trunk/libf/dynphy_lonlat/phylmd/etat0dyn_netcdf.F90 @ 4467

Last change on this file since 4467 was 4050, checked in by dcugnet, 3 years ago

Second commit for new tracers.

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