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

Last change on this file since 186 was 186, checked in by ymipsl, 10 years ago

Add new openMP parallelism based on distribution of domains on threads. There is no more limitation of number of threads by MPI process.

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