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

Last change on this file since 250 was 250, checked in by dubos, 10 years ago

Ignore dissipation if tau_XX all 0
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.