source: CONFIG/UNIFORM/v7/IPSLCM7/SOURCES/DYNAMICO/trunk/dissip/sponge.f90 @ 6362

Last change on this file since 6362 was 6362, checked in by aclsce, 15 months ago

Modified sponge parameters ICOLMDZ to use same treatment as LMDZ (iflag_sponge=2, mode_sponge=3).
Remind :
iflag_sponge !0 --> for no sponge

!1 --> for sponge over 4 topmost layers
!2 --> for sponge from top to ~1% of top layer pressure

mode_sponge !1 --> u and v relax towards 0

!2 --> u and v relax towards their zonal mean
!3 --> u,v and pot. temp. relax towards their zonal mean
File size: 15.4 KB
Line 
1MODULE sponge_mod
2  USE icosa
3
4  PRIVATE
5 
6  REAL,SAVE :: tau_sponge !inverse of charactericstic relaxation time scale at the topmost layer (Hz)
7  INTEGER,SAVE :: iflag_sponge !0 --> for no sponge
8                               !1 --> for sponge over 4 topmost layers
9                               !2 --> for sponge from top to ~1% of top layer pressure
10  INTEGER,SAVE :: mode_sponge !1 --> u and v relax towards 0
11                              !2 --> u and v relax towards their zonal mean
12                              !3 --> u,v and pot. temp. relax towards their zonal mean
13 
14  !$OMP THREADPRIVATE(tau_sponge,iflag_sponge,mode_sponge)
15  REAL,ALLOCATABLE,SAVE :: rdamp(:) ! quenching coefficient
16  !$OMP THREADPRIVATE(rdamp)
17  REAL,ALLOCATABLE,SAVE:: lambda(:) ! inverse or quenching time scale (Hz)
18  !$OMP THREADPRIVATE(lambda)
19
20  PUBLIC sponge, init_sponge, iflag_sponge, pre_init_sponge
21
22  INTEGER,SAVE :: ni_glo
23  !$OMP THREADPRIVATE(ni_glo)
24  INTEGER,SAVE :: nj_glo
25  !$OMP THREADPRIVATE(nj_glo)
26  INTEGER,SAVE :: ni
27  !$OMP THREADPRIVATE(ni)
28  INTEGER,SAVE :: nj
29  !$OMP THREADPRIVATE(nj)
30  INTEGER,SAVE :: ibegin
31  !$OMP THREADPRIVATE(ibegin)
32  INTEGER,SAVE :: jbegin
33  !$OMP THREADPRIVATE(jbegin)
34 
35  TYPE(t_field),POINTER,SAVE :: f_hflux(:)
36  TYPE(t_field),POINTER,SAVE :: f_ulon(:)
37  TYPE(t_field),POINTER,SAVE :: f_ulat(:)
38  TYPE(t_field),POINTER,SAVE :: f_ulon_sponge(:)
39  TYPE(t_field),POINTER,SAVE :: f_ulat_sponge(:)
40  TYPE(t_field),POINTER,SAVE :: f_rhodz_sponge(:)
41  TYPE(t_field),POINTER,SAVE :: f_theta_sponge(:)
42
43  TYPE(t_field),POINTER,SAVE :: f_dulon(:)
44  TYPE(t_field),POINTER,SAVE :: f_dulat(:)
45 
46  INTEGER,SAVE  :: n_sponge
47  !$OMP THREADPRIVATE(n_sponge)
48  INTEGER,SAVE  :: begin_sponge
49  !$OMP THREADPRIVATE(begin_sponge)
50 
51CONTAINS
52
53  SUBROUTINE pre_init_sponge
54  USE icosa
55  USE disvert_mod
56  IMPLICIT NONE
57    INTEGER :: l
58
59    tau_sponge = 1.e-5
60    CALL getin("tau_sponge",tau_sponge)
61    PRINT*,'tau_sponge = ',tau_sponge
62
63    iflag_sponge = 0
64    CALL getin("iflag_sponge",iflag_sponge)
65    PRINT*,'iflag_sponge = ',iflag_sponge
66   
67    mode_sponge = 1
68    CALL getin("mode_sponge",mode_sponge)
69    PRINT*,'mode_sponge = ',mode_sponge
70
71    IF (iflag_sponge == 0) THEN
72      PRINT*,'init_sponge: no sponge'
73      RETURN
74    ENDIF
75   
76    IF (.NOT. (iflag_sponge == 1 .OR. iflag_sponge == 2 .OR. iflag_sponge == 3)) THEN
77      PRINT*,'Bad selector for variable iflag_sponge : <',iflag_sponge,'> options are 0,1,2'
78      STOP
79    ELSE
80      IF (.NOT. (mode_sponge == 1 .OR. mode_sponge == 2 .OR. mode_sponge == 3)) THEN
81        PRINT*,'Bad selector for variable mode_sponge : <',mode_sponge,'> options are 0,1,2'
82        STOP
83      ENDIF
84    ENDIF
85
86
87    IF (iflag_sponge == 0) RETURN
88   
89    ALLOCATE(rdamp(llm))
90    ALLOCATE(lambda(llm))
91    !$acc enter data create(rdamp(:),lambda(:)) async
92   
93    IF (iflag_sponge == 1) THEN
94! sponge quenching over the topmost 4 atmospheric layers
95        lambda(:)=0.
96        lambda(llm)=tau_sponge
97        lambda(llm-1)=tau_sponge/2.
98        lambda(llm-2)=tau_sponge/4.
99        lambda(llm-3)=tau_sponge/8.
100        n_sponge=4
101        begin_sponge=llm-3
102    ELSE IF (iflag_sponge == 2 ) THEN
103! sponge quenching over topmost layers down to pressures which are
104! higher than 100 times the topmost layer pressure
105        lambda(:)=tau_sponge*max(presnivs(llm)/presnivs(:)-0.01,0.)
106        DO l=llm,1,-1
107          IF (lambda(l)/=0.) begin_sponge=l
108        ENDDO
109        n_sponge=llm-begin_sponge+1
110    ENDIF
111     
112    IF (iflag_sponge>=1 .AND. (mode_sponge==2 .OR. mode_sponge==3)) THEN
113     CALL pre_init_sponge_zonal_mean_u
114    ENDIF
115   
116  END SUBROUTINE pre_init_sponge
117 
118 
119  SUBROUTINE init_sponge
120  USE icosa
121  USE disvert_mod
122  USE omp_para
123  USE mpipara, ONLY: is_mpi_master
124  IMPLICIT NONE
125    INTEGER :: l
126 
127    IF (iflag_sponge == 0) RETURN
128
129    rdamp(:)=itau_dissip*lambda(:)
130
131
132    IF (is_mpi_master) THEN
133      PRINT*,'lambda = '
134      DO l=ll_begin,ll_end
135        PRINT*,l,lambda(l)
136      ENDDO
137      PRINT*,'rdamp = '
138      DO l=ll_begin,ll_end
139        PRINT*,l,rdamp(l)
140      ENDDO
141    ENDIF
142    !$acc update device(rdamp(:), lambda(:)) async
143 
144   IF (iflag_sponge>=1 .AND. (mode_sponge == 2 .OR. mode_sponge == 3)) THEN
145     CALL init_sponge_zonal_mean_u
146   ENDIF
147   
148  END SUBROUTINE init_sponge
149
150 
151  SUBROUTINE sponge(f_mass, f_ue,f_due,f_theta_rhodz,f_dtheta_rhodz)
152  USE icosa
153  USE theta2theta_rhodz_mod
154  USE exner_mod
155  USE geopotential_mod
156  USE trace
157  USE time_mod
158  USE omp_para
159  IMPLICIT NONE
160    TYPE(t_field),POINTER :: f_ue(:)
161    TYPE(t_field),POINTER :: f_mass(:)
162    TYPE(t_field),POINTER :: f_due(:)
163    TYPE(t_field),POINTER :: f_theta_rhodz(:)
164    TYPE(t_field),POINTER :: f_dtheta_rhodz(:)
165
166    REAL(rstd),POINTER         :: due(:,:)!,due_sponge(:,:)
167    REAL(rstd),POINTER         :: ue(:,:)
168    REAL(rstd),POINTER         :: theta_rhodz(:,:,:)
169    REAL(rstd),POINTER         :: dtheta_rhodz(:,:,:)!,dtheta_sponge(:,:)
170
171    INTEGER :: ind
172    INTEGER :: l,ij
173
174!$OMP BARRIER
175   
176    CALL trace_start("sponge")
177
178    IF (mode_sponge == 1) THEN
179      DO ind=1,ndomain
180        IF (.NOT. assigned_domain(ind)) CYCLE
181        CALL swap_dimensions(ind)
182        CALL swap_geometry(ind)
183        ue=f_ue(ind) 
184        due=f_due(ind) 
185        theta_rhodz=f_theta_rhodz(ind)
186        dtheta_rhodz=f_dtheta_rhodz(ind)
187
188       !$acc parallel loop collapse(2) present(ue(:,:), due(:,:), dtheta_rhodz)  async
189        DO l=ll_begin,ll_end
190         !$SIMD
191          DO ij=ij_begin,ij_end
192
193            due(ij+u_right,l) = ue(ij+u_right,l)
194            due(ij+u_lup,l)   = ue(ij+u_lup,l)
195            due(ij+u_ldown,l) = ue(ij+u_ldown,l)
196
197            dtheta_rhodz(ij,l,:) = 0.0
198          ENDDO
199        ENDDO
200       !$acc end parallel loop
201      END DO
202    ELSE IF (mode_sponge == 2 .OR. mode_sponge == 3) THEN
203      CALL sponge_zonal_mean_u(mode_sponge, f_mass, f_ue,f_due,f_theta_rhodz,f_dtheta_rhodz) 
204    ELSE
205      PRINT*,'Arnaud Bad selector for variable mode_sponge : <',mode_sponge,             &
206            '> options 2 and 3 not available for the moment!'
207      STOP
208    ENDIF
209
210    DO ind=1,ndomain
211      IF (.NOT. assigned_domain(ind)) CYCLE
212      CALL swap_dimensions(ind)
213      CALL swap_geometry(ind)
214      due=f_due(ind) 
215      dtheta_rhodz=f_dtheta_rhodz(ind)
216
217      !$acc parallel loop collapse(2) present(ue(:,:), due(:,:), rdamp(:), dtheta_rhodz)  async
218      DO l=ll_begin,ll_end
219      !$SIMD
220        DO ij=ij_begin,ij_end
221          due(ij+u_right,l)        = -rdamp(l)*due(ij+u_right,l)
222          due(ij+u_lup,l)          = -rdamp(l)*due(ij+u_lup,l)
223          due(ij+u_ldown,l)        = -rdamp(l)*due(ij+u_ldown,l)
224          dtheta_rhodz(ij,l,:)     = -rdamp(l)*dtheta_rhodz(ij,l,:)
225        ENDDO
226      ENDDO
227     !$acc end parallel loop
228    END DO
229
230   CALL trace_end("sponge")
231
232!$OMP BARRIER
233  END SUBROUTINE sponge
234
235  SUBROUTINE pre_init_sponge_zonal_mean_u()
236  USE xios_mod
237  USE grid_param, ONLY : iim_glo
238  USE disvert_mod
239  IMPLICIT NONE
240  INTEGER :: nlon
241  INTEGER :: nlat 
242
243!$OMP MASTER
244    nlon=360.*iim_glo/80.
245    nlat=180.*iim_glo/80.
246
247    CALL xios_set_axis_attr("lev_sponge", n_glo=n_sponge, value=presnivs(begin_sponge:))
248    CALL xios_set_domain_attr("to_sponge", ni_glo=nlon, nj_glo=nlat)
249    CALL xios_set_axis_attr("reduced_sponge", n_glo=nlat)
250    CALL xios_set_field_attr("ulon_sponge_reduced", read_access=.TRUE.)
251    CALL xios_set_field_attr("ulon_from_sponge", read_access=.TRUE.)
252    CALL xios_set_field_attr("ulat_sponge_reduced", read_access=.TRUE.)
253    CALL xios_set_field_attr("ulat_from_sponge", read_access=.TRUE.)
254    CALL xios_set_field_attr("rhodz_sponge_reduced", read_access=.TRUE.)
255    IF (mode_sponge==3) CALL xios_set_field_attr("theta_sponge_reduced", read_access=.TRUE.)
256    IF (mode_sponge==3) CALL xios_set_field_attr("theta_from_sponge", read_access=.TRUE.)
257!$OMP END MASTER
258
259  END SUBROUTINE  pre_init_sponge_zonal_mean_u
260
261
262  SUBROUTINE init_sponge_zonal_mean_u()
263  USE xios_mod
264  IMPLICIT NONE
265   
266!$OMP MASTER
267    CALL xios_get_domain_attr("to_sponge", ni_glo=ni_glo, nj_glo=nj_glo,ni=ni, nj=nj, ibegin=ibegin,jbegin=jbegin)
268!$OMP END MASTER
269    CALL allocate_field(f_hflux, field_u, type_real, llm, name='hflux') 
270    CALL allocate_field(f_ulon, field_t, type_real, llm, name='ulon') 
271    CALL allocate_field(f_ulat, field_t, type_real, llm, name='ulat') 
272    CALL allocate_field(f_ulon_sponge, field_t, type_real, n_sponge, name='ulon_sponge') 
273    CALL allocate_field(f_ulat_sponge, field_t, type_real, n_sponge, name='ulat_sponge') 
274    CALL allocate_field(f_rhodz_sponge, field_t, type_real, n_sponge, name='rhodz_sponge') 
275    CALL allocate_field(f_theta_sponge, field_t, type_real, n_sponge, name='theta_sponge') 
276    CALL allocate_field(f_dulon, field_t, type_real, llm,name='dulon_sponge') 
277    CALL allocate_field(f_dulat, field_t, type_real, llm,name='dulat_sponge') 
278
279  END SUBROUTINE init_sponge_zonal_mean_u
280
281 
282  SUBROUTINE sponge_zonal_mean_u(mode_sponge, f_rhodz, f_ue,f_due,f_theta_rhodz,f_dtheta_rhodz)
283  USE icosa
284  USE theta2theta_rhodz_mod
285  USE exner_mod
286  USE geopotential_mod
287  USE trace
288  USE time_mod
289  USE omp_para
290  USE wind_mod
291  USE xios_mod
292  IMPLICIT NONE
293    INTEGER,INTENT(IN)    :: mode_sponge
294    TYPE(t_field),POINTER :: f_rhodz(:)
295    TYPE(t_field),POINTER :: f_ue(:)
296    TYPE(t_field),POINTER :: f_due(:)
297    TYPE(t_field),POINTER :: f_theta_rhodz(:)
298    TYPE(t_field),POINTER :: f_dtheta_rhodz(:)
299
300    REAL(rstd),POINTER         :: due(:,:)!,due_sponge(:,:)
301    REAL(rstd),POINTER         :: ue(:,:)
302    REAL(rstd),POINTER         :: rhodz(:,:)
303    REAL(rstd),POINTER         :: hflux(:,:)
304    REAL(rstd),POINTER         :: theta_rhodz(:,:,:)
305    REAL(rstd),POINTER         :: dtheta_rhodz(:,:,:)!,dtheta_sponge(:,:)
306    REAL(rstd),POINTER         :: ulon(:,:)
307    REAL(rstd),POINTER         :: ulat(:,:)
308    REAL(rstd),POINTER         :: ulon_sponge(:,:)
309    REAL(rstd),POINTER         :: ulat_sponge(:,:)
310    REAL(rstd),POINTER         :: theta_sponge(:,:)
311    REAL(rstd),POINTER         :: rhodz_sponge(:,:)
312   
313    REAL(rstd) :: ulon_zonal_mean(nj_glo,n_sponge)
314    REAL(rstd) :: ulat_zonal_mean(nj_glo,n_sponge)
315    REAL(rstd) :: theta_zonal_mean(nj_glo,n_sponge)
316    REAL(rstd) :: rhodz_zonal_mean(nj_glo,n_sponge)
317    REAL(rstd) :: domain_zonal_mean(ni,nj,n_sponge)
318   
319    INTEGER :: i,j,ind,ij,l,ll 
320
321    DO ind = 1 , ndomain
322      CALL swap_dimensions(ind)
323      CALL swap_geometry(ind)
324      ue=f_ue(ind)
325      rhodz=f_rhodz(ind)
326      hflux = f_hflux(ind)
327         
328      DO l = ll_begin, ll_end
329         !DIR$ IVDEP
330         DO ij=ij_begin_ext,ij_end_ext         
331            hflux(ij+u_right,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_right,l))*ue(ij+u_right,l)
332            hflux(ij+u_lup,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_lup,l))*ue(ij+u_lup,l)
333            hflux(ij+u_ldown,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_ldown,l))*ue(ij+u_ldown,l)
334         ENDDO
335       ENDDO
336    ENDDO
337
338             
339    CALL transfert_request(f_hflux,req_e1_vect)
340    CALL un2ulonlat(f_hflux, f_ulon, f_ulat)
341
342    DO ind = 1 , ndomain
343      CALL swap_dimensions(ind)
344      CALL swap_geometry(ind)
345     
346      ulon=f_ulon(ind)
347      ulat=f_ulat(ind)
348      ulon_sponge = f_ulon_sponge(ind)
349      ulat_sponge = f_ulat_sponge(ind)
350      rhodz_sponge = f_rhodz_sponge(ind)
351      theta_sponge = f_theta_sponge(ind)
352      theta_rhodz = f_theta_rhodz(ind)
353             
354      DO l = 1, llm
355         IF (l<begin_sponge) CYCLE
356         ll=l-begin_sponge+1
357         !DIR$ IVDEP
358         DO ij=ij_begin_ext,ij_end_ext         
359            ulon_sponge(ij,ll)= ulon(ij,l) 
360            ulat_sponge(ij,ll)= ulat(ij,l) 
361            rhodz_sponge(ij,ll) = rhodz(ij,l)
362            IF (mode_sponge==3) theta_sponge(ij,ll) = theta_rhodz(ij,l,1)
363         ENDDO
364       ENDDO
365    ENDDO
366
367!$OMP BARRIER
368   
369    CALL xios_write_field("ulon_sponge",f_ulon_sponge)
370    CALL xios_write_field("ulat_sponge",f_ulat_sponge)
371    CALL xios_write_field("rhodz_sponge",f_rhodz_sponge)
372!$OMP MASTER
373    CALL xios_recv_field("ulon_sponge_reduced",ulon_zonal_mean)
374    CALL xios_recv_field("ulat_sponge_reduced",ulat_zonal_mean)
375    CALL xios_recv_field("rhodz_sponge_reduced",rhodz_zonal_mean)
376!$OMP END MASTER
377
378    IF (mode_sponge==3) THEN
379      CALL xios_write_field("theta_sponge",f_theta_sponge)
380!$OMP MASTER
381      CALL xios_recv_field("theta_sponge_reduced",theta_zonal_mean)
382!$OMP END MASTER
383    ENDIF
384   
385!$OMP MASTER
386    DO j=1,nj
387      DO i=1,ni
388        domain_zonal_mean(i,j,:)=ulon_zonal_mean(jbegin+j,:)/rhodz_zonal_mean(jbegin+j,:)
389      ENDDO
390    ENDDO
391   
392    CALL xios_send_field("ulon_sponge_mean",domain_zonal_mean)
393!$OMP END MASTER
394    CALL xios_read_field("ulon_from_sponge",f_ulon_sponge)
395 
396!$OMP MASTER
397    DO j=1,nj
398      DO i=1,ni
399        domain_zonal_mean(i,j,:)=ulat_zonal_mean(jbegin+j,:)/rhodz_zonal_mean(jbegin+j,:)
400      ENDDO
401    ENDDO
402   
403    CALL xios_send_field("ulat_sponge_mean",domain_zonal_mean)
404!$OMP END MASTER
405    CALL xios_read_field("ulat_from_sponge",f_ulat_sponge)
406   
407    IF (mode_sponge==3) THEN
408!$OMP MASTER
409      DO j=1,nj
410        DO i=1,ni
411          domain_zonal_mean(i,j,:)=theta_zonal_mean(jbegin+j,:)
412        ENDDO
413      ENDDO
414      CALL xios_send_field("theta_sponge_mean",domain_zonal_mean)
415!$OMP END MASTER
416      CALL xios_read_field("theta_from_sponge",f_theta_sponge)
417    ENDIF
418   
419
420    DO ind = 1 , ndomain
421      CALL swap_dimensions(ind)
422      CALL swap_geometry(ind)
423     
424      ulon=f_ulon(ind)
425      ulat=f_ulat(ind)
426      ulon_sponge = f_ulon_sponge(ind)
427      ulat_sponge = f_ulat_sponge(ind)
428      theta_sponge = f_theta_sponge(ind)
429      dtheta_rhodz = f_dtheta_rhodz(ind)   
430      theta_rhodz = f_theta_rhodz(ind)
431      DO l = ll_begin, ll_end
432        IF (l<begin_sponge) CYCLE
433        ll=l-begin_sponge+1
434        IF (mode_sponge==2) THEN
435          DO ij=ij_begin_ext,ij_end_ext         
436            ulon(ij,l) = ulon_sponge(ij,ll) 
437            ulat(ij,l) = ulat_sponge(ij,ll) 
438            dtheta_rhodz(ij,l,1)=theta_rhodz(ij,ll,1) 
439          ENDDO
440        ELSE  IF (mode_sponge==3) THEN
441          DO ij=ij_begin_ext,ij_end_ext         
442            ulon(ij,l) = ulon_sponge(ij,ll) 
443            ulat(ij,l) = ulat_sponge(ij,ll) 
444            dtheta_rhodz(ij,l,1) = theta_sponge(ij,ll)
445          ENDDO
446        ENDIF
447      ENDDO
448    ENDDO
449   
450    CALL transfert_request(f_ulon,req_i0)
451    CALL transfert_request(f_ulat,req_i0)
452    IF (mode_sponge==3) THEN
453      CALL transfert_request(f_dtheta_rhodz,req_i0)
454    ENDIF
455    CALL ulonlat2un(f_ulon,f_ulat,f_due)
456!$OMP BARRIER
457   
458    DO ind = 1 , ndomain
459      CALL swap_dimensions(ind)
460      CALL swap_geometry(ind)
461      due=f_due(ind)
462      ue=f_ue(ind)
463      dtheta_rhodz=f_dtheta_rhodz(ind)
464      theta_rhodz=f_theta_rhodz(ind)
465 
466       DO l = ll_begin, ll_end
467         IF (l<begin_sponge) THEN
468           !DIR$ IVDEP
469           DO ij=ij_begin,ij_end         
470              due(ij+u_right,l) = 0.
471              due(ij+u_lup,l)   = 0.
472              due(ij+u_ldown,l) = 0.
473              dtheta_rhodz(ij,l,:) = 0.
474           ENDDO         
475         ELSE
476           !DIR$ IVDEP
477           DO ij=ij_begin,ij_end         
478             due(ij+u_right,l) = ue(ij+u_right,l) - due(ij+u_right,l)
479             due(ij+u_lup,l)   = ue(ij+u_lup,l)   - due(ij+u_lup,l)
480             due(ij+u_ldown,l) = ue(ij+u_ldown,l) - due(ij+u_ldown,l)
481             dtheta_rhodz(ij,l,1) = theta_rhodz(ij,l,1) - dtheta_rhodz(ij,l,1)
482          ENDDO
483        ENDIF
484      ENDDO
485    ENDDO 
486   
487    CALL transfert_request(f_due,req_e1_vect)
488!ym to check results
489!    CALL un2ulonlat(f_due, f_dulon, f_dulat)
490!    CALL xios_write_field("dulon_sponge",f_dulon)
491!    CALL xios_write_field("dulat_sponge",f_dulat)
492!    CALL xios_write_field("dtheta_sponge",f_dtheta_rhodz)
493   
494  END SUBROUTINE sponge_zonal_mean_u
495   
496
497END MODULE sponge_mod
Note: See TracBrowser for help on using the repository browser.