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

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

trunk : GPU implementation with OpenACC ( merge from glcp.idris.fr )
File size: 34.2 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,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,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,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,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,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,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,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,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
892    gradivu_e = ue
893    CALL compute_gradiv_inplace(gradivu_e,llb,lle)
894
895  END SUBROUTINE compute_gradiv
896 
897  SUBROUTINE compute_gradiv_inplace(ue_gradivu_e,llb,lle)
898    USE geometry, ONLY : Ai, ne, le, de
899    INTEGER,INTENT(IN)     :: llb
900    INTEGER,INTENT(IN)     :: lle
901    REAL(rstd),INTENT(INOUT) :: ue_gradivu_e(iim*3*jjm,llm)
902    REAL(rstd) :: divu_i(iim*jjm,llb:lle)
903   
904    INTEGER :: ij,l
905
906    ! ue and gradivu_e second dimension is not always llm so use the bounds explicitly
907    !$acc data present( ue_gradivu_e(:,llb:lle), Ai(:), ne(:,:), le(:), de(:)) create(divu_i(:,:)) async
908
909    !$acc parallel loop collapse(2) async
910    DO l=llb,lle
911      !DIR$ SIMD
912      DO ij=ij_begin,ij_end
913          divu_i(ij,l)=1./Ai(ij)*(ne(ij,right)*ue_gradivu_e(ij+u_right,l)*le(ij+u_right)  +  &
914                             ne(ij,rup)*ue_gradivu_e(ij+u_rup,l)*le(ij+u_rup)        +  & 
915                             ne(ij,lup)*ue_gradivu_e(ij+u_lup,l)*le(ij+u_lup)        +  & 
916                             ne(ij,left)*ue_gradivu_e(ij+u_left,l)*le(ij+u_left)     +  & 
917                             ne(ij,ldown)*ue_gradivu_e(ij+u_ldown,l)*le(ij+u_ldown)  +  & 
918                             ne(ij,rdown)*ue_gradivu_e(ij+u_rdown,l)*le(ij+u_rdown))
919      ENDDO
920    ENDDO
921    !$acc end parallel loop
922   
923    !$acc parallel loop collapse(2) async
924    DO l=llb,lle
925      !DIR$ SIMD
926      DO ij=ij_begin,ij_end
927          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) ) 
928          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))   
929          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) )
930      ENDDO
931    ENDDO
932    !$acc end parallel loop
933
934    !$acc parallel loop collapse(2) async
935    DO l=llb,lle
936      !DIR$ SIMD
937      DO ij=ij_begin,ij_end
938          ue_gradivu_e(ij+u_right,l)=-ue_gradivu_e(ij+u_right,l)*cgraddiv
939          ue_gradivu_e(ij+u_lup,l)=-ue_gradivu_e(ij+u_lup,l)*cgraddiv
940          ue_gradivu_e(ij+u_ldown,l)=-ue_gradivu_e(ij+u_ldown,l)*cgraddiv
941      ENDDO
942    ENDDO
943    !$acc end parallel loop
944    !$acc end data   
945  END SUBROUTINE compute_gradiv_inplace
946
947  SUBROUTINE compute_divgrad(theta,divgrad_i,llb,lle)
948    INTEGER,INTENT(IN)     :: llb
949    INTEGER,INTENT(IN)     :: lle
950    REAL(rstd),INTENT(IN) :: theta(iim*jjm,1:lle)
951    REAL(rstd),INTENT(OUT) :: divgrad_i(iim*jjm,1:lle)
952
953    divgrad_i = theta
954    CALL compute_divgrad_inplace(divgrad_i,llb,lle)
955  END SUBROUTINE compute_divgrad
956 
957  SUBROUTINE compute_divgrad_inplace(theta_divgrad_i,llb,lle)
958    USE geometry, ONLY : Ai, ne, le, de
959    INTEGER,INTENT(IN)     :: llb
960    INTEGER,INTENT(IN)     :: lle
961    REAL(rstd),INTENT(INOUT) :: theta_divgrad_i(iim*jjm,1:lle)
962    REAL(rstd)  :: grad_e(3*iim*jjm,llb:lle)
963
964    INTEGER :: ij,l
965
966    ! theta and divgrad_i second dimension is not always llm so use the bounds explicitly
967    !$acc data present(theta_divgrad_i(:,llb:lle), Ai(:), de(:), ne(:,:), le(:)) create(grad_e(:,:)) async
968       
969    !$acc parallel loop collapse(2) async
970    DO l=llb,lle
971      !DIR$ SIMD
972      DO ij=ij_begin_ext,ij_end_ext 
973          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) )
974          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 ))
975          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) )
976      ENDDO
977    ENDDO
978    !$acc end parallel loop
979   
980   
981    !$acc parallel loop collapse(2) async
982    DO l=llb,lle
983      !DIR$ SIMD
984      DO ij=ij_begin,ij_end
985          theta_divgrad_i(ij,l)=1./Ai(ij)*(ne(ij,right)*grad_e(ij+u_right,l)*le(ij+u_right)  +  &
986                             ne(ij,rup)*grad_e(ij+u_rup,l)*le(ij+u_rup)              +  &
987                             ne(ij,lup)*grad_e(ij+u_lup,l)*le(ij+u_lup)              +  &
988                             ne(ij,left)*grad_e(ij+u_left,l)*le(ij+u_left)           +  &
989                             ne(ij,ldown)*grad_e(ij+u_ldown,l)*le(ij+u_ldown)        +  & 
990                             ne(ij,rdown)*grad_e(ij+u_rdown,l)*le(ij+u_rdown))
991      ENDDO
992    ENDDO
993    !$acc end parallel loop
994   
995    !$acc parallel loop collapse(2) async
996    DO l=llb,lle
997      DO ij=ij_begin,ij_end
998          theta_divgrad_i(ij,l)=-theta_divgrad_i(ij,l)*cdivgrad
999      ENDDO
1000    ENDDO
1001    !$acc end parallel loop
1002    !$acc end data
1003  END SUBROUTINE compute_divgrad_inplace
1004
1005  SUBROUTINE compute_gradrot(ue,gradrot_e,llb,lle)
1006    INTEGER,INTENT(IN)     :: llb
1007    INTEGER,INTENT(IN)     :: lle
1008    REAL(rstd),INTENT(IN) :: ue(iim*3*jjm,lle)
1009    REAL(rstd),INTENT(OUT) :: gradrot_e(iim*3*jjm,lle)
1010
1011    gradrot_e = ue
1012    CALL compute_gradrot_inplace(gradrot_e,llb,lle)
1013  END SUBROUTINE compute_gradrot
1014
1015  SUBROUTINE compute_gradrot_inplace(ue_gradrot_e,llb,lle)
1016    USE geometry, ONLY : Av, ne, le, de
1017    INTEGER,INTENT(IN)     :: llb
1018    INTEGER,INTENT(IN)     :: lle
1019    REAL(rstd),INTENT(INOUT) :: ue_gradrot_e(iim*3*jjm,lle)
1020    REAL(rstd) :: rot_v(2*iim*jjm,llb:lle)
1021
1022    INTEGER :: ij,l
1023
1024    ! ue and gradrot_e second dimension is not always llm so use the bounds explicitly
1025    ! gradrot_e should be copyout but using copy instead allows to compare the output
1026    ! more easily as the code sometimes uses unintialed values
1027    !$acc data present(ue_gradrot_e(:,llb:lle), Av(:), ne(:,:), de(:), le(:)) create(rot_v(:,:)) async
1028     
1029    !$acc parallel loop collapse(2) async
1030    DO l=llb,lle
1031!DIR$ SIMD
1032      DO ij=ij_begin_ext,ij_end_ext       
1033          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) &
1034                                + ne(ij+t_rup,left)*ue_gradrot_e(ij+t_rup+u_left,l)*de(ij+t_rup+u_left) &
1035                                - ne(ij,lup)*ue_gradrot_e(ij+u_lup,l)*de(ij+u_lup) )                               
1036          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) &
1037                                     + ne(ij+t_ldown,right)*ue_gradrot_e(ij+t_ldown+u_right,l)*de(ij+t_ldown+u_right) &
1038                                     - ne(ij,rdown)*ue_gradrot_e(ij+u_rdown,l)*de(ij+u_rdown) )
1039      ENDDO
1040    ENDDO                             
1041    !$acc end parallel loop
1042 
1043    !$acc parallel loop collapse(2) async
1044    DO l=llb,lle
1045!DIR$ SIMD
1046      DO ij=ij_begin,ij_end
1047          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))       
1048          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))       
1049          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))
1050      ENDDO
1051    ENDDO
1052    !$acc end parallel loop
1053
1054    !$acc parallel loop collapse(2) async
1055    DO l=llb,lle
1056!DIR$ SIMD
1057      DO ij=ij_begin,ij_end   
1058          ue_gradrot_e(ij+u_right,l)=-ue_gradrot_e(ij+u_right,l)*cgradrot       
1059          ue_gradrot_e(ij+u_lup,l)=-ue_gradrot_e(ij+u_lup,l)*cgradrot
1060          ue_gradrot_e(ij+u_ldown,l)=-ue_gradrot_e(ij+u_ldown,l)*cgradrot       
1061      ENDDO
1062    ENDDO
1063    !$acc end parallel loop
1064    !$acc end data   
1065  END SUBROUTINE compute_gradrot_inplace
1066
1067
1068END MODULE dissip_gcm_mod
Note: See TracBrowser for help on using the repository browser.