source: ICOSA_LMDZ/src/phylmd/icolmdz_param_gravity_wave.f90 @ 5874

Last change on this file since 5874 was 5724, checked in by yann meurdesoif, 5 months ago

Due to dynamico module reorganization, some module imports were missing.
YM
  • Property svn:executable set to *
File size: 18.4 KB
Line 
1MODULE icolmdz_param_gravity_wave
2 
3   INTEGER, PARAMETER :: gw_legacy=1
4   INTEGER, PARAMETER :: gw_sso=2
5   INTEGER,SAVE       :: gw_method
6   !$OMP THREADPRIVATE(gw_method)
7
8CONTAINS
9 
10  SUBROUTINE init_param_gravity_wave
11  USE getin_mod  ! from dynamico
12  USE xios_mod   ! from dynamico
13  IMPLICIT NONE
14    CHARACTER(LEN=255) :: param_gw_method
15    LOGICAL :: create_etat0_limit
16
17    create_etat0_limit = .FALSE.
18    CALL getin('create_etat0_limit', create_etat0_limit)
19    IF (.NOT. create_etat0_limit) RETURN
20
21    param_gw_method='sso'
22    CALL getin('param_gw_method', param_gw_method)
23    SELECT CASE (TRIM(param_gw_method))
24      CASE ('legacy')
25        gw_method = gw_legacy
26      CASE ('sso')
27        gw_method = gw_sso
28      CASE default
29         PRINT *, 'Bad selector for param_gw_method : <', TRIM(param_gw_method), &
30            ' > options are <legacy>, <sso>'
31         STOP
32    END SELECT   
33
34    SELECT CASE (gw_method)
35      CASE (gw_legacy)
36      !$OMP MASTER
37        CALL xios_set_file_attr("relief_gw",enabled=.TRUE.)
38        CALL xios_set_fieldgroup_attr("gw_read_access",read_access=.TRUE.)
39      !$OMP END MASTER
40      CASE (gw_sso)
41      !$OMP MASTER
42        CALL xios_set_file_attr("orography",enabled=.TRUE.)
43        CALL xios_set_fieldgroup_attr("gwsso_read_access",read_access=.TRUE.)
44      !$OMP END MASTER
45    END SELECT 
46
47  END SUBROUTINE init_param_gravity_wave
48
49  SUBROUTINE param_gravity_wave
50  USE getin_mod
51  IMPLICIT NONE
52   LOGICAL :: create_etat0_limit
53
54    create_etat0_limit = .FALSE.
55    CALL getin('create_etat0_limit', create_etat0_limit)
56    IF (.NOT. create_etat0_limit) RETURN
57     
58    SELECT CASE (gw_method)
59      CASE (gw_legacy)
60        CALL param_gravity_wave_legacy
61      CASE (gw_sso)
62        CALL param_gravity_wave_sso
63    END SELECT
64
65  END SUBROUTINE param_gravity_wave
66
67  SUBROUTINE param_gravity_wave_sso
68    ! from dynamico
69  USE icosa
70  USE xios_mod
71  USE mpipara
72  USE earth_const
73  USE transfert_mod
74  ! from lmdz physic
75  USE create_etat0_unstruct_mod, ONLY : init_param_gw_phy => init_param_gw
76  USE dimphy, ONLY : klon
77  !from icosa_lmdz
78  USE distrib_icosa_lmdz_mod
79  IMPLICIT NONE
80    TYPE(t_field), POINTER, SAVE :: f_zmea(:), f_zpic(:), f_zval(:), f_zstd(:), f_zsig(:), f_zgam(:), f_zthe(:)
81     REAL, ALLOCATABLE :: zmea_phy(:), zpic_phy(:), zval_phy(:), zstd_phy(:), zsig_phy(:), zgam_phy(:), zthe_phy(:)
82
83    CALL allocate_field(f_zmea, field_t, type_real)
84    CALL allocate_field(f_zpic, field_t, type_real)
85    CALL allocate_field(f_zval, field_t, type_real)
86    CALL allocate_field(f_zstd, field_t, type_real)
87    CALL allocate_field(f_zsig, field_t, type_real)
88    CALL allocate_field(f_zgam, field_t, type_real)
89    CALL allocate_field(f_zthe, field_t, type_real)
90
91    CALL xios_read_field("zmea_sso",f_zmea)
92    CALL transfert_request(f_zmea, req_i0)
93    CALL transfert_request(f_zmea, req_i1)
94
95    CALL xios_read_field("zstd_sso",f_zstd)
96    CALL transfert_request(f_zstd, req_i0)
97    CALL transfert_request(f_zstd, req_i1)
98
99    CALL xios_read_field("zsig_sso",f_zsig)
100    CALL transfert_request(f_zsig, req_i0)
101    CALL transfert_request(f_zsig, req_i1)
102
103    CALL xios_read_field("zgam_sso",f_zgam)
104    CALL transfert_request(f_zgam, req_i0)
105    CALL transfert_request(f_zgam, req_i1)
106
107    CALL xios_read_field("zthe_sso",f_zthe)
108    CALL transfert_request(f_zthe, req_i0)
109    CALL transfert_request(f_zthe, req_i1)
110
111    ALLOCATE(zmea_phy(klon), zpic_phy(klon), zval_phy(klon), zstd_phy(klon), zsig_phy(klon), zgam_phy(klon), zthe_phy(klon))
112    CALL transfer_icosa_to_lmdz(f_zmea  , zmea_phy)
113    CALL transfer_icosa_to_lmdz(f_zstd  , zstd_phy)
114    CALL transfer_icosa_to_lmdz(f_zsig  , zsig_phy)
115    CALL transfer_icosa_to_lmdz(f_zgam  , zgam_phy)
116    CALL transfer_icosa_to_lmdz(f_zthe  , zthe_phy)
117    zpic_phy(:) = zmea_phy(:) + zstd_phy(:)
118    zval_phy(:) = MAX(zmea_phy(:) - zstd_phy(:), 0.)
119    CALL init_param_gw_phy(zmea_phy, zpic_phy, zval_phy, zstd_phy, zsig_phy, zgam_phy, zthe_phy)
120
121  END SUBROUTINE param_gravity_wave_sso
122
123  SUBROUTINE param_gravity_wave_legacy
124  ! from dynamico
125  USE icosa
126  USE xios_mod
127  USE mpipara
128  USE earth_const
129  USE transfert_mod
130  ! from lmdz physic
131  USE create_etat0_unstruct_mod, ONLY : init_param_gw_phy => init_param_gw
132  USE dimphy, ONLY : klon
133  !from icosa_lmdz
134  USE distrib_icosa_lmdz_mod
135  USE math_const, ONLY: Pi
136  IMPLICIT NONE
137    INTEGER :: ibegin, jbegin, ni, nj, ni_glo,nj_glo
138    REAL, ALLOCATABLE :: z(:,:)
139    REAL, ALLOCATABLE :: lon1d(:), lon2d_in(:,:), lon2d_out(:,:)
140    REAL, ALLOCATABLE :: lat1d(:), lat2d_in(:,:), lat2d_out(:,:)
141    REAL,ALLOCATABLE :: delta_x(:), delta_y(:)
142    REAL, ALLOCATABLE :: mask_in(:,:)
143    REAL, ALLOCATABLE :: zmea_in(:,:), zpic_in(:,:), zval_in(:,:), ztz_in(:,:), zytzy_in(:,:), zxtzx_in(:,:), zxtzy_in(:,:)
144
145    TYPE(t_field), POINTER, SAVE :: f_mask(:), f_zmea(:), f_zpic(:), f_zval(:), f_ztz(:), f_zytzy(:), f_zxtzx(:), f_zxtzy(:)
146    TYPE(t_field), POINTER, SAVE :: f_zstd(:), f_zsig(:), f_zgam(:), f_zthe(:)
147    REAL(rstd),POINTER :: mask(:), zmea(:), zpic(:), zval(:), ztz(:), zytzy(:), zxtzx(:), zxtzy(:), zstd(:), zsig(:), zgam(:), zthe(:)
148    REAL, ALLOCATABLE :: zmea_phy(:), zpic_phy(:), zval_phy(:), zstd_phy(:), zsig_phy(:), zgam_phy(:), zthe_phy(:)
149
150    INTEGER :: i, j, ind
151
152    CALL allocate_field(f_mask, field_t, type_real)
153    CALL allocate_field(f_zmea, field_t, type_real)
154    CALL allocate_field(f_zpic, field_t, type_real)
155    CALL allocate_field(f_zval, field_t, type_real)
156    CALL allocate_field(f_ztz, field_t, type_real)
157    CALL allocate_field(f_zytzy, field_t, type_real)
158    CALL allocate_field(f_zxtzx, field_t, type_real)
159    CALL allocate_field(f_zxtzy, field_t, type_real)
160    CALL allocate_field(f_zstd, field_t, type_real)
161    CALL allocate_field(f_zsig, field_t, type_real)
162    CALL allocate_field(f_zgam, field_t, type_real)
163    CALL allocate_field(f_zthe, field_t, type_real)
164
165!$OMP MASTER
166    CALL xios_get_domain_attr("domain_relief_gw",ibegin=ibegin, ni=ni, ni_glo=ni_glo, jbegin=jbegin, nj=nj,  nj_glo=nj_glo)
167    PRINT*, ni,nj,ni_glo,nj_glo
168   
169    ALLOCATE( lon1d(ni), lat1d(nj))
170    CALL xios_get_domain_attr("domain_relief_gw",lonvalue_1d=lon1d, latvalue_1d=lat1d)
171       
172    ALLOCATE(z(0:ni+1,0:nj+1))
173    CALL xios_recv_field("relief_exp", z)
174
175    ALLOCATE(lon2d_in(ni,nj), lat2d_in(ni,nj) )
176    ALLOCATE(lon2d_out(0:ni+1,0:nj+1), lat2d_out(0:ni+1,0:nj+1) )
177   
178    DO i=1,ni
179     lon2d_in(i,:) = lon1d(i)
180    END DO
181
182    DO j=1,nj
183     lat2d_in(:,j) = lat1d(j)
184    END DO
185
186    CALL xios_send_field("lon2d_in", lon2d_in)
187    CALL xios_recv_field("lon2d_out", lon2d_out)
188    CALL xios_send_field("lat2d_in", lat2d_in)
189    CALL xios_recv_field("lat2d_out", lat2d_out)
190   
191    ALLOCATE(delta_x(nj))
192    ALLOCATE(delta_y(nj))
193   
194    delta_y(:) = Pi*radius / nj_glo  !! difference with  lmdz reference : delta_y=2*Pi*radius/nj_glo ??
195   
196    DO j=1,nj
197      delta_x(j) = 2*Pi*radius / ni_glo * cos( lat2d_out(1,j) * pi / 180.)
198    ENDDO
199   
200    !north pole
201    IF (jbegin==0) THEN
202      z(:,0) = z(:,1)
203      delta_y(1) = delta_y(1)*0.5
204    ENDIF
205   
206    !south pole
207    IF (jbegin+nj==nj_glo) THEN
208      z(:,nj+1) = z(:,nj)
209      delta_y(nj) = delta_y(nj)*0.5
210    ENDIF
211   
212    ALLOCATE(mask_in(ni,nj))
213    ALLOCATE(zmea_in(ni,nj), zpic_in(ni,nj), zval_in(ni,nj), ztz_in(ni,nj), zytzy_in(ni,nj), zxtzx_in(ni,nj), zxtzy_in(ni,nj) )
214   
215    mask_in(:,:)=0
216    DO j=1,nj
217      DO i=1,ni
218        IF ( z(i,j)>1. ) mask_in(i,j) = 1.
219        zmea_in(i,j)  = z(i,j)
220        zpic_in(i,j)  = z(i,j)
221        zval_in(i,j)  = z(i,j)
222        ztz_in(i,j)   = z(i,j)*z(i,j)
223        zytzy_in(i,j) = (z(i,j+1)-z(i,j-1))**2/(2*delta_y(j))**2
224        zxtzx_in(i,j) = (z(i+1,j)-z(i-1,j))**2/(2*delta_x(j))**2
225        zxtzy_in(i,j) = (z(i,j+1)-z(i,j-1)) / (2*delta_y(j)) *(z(i+1,j)-z(i-1,j)) / (2*delta_x(j))
226      ENDDO
227    ENDDO
228
229    CALL xios_send_field("mask_in",mask_in)
230    CALL xios_send_field("zmea_in",zmea_in)
231    CALL xios_send_field("zpic_in",zpic_in)
232    CALL xios_send_field("zval_in",zval_in)
233    CALL xios_send_field("ztz_in",ztz_in)
234    CALL xios_send_field("zytzy_in",zytzy_in)
235    CALL xios_send_field("zxtzx_in",zxtzx_in)
236    CALL xios_send_field("zxtzy_in",zxtzy_in)
237   
238    !! only for test, to compare projection on ico grid Vs lmdz_reg grid
239    !! CALL compute_regular_param
240
241    !$OMP END MASTER
242   
243    CALL xios_read_field("mask_out", f_mask)
244    CALL transfert_request(f_mask, req_i0)
245    CALL transfert_request(f_mask, req_i1)
246
247    CALL xios_read_field("zmea_out",f_zmea)
248    CALL transfert_request(f_zmea, req_i0)
249    CALL transfert_request(f_zmea, req_i1)
250
251    CALL xios_read_field("zpic_out",f_zpic)
252    CALL transfert_request(f_zpic, req_i0)
253    CALL transfert_request(f_zpic, req_i1)
254
255    CALL xios_read_field("zval_out",f_zval)
256    CALL transfert_request(f_zval, req_i0)
257    CALL transfert_request(f_zval, req_i1)
258
259    CALL xios_read_field("ztz_out", f_ztz)
260    CALL transfert_request(f_ztz, req_i0)
261    CALL transfert_request(f_ztz, req_i1)
262
263    CALL xios_read_field("zytzy_out", f_zytzy)
264    CALL transfert_request(f_zytzy, req_i0)
265    CALL transfert_request(f_zytzy, req_i1)
266
267    CALL xios_read_field("zxtzx_out", f_zxtzx)
268    CALL transfert_request(f_zxtzx, req_i0)
269    CALL transfert_request(f_zxtzx, req_i1)
270
271    CALL xios_read_field("zxtzy_out", f_zxtzy)
272    CALL transfert_request(f_zxtzy, req_i0)
273    CALL transfert_request(f_zxtzy, req_i1)
274
275   !$OMP BARRIER
276   !$OMP MASTER
277
278   !! compute zstd
279
280    DO ind=1,ndomain
281      CALL swap_dimensions(ind)
282      CALL swap_geometry(ind)
283      zmea =f_zmea(ind)
284      zstd = f_zstd(ind)
285      ztz  = f_ztz(ind)
286      CALL compute_zstd(zmea, ztz, zstd)
287    ENDDO
288
289   !! smmothing ?
290
291    DO ind=1,ndomain
292      CALL swap_dimensions(ind)
293      CALL swap_geometry(ind)
294     
295      zmea=f_zmea(ind)
296      zpic=f_zpic(ind)
297      zval=f_zval(ind)
298      zstd=f_zstd(ind)
299      zytzy=f_zytzy(ind)
300      zxtzx=f_zxtzx(ind)
301      zxtzy=f_zxtzy(ind)
302 
303      CALL compute_smoothing(zmea)
304      CALL compute_smoothing(zpic)
305      CALL compute_smoothing(zval)
306      CALL compute_smoothing(zstd)
307      CALL compute_smoothing(zytzy)
308      CALL compute_smoothing(zxtzx)
309      CALL compute_smoothing(zxtzy)
310 
311    ENDDO
312
313    !! SLOPE, ANISOTROPY AND THETA ANGLE
314   
315    DO ind=1,ndomain
316      CALL swap_dimensions(ind)
317      CALL swap_geometry(ind)
318
319      zytzy=f_zytzy(ind)
320      zxtzx=f_zxtzx(ind)
321      zxtzy=f_zxtzy(ind)
322      zsig=f_zsig(ind)
323      zgam=f_zgam(ind)
324      zthe=f_zthe(ind)
325      CALL compute_sigma_gamma_theta(zxtzx, zytzy, zxtzy, zsig, zgam, zthe)
326   
327    ENDDO
328
329    !! apply mask
330   
331    DO ind=1,ndomain
332      CALL swap_dimensions(ind)
333      CALL swap_geometry(ind)
334
335      mask=f_mask(ind)
336      zmea=f_zmea(ind)
337      zpic=f_zpic(ind)
338      zval=f_zval(ind)
339      zstd=f_zstd(ind)
340      zsig=f_zsig(ind)
341      zgam=f_zgam(ind)
342      zthe=f_zthe(ind)
343
344      CALL compute_masking(mask, zmea, zpic, zval, zstd, zsig, zgam, zthe)
345
346    ENDDO
347
348!$OMP END MASTER
349!$OMP BARRIER
350 
351  !! for debugging
352  !!  CALL xios_write_field("zmask", f_mask)
353  !!  CALL xios_write_field("zmea", f_zmea)
354  !!  CALL xios_write_field("zpic", f_zpic)
355  !!  CALL xios_write_field("zval", f_zval)
356  !!  CALL xios_write_field("zsig", f_zsig)
357  !!  CALL xios_write_field("zgam", f_zgam)
358  !!  CALL xios_write_field("zthe", f_zthe)
359   
360    ALLOCATE(zmea_phy(klon), zpic_phy(klon), zval_phy(klon), zstd_phy(klon), zsig_phy(klon), zgam_phy(klon), zthe_phy(klon))
361    CALL transfer_icosa_to_lmdz(f_zmea  , zmea_phy)
362    CALL transfer_icosa_to_lmdz(f_zpic  , zpic_phy)
363    CALL transfer_icosa_to_lmdz(f_zval  , zval_phy)
364    CALL transfer_icosa_to_lmdz(f_zstd  , zstd_phy)
365    CALL transfer_icosa_to_lmdz(f_zsig  , zsig_phy)
366    CALL transfer_icosa_to_lmdz(f_zgam  , zgam_phy)
367    CALL transfer_icosa_to_lmdz(f_zthe  , zthe_phy)
368    CALL init_param_gw_phy(zmea_phy, zpic_phy, zval_phy, zstd_phy, zsig_phy, zgam_phy, zthe_phy)
369
370  END SUBROUTINE param_gravity_wave_legacy
371
372  SUBROUTINE compute_smoothing(X)
373  USE icosa
374  IMPLICIT NONE
375    REAL(rstd),INTENT(INOUT)  :: X(iim*jjm)
376    REAL(rstd)                :: tmp(iim*jjm)
377    INTEGER :: i,j,ij
378
379    DO j=jj_begin,jj_end
380      DO i=ii_begin,ii_end
381        ij=(j-1)*iim+i
382        tmp(ij) = (X(ij) + (X(ij+t_right) + X(ij+t_rup) + X(ij+t_lup) + X(ij+t_left) + X(ij+t_ldown) + X(ij+t_rdown))/3.)/3.
383      ENDDO
384    ENDDO
385 
386    DO j=jj_begin,jj_end
387      DO i=ii_begin,ii_end
388        ij=(j-1)*iim+i
389        X(ij) = tmp(ij)
390      ENDDO
391    ENDDO
392 
393  END SUBROUTINE
394
395  SUBROUTINE compute_zstd(zmea, ztz, zstd)
396  USE icosa
397  IMPLICIT NONE
398    REAL(rstd),INTENT(IN)   :: zmea(iim*jjm)
399    REAL(rstd),INTENT(IN)   :: ztz(iim*jjm)
400    REAL(rstd),INTENT(OUT) :: zstd(iim*jjm)
401
402    INTEGER :: i,j,ij
403
404    DO j=jj_begin-1,jj_end+1
405      DO i=ii_begin-1,ii_end+1
406        ij=(j-1)*iim+i
407        zstd(ij) = ztz(ij) - zmea(ij)*zmea(ij)
408        IF (zstd(ij)<=0) zstd(ij)=0
409        zstd(ij)=SQRT(zstd(ij))
410      ENDDO
411    ENDDO
412 
413  END SUBROUTINE compute_zstd
414
415
416  SUBROUTINE compute_sigma_gamma_theta(zxtzx, zytzy, zxtzy, zsig, zgam, zthe)
417  USE icosa
418  USE math_const, ONLY : Pi
419  IMPLICIT NONE
420    REAL(rstd),INTENT(IN)   :: zxtzx(iim*jjm)
421    REAL(rstd),INTENT(IN)   :: zytzy(iim*jjm)
422    REAL(rstd),INTENT(IN)   :: zxtzy(iim*jjm)
423    REAL(rstd),INTENT(OUT)  :: zsig(iim*jjm)
424    REAL(rstd),INTENT(OUT)  :: zgam(iim*jjm)
425    REAL(rstd),INTENT(OUT)  :: zthe(iim*jjm)
426   
427    REAL(rstd) :: xk,xl,xm,xp,xq,xw
428
429    INTEGER :: i,j,ij
430
431    DO j=jj_begin,jj_end
432      DO i=ii_begin,ii_end
433     
434        ij=(j-1)*iim+i
435        xk=(zxtzx(ij)+zytzy(ij))/2.
436        xl=(zxtzx(ij)-zytzy(ij))/2.
437        xm=zxtzy(ij)
438        xp=xk-SQRT(xl**2+xm**2)
439        xq=xk+SQRT(xl**2+xm**2)
440        xw=1.e-8
441        IF(xp<=xw) xp=0.
442        IF(xq<=xw) xq=xw
443        IF(ABS(xm)<=xw) xm=xw*SIGN(1.,xm)
444       
445        !--- SLOPE, ANISOTROPY AND THETA ANGLE
446     
447        zsig(ij)=SQRT(xq)
448        zgam(ij)=xp/xq
449        zthe(ij)=90.*ATAN2(xm,xl)/Pi
450
451      ENDDO
452    ENDDO
453
454  END SUBROUTINE compute_sigma_gamma_theta
455
456  SUBROUTINE compute_masking(mask, zmea, zpic, zval, zstd, zsig, zgam, zthe)
457  USE icosa
458  IMPLICIT NONE
459    REAL(rstd),INTENT(IN)   :: mask(iim*jjm)
460    REAL(rstd),INTENT(OUT)   :: zmea(iim*jjm)
461    REAL(rstd),INTENT(OUT)   :: zpic(iim*jjm)
462    REAL(rstd),INTENT(OUT)   :: zval(iim*jjm)
463    REAL(rstd),INTENT(OUT)   :: zstd(iim*jjm)
464    REAL(rstd),INTENT(OUT)   :: zsig(iim*jjm)
465    REAL(rstd),INTENT(OUT)   :: zgam(iim*jjm)
466    REAL(rstd),INTENT(OUT)   :: zthe(iim*jjm)
467   
468
469    INTEGER :: i,j,ij
470
471    DO j=jj_begin,jj_end
472      DO i=ii_begin,ii_end
473        ij=(j-1)*iim+i
474        IF (mask(ij)<0.1) THEN
475          zmea(ij)=0.
476          zpic(ij)=0.
477          zval(ij)=0.
478          zstd(ij)=0.
479          zsig(ij)=0.
480          zgam(ij)=0.
481          zthe(ij)=0.
482        ENDIF
483      ENDDO
484    ENDDO
485 
486  END SUBROUTINE compute_masking
487
488  SUBROUTINE compute_regular_param
489  USE xios_mod
490  USE icosa
491  USE math_const, ONLY : Pi
492  IMPLICIT NONE
493    REAL, ALLOCATABLE :: mask(:,:) , zmea(:,:), zpic(:,:), zval(:,:), ztz(:,:), zstd(:,:)
494    REAL, ALLOCATABLE :: zytzy(:,:) , zxtzx(:,:), zxtzy(:,:), zsig(:,:) , zgam(:,:), zthe(:,:)
495    REAL :: xk, xl, xm, xp, xq, xw
496    INTEGER :: ni, nj, ibegin, jbegin, ni_glo, nj_glo
497    INTEGER :: i,j
498
499    CALL xios_get_domain_attr("regular_gw",ibegin=ibegin, ni=ni, ni_glo=ni_glo, jbegin=jbegin, nj=nj,  nj_glo=nj_glo)
500
501    ALLOCATE(mask(0:ni+1,0:nj+1), zmea(0:ni+1,0:nj+1), zpic(0:ni+1,0:nj+1), zval(0:ni+1,0:nj+1), ztz(0:ni+1,0:nj+1))
502    ALLOCATE(zstd(0:ni+1,0:nj+1), zytzy(0:ni+1,0:nj+1), zxtzx(0:ni+1,0:nj+1), zxtzy(0:ni+1,0:nj+1))
503    ALLOCATE(zsig(0:ni+1,0:nj+1), zgam(0:ni+1,0:nj+1), zthe(0:ni+1,0:nj+1))
504
505    CALL xios_recv_field("mask_reg_exp",mask)
506    CALL xios_recv_field("zmea_reg_exp",zmea)
507    CALL xios_recv_field("zpic_reg_exp",zpic)
508    CALL xios_recv_field("zval_reg_exp",zval)
509    CALL xios_recv_field("ztz_reg_exp",ztz)
510    CALL xios_recv_field("zytzy_reg_exp",zytzy)
511    CALL xios_recv_field("zxtzx_reg_exp",zxtzx)
512    CALL xios_recv_field("zxtzy_reg_exp",zxtzy)
513   
514    DO j=0,nj+1
515      DO i=0, ni+1
516        zstd(i,j) = ztz(i,j)-zmea(i,j)*zmea(i,j)
517      ENDDO
518    ENDDO
519
520    ! smoothing
521    CALL smoothing_regular(ni,nj, zmea)
522    CALL smoothing_regular(ni,nj, zpic)
523    CALL smoothing_regular(ni,nj, zval)
524    CALL smoothing_regular(ni,nj, zstd)
525    CALL smoothing_regular(ni,nj, zytzy)
526    CALL smoothing_regular(ni,nj, zxtzx)
527    CALL smoothing_regular(ni,nj, zxtzy)
528
529    DO j=1,nj
530      DO i=1, ni
531        xk=(zxtzx(i,j)+zytzy(i,j))/2.
532        xl=(zxtzx(i,j)-zytzy(i,j))/2.
533        xm=zxtzy(i,j)
534        xp=xk-SQRT(xl**2+xm**2)
535        xq=xk+SQRT(xl**2+xm**2)
536        xw=1.e-8
537        IF(xp<=xw) xp=0.
538        IF(xq<=xw) xq=xw
539        IF(ABS(xm)<=xw) xm=xw*SIGN(1.,xm)
540       
541        !--- SLOPE, ANISOTROPY AND THETA ANGLE
542     
543        zsig(i,j)=SQRT(xq)
544        zgam(i,j)=xp/xq
545        zthe(i,j)=90.*ATAN2(xm,xl)/Pi
546      ENDDO
547    ENDDO
548
549    DO j=1,nj
550      DO i=1,ni
551        IF (mask(i,j)<0.1) THEN
552          zmea(i,j)=0.
553          zpic(i,j)=0.
554          zval(i,j)=0.
555          zstd(i,j)=0.
556          zsig(i,j)=0.
557          zgam(i,j)=0.
558          zthe(i,j)=0.
559        ENDIF
560      ENDDO
561    ENDDO
562
563    CALL xios_send_field("zmask_regout", mask(1:ni,1:nj))
564    CALL xios_send_field("zmea_regout", zmea(1:ni,1:nj))
565    CALL xios_send_field("zpic_regout", zpic(1:ni,1:nj))
566    CALL xios_send_field("zval_regout", zval(1:ni,1:nj))
567    CALL xios_send_field("zsig_regout", zsig(1:ni,1:nj))
568    CALL xios_send_field("zgam_regout", zgam(1:ni,1:nj))
569    CALL xios_send_field("zthe_regout", zthe(1:ni,1:nj))
570
571  END SUBROUTINE compute_regular_param
572
573  SUBROUTINE smoothing_regular(ni,nj,var)
574  IMPLICIT NONE
575    INTEGER, INTENT(IN) :: ni, nj
576    REAL, INTENT(INOUT) :: var(0:ni+1,0:nj+1)
577    REAL :: tmp(0:ni+1,0:nj+1)
578    INTEGER :: i,j
579
580    DO j=1,nj
581      DO i=1,ni
582        tmp(i,j) = (var(i,j) + 0.5*(var(i+1,j) + var(i-1,j) + var(i,j+1) +var(i,j-1))              &
583                             + 0.25*(var(i+1,j+1) + var(i-1,j+1) + var(i+1,j-1) + var(i-1,j-1)))/4.
584      ENDDO
585    ENDDO
586   
587    DO j=1,nj
588      DO i=1,ni
589        var(i,j)=tmp(i,j)
590      ENDDO
591    ENDDO
592
593  END SUBROUTINE smoothing_regular
594
595
596END MODULE icolmdz_param_gravity_wave
Note: See TracBrowser for help on using the repository browser.