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

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

Synchronize trunk and Saturn branch.
Merge modification from trunk branch to Saturn branch.

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