source: codes/icosagcm/devel/src/dissip/dissip_gcm.f90 @ 714

Last change on this file since 714 was 714, checked in by dubos, 6 years ago

devel : backported from trunk commits r607,r648,r649,r667,r668,r669,r706
File size: 28.1 KB
Line 
1MODULE dissip_gcm_mod
2  USE icosa
3
4  PRIVATE
5
6  TYPE(t_field),POINTER,SAVE :: f_due_diss1(:)
7  TYPE(t_field),POINTER,SAVE :: f_due_diss2(:)
8
9  TYPE(t_field),POINTER,SAVE :: f_dtheta_diss(:)
10  TYPE(t_field),POINTER,SAVE :: f_dtheta_rhodz_diss(:)
11  TYPE(t_message),SAVE :: req_due, req_dtheta 
12 
13  INTEGER,SAVE :: nitergdiv=1
14!$OMP THREADPRIVATE(nitergdiv)
15  INTEGER,SAVE :: nitergrot=1
16!$OMP THREADPRIVATE(nitergrot)
17  INTEGER,SAVE :: niterdivgrad=1
18!$OMP THREADPRIVATE(niterdivgrad)
19
20  REAL,ALLOCATABLE,SAVE :: tau_graddiv(:)
21!$OMP THREADPRIVATE(tau_graddiv)
22  REAL,ALLOCATABLE,SAVE :: tau_gradrot(:)
23!$OMP THREADPRIVATE(tau_gradrot)
24  REAL,ALLOCATABLE,SAVE :: tau_divgrad(:)
25!$OMP THREADPRIVATE(tau_divgrad)
26 
27  REAL,SAVE :: cgraddiv
28!$OMP THREADPRIVATE(cgraddiv)
29  REAL,SAVE :: cgradrot
30!$OMP THREADPRIVATE(cgradrot)
31  REAL,SAVE :: cdivgrad
32!$OMP THREADPRIVATE(cdivgrad)
33
34  INTEGER, SAVE :: rayleigh_friction_type, rayleigh_shear
35!$OMP THREADPRIVATE(rayleigh_friction_type)
36!$OMP THREADPRIVATE(rayleigh_shear)
37  REAL, SAVE    :: rayleigh_tau, rayleigh_limlat
38!$OMP THREADPRIVATE(rayleigh_tau)
39!$OMP THREADPRIVATE(rayleigh_limlat)
40 
41  REAL,SAVE    :: dtdissip
42!$OMP THREADPRIVATE(dtdissip)
43 
44  PUBLIC init_dissip, dissip
45CONTAINS
46
47  SUBROUTINE allocate_dissip
48  USE icosa
49  IMPLICIT NONE 
50    CALL allocate_field(f_due_diss1,field_u,type_real,llm)
51    CALL allocate_field(f_due_diss2,field_u,type_real,llm)
52    CALL allocate_field(f_dtheta_diss,field_t,type_real,llm)
53    CALL allocate_field(f_dtheta_rhodz_diss,field_t,type_real,llm)
54    ALLOCATE(tau_graddiv(llm))
55    ALLOCATE(tau_gradrot(llm))   
56    ALLOCATE(tau_divgrad(llm))
57  END SUBROUTINE allocate_dissip
58 
59  SUBROUTINE init_dissip
60  USE icosa
61  USE disvert_mod
62  USE mpi_mod
63  USE mpipara
64  USE transfert_mod
65  USE time_mod
66  USE transfert_omp_mod
67  USE omp_para
68  IMPLICIT NONE
69 
70  TYPE(t_field),POINTER,SAVE  :: f_u(:)
71  TYPE(t_field),POINTER,SAVE  :: f_du(:)
72  REAL(rstd),POINTER     :: u(:)
73  REAL(rstd),POINTER     :: du(:)
74  TYPE(t_field),POINTER,SAVE  :: f_theta(:)
75  TYPE(t_field),POINTER ,SAVE :: f_dtheta(:)
76  REAL(rstd),POINTER    :: theta(:)
77  REAL(rstd),POINTER    :: dtheta(:)
78  REAL(rstd)            :: dumax,dumax1
79  REAL(rstd)            :: dthetamax,dthetamax1
80  REAL(rstd)            :: r
81  REAL(rstd)            :: tau
82  REAL(rstd)            :: zz, zvert, fact
83  INTEGER               :: l
84  CHARACTER(len=255)    :: rayleigh_friction_key
85  REAL(rstd)            :: mintau
86  INTEGER               :: seed_size
87  INTEGER,ALLOCATABLE   :: seed(:)
88 
89           
90  INTEGER :: i,j,ij,ind,it,iter,M
91
92   rayleigh_friction_key='none'
93   CALL getin("rayleigh_friction_type",rayleigh_friction_key)
94   SELECT CASE(TRIM(rayleigh_friction_key))
95   CASE('none')
96      rayleigh_friction_type=0
97      IF (is_master) PRINT *, 'No Rayleigh friction'
98   CASE('dcmip2_schaer_noshear')
99      rayleigh_friction_type=1
100      rayleigh_shear=0
101      IF (is_master) PRINT *, 'Rayleigh friction : Schaer-like mountain without shear DCMIP2.1'
102   CASE('dcmip2_schaer_shear')
103      rayleigh_shear=1
104      rayleigh_friction_type=2
105      IF (is_master) PRINT *, 'Rayleigh friction : Schaer-like mountain with shear DCMIP2.2'
106   CASE('giant_liu_schneider') 
107      rayleigh_friction_type=99
108      IF (is_master) PRINT *, 'Rayleigh friction : giant planets Liu Schneider 2010'
109   CASE DEFAULT
110      IF (is_master) PRINT *, 'Bad selector : rayleigh_friction_type =', TRIM(rayleigh_friction_key), &
111        ' in dissip_gcm.f90/init_dissip'
112      STOP
113   END SELECT
114
115   IF(rayleigh_friction_type>0) THEN
116      rayleigh_tau=0.
117      CALL getin("rayleigh_friction_tau",rayleigh_tau)
118      rayleigh_tau = rayleigh_tau / scale_factor
119      IF(rayleigh_tau<=0) THEN
120         IF (is_master) PRINT *, 'Forbidden : negative value for rayleigh_friction_tau =',rayleigh_tau
121         STOP
122      END IF
123      IF(rayleigh_friction_type == 99) THEN
124         rayleigh_limlat=0.
125         CALL getin("rayleigh_limlat",rayleigh_limlat)
126         rayleigh_limlat = rayleigh_limlat*3.14159/180.
127      ENDIF
128   END IF
129
130    CALL allocate_dissip
131    CALL allocate_field(f_u,field_u,type_real)
132    CALL allocate_field(f_du,field_u,type_real)
133    CALL allocate_field(f_theta,field_t,type_real)
134    CALL allocate_field(f_dtheta,field_t,type_real)
135   
136    CALL init_message(f_due_diss1,req_e1_vect,req_due)
137    CALL init_message(f_dtheta_diss,req_i1,req_dtheta)
138
139    tau_graddiv(:)=5000
140    CALL getin("tau_graddiv",tau)
141    tau_graddiv(:)=tau/scale_factor
142
143    CALL getin("nitergdiv",nitergdiv)
144   
145    tau_gradrot(:)=5000
146    CALL getin("tau_gradrot",tau)
147    tau_gradrot(:)=tau/scale_factor
148
149    CALL getin("nitergrot",nitergrot)
150   
151
152    tau_divgrad(:)=5000
153    CALL getin("tau_divgrad",tau)
154    tau_divgrad(:)=tau/scale_factor
155
156    CALL getin("niterdivgrad",niterdivgrad)
157
158
159    cgraddiv=-1
160    cdivgrad=-1
161    cgradrot=-1
162
163!$OMP BARRIER
164!$OMP MASTER
165    DO ind=1,ndomain
166      CALL swap_dimensions(ind)
167      CALL swap_geometry(ind)
168      u=f_u(ind)
169
170! set random seed to get reproductibility when using a different number of process
171      CALL RANDOM_SEED(size=M)
172      CALL RANDOM_SEED(put=(/(i,i=1,M)/))
173
174      DO j=jj_begin,jj_end
175        DO i=ii_begin,ii_end
176          ij=(j-1)*iim+i   
177          CALL RANDOM_NUMBER(r)
178          u(ij+u_right)=r-0.5
179          CALL RANDOM_NUMBER(r)
180          u(ij+u_lup)=r-0.5
181          CALL RANDOM_NUMBER(r)
182          u(ij+u_ldown)=r-0.5
183        ENDDO
184      ENDDO       
185    ENDDO
186!$OMP END MASTER 
187!$OMP BARRIER
188
189    DO it=1,20
190     
191      dumax=0
192      DO iter=1,nitergdiv
193        CALL transfert_request(f_u,req_e1_vect)
194        DO ind=1,ndomain
195          IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE
196          CALL swap_dimensions(ind)
197          CALL swap_geometry(ind)
198          u=f_u(ind)
199          du=f_du(ind)
200          CALL compute_gradiv(u,du,1,1)
201          u=du
202        ENDDO
203      ENDDO
204
205      CALL transfert_request(f_du,req_e1_vect)
206     
207      DO ind=1,ndomain
208        IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE
209        CALL swap_dimensions(ind)
210        CALL swap_geometry(ind)
211        u=f_u(ind)
212        du=f_du(ind)
213         
214        DO j=jj_begin,jj_end
215          DO i=ii_begin,ii_end
216            ij=(j-1)*iim+i   
217            if (le(ij+u_right)>1e-100) dumax=MAX(dumax,ABS(du(ij+u_right)))
218            if (le(ij+u_lup)>1e-100)   dumax=MAX(dumax,ABS(du(ij+u_lup)))
219            if (le(ij+u_ldown)>1e-100) dumax=MAX(dumax,ABS(du(ij+u_ldown)))
220          ENDDO
221        ENDDO
222      ENDDO
223
224      IF (using_mpi) THEN
225        CALL reduce_max_omp(dumax,dumax1)
226!$OMP MASTER       
227        CALL MPI_ALLREDUCE(dumax1,dumax,1,MPI_REAL8,MPI_MAX,comm_icosa,ierr)
228!$OMP END MASTER     
229        CALL bcast_omp(dumax) 
230      ELSE
231        CALL allreduce_max_omp(dumax,dumax1)
232        dumax=dumax1
233      ENDIF 
234                       
235      DO ind=1,ndomain
236        IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE
237        CALL swap_dimensions(ind)
238        CALL swap_geometry(ind)
239        u=f_u(ind)
240        du=f_du(ind)
241        u=du/dumax
242      ENDDO
243      IF (is_master) PRINT *,"gradiv : it :",it ,": dumax",dumax
244
245    ENDDO 
246    IF (is_master) PRINT *,"gradiv : dumax",dumax
247    IF (is_master) PRINT *, 'mean T-cell edge size (km)', 1.45*radius/iim_glo/1000., &
248                              'effective T-cell half-edge size (km)', dumax**(-.5/nitergdiv)/1000
249    IF (is_master)  PRINT *, 'Max. time step assuming c=340 m/s and Courant number=3 (ARK2.3) :', &
250                               3./340.*dumax**(-.5/nitergdiv)
251
252    cgraddiv=dumax**(-1./nitergdiv)
253    IF (is_master) PRINT *,"cgraddiv : ",cgraddiv
254
255!$OMP BARRIER
256!$OMP MASTER
257    DO ind=1,ndomain
258      CALL swap_dimensions(ind)
259      CALL swap_geometry(ind)
260      u=f_u(ind)
261
262! set random seed to get reproductibility when using a different number of process
263      CALL RANDOM_SEED(size=M)
264      CALL RANDOM_SEED(put=(/(i,i=1,M)/))
265 
266       DO j=jj_begin,jj_end
267        DO i=ii_begin,ii_end
268          ij=(j-1)*iim+i   
269          CALL RANDOM_NUMBER(r)
270          u(ij+u_right)=r-0.5
271          CALL RANDOM_NUMBER(r)
272          u(ij+u_lup)=r-0.5
273          CALL RANDOM_NUMBER(r)
274          u(ij+u_ldown)=r-0.5
275        ENDDO
276      ENDDO       
277    ENDDO
278!$OMP END MASTER 
279!$OMP BARRIER
280
281
282    DO it=1,20
283 
284      dumax=0
285      DO iter=1,nitergrot
286        CALL transfert_request(f_u,req_e1_vect)
287        DO ind=1,ndomain
288          IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE
289          CALL swap_dimensions(ind)
290          CALL swap_geometry(ind)
291          u=f_u(ind)
292          du=f_du(ind)
293          CALL compute_gradrot(u,du,1,1)
294          u=du
295        ENDDO
296      ENDDO
297
298      CALL transfert_request(f_du,req_e1_vect)
299     
300      DO ind=1,ndomain
301        IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE
302        CALL swap_dimensions(ind)
303        CALL swap_geometry(ind)
304        u=f_u(ind)
305        du=f_du(ind)
306       
307        DO j=jj_begin,jj_end
308          DO i=ii_begin,ii_end
309            ij=(j-1)*iim+i   
310            if (le(ij+u_right)>1e-100) dumax=MAX(dumax,ABS(du(ij+u_right)))
311            if (le(ij+u_lup)>1e-100) dumax=MAX(dumax,ABS(du(ij+u_lup)))
312            if (le(ij+u_ldown)>1e-100) dumax=MAX(dumax,ABS(du(ij+u_ldown)))
313          ENDDO
314        ENDDO
315      ENDDO
316
317      IF (using_mpi) THEN
318        CALL reduce_max_omp(dumax,dumax1)
319!$OMP MASTER       
320        CALL MPI_ALLREDUCE(dumax1,dumax,1,MPI_REAL8,MPI_MAX,comm_icosa,ierr)
321!$OMP END MASTER     
322        CALL bcast_omp(dumax) 
323      ELSE
324        CALL allreduce_max_omp(dumax,dumax1)
325        dumax=dumax1
326      ENDIF 
327
328     
329      DO ind=1,ndomain
330        IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE
331        CALL swap_dimensions(ind)
332        CALL swap_geometry(ind)
333        u=f_u(ind)
334        du=f_du(ind)
335        u=du/dumax
336      ENDDO
337     
338      IF (is_master) PRINT *,"gradrot : it :",it ,": dumax",dumax
339
340    ENDDO 
341    IF (is_master) PRINT *,"gradrot : dumax",dumax
342 
343    cgradrot=dumax**(-1./nitergrot) 
344    IF (is_master) PRINT *,"cgradrot : ",cgradrot
345   
346
347
348!$OMP BARRIER
349!$OMP MASTER
350    DO ind=1,ndomain
351      CALL swap_dimensions(ind)
352      CALL swap_geometry(ind)
353      theta=f_theta(ind)
354
355! set random seed to get reproductibility when using a different number of process
356      CALL RANDOM_SEED(size=M)
357      CALL RANDOM_SEED(put=(/(i,i=1,M)/))
358
359      DO j=jj_begin,jj_end
360        DO i=ii_begin,ii_end
361          ij=(j-1)*iim+i   
362          CALL RANDOM_NUMBER(r)
363          theta(ij)=r-0.5
364        ENDDO
365      ENDDO 
366    ENDDO
367!$OMP END MASTER
368!$OMP BARRIER
369
370    DO it=1,20
371
372      dthetamax=0
373      DO iter=1,niterdivgrad
374        CALL transfert_request(f_theta,req_i1)
375        DO ind=1,ndomain
376          IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE
377          CALL swap_dimensions(ind)
378          CALL swap_geometry(ind)
379          theta=f_theta(ind)
380          dtheta=f_dtheta(ind)
381          CALL compute_divgrad(theta,dtheta,1,1)
382          theta=dtheta
383        ENDDO
384      ENDDO
385
386      CALL transfert_request(f_dtheta,req_i1)
387     
388      DO ind=1,ndomain
389        IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE
390        CALL swap_dimensions(ind)
391        CALL swap_geometry(ind)
392        theta=f_theta(ind)
393        dtheta=f_dtheta(ind)
394       
395        DO j=jj_begin,jj_end
396          DO i=ii_begin,ii_end
397            ij=(j-1)*iim+i   
398            dthetamax=MAX(dthetamax,ABS(dtheta(ij)))
399          ENDDO
400        ENDDO
401      ENDDO
402
403      IF (using_mpi) THEN
404        CALL reduce_max_omp(dthetamax ,dthetamax1)
405!$OMP MASTER       
406        CALL MPI_ALLREDUCE(dthetamax1,dthetamax,1,MPI_REAL8,MPI_MAX,comm_icosa,ierr)
407!$OMP END MASTER     
408        CALL bcast_omp(dthetamax) 
409      ELSE
410        CALL allreduce_max_omp(dthetamax,dthetamax1)
411        dumax=dumax1
412      ENDIF 
413     
414      IF (is_master) PRINT *,"divgrad : it :",it ,": dthetamax",dthetamax
415
416      DO ind=1,ndomain
417        IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE
418        CALL swap_dimensions(ind)
419        CALL swap_geometry(ind)
420        theta=f_theta(ind)
421        dtheta=f_dtheta(ind)
422        theta=dtheta/dthetamax
423      ENDDO
424    ENDDO 
425
426    IF (is_master) PRINT *,"divgrad : divgrad",dthetamax
427 
428    cdivgrad=dthetamax**(-1./niterdivgrad) 
429    IF (is_master) PRINT *,"cdivgrad : ",cdivgrad
430
431     
432    fact=2
433    DO l=1,llm
434       IF(ap_bp_present) THEN ! height-dependent dissipation
435          zz= 1. - preff/presnivs(l)
436       ELSE
437          zz = 0.
438       END IF
439       zvert=fact-(fact-1)/(1+zz*zz)
440       tau_graddiv(l) = tau_graddiv(l)/zvert
441       tau_gradrot(l) = tau_gradrot(l)/zvert
442       tau_divgrad(l) = tau_divgrad(l)/zvert
443    ENDDO
444
445    mintau=tau_graddiv(1)
446    DO l=1,llm
447      mintau=MIN(mintau,tau_graddiv(l))
448      mintau=MIN(mintau,tau_gradrot(l))
449      mintau=MIN(mintau,tau_divgrad(l))
450    ENDDO
451
452    IF(rayleigh_friction_type>0) mintau=MIN(mintau,rayleigh_tau)
453
454    IF(mintau>0) THEN
455       IF (itau_dissip==0) THEN
456         IF (is_master) PRINT *,"init_dissip: Automatic computation of itau_dissip..."
457         itau_dissip=INT(mintau/dt)
458       ENDIF
459    ELSE
460       IF (is_master) PRINT *,"init_dissip: No dissipation time set, setting itau_dissip to 1000000000"
461       itau_dissip=100000000
462    END IF
463    itau_dissip=MAX(1,itau_dissip)
464    dtdissip=itau_dissip*dt
465    IF (is_master) THEN
466      PRINT *,"init_dissip: rayleigh_tau",rayleigh_tau, "mintau ",mintau
467      PRINT *,"init_dissip: itau_dissip",itau_dissip," dtdissip ",dtdissip
468    ENDIF
469
470  END SUBROUTINE init_dissip 
471 
472 
473  SUBROUTINE dissip(f_ps,f_mass,f_phis,f_geopot,f_theta_rhodz,f_ue, f_dtheta_rhodz,f_due)
474  USE icosa
475  USE theta2theta_rhodz_mod
476  USE pression_mod
477  USE exner_mod
478  USE geopotential_mod
479  USE trace
480  USE time_mod
481  USE omp_para
482  IMPLICIT NONE
483    TYPE(t_field),POINTER :: f_ps(:), f_mass(:), f_phis(:), f_geopot(:)
484    TYPE(t_field),POINTER :: f_theta_rhodz(:), f_dtheta_rhodz(:)
485    TYPE(t_field),POINTER :: f_ue(:), f_due(:)
486
487    REAL(rstd),POINTER         :: due(:,:)
488    REAL(rstd),POINTER         :: phi(:,:), ue(:,:)
489    REAL(rstd),POINTER         :: due_diss1(:,:)
490    REAL(rstd),POINTER         :: due_diss2(:,:)
491    REAL(rstd),POINTER         :: dtheta_rhodz(:,:,:)
492    REAL(rstd),POINTER         :: dtheta_rhodz_diss(:,:)
493
494    INTEGER :: ind, shear
495    INTEGER :: l,ij,nn
496
497!$OMP BARRIER
498   
499    CALL trace_start("dissip")
500    CALL gradiv(f_ue,f_due_diss1)
501    CALL gradrot(f_ue,f_due_diss2)
502
503    CALL divgrad_theta_rhodz(f_mass,f_theta_rhodz,f_dtheta_rhodz_diss)
504   
505    DO ind=1,ndomain
506      IF (.NOT. assigned_domain(ind)) CYCLE
507      CALL swap_dimensions(ind)
508      CALL swap_geometry(ind)
509      due=f_due(ind) 
510      due_diss1=f_due_diss1(ind)
511      due_diss2=f_due_diss2(ind)
512      dtheta_rhodz=f_dtheta_rhodz(ind)
513      dtheta_rhodz_diss=f_dtheta_rhodz_diss(ind)
514
515      DO l=ll_begin,ll_end
516!DIR$ SIMD
517        DO ij=ij_begin,ij_end
518
519            due(ij+u_right,l) = -0.5*( due_diss1(ij+u_right,l)/tau_graddiv(l) + due_diss2(ij+u_right,l)/tau_gradrot(l))*itau_dissip 
520            due(ij+u_lup,l)   = -0.5*( due_diss1(ij+u_lup,l)  /tau_graddiv(l) + due_diss2(ij+u_lup,l)  /tau_gradrot(l))*itau_dissip
521            due(ij+u_ldown,l) = -0.5*( due_diss1(ij+u_ldown,l)/tau_graddiv(l) + due_diss2(ij+u_ldown,l)/tau_gradrot(l))*itau_dissip
522
523            dtheta_rhodz(ij,l,1) = -0.5*dtheta_rhodz_diss(ij,l)/tau_divgrad(l)*itau_dissip
524        ENDDO
525      ENDDO
526
527!      dtheta_rhodz=0
528!      due=0
529
530      IF(rayleigh_friction_type>0) THEN
531       IF(rayleigh_friction_type<99) THEN
532         phi=f_geopot(ind)
533         ue=f_ue(ind)
534         DO l=ll_begin,ll_end
535            DO ij=ij_begin,ij_end
536               CALL relax(t_right, u_right)
537               CALL relax(t_lup,   u_lup)
538               CALL relax(t_ldown, u_ldown)
539            ENDDO
540         END DO
541       ELSE
542         ue=f_ue(ind)
543            DO ij=ij_begin,ij_end
544              nn = ij+u_right
545              IF (ABS(lat_e(nn)) .gt. rayleigh_limlat) THEN
546                !print*, "latitude", lat_e(nn)*180./3.14159
547                due(nn,ll_begin:ll_begin+1) = due(nn,ll_begin:ll_begin+1) - (ue(nn,ll_begin:ll_begin+1)/rayleigh_tau)
548              ENDIF
549              nn = ij+u_lup
550              IF (ABS(lat_e(nn)) .gt. rayleigh_limlat) THEN
551                due(nn,ll_begin:ll_begin+1) = due(nn,ll_begin:ll_begin+1) - (ue(nn,ll_begin:ll_begin+1)/rayleigh_tau)
552              ENDIF
553              nn = ij+u_ldown
554              IF (ABS(lat_e(nn)) .gt. rayleigh_limlat) THEN
555                due(nn,ll_begin:ll_begin+1) = due(nn,ll_begin:ll_begin+1) - (ue(nn,ll_begin:ll_begin+1)/rayleigh_tau)
556              ENDIF
557            ENDDO
558       ENDIF
559      END IF
560   END DO
561
562   CALL trace_end("dissip")
563
564   CALL write_dissip_tendencies
565!$OMP BARRIER
566
567    CONTAINS
568
569      SUBROUTINE relax(shift_t, shift_u)
570        USE dcmip_initial_conditions_test_1_2_3
571        REAL(rstd) :: z, ulon,ulat, lon,lat, & ! input to test2_schaer_mountain
572             p,hyam,hybm,w,t,phis,ps,rho,q, &   ! unused input/output to test2_schaer_mountain
573             fz, u3d(3), uref
574        REAL(rstd), PARAMETER :: zh=2e4,ztop=3e4  ! DCMIP values
575        LOGICAL :: hybrid_eta
576        INTEGER :: shift_u, shift_t, zcoords, nn
577        z = (phi(ij,l)+phi(ij+shift_t,l))/(2.*g)
578        IF(z>zh) THEN  ! relax only in the sponge layer z>zh
579
580           nn = ij+shift_u
581           zcoords = 1 ; hybrid_eta = .FALSE. ! use z instead of p or hyam/hybm
582           CALL test2_schaer_mountain(lon_e(nn),lat_e(nn),p,z,zcoords,hybrid_eta, &
583                hyam,hybm,shear,ulon,ulat,w,t,phis,ps,rho,q)
584!           u3d = ulon*elon_e(nn,:) + ulat*elat_e(nn,:)
585           u3d = ulon*elon_e(nn,:) ! ulat=0
586           uref = sum(u3d*ep_e(nn,:))
587
588           fz = sin((pi/2)*(z-zh)/(ztop-zh))
589           fz = fz*fz/rayleigh_tau
590           due(nn,l) = due(nn,l) - itau_dissip*fz*(ue(nn,l)-uref)
591!           fz = 1./rayleigh_tau
592!           due(nn,l) = due(nn,l) - fz*ue(nn,l)
593         END IF
594       END SUBROUTINE relax
595     
596       SUBROUTINE write_dissip_tendencies
597         USE observable_mod, ONLY : f_buf_ulon, f_buf_ulat
598         USE wind_mod
599         USE output_field_mod
600
601         CALL transfert_request(f_due_diss1,req_e1_vect)
602         CALL un2ulonlat(f_due_diss1, f_buf_ulon, f_buf_ulat, (1./(tau_graddiv(1))))
603         CALL output_field("dulon_diss1",f_buf_ulon)
604         CALL output_field("dulat_diss1",f_buf_ulat)
605!
606         CALL transfert_request(f_due_diss2,req_e1_vect)
607         CALL un2ulonlat(f_due_diss2, f_buf_ulon, f_buf_ulat, (1./(tau_graddiv(1))))
608         CALL output_field("dulon_diss2",f_buf_ulon)
609         CALL output_field("dulat_diss2",f_buf_ulat)
610       END SUBROUTINE write_dissip_tendencies
611
612  END SUBROUTINE dissip
613
614  SUBROUTINE gradiv(f_ue,f_due)
615  USE icosa
616  USE trace
617  USE omp_para
618  IMPLICIT NONE
619    TYPE(t_field),POINTER :: f_ue(:)
620    TYPE(t_field),POINTER :: f_due(:)
621    REAL(rstd),POINTER    :: ue(:,:)
622    REAL(rstd),POINTER    :: due(:,:)
623    INTEGER :: ind
624    INTEGER :: it,l,ij
625       
626    CALL trace_start("gradiv")
627
628    DO ind=1,ndomain
629      IF (.NOT. assigned_domain(ind)) CYCLE
630      CALL swap_dimensions(ind)
631      CALL swap_geometry(ind)
632      ue=f_ue(ind)
633      due=f_due(ind) 
634      DO  l = ll_begin, ll_end
635!DIR$ SIMD
636        DO ij=ij_begin,ij_end
637             due(ij+u_right,l)=ue(ij+u_right,l)
638             due(ij+u_lup,l)=ue(ij+u_lup,l)
639             due(ij+u_ldown,l)=ue(ij+u_ldown,l)
640        ENDDO
641      ENDDO
642    ENDDO
643
644    DO it=1,nitergdiv
645       
646      CALL send_message(f_due,req_due)
647      CALL wait_message(req_due)
648       
649      DO ind=1,ndomain
650        IF (.NOT. assigned_domain(ind)) CYCLE
651        CALL swap_dimensions(ind)
652        CALL swap_geometry(ind)
653        due=f_due(ind) 
654        CALL compute_gradiv(due,due,ll_begin,ll_end)
655      ENDDO
656    ENDDO
657
658   CALL trace_end("gradiv")
659
660  END SUBROUTINE gradiv
661 
662
663  SUBROUTINE gradrot(f_ue,f_due)
664  USE icosa
665  USE trace
666  USE omp_para
667  IMPLICIT NONE
668    TYPE(t_field),POINTER :: f_ue(:)
669    TYPE(t_field),POINTER :: f_due(:)
670    REAL(rstd),POINTER    :: ue(:,:)
671    REAL(rstd),POINTER    :: due(:,:)
672    INTEGER :: ind
673    INTEGER :: it,l,ij
674       
675    CALL trace_start("gradrot")
676
677    DO ind=1,ndomain
678      IF (.NOT. assigned_domain(ind)) CYCLE
679      CALL swap_dimensions(ind)
680      CALL swap_geometry(ind)
681      ue=f_ue(ind)
682      due=f_due(ind) 
683      DO  l = ll_begin, ll_end
684!DIR$ SIMD
685        DO ij=ij_begin,ij_end
686             due(ij+u_right,l)=ue(ij+u_right,l)
687             due(ij+u_lup,l)=ue(ij+u_lup,l)
688             due(ij+u_ldown,l)=ue(ij+u_ldown,l)
689        ENDDO
690      ENDDO
691    ENDDO
692
693    DO it=1,nitergrot
694       
695      CALL send_message(f_due,req_due)
696      CALL wait_message(req_due)
697       
698      DO ind=1,ndomain
699        IF (.NOT. assigned_domain(ind)) CYCLE
700        CALL swap_dimensions(ind)
701        CALL swap_geometry(ind)
702        due=f_due(ind) 
703        CALL compute_gradrot(due,due,ll_begin,ll_end)
704      ENDDO
705
706    ENDDO
707
708    CALL trace_end("gradrot")
709
710  END SUBROUTINE gradrot
711 
712  SUBROUTINE divgrad(f_theta,f_dtheta)
713  USE icosa
714  USE trace
715  USE omp_para
716  IMPLICIT NONE
717    TYPE(t_field),POINTER :: f_theta(:)
718    TYPE(t_field),POINTER :: f_dtheta(:)
719    REAL(rstd),POINTER    :: theta(:,:)
720    REAL(rstd),POINTER    :: dtheta(:,:)
721    INTEGER :: ind
722    INTEGER :: it
723
724    CALL trace_start("divgrad")
725       
726    DO ind=1,ndomain
727      IF (.NOT. assigned_domain(ind)) CYCLE
728      CALL swap_dimensions(ind)
729      CALL swap_geometry(ind)
730      theta=f_theta(ind)
731      dtheta=f_dtheta(ind) 
732      dtheta=theta
733    ENDDO
734
735    DO it=1,niterdivgrad
736       
737      CALL transfert_request(f_dtheta,req_i1)
738       
739      DO ind=1,ndomain
740        IF (.NOT. assigned_domain(ind)) CYCLE
741        CALL swap_dimensions(ind)
742        CALL swap_geometry(ind)
743        dtheta=f_dtheta(ind) 
744        CALL compute_divgrad(dtheta,dtheta,ll_begin,ll_end)
745      ENDDO
746
747    ENDDO
748
749    CALL trace_end("divgrad")
750
751  END SUBROUTINE divgrad
752   
753  SUBROUTINE divgrad_theta_rhodz(f_mass,f_theta_rhodz,f_dtheta_rhodz)
754  USE icosa
755  USE trace
756  USE omp_para
757  IMPLICIT NONE
758    TYPE(t_field),POINTER :: f_mass(:)
759    TYPE(t_field),POINTER :: f_theta_rhodz(:)
760    TYPE(t_field),POINTER :: f_dtheta_rhodz(:)
761   
762    REAL(rstd),POINTER :: mass(:,:)
763    REAL(rstd),POINTER :: theta_rhodz(:,:,:)
764    REAL(rstd),POINTER :: dtheta_rhodz(:,:)
765
766    INTEGER :: ind
767    INTEGER :: it,l,ij
768
769    CALL trace_start("divgrad")
770
771    DO ind=1,ndomain
772      IF (.NOT. assigned_domain(ind)) CYCLE
773      CALL swap_dimensions(ind)
774      CALL swap_geometry(ind)
775      mass=f_mass(ind)
776      theta_rhodz=f_theta_rhodz(ind)
777      dtheta_rhodz=f_dtheta_rhodz(ind) 
778      DO  l = ll_begin, ll_end
779!DIR$ SIMD
780        DO ij=ij_begin,ij_end
781            dtheta_rhodz(ij,l) = theta_rhodz(ij,l,1) / mass(ij,l)
782        ENDDO
783      ENDDO
784    ENDDO
785
786    DO it=1,niterdivgrad
787       
788      CALL send_message(f_dtheta_rhodz,req_dtheta)
789      CALL wait_message(req_dtheta)
790      DO ind=1,ndomain
791        IF (.NOT. assigned_domain(ind)) CYCLE
792        CALL swap_dimensions(ind)
793        CALL swap_geometry(ind)
794        dtheta_rhodz=f_dtheta_rhodz(ind) 
795        CALL compute_divgrad(dtheta_rhodz,dtheta_rhodz,ll_begin,ll_end)
796      ENDDO
797
798    ENDDO
799
800    DO ind=1,ndomain
801      IF (.NOT. assigned_domain(ind)) CYCLE
802      CALL swap_dimensions(ind)
803      CALL swap_geometry(ind)
804      dtheta_rhodz=f_dtheta_rhodz(ind) 
805      mass=f_mass(ind)
806   
807      DO  l = ll_begin, ll_end
808!DIR$ SIMD
809        DO ij=ij_begin,ij_end
810            dtheta_rhodz(ij,l) = dtheta_rhodz(ij,l) * mass(ij,l)
811        ENDDO
812      ENDDO
813    ENDDO
814
815
816    CALL trace_end("divgrad")
817
818  END SUBROUTINE divgrad_theta_rhodz 
819 
820  SUBROUTINE compute_gradiv(ue,gradivu_e,llb,lle)
821  USE icosa
822  IMPLICIT NONE
823    INTEGER,INTENT(IN)     :: llb
824    INTEGER,INTENT(IN)     :: lle
825    REAL(rstd),INTENT(IN)  :: ue(iim*3*jjm,llm)
826    REAL(rstd),INTENT(OUT) :: gradivu_e(iim*3*jjm,llm)
827    REAL(rstd) :: divu_i(iim*jjm,llb:lle)
828   
829    INTEGER :: ij,l
830
831    DO l=llb,lle
832!DIR$ SIMD
833      DO ij=ij_begin,ij_end
834          divu_i(ij,l)=1./Ai(ij)*(ne(ij,right)*ue(ij+u_right,l)*le(ij+u_right)  +  &
835                             ne(ij,rup)*ue(ij+u_rup,l)*le(ij+u_rup)        +  & 
836                             ne(ij,lup)*ue(ij+u_lup,l)*le(ij+u_lup)        +  & 
837                             ne(ij,left)*ue(ij+u_left,l)*le(ij+u_left)     +  & 
838                             ne(ij,ldown)*ue(ij+u_ldown,l)*le(ij+u_ldown)  +  & 
839                             ne(ij,rdown)*ue(ij+u_rdown,l)*le(ij+u_rdown))
840      ENDDO
841    ENDDO
842   
843    DO l=llb,lle
844!DIR$ SIMD
845      DO ij=ij_begin,ij_end
846 
847          gradivu_e(ij+u_right,l)=-1/de(ij+u_right)*(ne(ij,right)*divu_i(ij,l)+ ne(ij+t_right,left)*divu_i(ij+t_right,l) )       
848
849          gradivu_e(ij+u_lup,l)=-1/de(ij+u_lup)*(ne(ij,lup)*divu_i(ij,l)+ ne(ij+t_lup,rdown)*divu_i(ij+t_lup,l))       
850   
851          gradivu_e(ij+u_ldown,l)=-1/de(ij+u_ldown)*(ne(ij,ldown)*divu_i(ij,l)+ne(ij+t_ldown,rup)*divu_i(ij+t_ldown,l) )
852
853      ENDDO
854    ENDDO
855
856    DO l=llb,lle
857!DIR$ SIMD
858      DO ij=ij_begin,ij_end
859          gradivu_e(ij+u_right,l)=-gradivu_e(ij+u_right,l)*cgraddiv
860          gradivu_e(ij+u_lup,l)=-gradivu_e(ij+u_lup,l)*cgraddiv
861          gradivu_e(ij+u_ldown,l)=-gradivu_e(ij+u_ldown,l)*cgraddiv
862      ENDDO
863    ENDDO
864
865   
866  END SUBROUTINE compute_gradiv
867 
868  SUBROUTINE compute_divgrad(theta,divgrad_i,llb,lle)
869  USE icosa
870  IMPLICIT NONE
871    INTEGER,INTENT(IN)     :: llb
872    INTEGER,INTENT(IN)     :: lle
873    REAL(rstd),INTENT(IN)  :: theta(iim*jjm,llm)
874    REAL(rstd),INTENT(OUT) :: divgrad_i(iim*jjm,llm)
875    REAL(rstd)  :: grad_e(3*iim*jjm,llb:lle)
876
877    INTEGER :: ij,l
878
879       
880    DO l=llb,lle
881!DIR$ SIMD
882      DO ij=ij_begin_ext,ij_end_ext
883 
884          grad_e(ij+u_right,l)=-1/de(ij+u_right)*(ne(ij,right)*theta(ij,l)+ ne(ij+t_right,left)*theta(ij+t_right,l) )       
885
886          grad_e(ij+u_lup,l)=-1/de(ij+u_lup)*(ne(ij,lup)*theta(ij,l)+ ne(ij+t_lup,rdown)*theta(ij+t_lup,l ))       
887   
888          grad_e(ij+u_ldown,l)=-1/de(ij+u_ldown)*(ne(ij,ldown)*theta(ij,l)+ne(ij+t_ldown,rup)*theta(ij+t_ldown,l) )
889
890      ENDDO
891    ENDDO
892   
893   
894    DO l=llb,lle
895!DIR$ SIMD
896      DO ij=ij_begin,ij_end
897
898          divgrad_i(ij,l)=1./Ai(ij)*(ne(ij,right)*grad_e(ij+u_right,l)*le(ij+u_right)  +  &
899                             ne(ij,rup)*grad_e(ij+u_rup,l)*le(ij+u_rup)              +  & 
900                             ne(ij,lup)*grad_e(ij+u_lup,l)*le(ij+u_lup)              +  & 
901                             ne(ij,left)*grad_e(ij+u_left,l)*le(ij+u_left)           +  & 
902                             ne(ij,ldown)*grad_e(ij+u_ldown,l)*le(ij+u_ldown)        +  & 
903                             ne(ij,rdown)*grad_e(ij+u_rdown,l)*le(ij+u_rdown))
904      ENDDO
905    ENDDO
906   
907    DO l=llb,lle
908      DO ij=ij_begin,ij_end
909          divgrad_i(ij,l)=-divgrad_i(ij,l)*cdivgrad
910      ENDDO
911    ENDDO
912
913  END SUBROUTINE compute_divgrad
914
915   
916  SUBROUTINE compute_gradrot(ue,gradrot_e,llb,lle)
917  USE icosa
918  IMPLICIT NONE
919    INTEGER,INTENT(IN)     :: llb
920    INTEGER,INTENT(IN)     :: lle
921    REAL(rstd),INTENT(IN)  :: ue(iim*3*jjm,llm)
922    REAL(rstd),INTENT(OUT) :: gradrot_e(iim*3*jjm,llm)
923    REAL(rstd) :: rot_v(2*iim*jjm,llb:lle)
924
925    INTEGER :: ij,l
926     
927    DO l=llb,lle
928!DIR$ SIMD
929      DO ij=ij_begin_ext,ij_end_ext
930       
931          rot_v(ij+z_up,l)= 1./Av(ij+z_up)*(  ne(ij,rup)*ue(ij+u_rup,l)*de(ij+u_rup)                     &
932                                + ne(ij+t_rup,left)*ue(ij+t_rup+u_left,l)*de(ij+t_rup+u_left)          &
933                                - ne(ij,lup)*ue(ij+u_lup,l)*de(ij+u_lup) ) 
934                             
935          rot_v(ij+z_down,l) = 1./Av(ij+z_down)*( ne(ij,ldown)*ue(ij+u_ldown,l)*de(ij+u_ldown)                 &
936                                     + ne(ij+t_ldown,right)*ue(ij+t_ldown+u_right,l)*de(ij+t_ldown+u_right)  &
937                                     - ne(ij,rdown)*ue(ij+u_rdown,l)*de(ij+u_rdown) )
938 
939      ENDDO
940    ENDDO                             
941 
942    DO l=llb,lle
943!DIR$ SIMD
944      DO ij=ij_begin,ij_end
945 
946          gradrot_e(ij+u_right,l)=1/le(ij+u_right)*ne(ij,right)*(rot_v(ij+z_rdown,l)-rot_v(ij+z_rup,l)) 
947         
948          gradrot_e(ij+u_lup,l)=1/le(ij+u_lup)*ne(ij,lup)*(rot_v(ij+z_up,l)-rot_v(ij+z_lup,l)) 
949       
950          gradrot_e(ij+u_ldown,l)=1/le(ij+u_ldown)*ne(ij,ldown)*(rot_v(ij+z_ldown,l)-rot_v(ij+z_down,l))
951       
952      ENDDO
953    ENDDO
954
955    DO l=llb,lle
956!DIR$ SIMD
957      DO ij=ij_begin,ij_end
958   
959          gradrot_e(ij+u_right,l)=-gradrot_e(ij+u_right,l)*cgradrot       
960          gradrot_e(ij+u_lup,l)=-gradrot_e(ij+u_lup,l)*cgradrot
961          gradrot_e(ij+u_ldown,l)=-gradrot_e(ij+u_ldown,l)*cgradrot
962       
963      ENDDO
964    ENDDO 
965   
966  END SUBROUTINE compute_gradrot
967
968
969END MODULE dissip_gcm_mod
970       
Note: See TracBrowser for help on using the repository browser.