source: codes/icosagcm/devel/src/dissip/dissip_gcm.f90 @ 550

Last change on this file since 550 was 550, checked in by dubos, 7 years ago

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