source: LMDZ5/trunk/libf/dynphy_lonlat/phylmd/etat0phys_netcdf.F90 @ 2457

Last change on this file since 2457 was 2453, checked in by fhourdin, 10 years ago

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