source: codes/icosagcm/trunk/src/dissip_gcm.f90 @ 532

Last change on this file since 532 was 525, checked in by dubos, 7 years ago

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