source: codes/icosagcm/trunk/src/dissip/dissip_gcm.F90 @ 954

Last change on this file since 954 was 954, checked in by adurocher, 5 years ago

trunk : Added metric terms to kernels parameters to avoid Host/GPU transferts

Metric terms are now subroutine parameters instead of module variables in kernel subroutines. Dummy arguments for metric terms are now defined as fixed-size arrays, and arrays dimensions are well known when entering an 'acc data' region. Array descriptors are no longer transferred form host to device each time the data region is executed.
File size: 35.3 KB
Line 
1MODULE dissip_gcm_mod
2  USE icosa
3  USE abort_mod
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, rayleigh_limlat
38!$OMP THREADPRIVATE(rayleigh_tau)
39!$OMP THREADPRIVATE(rayleigh_limlat)
40 
41  REAL,SAVE    :: dtdissip
42!$OMP THREADPRIVATE(dtdissip)
43 
44  PUBLIC init_dissip, dissip
45CONTAINS
46
47  SUBROUTINE allocate_dissip
48  USE icosa
49  IMPLICIT NONE 
50    CALL allocate_field(f_due_diss1,field_u,type_real,llm,ondevice=.TRUE.)
51    CALL allocate_field(f_due_diss2,field_u,type_real,llm,ondevice=.TRUE.)
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,ondevice=.TRUE.)
54    ALLOCATE(tau_graddiv(llm))
55    ALLOCATE(tau_gradrot(llm))   
56    ALLOCATE(tau_divgrad(llm))
57    !$acc enter data create(tau_graddiv(:),tau_gradrot(:),tau_divgrad(:)) async
58  END SUBROUTINE allocate_dissip
59 
60  SUBROUTINE init_dissip
61  USE icosa
62  USE disvert_mod
63  USE mpi_mod
64  USE mpipara
65  USE transfert_mod
66  USE time_mod
67  USE transfert_omp_mod
68  USE omp_para
69  USE abort_mod
70  IMPLICIT NONE
71 
72  TYPE(t_field),POINTER,SAVE  :: f_u(:)
73  TYPE(t_field),POINTER,SAVE  :: f_du(:)
74  REAL(rstd),POINTER     :: u(:)
75  REAL(rstd),POINTER     :: du(:)
76  TYPE(t_field),POINTER,SAVE  :: f_theta(:)
77  TYPE(t_field),POINTER ,SAVE :: f_dtheta(:)
78  REAL(rstd),POINTER    :: theta(:)
79  REAL(rstd),POINTER    :: dtheta(:)
80  REAL(rstd)            :: dumax,dumax1
81  REAL(rstd)            :: dthetamax,dthetamax1
82  REAL(rstd)            :: r
83  REAL(rstd)            :: tau
84  REAL(rstd)            :: zz, zvert, fact
85  INTEGER               :: l
86  CHARACTER(len=255)    :: rayleigh_friction_key
87  REAL(rstd)            :: mintau
88 
89  ! New variables added for dissipation vertical profile (SF, 19/09/18)
90  INTEGER    :: vert_prof_dissip
91  REAL(rstd) :: dissip_factz,dissip_deltaz,dissip_zref,pseudoz
92                       
93  INTEGER :: i,j,ij,ind,it,iter,M
94
95   rayleigh_friction_key='none'
96   CALL getin("rayleigh_friction_type",rayleigh_friction_key)
97   SELECT CASE(TRIM(rayleigh_friction_key))
98   CASE('none')
99      rayleigh_friction_type=0
100      IF (is_master) PRINT *, 'No Rayleigh friction'
101   CASE('dcmip2_schaer_noshear')
102      CALL abort_acc("rayleigh_friction_type /= 'none'")
103      rayleigh_friction_type=1
104      rayleigh_shear=0
105      IF (is_master) PRINT *, 'Rayleigh friction : Schaer-like mountain without shear DCMIP2.1'
106   CASE('dcmip2_schaer_shear')
107      CALL abort_acc("rayleigh_friction_type /= 'none'")
108      rayleigh_shear=1
109      rayleigh_friction_type=2
110      IF (is_master) PRINT *, 'Rayleigh friction : Schaer-like mountain with shear DCMIP2.2'
111   CASE('giant_liu_schneider') 
112      CALL abort_acc("rayleigh_friction_type /= 'none'")
113      rayleigh_friction_type=99
114      IF (is_master) PRINT *, 'Rayleigh friction : giant planets Liu Schneider 2010'
115   CASE DEFAULT
116      IF (is_master) PRINT *, 'Bad selector : rayleigh_friction_type =', TRIM(rayleigh_friction_key), &
117        ' in dissip_gcm.f90/init_dissip'
118      STOP
119   END SELECT
120
121   IF(rayleigh_friction_type>0) THEN
122      rayleigh_tau=0.
123      CALL getin("rayleigh_friction_tau",rayleigh_tau)
124      rayleigh_tau = rayleigh_tau / scale_factor
125      IF(rayleigh_tau<=0) THEN
126         IF (is_master) PRINT *, 'Forbidden : negative value for rayleigh_friction_tau =',rayleigh_tau
127         STOP
128      END IF
129      IF(rayleigh_friction_type == 99) THEN
130         rayleigh_limlat=0.
131         CALL getin("rayleigh_limlat",rayleigh_limlat)
132         rayleigh_limlat = rayleigh_limlat*3.14159/180.
133      ENDIF
134   END IF
135
136    CALL allocate_dissip
137    CALL allocate_field(f_u,field_u,type_real,ondevice=.TRUE.)
138    CALL allocate_field(f_du,field_u,type_real,ondevice=.TRUE.)
139    CALL allocate_field(f_theta,field_t,type_real,ondevice=.TRUE.)
140    CALL allocate_field(f_dtheta,field_t,type_real,ondevice=.TRUE.)
141   
142    CALL init_message(f_due_diss1,req_e1_vect,req_due)
143    CALL init_message(f_dtheta_diss,req_i1,req_dtheta)
144
145    tau_graddiv(:)=5000
146    CALL getin("tau_graddiv",tau)
147    tau_graddiv(:)=tau/scale_factor
148
149    CALL getin("nitergdiv",nitergdiv)
150   
151    tau_gradrot(:)=5000
152    CALL getin("tau_gradrot",tau)
153    tau_gradrot(:)=tau/scale_factor
154
155    CALL getin("nitergrot",nitergrot)
156   
157
158    tau_divgrad(:)=5000
159    CALL getin("tau_divgrad",tau)
160    tau_divgrad(:)=tau/scale_factor
161
162    CALL getin("niterdivgrad",niterdivgrad)
163
164
165    cgraddiv=-1
166    cdivgrad=-1
167    cgradrot=-1
168
169!$OMP BARRIER
170!$OMP MASTER
171    DO ind=1,ndomain
172      CALL swap_dimensions(ind)
173      CALL swap_geometry(ind)
174      u=f_u(ind)
175
176! set random seed to get reproductibility when using a different number of process
177      CALL RANDOM_SEED(size=M)
178      CALL RANDOM_SEED(put=(/(i,i=1,M)/))
179
180      ! This cannot be ported on GPU due to compiler limitations
181      DO j=jj_begin,jj_end
182        DO i=ii_begin,ii_end
183          ij=(j-1)*iim+i   
184          CALL RANDOM_NUMBER(r)
185          u(ij+u_right)=r-0.5
186          CALL RANDOM_NUMBER(r)
187          u(ij+u_lup)=r-0.5
188          CALL RANDOM_NUMBER(r)
189          u(ij+u_ldown)=r-0.5
190        ENDDO
191      ENDDO       
192    ENDDO
193!$OMP END MASTER 
194!$OMP BARRIER
195
196    DO it=1,20
197     
198      dumax=0
199      DO iter=1,nitergdiv
200        CALL update_device_field(f_u)
201        CALL transfert_request(f_u,req_e1_vect)
202         
203        DO ind=1,ndomain
204          IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE
205          CALL swap_dimensions(ind)
206          CALL swap_geometry(ind)
207          u=f_u(ind)
208          du=f_du(ind)
209          CALL compute_gradiv_inplace(u,Ai,ne,le,de,1,1)
210          ! This should be ported on GPU but we are running into compiler issues...
211          !$acc update host(u(:)) wait
212          du=u
213        ENDDO
214      ENDDO
215
216      CALL update_device_field(f_du)
217      CALL transfert_request(f_du,req_e1_vect)
218      CALL update_host_field(f_du)
219
220      DO ind=1,ndomain
221        IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE
222        CALL swap_dimensions(ind)
223        CALL swap_geometry(ind)
224        du=f_du(ind)
225         
226        ! Not ported on GPU because all the other kernels cannot be ported
227        DO j=jj_begin,jj_end
228          DO i=ii_begin,ii_end
229            ij=(j-1)*iim+i   
230            if (le(ij+u_right)>1e-100) dumax=MAX(dumax,ABS(du(ij+u_right)))
231            if (le(ij+u_lup)>1e-100)   dumax=MAX(dumax,ABS(du(ij+u_lup)))
232            if (le(ij+u_ldown)>1e-100) dumax=MAX(dumax,ABS(du(ij+u_ldown)))
233          ENDDO
234        ENDDO
235      ENDDO
236
237      IF (using_mpi) THEN
238        CALL reduce_max_omp(dumax,dumax1)
239!$OMP MASTER       
240        CALL MPI_ALLREDUCE(dumax1,dumax,1,MPI_REAL8,MPI_MAX,comm_icosa,ierr)
241!$OMP END MASTER     
242        CALL bcast_omp(dumax) 
243      ELSE
244        CALL allreduce_max_omp(dumax,dumax1)
245        dumax=dumax1
246      ENDIF 
247                       
248      DO ind=1,ndomain
249        IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE
250        CALL swap_dimensions(ind)
251        CALL swap_geometry(ind)
252        u=f_u(ind)
253        du=f_du(ind)
254        ! This should be ported on GPU but we are running into compiler issues...
255        u=du/dumax
256      ENDDO
257      IF (is_master) PRINT *,"gradiv : it :",it ,": dumax",dumax
258
259    ENDDO 
260    IF (is_master) PRINT *,"gradiv : dumax",dumax
261    IF (is_master) PRINT *, 'mean T-cell edge size (km)', 1.45*radius/iim_glo/1000., &
262                              'effective T-cell half-edge size (km)', dumax**(-.5/nitergdiv)/1000
263    IF (is_master)  PRINT *, 'Max. time step assuming c=340 m/s and Courant number=3 (ARK2.3) :', &
264                               3./340.*dumax**(-.5/nitergdiv)
265
266    cgraddiv=dumax**(-1./nitergdiv)
267    IF (is_master) PRINT *,"cgraddiv : ",cgraddiv
268
269!$OMP BARRIER
270!$OMP MASTER
271    DO ind=1,ndomain
272      CALL swap_dimensions(ind)
273      CALL swap_geometry(ind)
274      u=f_u(ind)
275
276! set random seed to get reproductibility when using a different number of process
277      CALL RANDOM_SEED(size=M)
278      CALL RANDOM_SEED(put=(/(i,i=1,M)/))
279 
280      ! This cannot be ported on GPU due to compiler limitations
281       DO j=jj_begin,jj_end
282        DO i=ii_begin,ii_end
283          ij=(j-1)*iim+i   
284          CALL RANDOM_NUMBER(r)
285          u(ij+u_right)=r-0.5
286          CALL RANDOM_NUMBER(r)
287          u(ij+u_lup)=r-0.5
288          CALL RANDOM_NUMBER(r)
289          u(ij+u_ldown)=r-0.5
290        ENDDO
291      ENDDO       
292    ENDDO
293!$OMP END MASTER 
294!$OMP BARRIER
295
296    DO it=1,20
297 
298      dumax=0
299      DO iter=1,nitergrot
300        CALL update_device_field(f_u)
301        CALL transfert_request(f_u,req_e1_vect)
302        DO ind=1,ndomain
303          IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE
304          CALL swap_dimensions(ind)
305          CALL swap_geometry(ind)
306          u=f_u(ind)
307          du=f_du(ind)
308          CALL compute_gradrot_inplace(u,Av,ne,le,de,1,1)
309          ! This should be ported on GPU but we are running into compiler issues...
310          !$acc update host(u(:)) wait
311          du=u
312        ENDDO
313      ENDDO
314     
315      CALL update_device_field(f_du)
316      CALL transfert_request(f_du,req_e1_vect)
317      CALL update_host_field(f_du)
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        du=f_du(ind)
324       
325        ! Not ported on GPU because all the other kernels cannot be ported
326        DO j=jj_begin,jj_end
327          DO i=ii_begin,ii_end
328            ij=(j-1)*iim+i   
329            if (le(ij+u_right)>1e-100) dumax=MAX(dumax,ABS(du(ij+u_right)))
330            if (le(ij+u_lup)>1e-100) dumax=MAX(dumax,ABS(du(ij+u_lup)))
331            if (le(ij+u_ldown)>1e-100) dumax=MAX(dumax,ABS(du(ij+u_ldown)))
332          ENDDO
333        ENDDO
334      ENDDO
335
336      IF (using_mpi) THEN
337        CALL reduce_max_omp(dumax,dumax1)
338!$OMP MASTER       
339        CALL MPI_ALLREDUCE(dumax1,dumax,1,MPI_REAL8,MPI_MAX,comm_icosa,ierr)
340!$OMP END MASTER     
341        CALL bcast_omp(dumax) 
342      ELSE
343        CALL allreduce_max_omp(dumax,dumax1)
344        dumax=dumax1
345      ENDIF 
346
347     
348      DO ind=1,ndomain
349        IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE
350        CALL swap_dimensions(ind)
351        CALL swap_geometry(ind)
352        u=f_u(ind)
353        du=f_du(ind)
354        ! This should be ported on GPU but we are running into compiler issues...
355        u=du/dumax
356      ENDDO
357     
358      IF (is_master) PRINT *,"gradrot : it :",it ,": dumax",dumax
359
360    ENDDO 
361    IF (is_master) PRINT *,"gradrot : dumax",dumax
362 
363    cgradrot=dumax**(-1./nitergrot) 
364    IF (is_master) PRINT *,"cgradrot : ",cgradrot
365   
366
367
368!$OMP BARRIER
369!$OMP MASTER
370    DO ind=1,ndomain
371      CALL swap_dimensions(ind)
372      CALL swap_geometry(ind)
373      theta=f_theta(ind)
374
375! set random seed to get reproductibility when using a different number of process
376      CALL RANDOM_SEED(size=M)
377      CALL RANDOM_SEED(put=(/(i,i=1,M)/))
378
379      ! This cannot be ported on GPU due to compiler limitations
380      DO j=jj_begin,jj_end
381        DO i=ii_begin,ii_end
382          ij=(j-1)*iim+i   
383          CALL RANDOM_NUMBER(r)
384          theta(ij)=r-0.5
385        ENDDO
386      ENDDO 
387    ENDDO
388!$OMP END MASTER
389!$OMP BARRIER
390
391    DO it=1,20
392
393      dthetamax=0
394      DO iter=1,niterdivgrad
395        CALL update_device_field(f_theta)
396        CALL transfert_request(f_theta,req_i1)
397        DO ind=1,ndomain
398          IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE
399          CALL swap_dimensions(ind)
400          CALL swap_geometry(ind)
401          theta=f_theta(ind)
402          dtheta=f_dtheta(ind)
403          CALL compute_divgrad_inplace(theta,Ai,ne,le,de,1,1)
404          ! This should be ported on GPU but we are running into compiler issues...
405          !$acc update host(theta(:)) wait
406          dtheta=theta
407        ENDDO
408      ENDDO
409
410      CALL update_device_field(f_dtheta)
411      CALL transfert_request(f_dtheta,req_i1)
412      CALL update_host_field(f_dtheta)
413     
414      DO ind=1,ndomain
415        IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE
416        CALL swap_dimensions(ind)
417        CALL swap_geometry(ind)
418        dtheta=f_dtheta(ind)
419       
420        ! Not ported on GPU because all the other kernels cannot be ported
421        DO j=jj_begin,jj_end
422          DO i=ii_begin,ii_end
423            ij=(j-1)*iim+i   
424            dthetamax=MAX(dthetamax,ABS(dtheta(ij)))
425          ENDDO
426        ENDDO
427      ENDDO
428
429      IF (using_mpi) THEN
430        CALL reduce_max_omp(dthetamax ,dthetamax1)
431!$OMP MASTER       
432        CALL MPI_ALLREDUCE(dthetamax1,dthetamax,1,MPI_REAL8,MPI_MAX,comm_icosa,ierr)
433!$OMP END MASTER     
434        CALL bcast_omp(dthetamax) 
435      ELSE
436        CALL allreduce_max_omp(dthetamax,dthetamax1)
437        dumax=dumax1
438      ENDIF 
439     
440      IF (is_master) PRINT *,"divgrad : it :",it ,": dthetamax",dthetamax
441
442      DO ind=1,ndomain
443        IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE
444        CALL swap_dimensions(ind)
445        CALL swap_geometry(ind)
446        theta=f_theta(ind)
447        dtheta=f_dtheta(ind)
448        ! This should be ported on GPU but we are running into compiler issues...
449        theta=dtheta/dthetamax
450      ENDDO
451    ENDDO 
452
453    IF (is_master) PRINT *,"divgrad : divgrad",dthetamax
454 
455    cdivgrad=dthetamax**(-1./niterdivgrad) 
456    IF (is_master) PRINT *,"cdivgrad : ",cdivgrad
457
458    ! Added CMIP5 dissipation vertical profile (SF, 19/09/18)
459    vert_prof_dissip = 1
460    CALL getin('vert_prof_dissip',vert_prof_dissip)
461    IF (vert_prof_dissip == 1) then
462       dissip_zref = 30.
463       dissip_deltaz = 10.
464       dissip_factz = 4.
465       CALL getin('dissip_factz',dissip_factz )
466       CALL getin('dissip_deltaz',dissip_deltaz )
467       CALL getin('dissip_zref',dissip_zref )
468       DO l=1,llm
469          pseudoz=8.*log(preff/presnivs(l))
470          zvert=1+ &
471               (tanh((pseudoz-dissip_zref)/dissip_deltaz)+1.)/2. &
472               *(dissip_factz-1.)
473          tau_graddiv(l) = tau_graddiv(l)/zvert
474          tau_gradrot(l) = tau_gradrot(l)/zvert
475          tau_divgrad(l) = tau_divgrad(l)/zvert
476       ENDDO
477    ELSE     
478       fact=2
479       DO l=1,llm
480          IF(ap_bp_present) THEN ! height-dependent dissipation
481             zz= 1. - preff/presnivs(l)
482          ELSE
483             zz = 0.
484          END IF
485          zvert=fact-(fact-1)/(1+zz*zz)
486          tau_graddiv(l) = tau_graddiv(l)/zvert
487          tau_gradrot(l) = tau_gradrot(l)/zvert
488          tau_divgrad(l) = tau_divgrad(l)/zvert
489       ENDDO
490    ENDIF
491
492    mintau=tau_graddiv(1)
493    DO l=1,llm
494      mintau=MIN(mintau,tau_graddiv(l))
495      mintau=MIN(mintau,tau_gradrot(l))
496      mintau=MIN(mintau,tau_divgrad(l))
497    ENDDO
498
499    IF(rayleigh_friction_type>0) mintau=MIN(mintau,rayleigh_tau)
500
501    IF(mintau>0) THEN
502       IF (itau_dissip==0) THEN
503         IF (is_master) PRINT *,"init_dissip: Automatic computation of itau_dissip..."
504         itau_dissip=INT(mintau/dt)
505       ENDIF
506    ELSE
507       IF (is_master) PRINT *,"init_dissip: No dissipation time set, setting itau_dissip to 1000000000"
508       itau_dissip=100000000
509    END IF
510    itau_dissip=MAX(1,itau_dissip)
511    dtdissip=itau_dissip*dt
512    IF (is_master) THEN
513      PRINT *,"init_dissip: rayleigh_tau",rayleigh_tau, "mintau ",mintau
514      PRINT *,"init_dissip: itau_dissip",itau_dissip," dtdissip ",dtdissip
515    ENDIF
516
517    !$acc update device(tau_graddiv(:),tau_gradrot(:),tau_divgrad(:)) async
518
519  END SUBROUTINE init_dissip 
520 
521 
522  SUBROUTINE dissip(f_ps,f_mass,f_phis,f_geopot,f_theta_rhodz,f_ue, f_dtheta_rhodz,f_due)
523  USE icosa
524  USE theta2theta_rhodz_mod
525  USE pression_mod
526  USE exner_mod
527  USE geopotential_mod
528  USE trace
529  USE time_mod
530  USE omp_para
531  USE abort_mod
532  IMPLICIT NONE
533    TYPE(t_field),POINTER :: f_ps(:), f_mass(:), f_phis(:), f_geopot(:)
534    TYPE(t_field),POINTER :: f_theta_rhodz(:), f_dtheta_rhodz(:)
535    TYPE(t_field),POINTER :: f_ue(:), f_due(:)
536
537    REAL(rstd),POINTER         :: due(:,:)
538    REAL(rstd),POINTER         :: phi(:,:), ue(:,:)
539    REAL(rstd),POINTER         :: due_diss1(:,:)
540    REAL(rstd),POINTER         :: due_diss2(:,:)
541    REAL(rstd),POINTER         :: dtheta_rhodz(:,:,:)
542    REAL(rstd),POINTER         :: dtheta_rhodz_diss(:,:)
543
544    INTEGER :: ind, shear
545    INTEGER :: l,ij,nn
546
547!$OMP BARRIER
548   
549    CALL trace_start("dissip")
550    CALL gradiv(f_ue,f_due_diss1)
551    CALL gradrot(f_ue,f_due_diss2)
552
553    CALL divgrad_theta_rhodz(f_mass,f_theta_rhodz,f_dtheta_rhodz_diss)
554   
555    DO ind=1,ndomain
556      IF (.NOT. assigned_domain(ind)) CYCLE
557      CALL swap_dimensions(ind)
558      CALL swap_geometry(ind)
559      due=f_due(ind) 
560      due_diss1=f_due_diss1(ind)
561      due_diss2=f_due_diss2(ind)
562      dtheta_rhodz=f_dtheta_rhodz(ind)
563      dtheta_rhodz_diss=f_dtheta_rhodz_diss(ind)
564
565      !$acc parallel loop collapse(2) present(due(:,:), dtheta_rhodz(:,:,:), due_diss1(:,:), due_diss2(:,:), dtheta_rhodz_diss(:,:), tau_graddiv(:), tau_gradrot(:), tau_divgrad(:)) async
566      DO l=ll_begin,ll_end
567!DIR$ SIMD
568        DO ij=ij_begin,ij_end
569
570            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 
571            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
572            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
573
574            dtheta_rhodz(ij,l,1) = -0.5*dtheta_rhodz_diss(ij,l)/tau_divgrad(l)*itau_dissip
575        ENDDO
576      ENDDO
577      !$acc end parallel loop
578
579!      dtheta_rhodz=0
580!      due=0
581
582      IF(rayleigh_friction_type>0) THEN
583        CALL abort_acc("dissip/rayleigh_friction_type>0")
584       IF(rayleigh_friction_type<99) THEN
585         phi=f_geopot(ind)
586         ue=f_ue(ind)
587         DO l=ll_begin,ll_end
588            DO ij=ij_begin,ij_end
589               CALL relax(t_right, u_right)
590               CALL relax(t_lup,   u_lup)
591               CALL relax(t_ldown, u_ldown)
592            ENDDO
593         END DO
594       ELSE
595         ue=f_ue(ind)
596            !$acc parallel loop present(ue(:,:), due(:,:), lat_e(:))
597            DO ij=ij_begin,ij_end
598              nn = ij+u_right
599              IF (ABS(lat_e(nn)) .gt. rayleigh_limlat) THEN
600                !print*, "latitude", lat_e(nn)*180./3.14159
601                due(nn,ll_begin:ll_begin+1) = due(nn,ll_begin:ll_begin+1) - (ue(nn,ll_begin:ll_begin+1)/rayleigh_tau)
602              ENDIF
603              nn = ij+u_lup
604              IF (ABS(lat_e(nn)) .gt. rayleigh_limlat) THEN
605                due(nn,ll_begin:ll_begin+1) = due(nn,ll_begin:ll_begin+1) - (ue(nn,ll_begin:ll_begin+1)/rayleigh_tau)
606              ENDIF
607              nn = ij+u_ldown
608              IF (ABS(lat_e(nn)) .gt. rayleigh_limlat) THEN
609                due(nn,ll_begin:ll_begin+1) = due(nn,ll_begin:ll_begin+1) - (ue(nn,ll_begin:ll_begin+1)/rayleigh_tau)
610              ENDIF
611            ENDDO
612            !$acc end parallel loop
613       ENDIF
614      END IF
615   END DO
616   CALL trace_end("dissip")
617
618   !CALL write_dissip_tendencies
619!$OMP BARRIER
620
621    CONTAINS
622
623      SUBROUTINE relax(shift_t, shift_u)
624        USE dcmip_initial_conditions_test_1_2_3
625        REAL(rstd) :: z, ulon,ulat, & ! input to test2_schaer_mountain
626             p,hyam,hybm,w,t,phis,ps,rho,q, &   ! unused input/output to test2_schaer_mountain
627             fz, u3d(3), uref
628        REAL(rstd), PARAMETER :: zh=2e4,ztop=3e4  ! DCMIP values
629        LOGICAL :: hybrid_eta
630        INTEGER :: shift_u, shift_t, zcoords, nn
631        z = (phi(ij,l)+phi(ij+shift_t,l))/(2.*g)
632        IF(z>zh) THEN  ! relax only in the sponge layer z>zh
633
634           nn = ij+shift_u
635           zcoords = 1 ; hybrid_eta = .FALSE. ! use z instead of p or hyam/hybm
636           CALL test2_schaer_mountain(lon_e(nn),lat_e(nn),p,z,zcoords,hybrid_eta, &
637                hyam,hybm,shear,ulon,ulat,w,t,phis,ps,rho,q)
638!           u3d = ulon*elon_e(nn,:) + ulat*elat_e(nn,:)
639           u3d = ulon*elon_e(nn,:) ! ulat=0
640           uref = sum(u3d*ep_e(nn,:))
641
642           fz = sin((pi/2)*(z-zh)/(ztop-zh))
643           fz = fz*fz/rayleigh_tau
644           due(nn,l) = due(nn,l) - itau_dissip*fz*(ue(nn,l)-uref)
645!           fz = 1./rayleigh_tau
646!           due(nn,l) = due(nn,l) - fz*ue(nn,l)
647         END IF
648       END SUBROUTINE relax
649     
650       SUBROUTINE write_dissip_tendencies
651         USE observable_mod, ONLY : f_buf_ulon, f_buf_ulat
652         USE wind_mod
653         USE output_field_mod
654
655         CALL transfert_request(f_due_diss1,req_e1_vect)
656         CALL un2ulonlat(f_due_diss1, f_buf_ulon, f_buf_ulat, (1./(tau_graddiv(1))))
657         CALL output_field("dulon_diss1",f_buf_ulon)
658         CALL output_field("dulat_diss1",f_buf_ulat)
659!
660         CALL transfert_request(f_due_diss2,req_e1_vect)
661         CALL un2ulonlat(f_due_diss2, f_buf_ulon, f_buf_ulat, (1./(tau_graddiv(1))))
662         CALL output_field("dulon_diss2",f_buf_ulon)
663         CALL output_field("dulat_diss2",f_buf_ulat)
664       END SUBROUTINE write_dissip_tendencies
665
666  END SUBROUTINE dissip
667
668  SUBROUTINE gradiv(f_ue,f_due)
669  USE icosa
670  USE trace
671  USE omp_para
672  IMPLICIT NONE
673    TYPE(t_field),POINTER :: f_ue(:)
674    TYPE(t_field),POINTER :: f_due(:)
675    REAL(rstd),POINTER    :: ue(:,:)
676    REAL(rstd),POINTER    :: due(:,:)
677    INTEGER :: ind
678    INTEGER :: it,l,ij
679       
680    CALL trace_start("gradiv")
681
682    DO ind=1,ndomain
683      IF (.NOT. assigned_domain(ind)) CYCLE
684      CALL swap_dimensions(ind)
685      CALL swap_geometry(ind)
686      ue=f_ue(ind)
687      due=f_due(ind) 
688      !$acc parallel loop present(ue(:,:), due(:,:)) async
689      DO  l = ll_begin, ll_end
690        !$acc loop
691!DIR$ SIMD
692        DO ij=ij_begin,ij_end
693             due(ij+u_right,l)=ue(ij+u_right,l)
694             due(ij+u_lup,l)=ue(ij+u_lup,l)
695             due(ij+u_ldown,l)=ue(ij+u_ldown,l)
696        ENDDO
697      ENDDO
698    ENDDO
699
700    DO it=1,nitergdiv
701      CALL send_message(f_due,req_due)
702      CALL wait_message(req_due)
703
704      DO ind=1,ndomain
705        IF (.NOT. assigned_domain(ind)) CYCLE
706        CALL swap_dimensions(ind)
707        CALL swap_geometry(ind)
708        due=f_due(ind)
709        CALL compute_gradiv_inplace(due,Ai,ne,le,de,ll_begin,ll_end)
710      ENDDO
711    ENDDO
712
713   CALL trace_end("gradiv")
714
715  END SUBROUTINE gradiv
716 
717
718  SUBROUTINE gradrot(f_ue,f_due)
719  USE icosa
720  USE trace
721  USE omp_para
722  IMPLICIT NONE
723    TYPE(t_field),POINTER :: f_ue(:)
724    TYPE(t_field),POINTER :: f_due(:)
725    REAL(rstd),POINTER    :: ue(:,:)
726    REAL(rstd),POINTER    :: due(:,:)
727    INTEGER :: ind
728    INTEGER :: it,l,ij
729       
730    CALL trace_start("gradrot")
731
732    DO ind=1,ndomain
733      IF (.NOT. assigned_domain(ind)) CYCLE
734      CALL swap_dimensions(ind)
735      CALL swap_geometry(ind)
736      ue=f_ue(ind)
737      due=f_due(ind) 
738      !$acc parallel loop present(ue(:,:), due(:,:)) async
739      DO  l = ll_begin, ll_end
740         !$acc loop
741!DIR$ SIMD
742        DO ij=ij_begin,ij_end
743             due(ij+u_right,l)=ue(ij+u_right,l)
744             due(ij+u_lup,l)=ue(ij+u_lup,l)
745             due(ij+u_ldown,l)=ue(ij+u_ldown,l)
746        ENDDO
747      ENDDO
748    ENDDO
749
750    DO it=1,nitergrot
751      CALL send_message(f_due,req_due)
752      CALL wait_message(req_due)
753       
754      DO ind=1,ndomain
755        IF (.NOT. assigned_domain(ind)) CYCLE
756        CALL swap_dimensions(ind)
757        CALL swap_geometry(ind)
758        due=f_due(ind) 
759        CALL compute_gradrot_inplace(due,Av,ne,le,de,ll_begin,ll_end)
760      ENDDO
761
762    ENDDO
763
764    CALL trace_end("gradrot")
765
766  END SUBROUTINE gradrot
767 
768  SUBROUTINE divgrad(f_theta,f_dtheta)
769  USE icosa
770  USE trace
771  USE omp_para
772  IMPLICIT NONE
773    TYPE(t_field),POINTER :: f_theta(:)
774    TYPE(t_field),POINTER :: f_dtheta(:)
775    REAL(rstd),POINTER    :: theta(:,:)
776    REAL(rstd),POINTER    :: dtheta(:,:)
777    INTEGER :: ind, l, ij
778    INTEGER :: it
779
780    CALL trace_start("divgrad")
781       
782    DO ind=1,ndomain
783      IF (.NOT. assigned_domain(ind)) CYCLE
784      CALL swap_dimensions(ind)
785      CALL swap_geometry(ind)
786      theta=f_theta(ind)
787      dtheta=f_dtheta(ind) 
788      ! Replace Fortran 90 construct dtheta=theta because it confuses PGI acc kernels...
789      !$acc parallel loop collapse(2) present(theta(:,:), dtheta(:,:)) async
790      DO l=ll_begin,ll_end
791!DIR$ SIMD
792        DO ij=1,iim*jjm
793          dtheta(ij,l)=theta(ij,l)
794        ENDDO
795      ENDDO
796      !$acc end parallel loop
797    ENDDO
798
799    DO it=1,niterdivgrad
800      CALL transfert_request(f_dtheta,req_i1)
801       
802      DO ind=1,ndomain
803        IF (.NOT. assigned_domain(ind)) CYCLE
804        CALL swap_dimensions(ind)
805        CALL swap_geometry(ind)
806        dtheta=f_dtheta(ind) 
807        CALL compute_divgrad_inplace(dtheta,Ai,ne,le,de,ll_begin,ll_end)
808      ENDDO
809
810    ENDDO
811
812    CALL trace_end("divgrad")
813
814  END SUBROUTINE divgrad
815   
816  SUBROUTINE divgrad_theta_rhodz(f_mass,f_theta_rhodz,f_dtheta_rhodz)
817  USE icosa
818  USE trace
819  USE omp_para
820  IMPLICIT NONE
821    TYPE(t_field),POINTER :: f_mass(:)
822    TYPE(t_field),POINTER :: f_theta_rhodz(:)
823    TYPE(t_field),POINTER :: f_dtheta_rhodz(:)
824   
825    REAL(rstd),POINTER :: mass(:,:)
826    REAL(rstd),POINTER :: theta_rhodz(:,:,:)
827    REAL(rstd),POINTER :: dtheta_rhodz(:,:)
828
829    INTEGER :: ind
830    INTEGER :: it,l,ij
831
832    CALL trace_start("divgrad")
833
834    DO ind=1,ndomain
835      IF (.NOT. assigned_domain(ind)) CYCLE
836      CALL swap_dimensions(ind)
837      CALL swap_geometry(ind)
838      mass=f_mass(ind)
839      theta_rhodz=f_theta_rhodz(ind)
840      dtheta_rhodz=f_dtheta_rhodz(ind) 
841      !$acc parallel loop present(mass(:,:), theta_rhodz(:,:,:), dtheta_rhodz(:,:)) async
842      DO  l = ll_begin, ll_end
843        !$acc loop
844!DIR$ SIMD
845        DO ij=ij_begin,ij_end
846            dtheta_rhodz(ij,l) = theta_rhodz(ij,l,1) / mass(ij,l)
847        ENDDO
848      ENDDO
849    ENDDO
850
851    DO it=1,niterdivgrad
852      CALL send_message(f_dtheta_rhodz,req_dtheta)
853      CALL wait_message(req_dtheta)
854      DO ind=1,ndomain
855        IF (.NOT. assigned_domain(ind)) CYCLE
856        CALL swap_dimensions(ind)
857        CALL swap_geometry(ind)
858        dtheta_rhodz=f_dtheta_rhodz(ind) 
859        CALL compute_divgrad_inplace(dtheta_rhodz,Ai,ne,le,de,ll_begin,ll_end)
860      ENDDO
861
862    ENDDO
863
864    DO ind=1,ndomain
865      IF (.NOT. assigned_domain(ind)) CYCLE
866      CALL swap_dimensions(ind)
867      CALL swap_geometry(ind)
868      dtheta_rhodz=f_dtheta_rhodz(ind) 
869      mass=f_mass(ind)
870   
871      !$acc parallel loop collapse(2) present(mass(:,:), dtheta_rhodz(:,:)) async
872      DO  l = ll_begin, ll_end
873!DIR$ SIMD
874        DO ij=ij_begin,ij_end
875            dtheta_rhodz(ij,l) = dtheta_rhodz(ij,l) * mass(ij,l)
876        ENDDO
877      ENDDO
878      !$acc end parallel loop
879    ENDDO
880
881
882    CALL trace_end("divgrad")
883
884  END SUBROUTINE divgrad_theta_rhodz
885
886  SUBROUTINE compute_gradiv(ue,gradivu_e,Ai,ne,le,de,llb,lle)
887    INTEGER,INTENT(IN)     :: llb
888    INTEGER,INTENT(IN)     :: lle
889    REAL(rstd),INTENT(OUT) :: gradivu_e(iim*3*jjm,llm)
890    REAL(rstd),INTENT(IN)  :: ue(iim*3*jjm,llm)
891    REAL(rstd),INTENT(IN)  :: Ai(iim*jjm)
892    INTEGER,INTENT(IN)     :: ne(iim*jjm,6)
893    REAL(rstd),INTENT(IN)  :: le(iim*3*jjm)
894    REAL(rstd),INTENT(IN)  :: de(iim*3*jjm)
895
896    gradivu_e = ue
897    CALL compute_gradiv_inplace(gradivu_e,Ai,ne,le,de,llb,lle)
898
899  END SUBROUTINE compute_gradiv
900 
901  SUBROUTINE compute_gradiv_inplace(ue_gradivu_e,Ai,ne,le,de,llb,lle)
902    INTEGER,INTENT(IN)     :: llb
903    INTEGER,INTENT(IN)     :: lle
904    REAL(rstd),INTENT(INOUT) :: ue_gradivu_e(iim*3*jjm,llm)
905    REAL(rstd),INTENT(IN)  :: Ai(iim*jjm)
906    INTEGER,INTENT(IN)     :: ne(iim*jjm,6)
907    REAL(rstd),INTENT(IN)  :: le(iim*3*jjm)
908    REAL(rstd),INTENT(IN)  :: de(iim*3*jjm)
909    REAL(rstd) :: divu_i(iim*jjm,llb:lle)
910   
911    INTEGER :: ij,l
912
913    ! ue and gradivu_e second dimension is not always llm so use the bounds explicitly
914    !$acc data present( ue_gradivu_e(:,llb:lle), Ai(:), ne(:,:), le(:), de(:)) create(divu_i(:,:)) async
915
916    !$acc parallel loop collapse(2) async
917    DO l=llb,lle
918      !DIR$ SIMD
919      DO ij=ij_begin,ij_end
920          divu_i(ij,l)=1./Ai(ij)*(ne(ij,right)*ue_gradivu_e(ij+u_right,l)*le(ij+u_right)  +  &
921                             ne(ij,rup)*ue_gradivu_e(ij+u_rup,l)*le(ij+u_rup)        +  & 
922                             ne(ij,lup)*ue_gradivu_e(ij+u_lup,l)*le(ij+u_lup)        +  & 
923                             ne(ij,left)*ue_gradivu_e(ij+u_left,l)*le(ij+u_left)     +  & 
924                             ne(ij,ldown)*ue_gradivu_e(ij+u_ldown,l)*le(ij+u_ldown)  +  & 
925                             ne(ij,rdown)*ue_gradivu_e(ij+u_rdown,l)*le(ij+u_rdown))
926      ENDDO
927    ENDDO
928    !$acc end parallel loop
929   
930    !$acc parallel loop collapse(2) async
931    DO l=llb,lle
932      !DIR$ SIMD
933      DO ij=ij_begin,ij_end
934          ue_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) ) 
935          ue_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))   
936          ue_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) )
937      ENDDO
938    ENDDO
939    !$acc end parallel loop
940
941    !$acc parallel loop collapse(2) async
942    DO l=llb,lle
943      !DIR$ SIMD
944      DO ij=ij_begin,ij_end
945          ue_gradivu_e(ij+u_right,l)=-ue_gradivu_e(ij+u_right,l)*cgraddiv
946          ue_gradivu_e(ij+u_lup,l)=-ue_gradivu_e(ij+u_lup,l)*cgraddiv
947          ue_gradivu_e(ij+u_ldown,l)=-ue_gradivu_e(ij+u_ldown,l)*cgraddiv
948      ENDDO
949    ENDDO
950    !$acc end parallel loop
951    !$acc end data   
952  END SUBROUTINE compute_gradiv_inplace
953
954  SUBROUTINE compute_divgrad(theta,divgrad_i,Ai,ne,le,de,llb,lle)
955    INTEGER,INTENT(IN)     :: llb
956    INTEGER,INTENT(IN)     :: lle
957    REAL(rstd),INTENT(IN) :: theta(iim*jjm,1:lle)
958    REAL(rstd),INTENT(OUT) :: divgrad_i(iim*jjm,1:lle)
959    REAL(rstd),INTENT(IN)  :: Ai(iim*jjm)
960    INTEGER,INTENT(IN)     :: ne(iim*jjm,6)
961    REAL(rstd),INTENT(IN)  :: le(iim*3*jjm)
962    REAL(rstd),INTENT(IN)  :: de(iim*3*jjm)
963
964    divgrad_i = theta
965    CALL compute_divgrad_inplace(divgrad_i,Ai,ne,le,de,llb,lle)
966  END SUBROUTINE compute_divgrad
967 
968  SUBROUTINE compute_divgrad_inplace(theta_divgrad_i,Ai,ne,le,de,llb,lle)
969    INTEGER,INTENT(IN)     :: llb
970    INTEGER,INTENT(IN)     :: lle
971    REAL(rstd),INTENT(INOUT) :: theta_divgrad_i(iim*jjm,1:lle)
972    REAL(rstd),INTENT(IN)  :: Ai(iim*jjm)
973    INTEGER,INTENT(IN)     :: ne(iim*jjm,6)
974    REAL(rstd),INTENT(IN)  :: le(iim*3*jjm)
975    REAL(rstd),INTENT(IN)  :: de(iim*3*jjm)
976    REAL(rstd)  :: grad_e(3*iim*jjm,llb:lle)
977
978    INTEGER :: ij,l
979
980    ! theta and divgrad_i second dimension is not always llm so use the bounds explicitly
981    !$acc data present(theta_divgrad_i(:,llb:lle), Ai(:), de(:), ne(:,:), le(:)) create(grad_e(:,:)) async
982       
983    !$acc parallel loop collapse(2) async
984    DO l=llb,lle
985      !DIR$ SIMD
986      DO ij=ij_begin_ext,ij_end_ext 
987          grad_e(ij+u_right,l)=-1/de(ij+u_right)*(ne(ij,right)*theta_divgrad_i(ij,l)+ ne(ij+t_right,left)*theta_divgrad_i(ij+t_right,l) )
988          grad_e(ij+u_lup,l)=-1/de(ij+u_lup)*(ne(ij,lup)*theta_divgrad_i(ij,l)+ ne(ij+t_lup,rdown)*theta_divgrad_i(ij+t_lup,l ))
989          grad_e(ij+u_ldown,l)=-1/de(ij+u_ldown)*(ne(ij,ldown)*theta_divgrad_i(ij,l)+ne(ij+t_ldown,rup)*theta_divgrad_i(ij+t_ldown,l) )
990      ENDDO
991    ENDDO
992    !$acc end parallel loop
993   
994   
995    !$acc parallel loop collapse(2) async
996    DO l=llb,lle
997      !DIR$ SIMD
998      DO ij=ij_begin,ij_end
999          theta_divgrad_i(ij,l)=1./Ai(ij)*(ne(ij,right)*grad_e(ij+u_right,l)*le(ij+u_right)  +  &
1000                             ne(ij,rup)*grad_e(ij+u_rup,l)*le(ij+u_rup)              +  &
1001                             ne(ij,lup)*grad_e(ij+u_lup,l)*le(ij+u_lup)              +  &
1002                             ne(ij,left)*grad_e(ij+u_left,l)*le(ij+u_left)           +  &
1003                             ne(ij,ldown)*grad_e(ij+u_ldown,l)*le(ij+u_ldown)        +  & 
1004                             ne(ij,rdown)*grad_e(ij+u_rdown,l)*le(ij+u_rdown))
1005      ENDDO
1006    ENDDO
1007    !$acc end parallel loop
1008   
1009    !$acc parallel loop collapse(2) async
1010    DO l=llb,lle
1011      DO ij=ij_begin,ij_end
1012          theta_divgrad_i(ij,l)=-theta_divgrad_i(ij,l)*cdivgrad
1013      ENDDO
1014    ENDDO
1015    !$acc end parallel loop
1016    !$acc end data
1017  END SUBROUTINE compute_divgrad_inplace
1018
1019  SUBROUTINE compute_gradrot(ue,gradrot_e,Av,ne,le,de,llb,lle)
1020    INTEGER,INTENT(IN)     :: llb
1021    INTEGER,INTENT(IN)     :: lle
1022    REAL(rstd),INTENT(IN) :: ue(iim*3*jjm,lle)
1023    REAL(rstd),INTENT(OUT) :: gradrot_e(iim*3*jjm,lle)
1024    REAL(rstd),INTENT(IN)  :: Av(2*iim*jjm)
1025    INTEGER,INTENT(IN)     :: ne(iim*jjm,6)
1026    REAL(rstd),INTENT(IN)  :: le(iim*3*jjm)
1027    REAL(rstd),INTENT(IN)  :: de(iim*3*jjm)
1028
1029    gradrot_e = ue
1030    CALL compute_gradrot_inplace(gradrot_e,Av,ne,le,de,llb,lle)
1031  END SUBROUTINE compute_gradrot
1032
1033  SUBROUTINE compute_gradrot_inplace(ue_gradrot_e,Av,ne,le,de,llb,lle)
1034    INTEGER,INTENT(IN)     :: llb
1035    INTEGER,INTENT(IN)     :: lle
1036    REAL(rstd),INTENT(INOUT) :: ue_gradrot_e(iim*3*jjm,lle)
1037    REAL(rstd),INTENT(IN)  :: Av(2*iim*jjm)
1038    INTEGER,INTENT(IN)     :: ne(iim*jjm,6)
1039    REAL(rstd),INTENT(IN)  :: le(iim*3*jjm)
1040    REAL(rstd),INTENT(IN)  :: de(iim*3*jjm)
1041    REAL(rstd) :: rot_v(2*iim*jjm,llb:lle)
1042
1043    INTEGER :: ij,l
1044
1045    ! ue and gradrot_e second dimension is not always llm so use the bounds explicitly
1046    ! gradrot_e should be copyout but using copy instead allows to compare the output
1047    ! more easily as the code sometimes uses unintialed values
1048    !$acc data present(ue_gradrot_e(:,llb:lle), Av(:), ne(:,:), de(:), le(:)) create(rot_v(:,:)) async
1049     
1050    !$acc parallel loop collapse(2) async
1051    DO l=llb,lle
1052!DIR$ SIMD
1053      DO ij=ij_begin_ext,ij_end_ext       
1054          rot_v(ij+z_up,l)= 1./Av(ij+z_up)*(  ne(ij,rup)*ue_gradrot_e(ij+u_rup,l)*de(ij+u_rup) &
1055                                + ne(ij+t_rup,left)*ue_gradrot_e(ij+t_rup+u_left,l)*de(ij+t_rup+u_left) &
1056                                - ne(ij,lup)*ue_gradrot_e(ij+u_lup,l)*de(ij+u_lup) )                               
1057          rot_v(ij+z_down,l) = 1./Av(ij+z_down)*( ne(ij,ldown)*ue_gradrot_e(ij+u_ldown,l)*de(ij+u_ldown) &
1058                                     + ne(ij+t_ldown,right)*ue_gradrot_e(ij+t_ldown+u_right,l)*de(ij+t_ldown+u_right) &
1059                                     - ne(ij,rdown)*ue_gradrot_e(ij+u_rdown,l)*de(ij+u_rdown) )
1060      ENDDO
1061    ENDDO                             
1062    !$acc end parallel loop
1063 
1064    !$acc parallel loop collapse(2) async
1065    DO l=llb,lle
1066!DIR$ SIMD
1067      DO ij=ij_begin,ij_end
1068          ue_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))       
1069          ue_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))       
1070          ue_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))
1071      ENDDO
1072    ENDDO
1073    !$acc end parallel loop
1074
1075    !$acc parallel loop collapse(2) async
1076    DO l=llb,lle
1077!DIR$ SIMD
1078      DO ij=ij_begin,ij_end   
1079          ue_gradrot_e(ij+u_right,l)=-ue_gradrot_e(ij+u_right,l)*cgradrot       
1080          ue_gradrot_e(ij+u_lup,l)=-ue_gradrot_e(ij+u_lup,l)*cgradrot
1081          ue_gradrot_e(ij+u_ldown,l)=-ue_gradrot_e(ij+u_ldown,l)*cgradrot       
1082      ENDDO
1083    ENDDO
1084    !$acc end parallel loop
1085    !$acc end data   
1086  END SUBROUTINE compute_gradrot_inplace
1087
1088
1089END MODULE dissip_gcm_mod
Note: See TracBrowser for help on using the repository browser.