source: codes/icosagcm/branches/SATURN_DYNAMICO/ICOSAGCM/src/dissip_gcm.f90 @ 316

Last change on this file since 316 was 298, checked in by milmd, 10 years ago

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