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

Last change on this file since 2748 was 2738, checked in by oboucher, 8 years ago

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