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

Last change on this file since 360 was 360, checked in by dubos, 9 years ago

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