source: LMDZ5/branches/IPSLCM5A2.1_ISO/libf/dynphy_lonlat/phyiso/etat0phys_netcdf.F90

Last change on this file was 3331, checked in by acozic, 7 years ago

Add modification for isotopes

  • Property svn:executable set to *
File size: 21.3 KB
Line 
1MODULE etat0phys
2!
3!*******************************************************************************
4! Purpose: Create physical 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!    *  etat0phys_netcdf routine can access to NetCDF data through subroutines:
12!         "start_init_phys" for variables contained in file "ECPHY.nc":
13!            'ST'     : Surface temperature
14!            'CDSW'   : Soil moisture
15!         "start_init_orog" for variables contained in file "Relief.nc":
16!            'RELIEF' : High resolution orography
17!
18!    * The land mask and corresponding weights can be:
19!      1) computed using the ocean mask from the ocean model (to ensure ocean
20!         fractions are the same for atmosphere and ocean) for coupled runs.
21!         File name: "o2a.nc"  ;  variable name: "OceMask"
22!      2) computed from topography file "Relief.nc" for forced runs.
23!
24!    * Allowed values for read_climoz flag are 0, 1 and 2:
25!      0: do not read an ozone climatology
26!      1: read a single ozone climatology that will be used day and night
27!      2: read two ozone climatologies, the average day and night climatology
28!         and the daylight climatology
29!-------------------------------------------------------------------------------
30!    * There is a big mess with the longitude size. Should it be iml or iml+1 ?
31!  I have chosen to use the iml+1 as an argument to this routine and we declare
32!  internaly smaller fields when needed. This needs to be cleared once and for
33!  all in LMDZ. A convention is required.
34!-------------------------------------------------------------------------------
35
36  USE ioipsl,             ONLY: flininfo, flinopen, flinget, flinclo
37  USE assert_eq_m,        ONLY: assert_eq
38  USE dimphy,             ONLY: klon
39  USE conf_dat_m,         ONLY: conf_dat2d
40  USE phys_state_var_mod, ONLY: zmea, zstd, zsig, zgam, zthe, zpic, zval, z0m, &
41          solsw, radsol, t_ancien, wake_deltat, wake_s,  rain_fall, qsol, z0h, &
42          sollw, rugoro, q_ancien, wake_deltaq, wake_pe, snow_fall, ratqs,w01, &
43    sig1, ftsol, clwcon, fm_therm, wake_Cstar,  pctsrf,  entr_therm,radpas, f0,&
44    zmax0,fevap, rnebcon,falb_dir, wake_fip,    agesno,  detr_therm, pbl_tke,  &
45#ifdef ISO
46    xt_ancien,wake_deltaxt,xtrain_fall,xtsnow_fall,fxtevap, &
47#endif
48    phys_state_var_init
49  USE comconst_mod, ONLY: pi, dtvr
50
51  PRIVATE
52  PUBLIC :: etat0phys_netcdf
53
54  include "iniprint.h"
55  include "dimensions.h"
56  include "paramet.h"
57  include "comgeom2.h"
58  include "dimsoil.h"
59  include "clesphys.h"
60  REAL, SAVE :: deg2rad
61  REAL, SAVE, ALLOCATABLE :: tsol(:)
62  INTEGER,            SAVE      :: iml_phys, jml_phys, llm_phys, ttm_phys, fid_phys
63  REAL, ALLOCATABLE,  SAVE      :: lon_phys(:,:), lat_phys(:,:), levphys_ini(:)
64  CHARACTER(LEN=256), PARAMETER :: orofname="Relief.nc", orogvar="RELIEF"
65  CHARACTER(LEN=256), PARAMETER :: phyfname="ECPHY.nc",  psrfvar="SP"
66  CHARACTER(LEN=256), PARAMETER :: qsolvar="CDSW",       tsrfvar="ST"
67
68
69CONTAINS
70
71
72!-------------------------------------------------------------------------------
73!
74SUBROUTINE etat0phys_netcdf(masque, phis)
75!
76!-------------------------------------------------------------------------------
77! Purpose: Creates initial states
78!-------------------------------------------------------------------------------
79! Notes:  1) This routine is designed to work for Earth
80!         2) If masque(:,:)/=-99999., masque and phis are already known.
81!         Otherwise: compute it.
82!-------------------------------------------------------------------------------
83  USE control_mod
84  USE fonte_neige_mod
85  USE pbl_surface_mod
86  USE regr_lat_time_climoz_m, ONLY: regr_lat_time_climoz
87  USE indice_sol_mod
88  USE conf_phys_m, ONLY: conf_phys
89  USE init_ssrf_m, ONLY: start_init_subsurf
90#ifdef ISO
91    use infotrac, ONLY: niso,iso_num,tnat
92#endif
93  !use ioipsl_getincom
94  IMPLICIT NONE
95!-------------------------------------------------------------------------------
96! Arguments:
97  REAL,    INTENT(INOUT) :: masque(:,:) !--- Land mask           dim(iip1,jjp1)
98  REAL,    INTENT(INOUT) :: phis  (:,:) !--- Ground geopotential dim(iip1,jjp1)
99!-------------------------------------------------------------------------------
100! Local variables:
101  CHARACTER(LEN=256) :: modname="etat0phys_netcdf", fmt
102  INTEGER            :: i, j, l, ji, iml, jml
103  LOGICAL            :: read_mask
104  REAL               :: phystep, dummy
105  REAL, DIMENSION(SIZE(masque,1),SIZE(masque,2)) :: masque_tmp,phiso
106  REAL, DIMENSION(klon)               :: sn, rugmer, run_off_lic_0, fder
107  REAL, DIMENSION(klon,nbsrf)         :: qsolsrf, snsrf
108  REAL, DIMENSION(klon,nsoilmx,nbsrf) :: tsoil
109#ifdef ISO
110  REAL, DIMENSION(niso,klon)               :: xtrun_off_lic_0
111  REAL, DIMENSION(niso,klon)               :: Rland_ice
112  REAL, DIMENSION(niso,klon,nbsrf)         :: xtsnsrf
113  INTEGER             :: ixt
114#endif
115
116!--- Arguments for conf_phys
117  LOGICAL :: ok_journe, ok_mensuel, ok_instan, ok_hf, ok_LES, callstats
118  REAL    :: solarlong0, seuil_inversion, fact_cldcon, facttemps
119  LOGICAL :: ok_newmicro
120  INTEGER :: iflag_radia, iflag_cldcon, iflag_ratqs
121  REAL    :: ratqsbas, ratqshaut, tau_ratqs
122  LOGICAL :: ok_ade, ok_aie, ok_cdnc, aerosol_couple
123  INTEGER :: flag_aerosol
124  INTEGER :: flag_aerosol_strat
125  LOGICAL :: new_aod
126  REAL    :: bl95_b0, bl95_b1
127  INTEGER :: read_climoz                        !--- Read ozone climatology
128  REAL    :: alp_offset
129  LOGICAL :: filtre_oro=.false.
130
131  deg2rad= pi/180.0
132  iml=assert_eq(SIZE(masque,1),SIZE(phis,1),TRIM(modname)//" iml")
133  jml=assert_eq(SIZE(masque,2),SIZE(phis,2),TRIM(modname)//" jml")
134
135! Physics configuration
136!*******************************************************************************
137  CALL conf_phys(  ok_journe, ok_mensuel, ok_instan, ok_hf, ok_LES,     &
138                   callstats,                                           &
139                   solarlong0,seuil_inversion,                          &
140                   fact_cldcon, facttemps,ok_newmicro,iflag_radia,      &
141                   iflag_cldcon,                                        &
142                   iflag_ratqs,ratqsbas,ratqshaut,tau_ratqs,            &
143                   ok_ade, ok_aie, ok_cdnc, aerosol_couple,             &
144                   flag_aerosol, flag_aerosol_strat, new_aod,           &
145                   bl95_b0, bl95_b1,                                    &
146                   read_climoz,                                         &
147                   alp_offset)
148  CALL phys_state_var_init(read_climoz)
149
150!--- Initial atmospheric CO2 conc. from .def file
151  co2_ppm0 = co2_ppm
152
153! Compute ground geopotential, sub-cells quantities and possibly the mask.
154!*******************************************************************************
155  read_mask=ANY(masque/=-99999.); masque_tmp=masque
156  CALL start_init_orog(rlonv, rlatu, phis, masque_tmp)
157
158  CALL getin('filtre_oro',filtre_oro)
159  IF (filtre_oro) CALL filtreoro(size(phis,1),size(phis,2),phis,masque_tmp,rlatu)
160
161  WRITE(fmt,"(i4,'i1)')")iml ; fmt='('//ADJUSTL(fmt)
162  IF(.NOT.read_mask) THEN                       !--- Keep mask form orography
163    masque=masque_tmp
164    IF(prt_level>=1) THEN
165      WRITE(lunout,*)'BUILT MASK :'
166      WRITE(lunout,fmt) NINT(masque)
167    END IF
168    WHERE(   masque(:,:)<EPSFRA) masque(:,:)=0.
169    WHERE(1.-masque(:,:)<EPSFRA) masque(:,:)=1.
170  END IF
171  CALL gr_dyn_fi(1,iml,jml,klon,masque,zmasq) !--- Land mask to physical grid
172
173! Compute tsol and qsol on physical grid, knowing phis on 2D grid.
174!*******************************************************************************
175  CALL start_init_phys(rlonu, rlatv, phis)
176
177! Some initializations.
178!*******************************************************************************
179  sn    (:) = 0.0                               !--- Snow
180  radsol(:) = 0.0                               !--- Net radiation at ground
181  rugmer(:) = 0.001                             !--- Ocean rugosity
182  IF(read_climoz>=1) &                          !--- Ozone climatology
183    CALL regr_lat_time_climoz(read_climoz)
184
185! Sub-surfaces initialization.
186!*******************************************************************************
187  CALL start_init_subsurf(read_mask)
188
189! Write physical initial state
190!*******************************************************************************
191  WRITE(lunout,*)'phystep ',dtvr,iphysiq,nbapp_rad
192  phystep = dtvr * FLOAT(iphysiq)
193  radpas  = NINT (86400./phystep/ FLOAT(nbapp_rad) )
194  WRITE(lunout,*)'phystep =', phystep, radpas
195
196! Init: ftsol, snsrf, qsolsrf, tsoil, rain_fall, snow_fall, solsw, sollw, z0
197!*******************************************************************************
198  DO i=1,nbsrf; ftsol(:,i) = tsol; END DO
199#ifdef ISO
200  DO i=1,nbsrf
201        snsrf(:,i) = sn
202        do ixt=1,niso
203           xtsnsrf(ixt,:,i) = snsrf(:,i)*tnat(iso_num(ixt))
204        enddo
205  END DO
206#else
207  DO i=1,nbsrf; snsrf(:,i) = sn;   END DO
208#endif
209  falb_dir(:, :, is_ter) = 0.08
210  falb_dir(:, :, is_lic) = 0.6
211  falb_dir(:, :, is_oce) = 0.5
212  falb_dir(:, :, is_sic) = 0.6
213  fevap(:,:) = 0.
214  DO i=1,nbsrf; qsolsrf(:,i)=150.; END DO
215  DO i=1,nbsrf; DO j=1,nsoilmx; tsoil(:,j,i) = tsol; END DO; END DO
216  rain_fall  = 0.
217  snow_fall  = 0.
218  solsw      = 165.
219  sollw      = -53.
220  t_ancien   = 273.15
221  q_ancien   = 0.
222  agesno     = 0.
223
224  z0m(:,is_oce) = rugmer(:)
225  z0m(:,is_ter) = MAX(1.0e-05,zstd(:)*zsig(:)/2.0)
226  z0m(:,is_lic) = MAX(1.0e-05,zstd(:)*zsig(:)/2.0)
227  z0m(:,is_sic) = 0.001
228  z0h(:,:)=z0m(:,:)
229
230  fder    = 0.0
231  clwcon  = 0.0
232  rnebcon = 0.0
233  ratqs   = 0.0
234  run_off_lic_0 = 0.0
235  rugoro  = 0.0
236#ifdef ISO
237  fxtevap(:,:,:) = 0.
238  xtrain_fall  = 0.
239  xtsnow_fall  = 0.
240  xtrun_off_lic_0 = 0.0
241  xt_ancien   = 0.
242#endif
243
244! Before phyredem calling, surface modules and values to be saved in startphy.nc
245! are initialized
246!*******************************************************************************
247  dummy            = 1.0
248  pbl_tke(:,:,:)   = 1.e-8
249  zmax0(:)         = 40.
250  f0(:)            = 1.e-5
251  sig1(:,:)        = 0.
252  w01(:,:)         = 0.
253  wake_deltat(:,:) = 0.
254  wake_deltaq(:,:) = 0.
255  wake_s(:)        = 0.
256  wake_cstar(:)    = 0.
257  wake_fip(:)      = 0.
258  wake_pe          = 0.
259  fm_therm         = 0.
260  entr_therm       = 0.
261  detr_therm       = 0.
262#ifdef ISO
263  wake_deltaxt(:,:,:) = 0.
264#endif
265
266
267  CALL fonte_neige_init(run_off_lic_0 &
268#ifdef ISO
269    &  ,xtrun_off_lic_0 &
270#endif     
271    &   )
272  CALL pbl_surface_init( fder, snsrf, qsolsrf, tsoil  &
273#ifdef ISO       
274     &   ,xtsnsrf,Rland_ice &
275#endif
276     &  )
277  CALL phyredem( "startphy.nc" )
278
279!  WRITE(lunout,*)'CCCCCCCCCCCCCCCCCC REACTIVER SORTIE VISU DANS ETAT0'
280!  WRITE(lunout,*)'entree histclo'
281  CALL histclo()
282
283END SUBROUTINE etat0phys_netcdf
284!
285!-------------------------------------------------------------------------------
286
287
288!-------------------------------------------------------------------------------
289!
290SUBROUTINE start_init_orog(lon_in,lat_in,phis,masque)
291!
292!===============================================================================
293! Comment:
294!   This routine launch grid_noro, which computes parameters for SSO scheme as
295!   described in LOTT & MILLER (1997) and LOTT(1999).
296!===============================================================================
297  USE grid_noro_m, ONLY: grid_noro
298  IMPLICIT NONE
299!-------------------------------------------------------------------------------
300! Arguments:
301  REAL,    INTENT(IN)    :: lon_in(:), lat_in(:)   ! dim (iml) (jml)
302  REAL,    INTENT(INOUT) :: phis(:,:), masque(:,:) ! dim (iml,jml)
303!-------------------------------------------------------------------------------
304! Local variables:
305  CHARACTER(LEN=256) :: modname
306  INTEGER            :: fid, llm_tmp,ttm_tmp, iml,jml, iml_rel,jml_rel, itau(1)
307  REAL               :: lev(1), date, dt
308  REAL, ALLOCATABLE  :: lon_rad(:), lon_ini(:), lon_rel(:,:), relief_hi(:,:)
309  REAL, ALLOCATABLE  :: lat_rad(:), lat_ini(:), lat_rel(:,:), tmp_var  (:,:)
310  REAL, ALLOCATABLE  :: zmea0(:,:), zstd0(:,:), zsig0(:,:)
311  REAL, ALLOCATABLE  :: zgam0(:,:), zthe0(:,:), zpic0(:,:), zval0(:,:)
312!-------------------------------------------------------------------------------
313  modname="start_init_orog"
314  iml=assert_eq(SIZE(lon_in),SIZE(phis,1),SIZE(masque,1),TRIM(modname)//" iml")
315  jml=assert_eq(SIZE(lat_in),SIZE(phis,2),SIZE(masque,2),TRIM(modname)//" jml")
316
317!--- HIGH RESOLUTION OROGRAPHY
318  CALL flininfo(orofname, iml_rel, jml_rel, llm_tmp, ttm_tmp, fid)
319
320  ALLOCATE(lat_rel(iml_rel,jml_rel),lon_rel(iml_rel,jml_rel))
321  CALL flinopen(orofname, .FALSE., iml_rel, jml_rel, llm_tmp, lon_rel, lat_rel,&
322                lev, ttm_tmp, itau, date, dt, fid)
323  ALLOCATE(relief_hi(iml_rel,jml_rel))
324  CALL flinget(fid, orogvar, iml_rel, jml_rel, llm_tmp, ttm_tmp, 1,1, relief_hi)
325  CALL flinclo(fid)
326
327!--- IF ANGLES ARE IN DEGREES, THEY ARE CONVERTED INTO RADIANS
328  ALLOCATE(lon_ini(iml_rel),lat_ini(jml_rel))
329  lon_ini(:)=lon_rel(:,1); IF(MAXVAL(lon_rel)>pi) lon_ini=lon_ini*deg2rad
330  lat_ini(:)=lat_rel(1,:); IF(MAXVAL(lat_rel)>pi) lat_ini=lat_ini*deg2rad
331
332!--- FIELDS ARE PROCESSED TO BE ON STANDARD ANGULAR DOMAINS
333  ALLOCATE(lon_rad(iml_rel),lat_rad(jml_rel))
334  CALL conf_dat2d(orogvar, lon_ini, lat_ini, lon_rad, lat_rad, relief_hi,.FALSE.)
335  DEALLOCATE(lon_ini,lat_ini)
336
337!--- COMPUTING THE REQUIRED FIELDS USING ROUTINE grid_noro
338  WRITE(lunout,*)
339  WRITE(lunout,*)'*** Compute parameters needed for gravity wave drag code ***'
340
341!--- ALLOCATIONS OF SUB-CELL SCALES QUANTITIES
342  ALLOCATE(zmea0(iml,jml),zstd0(iml,jml)) !--- Mean orography and std deviation
343  ALLOCATE(zsig0(iml,jml),zgam0(iml,jml)) !--- Slope and nisotropy
344  ALLOCATE(zthe0(iml,jml))                !--- Highest slope orientation
345  ALLOCATE(zpic0(iml,jml),zval0(iml,jml)) !--- Peaks and valley heights
346
347!--- CALL OROGRAPHY MODULE TO COMPUTE FIELDS
348  CALL grid_noro(lon_rad,lat_rad,relief_hi,lon_in,lat_in,phis,zmea0,zstd0,     &
349                                      zsig0,zgam0,zthe0,zpic0,zval0,masque)
350  phis = phis * 9.81
351  phis(iml,:) = phis(1,:)
352  DEALLOCATE(relief_hi,lon_rad,lat_rad)
353
354!--- PUT QUANTITIES TO PHYSICAL GRID
355  CALL gr_dyn_fi(1,iml,jml,klon,zmea0,zmea); DEALLOCATE(zmea0)
356  CALL gr_dyn_fi(1,iml,jml,klon,zstd0,zstd); DEALLOCATE(zstd0)
357  CALL gr_dyn_fi(1,iml,jml,klon,zsig0,zsig); DEALLOCATE(zsig0)
358  CALL gr_dyn_fi(1,iml,jml,klon,zgam0,zgam); DEALLOCATE(zgam0)
359  CALL gr_dyn_fi(1,iml,jml,klon,zthe0,zthe); DEALLOCATE(zthe0)
360  CALL gr_dyn_fi(1,iml,jml,klon,zpic0,zpic); DEALLOCATE(zpic0)
361  CALL gr_dyn_fi(1,iml,jml,klon,zval0,zval); DEALLOCATE(zval0)
362
363
364END SUBROUTINE start_init_orog
365!
366!-------------------------------------------------------------------------------
367
368
369!-------------------------------------------------------------------------------
370!
371SUBROUTINE start_init_phys(lon_in,lat_in,phis)
372!
373!===============================================================================
374! Purpose:   Compute tsol and qsol, knowing phis.
375!===============================================================================
376  IMPLICIT NONE
377!-------------------------------------------------------------------------------
378! Arguments:
379  REAL,    INTENT(IN) :: lon_in(:),  lat_in(:)       ! dim (iml) (jml2)
380  REAL,    INTENT(IN) :: phis(:,:)                   ! dim (iml,jml)
381!-------------------------------------------------------------------------------
382! Local variables:
383  CHARACTER(LEN=256) :: modname
384  REAL               :: date, dt
385  INTEGER            :: iml, jml, jml2, itau(1)
386  REAL, ALLOCATABLE  :: lon_rad(:), lon_ini(:), var_ana(:,:)
387  REAL, ALLOCATABLE  :: lat_rad(:), lat_ini(:)
388  REAL, ALLOCATABLE  :: ts(:,:), qs(:,:)
389!-------------------------------------------------------------------------------
390  modname="start_init_phys"
391  iml=assert_eq(SIZE(lon_in),SIZE(phis,1),TRIM(modname)//" iml")
392  jml=SIZE(phis,2); jml2=SIZE(lat_in)
393
394  WRITE(lunout,*)'Opening the surface analysis'
395  CALL flininfo(phyfname, iml_phys, jml_phys, llm_phys, ttm_phys, fid_phys)
396  WRITE(lunout,*) 'Values read: ',  iml_phys, jml_phys, llm_phys, ttm_phys
397
398  ALLOCATE(lat_phys(iml_phys,jml_phys),lon_phys(iml_phys,jml_phys))
399  ALLOCATE(levphys_ini(llm_phys))
400  CALL flinopen(phyfname, .FALSE., iml_phys, jml_phys, llm_phys,              &
401                lon_phys,lat_phys,levphys_ini,ttm_phys,itau,date,dt,fid_phys)
402
403!--- IF ANGLES ARE IN DEGREES, THEY ARE CONVERTED INTO RADIANS
404  ALLOCATE(lon_ini(iml_phys),lat_ini(jml_phys))
405  lon_ini(:)=lon_phys(:,1); IF(MAXVAL(lon_phys)>pi) lon_ini=lon_ini*deg2rad
406  lat_ini(:)=lat_phys(1,:); IF(MAXVAL(lat_phys)>pi) lat_ini=lat_ini*deg2rad
407
408  ALLOCATE(var_ana(iml_phys,jml_phys),lon_rad(iml_phys),lat_rad(jml_phys))
409  CALL get_var_phys(tsrfvar,ts)                   !--- SURFACE TEMPERATURE
410  CALL get_var_phys(qsolvar,qs)                   !--- SOIL MOISTURE
411  CALL flinclo(fid_phys)
412  DEALLOCATE(var_ana,lon_rad,lat_rad,lon_ini,lat_ini)
413
414!--- TSOL AND QSOL ON PHYSICAL GRID
415  ALLOCATE(tsol(klon))
416  CALL gr_dyn_fi(1,iml,jml,klon,ts,tsol)
417  CALL gr_dyn_fi(1,iml,jml,klon,qs,qsol)
418  DEALLOCATE(ts,qs)
419
420CONTAINS
421
422!-------------------------------------------------------------------------------
423!
424SUBROUTINE get_var_phys(title,field)
425!
426!-------------------------------------------------------------------------------
427  IMPLICIT NONE
428!-------------------------------------------------------------------------------
429! Arguments:
430  CHARACTER(LEN=*),  INTENT(IN)    :: title
431  REAL, ALLOCATABLE, INTENT(INOUT) :: field(:,:)
432!-------------------------------------------------------------------------------
433! Local variables:
434  INTEGER :: tllm
435!-------------------------------------------------------------------------------
436  SELECT CASE(title)
437    CASE(psrfvar);         tllm=0
438    CASE(tsrfvar,qsolvar); tllm=llm_phys
439  END SELECT
440  IF(ALLOCATED(field)) RETURN
441  ALLOCATE(field(iml,jml)); field(:,:)=0.
442  CALL flinget(fid_phys,title,iml_phys,jml_phys,tllm,ttm_phys,1,1,var_ana)
443  CALL conf_dat2d(title, lon_ini, lat_ini, lon_rad, lat_rad, var_ana, .TRUE.)
444  CALL interp_startvar(title, .TRUE., lon_rad, lat_rad, var_ana,               &
445                                      lon_in,  lat_in,  field)
446
447END SUBROUTINE get_var_phys
448!
449!-------------------------------------------------------------------------------
450!
451END SUBROUTINE start_init_phys
452!
453!-------------------------------------------------------------------------------
454
455
456!-------------------------------------------------------------------------------
457!
458SUBROUTINE interp_startvar(nam,ibeg,lon,lat,vari,lon2,lat2,varo)
459!
460!-------------------------------------------------------------------------------
461  USE inter_barxy_m, ONLY: inter_barxy
462  IMPLICIT NONE
463!-------------------------------------------------------------------------------
464! Arguments:
465  CHARACTER(LEN=*), INTENT(IN)  :: nam
466  LOGICAL,          INTENT(IN)  :: ibeg
467  REAL,             INTENT(IN)  :: lon(:), lat(:)   ! dim (ii) (jj)
468  REAL,             INTENT(IN)  :: vari(:,:)        ! dim (ii,jj)
469  REAL,             INTENT(IN)  :: lon2(:), lat2(:) ! dim (i1) (j2)
470  REAL,             INTENT(OUT) :: varo(:,:)        ! dim (i1) (j1)
471!-------------------------------------------------------------------------------
472! Local variables:
473  CHARACTER(LEN=256) :: modname
474  INTEGER            :: ii, jj, i1, j1, j2
475  REAL, ALLOCATABLE  :: vtmp(:,:)
476!-------------------------------------------------------------------------------
477  modname="interp_startvar"
478  ii=assert_eq(SIZE(lon), SIZE(vari,1),TRIM(modname)//" ii")
479  jj=assert_eq(SIZE(lat), SIZE(vari,2),TRIM(modname)//" jj")
480  i1=assert_eq(SIZE(lon2),SIZE(varo,1),TRIM(modname)//" i1")
481  j1=SIZE(varo,2); j2=SIZE(lat2)
482  ALLOCATE(vtmp(i1-1,j1))
483  IF(ibeg.AND.prt_level>1) THEN
484    WRITE(lunout,*)"--------------------------------------------------------"
485    WRITE(lunout,*)"$$$ Interpolation barycentrique pour "//TRIM(nam)//" $$$"
486    WRITE(lunout,*)"--------------------------------------------------------"
487  END IF
488  CALL inter_barxy(lon, lat(:jj-1), vari, lon2(:i1-1), lat2, vtmp)
489  CALL gr_int_dyn(vtmp, varo, i1-1, j1)
490
491END SUBROUTINE interp_startvar
492!
493!-------------------------------------------------------------------------------
494!
495!*******************************************************************************
496
497SUBROUTINE filtreoro(imp1,jmp1,phis,masque,rlatu)
498
499IMPLICIT NONE
500
501  INTEGER imp1,jmp1
502  REAL, DIMENSION(imp1,jmp1) :: phis,masque
503  REAL, DIMENSION(jmp1) :: rlatu
504  REAL, DIMENSION(imp1) :: wwf
505  REAL, DIMENSION(imp1,jmp1) :: phiso
506  INTEGER :: ifiltre,ifi,ii,i,j
507  REAL :: coslat0,ssz
508
509  coslat0=0.5
510  phiso=phis
511  do j=2,jmp1-1
512     print*,'avant if ',cos(rlatu(j)),coslat0
513     if (cos(rlatu(j))<coslat0) then
514         ! nb de pts affectes par le filtrage de part et d'autre du pt
515         ifiltre=(coslat0/cos(rlatu(j))-1.)/2.
516         wwf=0.
517         do i=1,ifiltre
518            wwf(i)=1.
519         enddo
520         wwf(ifiltre+1)=(coslat0/cos(rlatu(j))-1.)/2.-ifiltre
521         do i=1,imp1-1
522            if (masque(i,j)>0.9) then
523               ssz=phis(i,j)
524               do ifi=1,ifiltre+1
525                  ii=i+ifi
526                  if (ii>imp1-1) ii=ii-imp1+1
527                  ssz=ssz+wwf(ifi)*phis(ii,j)
528                  ii=i-ifi
529                  if (ii<1) ii=ii+imp1-1
530                  ssz=ssz+wwf(ifi)*phis(ii,j)
531               enddo
532               phis(i,j)=ssz*cos(rlatu(j))/coslat0
533            endif
534         enddo
535         print*,'j=',j,coslat0/cos(rlatu(j)), (1.+2.*sum(wwf))*cos(rlatu(j))/coslat0
536     endif
537  enddo
538  call dump2d(imp1,jmp1,phis,'phis ')
539  call dump2d(imp1,jmp1,masque,'masque ')
540  call dump2d(imp1,jmp1,phis-phiso,'dphis ')
541
542END SUBROUTINE filtreoro
543
544
545END MODULE etat0phys
Note: See TracBrowser for help on using the repository browser.