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

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

Merging OpenMP parallisme mode : by subdomain and on vertical level.
This feature is actually experimental but may be retro-compatible with the last method based only on subdomain

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