source: codes/icosagcm/trunk/src/dissip/dissip_gcm.f90 @ 648

Last change on this file since 648 was 648, checked in by aslmd, 6 years ago

added an option for simple Rayleigh friction at bottom (giant planets, but could be used for any case really)
File size: 27.2 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       itau_dissip=INT(mintau/dt)
456       dtdissip=itau_dissip*dt
457    ELSE
458       IF (is_master) PRINT *,"No dissipation time set, setting itau_dissip to 1000000000"
459       itau_dissip=100000000
460    END IF
461    itau_dissip=MAX(1,itau_dissip)
462    IF (is_master) PRINT *,"rayleigh_tau",rayleigh_tau, "mintau ",mintau, &
463         "itau_dissip",itau_dissip," dtdissip ",dtdissip
464
465  END SUBROUTINE init_dissip 
466 
467 
468  SUBROUTINE dissip(f_ps,f_mass,f_phis,f_geopot,f_theta_rhodz,f_ue, f_dtheta_rhodz,f_due)
469  USE icosa
470  USE theta2theta_rhodz_mod
471  USE pression_mod
472  USE exner_mod
473  USE geopotential_mod
474  USE trace
475  USE time_mod
476  USE omp_para
477  IMPLICIT NONE
478    TYPE(t_field),POINTER :: f_ps(:), f_mass(:), f_phis(:), f_geopot(:)
479    TYPE(t_field),POINTER :: f_theta_rhodz(:), f_dtheta_rhodz(:)
480    TYPE(t_field),POINTER :: f_ue(:), f_due(:)
481
482    REAL(rstd),POINTER         :: due(:,:)
483    REAL(rstd),POINTER         :: phi(:,:), ue(:,:)
484    REAL(rstd),POINTER         :: due_diss1(:,:)
485    REAL(rstd),POINTER         :: due_diss2(:,:)
486    REAL(rstd),POINTER         :: dtheta_rhodz(:,:,:)
487    REAL(rstd),POINTER         :: dtheta_rhodz_diss(:,:)
488
489    INTEGER :: ind, shear
490    INTEGER :: l,ij,nn
491
492!$OMP BARRIER
493   
494    CALL trace_start("dissip")
495    CALL gradiv(f_ue,f_due_diss1)
496    CALL gradrot(f_ue,f_due_diss2)
497
498    CALL divgrad_theta_rhodz(f_mass,f_theta_rhodz,f_dtheta_rhodz_diss)
499   
500    DO ind=1,ndomain
501      IF (.NOT. assigned_domain(ind)) CYCLE
502      CALL swap_dimensions(ind)
503      CALL swap_geometry(ind)
504      due=f_due(ind) 
505      due_diss1=f_due_diss1(ind)
506      due_diss2=f_due_diss2(ind)
507      dtheta_rhodz=f_dtheta_rhodz(ind)
508      dtheta_rhodz_diss=f_dtheta_rhodz_diss(ind)
509
510      DO l=ll_begin,ll_end
511!DIR$ SIMD
512        DO ij=ij_begin,ij_end
513
514            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 
515            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
516            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
517
518            dtheta_rhodz(ij,l,1) = -0.5*dtheta_rhodz_diss(ij,l)/tau_divgrad(l)*itau_dissip
519        ENDDO
520      ENDDO
521
522!      dtheta_rhodz=0
523!      due=0
524
525      IF(rayleigh_friction_type>0) THEN
526       IF(rayleigh_friction_type<99) THEN
527         phi=f_geopot(ind)
528         ue=f_ue(ind)
529         DO l=ll_begin,ll_end
530            DO ij=ij_begin,ij_end
531               CALL relax(t_right, u_right)
532               CALL relax(t_lup,   u_lup)
533               CALL relax(t_ldown, u_ldown)
534            ENDDO
535         END DO
536       ELSE
537         ue=f_ue(ind)
538            DO ij=ij_begin,ij_end
539              nn = ij+u_right
540              IF (ABS(lat_e(nn)) .gt. rayleigh_limlat) THEN
541                !print*, "latitude", lat_e(nn)*180./3.14159
542                due(nn,ll_begin:ll_begin+1) = due(nn,ll_begin:ll_begin+1) - (ue(nn,ll_begin:ll_begin+1)/rayleigh_tau)
543              ENDIF
544              nn = ij+u_lup
545              IF (ABS(lat_e(nn)) .gt. rayleigh_limlat) THEN
546                due(nn,ll_begin:ll_begin+1) = due(nn,ll_begin:ll_begin+1) - (ue(nn,ll_begin:ll_begin+1)/rayleigh_tau)
547              ENDIF
548              nn = ij+u_ldown
549              IF (ABS(lat_e(nn)) .gt. rayleigh_limlat) THEN
550                due(nn,ll_begin:ll_begin+1) = due(nn,ll_begin:ll_begin+1) - (ue(nn,ll_begin:ll_begin+1)/rayleigh_tau)
551              ENDIF
552            ENDDO
553       ENDIF
554      END IF
555   END DO
556
557   CALL trace_end("dissip")
558
559!$OMP BARRIER
560
561    CONTAINS
562
563      SUBROUTINE relax(shift_t, shift_u)
564        USE dcmip_initial_conditions_test_1_2_3
565        REAL(rstd) :: z, ulon,ulat, lon,lat, & ! input to test2_schaer_mountain
566             p,hyam,hybm,w,t,phis,ps,rho,q, &   ! unused input/output to test2_schaer_mountain
567             fz, u3d(3), uref
568        REAL(rstd), PARAMETER :: zh=2e4,ztop=3e4  ! DCMIP values
569        LOGICAL :: hybrid_eta
570        INTEGER :: shift_u, shift_t, zcoords, nn
571        z = (phi(ij,l)+phi(ij+shift_t,l))/(2.*g)
572        IF(z>zh) THEN  ! relax only in the sponge layer z>zh
573
574           nn = ij+shift_u
575           zcoords = 1 ; hybrid_eta = .FALSE. ! use z instead of p or hyam/hybm
576           CALL test2_schaer_mountain(lon_e(nn),lat_e(nn),p,z,zcoords,hybrid_eta, &
577                hyam,hybm,shear,ulon,ulat,w,t,phis,ps,rho,q)
578!           u3d = ulon*elon_e(nn,:) + ulat*elat_e(nn,:)
579           u3d = ulon*elon_e(nn,:) ! ulat=0
580           uref = sum(u3d*ep_e(nn,:))
581
582           fz = sin((pi/2)*(z-zh)/(ztop-zh))
583           fz = fz*fz/rayleigh_tau
584           due(nn,l) = due(nn,l) - itau_dissip*fz*(ue(nn,l)-uref)
585!           fz = 1./rayleigh_tau
586!           due(nn,l) = due(nn,l) - fz*ue(nn,l)
587         END IF
588       END SUBROUTINE relax
589     
590  END SUBROUTINE dissip
591
592  SUBROUTINE gradiv(f_ue,f_due)
593  USE icosa
594  USE trace
595  USE omp_para
596  IMPLICIT NONE
597    TYPE(t_field),POINTER :: f_ue(:)
598    TYPE(t_field),POINTER :: f_due(:)
599    REAL(rstd),POINTER    :: ue(:,:)
600    REAL(rstd),POINTER    :: due(:,:)
601    INTEGER :: ind
602    INTEGER :: it,l,ij
603       
604    CALL trace_start("gradiv")
605
606    DO ind=1,ndomain
607      IF (.NOT. assigned_domain(ind)) CYCLE
608      CALL swap_dimensions(ind)
609      CALL swap_geometry(ind)
610      ue=f_ue(ind)
611      due=f_due(ind) 
612      DO  l = ll_begin, ll_end
613!DIR$ SIMD
614        DO ij=ij_begin,ij_end
615             due(ij+u_right,l)=ue(ij+u_right,l)
616             due(ij+u_lup,l)=ue(ij+u_lup,l)
617             due(ij+u_ldown,l)=ue(ij+u_ldown,l)
618        ENDDO
619      ENDDO
620    ENDDO
621
622    DO it=1,nitergdiv
623       
624      CALL send_message(f_due,req_due)
625      CALL wait_message(req_due)
626       
627      DO ind=1,ndomain
628        IF (.NOT. assigned_domain(ind)) CYCLE
629        CALL swap_dimensions(ind)
630        CALL swap_geometry(ind)
631        due=f_due(ind) 
632        CALL compute_gradiv(due,due,ll_begin,ll_end)
633      ENDDO
634    ENDDO
635
636   CALL trace_end("gradiv")
637
638  END SUBROUTINE gradiv
639 
640
641  SUBROUTINE gradrot(f_ue,f_due)
642  USE icosa
643  USE trace
644  USE omp_para
645  IMPLICIT NONE
646    TYPE(t_field),POINTER :: f_ue(:)
647    TYPE(t_field),POINTER :: f_due(:)
648    REAL(rstd),POINTER    :: ue(:,:)
649    REAL(rstd),POINTER    :: due(:,:)
650    INTEGER :: ind
651    INTEGER :: it,l,ij
652       
653    CALL trace_start("gradrot")
654
655    DO ind=1,ndomain
656      IF (.NOT. assigned_domain(ind)) CYCLE
657      CALL swap_dimensions(ind)
658      CALL swap_geometry(ind)
659      ue=f_ue(ind)
660      due=f_due(ind) 
661      DO  l = ll_begin, ll_end
662!DIR$ SIMD
663        DO ij=ij_begin,ij_end
664             due(ij+u_right,l)=ue(ij+u_right,l)
665             due(ij+u_lup,l)=ue(ij+u_lup,l)
666             due(ij+u_ldown,l)=ue(ij+u_ldown,l)
667        ENDDO
668      ENDDO
669    ENDDO
670
671    DO it=1,nitergrot
672       
673      CALL send_message(f_due,req_due)
674      CALL wait_message(req_due)
675       
676      DO ind=1,ndomain
677        IF (.NOT. assigned_domain(ind)) CYCLE
678        CALL swap_dimensions(ind)
679        CALL swap_geometry(ind)
680        due=f_due(ind) 
681        CALL compute_gradrot(due,due,ll_begin,ll_end)
682      ENDDO
683
684    ENDDO
685
686    CALL trace_end("gradrot")
687
688  END SUBROUTINE gradrot
689 
690  SUBROUTINE divgrad(f_theta,f_dtheta)
691  USE icosa
692  USE trace
693  USE omp_para
694  IMPLICIT NONE
695    TYPE(t_field),POINTER :: f_theta(:)
696    TYPE(t_field),POINTER :: f_dtheta(:)
697    REAL(rstd),POINTER    :: theta(:,:)
698    REAL(rstd),POINTER    :: dtheta(:,:)
699    INTEGER :: ind
700    INTEGER :: it
701
702    CALL trace_start("divgrad")
703       
704    DO ind=1,ndomain
705      IF (.NOT. assigned_domain(ind)) CYCLE
706      CALL swap_dimensions(ind)
707      CALL swap_geometry(ind)
708      theta=f_theta(ind)
709      dtheta=f_dtheta(ind) 
710      dtheta=theta
711    ENDDO
712
713    DO it=1,niterdivgrad
714       
715      CALL transfert_request(f_dtheta,req_i1)
716       
717      DO ind=1,ndomain
718        IF (.NOT. assigned_domain(ind)) CYCLE
719        CALL swap_dimensions(ind)
720        CALL swap_geometry(ind)
721        dtheta=f_dtheta(ind) 
722        CALL compute_divgrad(dtheta,dtheta,ll_begin,ll_end)
723      ENDDO
724
725    ENDDO
726
727    CALL trace_end("divgrad")
728
729  END SUBROUTINE divgrad
730   
731  SUBROUTINE divgrad_theta_rhodz(f_mass,f_theta_rhodz,f_dtheta_rhodz)
732  USE icosa
733  USE trace
734  USE omp_para
735  IMPLICIT NONE
736    TYPE(t_field),POINTER :: f_mass(:)
737    TYPE(t_field),POINTER :: f_theta_rhodz(:)
738    TYPE(t_field),POINTER :: f_dtheta_rhodz(:)
739   
740    REAL(rstd),POINTER :: mass(:,:)
741    REAL(rstd),POINTER :: theta_rhodz(:,:,:)
742    REAL(rstd),POINTER :: dtheta_rhodz(:,:)
743
744    INTEGER :: ind
745    INTEGER :: it,l,ij
746
747    CALL trace_start("divgrad")
748
749    DO ind=1,ndomain
750      IF (.NOT. assigned_domain(ind)) CYCLE
751      CALL swap_dimensions(ind)
752      CALL swap_geometry(ind)
753      mass=f_mass(ind)
754      theta_rhodz=f_theta_rhodz(ind)
755      dtheta_rhodz=f_dtheta_rhodz(ind) 
756      DO  l = ll_begin, ll_end
757!DIR$ SIMD
758        DO ij=ij_begin,ij_end
759            dtheta_rhodz(ij,l) = theta_rhodz(ij,l,1) / mass(ij,l)
760        ENDDO
761      ENDDO
762    ENDDO
763
764    DO it=1,niterdivgrad
765       
766      CALL send_message(f_dtheta_rhodz,req_dtheta)
767      CALL wait_message(req_dtheta)
768      DO ind=1,ndomain
769        IF (.NOT. assigned_domain(ind)) CYCLE
770        CALL swap_dimensions(ind)
771        CALL swap_geometry(ind)
772        dtheta_rhodz=f_dtheta_rhodz(ind) 
773        CALL compute_divgrad(dtheta_rhodz,dtheta_rhodz,ll_begin,ll_end)
774      ENDDO
775
776    ENDDO
777
778    DO ind=1,ndomain
779      IF (.NOT. assigned_domain(ind)) CYCLE
780      CALL swap_dimensions(ind)
781      CALL swap_geometry(ind)
782      dtheta_rhodz=f_dtheta_rhodz(ind) 
783      mass=f_mass(ind)
784   
785      DO  l = ll_begin, ll_end
786!DIR$ SIMD
787        DO ij=ij_begin,ij_end
788            dtheta_rhodz(ij,l) = dtheta_rhodz(ij,l) * mass(ij,l)
789        ENDDO
790      ENDDO
791    ENDDO
792
793
794    CALL trace_end("divgrad")
795
796  END SUBROUTINE divgrad_theta_rhodz 
797 
798  SUBROUTINE compute_gradiv(ue,gradivu_e,llb,lle)
799  USE icosa
800  IMPLICIT NONE
801    INTEGER,INTENT(IN)     :: llb
802    INTEGER,INTENT(IN)     :: lle
803    REAL(rstd),INTENT(IN)  :: ue(iim*3*jjm,llm)
804    REAL(rstd),INTENT(OUT) :: gradivu_e(iim*3*jjm,llm)
805    REAL(rstd) :: divu_i(iim*jjm,llb:lle)
806   
807    INTEGER :: ij,l
808
809    DO l=llb,lle
810!DIR$ SIMD
811      DO ij=ij_begin,ij_end
812          divu_i(ij,l)=1./Ai(ij)*(ne(ij,right)*ue(ij+u_right,l)*le(ij+u_right)  +  &
813                             ne(ij,rup)*ue(ij+u_rup,l)*le(ij+u_rup)        +  & 
814                             ne(ij,lup)*ue(ij+u_lup,l)*le(ij+u_lup)        +  & 
815                             ne(ij,left)*ue(ij+u_left,l)*le(ij+u_left)     +  & 
816                             ne(ij,ldown)*ue(ij+u_ldown,l)*le(ij+u_ldown)  +  & 
817                             ne(ij,rdown)*ue(ij+u_rdown,l)*le(ij+u_rdown))
818      ENDDO
819    ENDDO
820   
821    DO l=llb,lle
822!DIR$ SIMD
823      DO ij=ij_begin,ij_end
824 
825          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) )       
826
827          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))       
828   
829          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) )
830
831      ENDDO
832    ENDDO
833
834    DO l=llb,lle
835!DIR$ SIMD
836      DO ij=ij_begin,ij_end
837          gradivu_e(ij+u_right,l)=-gradivu_e(ij+u_right,l)*cgraddiv
838          gradivu_e(ij+u_lup,l)=-gradivu_e(ij+u_lup,l)*cgraddiv
839          gradivu_e(ij+u_ldown,l)=-gradivu_e(ij+u_ldown,l)*cgraddiv
840      ENDDO
841    ENDDO
842
843   
844  END SUBROUTINE compute_gradiv
845 
846  SUBROUTINE compute_divgrad(theta,divgrad_i,llb,lle)
847  USE icosa
848  IMPLICIT NONE
849    INTEGER,INTENT(IN)     :: llb
850    INTEGER,INTENT(IN)     :: lle
851    REAL(rstd),INTENT(IN)  :: theta(iim*jjm,llm)
852    REAL(rstd),INTENT(OUT) :: divgrad_i(iim*jjm,llm)
853    REAL(rstd)  :: grad_e(3*iim*jjm,llb:lle)
854
855    INTEGER :: ij,l
856
857       
858    DO l=llb,lle
859!DIR$ SIMD
860      DO ij=ij_begin_ext,ij_end_ext
861 
862          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) )       
863
864          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 ))       
865   
866          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) )
867
868      ENDDO
869    ENDDO
870   
871   
872    DO l=llb,lle
873!DIR$ SIMD
874      DO ij=ij_begin,ij_end
875
876          divgrad_i(ij,l)=1./Ai(ij)*(ne(ij,right)*grad_e(ij+u_right,l)*le(ij+u_right)  +  &
877                             ne(ij,rup)*grad_e(ij+u_rup,l)*le(ij+u_rup)              +  & 
878                             ne(ij,lup)*grad_e(ij+u_lup,l)*le(ij+u_lup)              +  & 
879                             ne(ij,left)*grad_e(ij+u_left,l)*le(ij+u_left)           +  & 
880                             ne(ij,ldown)*grad_e(ij+u_ldown,l)*le(ij+u_ldown)        +  & 
881                             ne(ij,rdown)*grad_e(ij+u_rdown,l)*le(ij+u_rdown))
882      ENDDO
883    ENDDO
884   
885    DO l=llb,lle
886      DO ij=ij_begin,ij_end
887          divgrad_i(ij,l)=-divgrad_i(ij,l)*cdivgrad
888      ENDDO
889    ENDDO
890
891  END SUBROUTINE compute_divgrad
892
893   
894  SUBROUTINE compute_gradrot(ue,gradrot_e,llb,lle)
895  USE icosa
896  IMPLICIT NONE
897    INTEGER,INTENT(IN)     :: llb
898    INTEGER,INTENT(IN)     :: lle
899    REAL(rstd),INTENT(IN)  :: ue(iim*3*jjm,llm)
900    REAL(rstd),INTENT(OUT) :: gradrot_e(iim*3*jjm,llm)
901    REAL(rstd) :: rot_v(2*iim*jjm,llb:lle)
902
903    INTEGER :: ij,l
904     
905    DO l=llb,lle
906!DIR$ SIMD
907      DO ij=ij_begin_ext,ij_end_ext
908       
909          rot_v(ij+z_up,l)= 1./Av(ij+z_up)*(  ne(ij,rup)*ue(ij+u_rup,l)*de(ij+u_rup)                     &
910                                + ne(ij+t_rup,left)*ue(ij+t_rup+u_left,l)*de(ij+t_rup+u_left)          &
911                                - ne(ij,lup)*ue(ij+u_lup,l)*de(ij+u_lup) ) 
912                             
913          rot_v(ij+z_down,l) = 1./Av(ij+z_down)*( ne(ij,ldown)*ue(ij+u_ldown,l)*de(ij+u_ldown)                 &
914                                     + ne(ij+t_ldown,right)*ue(ij+t_ldown+u_right,l)*de(ij+t_ldown+u_right)  &
915                                     - ne(ij,rdown)*ue(ij+u_rdown,l)*de(ij+u_rdown) )
916 
917      ENDDO
918    ENDDO                             
919 
920    DO l=llb,lle
921!DIR$ SIMD
922      DO ij=ij_begin,ij_end
923 
924          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)) 
925         
926          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)) 
927       
928          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))
929       
930      ENDDO
931    ENDDO
932
933    DO l=llb,lle
934!DIR$ SIMD
935      DO ij=ij_begin,ij_end
936   
937          gradrot_e(ij+u_right,l)=-gradrot_e(ij+u_right,l)*cgradrot       
938          gradrot_e(ij+u_lup,l)=-gradrot_e(ij+u_lup,l)*cgradrot
939          gradrot_e(ij+u_ldown,l)=-gradrot_e(ij+u_ldown,l)*cgradrot
940       
941      ENDDO
942    ENDDO 
943   
944  END SUBROUTINE compute_gradrot
945
946
947END MODULE dissip_gcm_mod
948       
Note: See TracBrowser for help on using the repository browser.