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

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

Fix halo problem when computing parameters of legacy gravity wave parametrization.
YM
  • Property svn:executable set to *
File size: 18.3 KB
RevLine 
[4855]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
12  USE xios_mod
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
[4858]21    param_gw_method='sso'
[4855]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  IMPLICIT NONE
136    INTEGER :: ibegin, jbegin, ni, nj, ni_glo,nj_glo
137    REAL, ALLOCATABLE :: z(:,:)
138    REAL, ALLOCATABLE :: lon1d(:), lon2d_in(:,:), lon2d_out(:,:)
139    REAL, ALLOCATABLE :: lat1d(:), lat2d_in(:,:), lat2d_out(:,:)
140    REAL,ALLOCATABLE :: delta_x(:), delta_y(:)
141    REAL, ALLOCATABLE :: mask_in(:,:)
142    REAL, ALLOCATABLE :: zmea_in(:,:), zpic_in(:,:), zval_in(:,:), ztz_in(:,:), zytzy_in(:,:), zxtzx_in(:,:), zxtzy_in(:,:)
143
144    TYPE(t_field), POINTER, SAVE :: f_mask(:), f_zmea(:), f_zpic(:), f_zval(:), f_ztz(:), f_zytzy(:), f_zxtzx(:), f_zxtzy(:)
145    TYPE(t_field), POINTER, SAVE :: f_zstd(:), f_zsig(:), f_zgam(:), f_zthe(:)
146    REAL(rstd),POINTER :: mask(:), zmea(:), zpic(:), zval(:), ztz(:), zytzy(:), zxtzx(:), zxtzy(:), zstd(:), zsig(:), zgam(:), zthe(:)
147    REAL, ALLOCATABLE :: zmea_phy(:), zpic_phy(:), zval_phy(:), zstd_phy(:), zsig_phy(:), zgam_phy(:), zthe_phy(:)
148
149    INTEGER :: i, j, ind
150
151    CALL allocate_field(f_mask, field_t, type_real)
152    CALL allocate_field(f_zmea, field_t, type_real)
153    CALL allocate_field(f_zpic, field_t, type_real)
154    CALL allocate_field(f_zval, field_t, type_real)
155    CALL allocate_field(f_ztz, field_t, type_real)
156    CALL allocate_field(f_zytzy, field_t, type_real)
157    CALL allocate_field(f_zxtzx, field_t, type_real)
158    CALL allocate_field(f_zxtzy, field_t, type_real)
159    CALL allocate_field(f_zstd, field_t, type_real)
160    CALL allocate_field(f_zsig, field_t, type_real)
161    CALL allocate_field(f_zgam, field_t, type_real)
162    CALL allocate_field(f_zthe, field_t, type_real)
163
164!$OMP MASTER
165    CALL xios_get_domain_attr("domain_relief_gw",ibegin=ibegin, ni=ni, ni_glo=ni_glo, jbegin=jbegin, nj=nj,  nj_glo=nj_glo)
166    PRINT*, ni,nj,ni_glo,nj_glo
167   
168    ALLOCATE( lon1d(ni), lat1d(nj))
169    CALL xios_get_domain_attr("domain_relief_gw",lonvalue_1d=lon1d, latvalue_1d=lat1d)
170       
171    ALLOCATE(z(0:ni+1,0:nj+1))
172    CALL xios_recv_field("relief_exp", z)
173
174    ALLOCATE(lon2d_in(ni,nj), lat2d_in(ni,nj) )
175    ALLOCATE(lon2d_out(0:ni+1,0:nj+1), lat2d_out(0:ni+1,0:nj+1) )
176   
177    DO i=1,ni
178     lon2d_in(i,:) = lon1d(i)
179    END DO
180
181    DO j=1,nj
182     lat2d_in(:,j) = lat1d(j)
183    END DO
184
185    CALL xios_send_field("lon2d_in", lon2d_in)
186    CALL xios_recv_field("lon2d_out", lon2d_out)
187    CALL xios_send_field("lat2d_in", lat2d_in)
188    CALL xios_recv_field("lat2d_out", lat2d_out)
189   
190    ALLOCATE(delta_x(nj))
191    ALLOCATE(delta_y(nj))
192   
193    delta_y(:) = Pi*radius / nj_glo  !! difference with  lmdz reference : delta_y=2*Pi*radius/nj_glo ??
194   
195    DO j=1,nj
196      delta_x(j) = 2*Pi*radius / ni_glo * cos( lat2d_out(1,j) * pi / 180.)
197    ENDDO
198   
199    !north pole
200    IF (jbegin==0) THEN
201      z(:,0) = z(:,1)
202      delta_y(1) = delta_y(1)*0.5
203    ENDIF
204   
205    !south pole
206    IF (jbegin+nj==nj_glo) THEN
207      z(:,nj+1) = z(:,nj)
208      delta_y(nj) = delta_y(nj)*0.5
209    ENDIF
210   
211    ALLOCATE(mask_in(ni,nj))
212    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) )
213   
214    mask_in(:,:)=0
215    DO j=1,nj
216      DO i=1,ni
217        IF ( z(i,j)>1. ) mask_in(i,j) = 1.
218        zmea_in(i,j)  = z(i,j)
219        zpic_in(i,j)  = z(i,j)
220        zval_in(i,j)  = z(i,j)
221        ztz_in(i,j)   = z(i,j)*z(i,j)
222        zytzy_in(i,j) = (z(i,j+1)-z(i,j-1))**2/(2*delta_y(j))**2
223        zxtzx_in(i,j) = (z(i+1,j)-z(i-1,j))**2/(2*delta_x(j))**2
224        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))
225      ENDDO
226    ENDDO
227
228    CALL xios_send_field("mask_in",mask_in)
229    CALL xios_send_field("zmea_in",zmea_in)
230    CALL xios_send_field("zpic_in",zpic_in)
231    CALL xios_send_field("zval_in",zval_in)
232    CALL xios_send_field("ztz_in",ztz_in)
233    CALL xios_send_field("zytzy_in",zytzy_in)
234    CALL xios_send_field("zxtzx_in",zxtzx_in)
235    CALL xios_send_field("zxtzy_in",zxtzy_in)
236   
237    !! only for test, to compare projection on ico grid Vs lmdz_reg grid
238    !! CALL compute_regular_param
239
240    !$OMP END MASTER
241   
242    CALL xios_read_field("mask_out", f_mask)
243    CALL transfert_request(f_mask, req_i0)
244    CALL transfert_request(f_mask, req_i1)
245
246    CALL xios_read_field("zmea_out",f_zmea)
247    CALL transfert_request(f_zmea, req_i0)
248    CALL transfert_request(f_zmea, req_i1)
249
250    CALL xios_read_field("zpic_out",f_zpic)
251    CALL transfert_request(f_zpic, req_i0)
252    CALL transfert_request(f_zpic, req_i1)
253
254    CALL xios_read_field("zval_out",f_zval)
255    CALL transfert_request(f_zval, req_i0)
256    CALL transfert_request(f_zval, req_i1)
257
258    CALL xios_read_field("ztz_out", f_ztz)
259    CALL transfert_request(f_ztz, req_i0)
260    CALL transfert_request(f_ztz, req_i1)
261
262    CALL xios_read_field("zytzy_out", f_zytzy)
263    CALL transfert_request(f_zytzy, req_i0)
264    CALL transfert_request(f_zytzy, req_i1)
265
266    CALL xios_read_field("zxtzx_out", f_zxtzx)
267    CALL transfert_request(f_zxtzx, req_i0)
268    CALL transfert_request(f_zxtzx, req_i1)
269
270    CALL xios_read_field("zxtzy_out", f_zxtzy)
271    CALL transfert_request(f_zxtzy, req_i0)
272    CALL transfert_request(f_zxtzy, req_i1)
273
274   !$OMP BARRIER
275   !$OMP MASTER
276
277   !! compute zstd
278
279    DO ind=1,ndomain
280      CALL swap_dimensions(ind)
281      CALL swap_geometry(ind)
282      zmea =f_zmea(ind)
283      zstd = f_zstd(ind)
284      ztz  = f_ztz(ind)
285      CALL compute_zstd(zmea, ztz, zstd)
286    ENDDO
287
288   !! smmothing ?
289
290    DO ind=1,ndomain
291      CALL swap_dimensions(ind)
292      CALL swap_geometry(ind)
293     
294      zmea=f_zmea(ind)
295      zpic=f_zpic(ind)
296      zval=f_zval(ind)
297      zstd=f_zstd(ind)
298      zytzy=f_zytzy(ind)
299      zxtzx=f_zxtzx(ind)
300      zxtzy=f_zxtzy(ind)
301 
302      CALL compute_smoothing(zmea)
303      CALL compute_smoothing(zpic)
304      CALL compute_smoothing(zval)
305      CALL compute_smoothing(zstd)
306      CALL compute_smoothing(zytzy)
307      CALL compute_smoothing(zxtzx)
308      CALL compute_smoothing(zxtzy)
309 
310    ENDDO
311
312    !! SLOPE, ANISOTROPY AND THETA ANGLE
313   
314    DO ind=1,ndomain
315      CALL swap_dimensions(ind)
316      CALL swap_geometry(ind)
317
318      zytzy=f_zytzy(ind)
319      zxtzx=f_zxtzx(ind)
320      zxtzy=f_zxtzy(ind)
321      zsig=f_zsig(ind)
322      zgam=f_zgam(ind)
323      zthe=f_zthe(ind)
324      CALL compute_sigma_gamma_theta(zxtzx, zytzy, zxtzy, zsig, zgam, zthe)
325   
326    ENDDO
327
328    !! apply mask
329   
330    DO ind=1,ndomain
331      CALL swap_dimensions(ind)
332      CALL swap_geometry(ind)
333
334      mask=f_mask(ind)
335      zmea=f_zmea(ind)
336      zpic=f_zpic(ind)
337      zval=f_zval(ind)
338      zstd=f_zstd(ind)
339      zsig=f_zsig(ind)
340      zgam=f_zgam(ind)
341      zthe=f_zthe(ind)
342
343      CALL compute_masking(mask, zmea, zpic, zval, zstd, zsig, zgam, zthe)
344
345    ENDDO
346
347!$OMP END MASTER
348!$OMP BARRIER
349 
350  !! for debugging
351  !!  CALL xios_write_field("zmask", f_mask)
352  !!  CALL xios_write_field("zmea", f_zmea)
353  !!  CALL xios_write_field("zpic", f_zpic)
354  !!  CALL xios_write_field("zval", f_zval)
355  !!  CALL xios_write_field("zsig", f_zsig)
356  !!  CALL xios_write_field("zgam", f_zgam)
357  !!  CALL xios_write_field("zthe", f_zthe)
358   
359    ALLOCATE(zmea_phy(klon), zpic_phy(klon), zval_phy(klon), zstd_phy(klon), zsig_phy(klon), zgam_phy(klon), zthe_phy(klon))
360    CALL transfer_icosa_to_lmdz(f_zmea  , zmea_phy)
361    CALL transfer_icosa_to_lmdz(f_zpic  , zpic_phy)
362    CALL transfer_icosa_to_lmdz(f_zval  , zval_phy)
363    CALL transfer_icosa_to_lmdz(f_zstd  , zstd_phy)
364    CALL transfer_icosa_to_lmdz(f_zsig  , zsig_phy)
365    CALL transfer_icosa_to_lmdz(f_zgam  , zgam_phy)
366    CALL transfer_icosa_to_lmdz(f_zthe  , zthe_phy)
367    CALL init_param_gw_phy(zmea_phy, zpic_phy, zval_phy, zstd_phy, zsig_phy, zgam_phy, zthe_phy)
368
369  END SUBROUTINE param_gravity_wave_legacy
370
371  SUBROUTINE compute_smoothing(X)
372  USE icosa
373  IMPLICIT NONE
374    REAL(rstd),INTENT(INOUT)  :: X(iim*jjm)
375    REAL(rstd)                :: tmp(iim*jjm)
376    INTEGER :: i,j,ij
377
378    DO j=jj_begin,jj_end
379      DO i=ii_begin,ii_end
380        ij=(j-1)*iim+i
381        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.
382      ENDDO
383    ENDDO
384 
385    DO j=jj_begin,jj_end
386      DO i=ii_begin,ii_end
387        ij=(j-1)*iim+i
388        X(ij) = tmp(ij)
389      ENDDO
390    ENDDO
391 
392  END SUBROUTINE
393
394  SUBROUTINE compute_zstd(zmea, ztz, zstd)
395  USE icosa
396  IMPLICIT NONE
397    REAL(rstd),INTENT(IN)   :: zmea(iim*jjm)
398    REAL(rstd),INTENT(IN)   :: ztz(iim*jjm)
399    REAL(rstd),INTENT(OUT) :: zstd(iim*jjm)
400
401    INTEGER :: i,j,ij
402
[4966]403    DO j=jj_begin-1,jj_end+1
404      DO i=ii_begin-1,ii_end+1
[4855]405        ij=(j-1)*iim+i
406        zstd(ij) = ztz(ij) - zmea(ij)*zmea(ij)
[4966]407        IF (zstd(ij)<=0) zstd(ij)=0
[4855]408        zstd(ij)=SQRT(zstd(ij))
409      ENDDO
410    ENDDO
411 
412  END SUBROUTINE compute_zstd
413
414
415  SUBROUTINE compute_sigma_gamma_theta(zxtzx, zytzy, zxtzy, zsig, zgam, zthe)
416  USE icosa
417  IMPLICIT NONE
418    REAL(rstd),INTENT(IN)   :: zxtzx(iim*jjm)
419    REAL(rstd),INTENT(IN)   :: zytzy(iim*jjm)
420    REAL(rstd),INTENT(IN)   :: zxtzy(iim*jjm)
421    REAL(rstd),INTENT(OUT)  :: zsig(iim*jjm)
422    REAL(rstd),INTENT(OUT)  :: zgam(iim*jjm)
423    REAL(rstd),INTENT(OUT)  :: zthe(iim*jjm)
424   
425    REAL(rstd) :: xk,xl,xm,xp,xq,xw
426
427    INTEGER :: i,j,ij
428
429    DO j=jj_begin,jj_end
430      DO i=ii_begin,ii_end
431     
432        ij=(j-1)*iim+i
433        xk=(zxtzx(ij)+zytzy(ij))/2.
434        xl=(zxtzx(ij)-zytzy(ij))/2.
435        xm=zxtzy(ij)
436        xp=xk-SQRT(xl**2+xm**2)
437        xq=xk+SQRT(xl**2+xm**2)
438        xw=1.e-8
439        IF(xp<=xw) xp=0.
440        IF(xq<=xw) xq=xw
441        IF(ABS(xm)<=xw) xm=xw*SIGN(1.,xm)
442       
443        !--- SLOPE, ANISOTROPY AND THETA ANGLE
444     
445        zsig(ij)=SQRT(xq)
446        zgam(ij)=xp/xq
447        zthe(ij)=90.*ATAN2(xm,xl)/Pi
448
449      ENDDO
450    ENDDO
451
452  END SUBROUTINE compute_sigma_gamma_theta
453
454  SUBROUTINE compute_masking(mask, zmea, zpic, zval, zstd, zsig, zgam, zthe)
455  USE icosa
456  IMPLICIT NONE
457    REAL(rstd),INTENT(IN)   :: mask(iim*jjm)
458    REAL(rstd),INTENT(OUT)   :: zmea(iim*jjm)
459    REAL(rstd),INTENT(OUT)   :: zpic(iim*jjm)
460    REAL(rstd),INTENT(OUT)   :: zval(iim*jjm)
461    REAL(rstd),INTENT(OUT)   :: zstd(iim*jjm)
462    REAL(rstd),INTENT(OUT)   :: zsig(iim*jjm)
463    REAL(rstd),INTENT(OUT)   :: zgam(iim*jjm)
464    REAL(rstd),INTENT(OUT)   :: zthe(iim*jjm)
465   
466
467    INTEGER :: i,j,ij
468
469    DO j=jj_begin,jj_end
470      DO i=ii_begin,ii_end
471        ij=(j-1)*iim+i
472        IF (mask(ij)<0.1) THEN
473          zmea(ij)=0.
474          zpic(ij)=0.
475          zval(ij)=0.
476          zstd(ij)=0.
477          zsig(ij)=0.
478          zgam(ij)=0.
479          zthe(ij)=0.
480        ENDIF
481      ENDDO
482    ENDDO
483 
484  END SUBROUTINE compute_masking
485
486  SUBROUTINE compute_regular_param
487  USE xios_mod
488  USE icosa
489  IMPLICIT NONE
490    REAL, ALLOCATABLE :: mask(:,:) , zmea(:,:), zpic(:,:), zval(:,:), ztz(:,:), zstd(:,:)
491    REAL, ALLOCATABLE :: zytzy(:,:) , zxtzx(:,:), zxtzy(:,:), zsig(:,:) , zgam(:,:), zthe(:,:)
492    REAL :: xk, xl, xm, xp, xq, xw
493    INTEGER :: ni, nj, ibegin, jbegin, ni_glo, nj_glo
494    INTEGER :: i,j
495
496    CALL xios_get_domain_attr("regular_gw",ibegin=ibegin, ni=ni, ni_glo=ni_glo, jbegin=jbegin, nj=nj,  nj_glo=nj_glo)
497
498    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))
499    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))
500    ALLOCATE(zsig(0:ni+1,0:nj+1), zgam(0:ni+1,0:nj+1), zthe(0:ni+1,0:nj+1))
501
502    CALL xios_recv_field("mask_reg_exp",mask)
503    CALL xios_recv_field("zmea_reg_exp",zmea)
504    CALL xios_recv_field("zpic_reg_exp",zpic)
505    CALL xios_recv_field("zval_reg_exp",zval)
506    CALL xios_recv_field("ztz_reg_exp",ztz)
507    CALL xios_recv_field("zytzy_reg_exp",zytzy)
508    CALL xios_recv_field("zxtzx_reg_exp",zxtzx)
509    CALL xios_recv_field("zxtzy_reg_exp",zxtzy)
510   
511    DO j=0,nj+1
512      DO i=0, ni+1
513        zstd(i,j) = ztz(i,j)-zmea(i,j)*zmea(i,j)
514      ENDDO
515    ENDDO
516
517    ! smoothing
518    CALL smoothing_regular(ni,nj, zmea)
519    CALL smoothing_regular(ni,nj, zpic)
520    CALL smoothing_regular(ni,nj, zval)
521    CALL smoothing_regular(ni,nj, zstd)
522    CALL smoothing_regular(ni,nj, zytzy)
523    CALL smoothing_regular(ni,nj, zxtzx)
524    CALL smoothing_regular(ni,nj, zxtzy)
525
526    DO j=1,nj
527      DO i=1, ni
528        xk=(zxtzx(i,j)+zytzy(i,j))/2.
529        xl=(zxtzx(i,j)-zytzy(i,j))/2.
530        xm=zxtzy(i,j)
531        xp=xk-SQRT(xl**2+xm**2)
532        xq=xk+SQRT(xl**2+xm**2)
533        xw=1.e-8
534        IF(xp<=xw) xp=0.
535        IF(xq<=xw) xq=xw
536        IF(ABS(xm)<=xw) xm=xw*SIGN(1.,xm)
537       
538        !--- SLOPE, ANISOTROPY AND THETA ANGLE
539     
540        zsig(i,j)=SQRT(xq)
541        zgam(i,j)=xp/xq
542        zthe(i,j)=90.*ATAN2(xm,xl)/Pi
543      ENDDO
544    ENDDO
545
546    DO j=1,nj
547      DO i=1,ni
548        IF (mask(i,j)<0.1) THEN
549          zmea(i,j)=0.
550          zpic(i,j)=0.
551          zval(i,j)=0.
552          zstd(i,j)=0.
553          zsig(i,j)=0.
554          zgam(i,j)=0.
555          zthe(i,j)=0.
556        ENDIF
557      ENDDO
558    ENDDO
559
560    CALL xios_send_field("zmask_regout", mask(1:ni,1:nj))
561    CALL xios_send_field("zmea_regout", zmea(1:ni,1:nj))
562    CALL xios_send_field("zpic_regout", zpic(1:ni,1:nj))
563    CALL xios_send_field("zval_regout", zval(1:ni,1:nj))
564    CALL xios_send_field("zsig_regout", zsig(1:ni,1:nj))
565    CALL xios_send_field("zgam_regout", zgam(1:ni,1:nj))
566    CALL xios_send_field("zthe_regout", zthe(1:ni,1:nj))
567
568  END SUBROUTINE compute_regular_param
569
570  SUBROUTINE smoothing_regular(ni,nj,var)
571  IMPLICIT NONE
572    INTEGER, INTENT(IN) :: ni, nj
573    REAL, INTENT(INOUT) :: var(0:ni+1,0:nj+1)
574    REAL :: tmp(0:ni+1,0:nj+1)
575    INTEGER :: i,j
576
577    DO j=1,nj
578      DO i=1,ni
579        tmp(i,j) = (var(i,j) + 0.5*(var(i+1,j) + var(i-1,j) + var(i,j+1) +var(i,j-1))              &
580                             + 0.25*(var(i+1,j+1) + var(i-1,j+1) + var(i+1,j-1) + var(i-1,j-1)))/4.
581      ENDDO
582    ENDDO
583   
584    DO j=1,nj
585      DO i=1,ni
586        var(i,j)=tmp(i,j)
587      ENDDO
588    ENDDO
589
590  END SUBROUTINE smoothing_regular
591
592
[4858]593END MODULE icolmdz_param_gravity_wave
Note: See TracBrowser for help on using the repository browser.