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

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

Merge 'mpi_rewrite' into trunk
File size: 35.5 KB
Line 
1MODULE dissip_gcm_mod
2  USE icosa
3  USE abort_mod
4  PRIVATE
5
6  TYPE(t_field),POINTER,SAVE :: f_due_diss_gradiv(:)
7  TYPE(t_field),POINTER,SAVE :: f_due_diss_gradrot(:)
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_gradiv, req_due_gradrot, 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_diss_gradiv,field_u,type_real,llm,ondevice=.TRUE.)
51    CALL allocate_field(f_due_diss_gradrot,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_diss_gradiv,req_e1_vect,req_due_gradiv)
143    CALL init_message(f_due_diss_gradrot,req_e1_vect,req_due_gradrot)
144    CALL init_message(f_dtheta_rhodz_diss,req_i1,req_dtheta)
145
146    tau_graddiv(:)=5000
147    CALL getin("tau_graddiv",tau)
148    tau_graddiv(:)=tau/scale_factor
149
150    CALL getin("nitergdiv",nitergdiv)
151   
152    tau_gradrot(:)=5000
153    CALL getin("tau_gradrot",tau)
154    tau_gradrot(:)=tau/scale_factor
155
156    CALL getin("nitergrot",nitergrot)
157   
158
159    tau_divgrad(:)=5000
160    CALL getin("tau_divgrad",tau)
161    tau_divgrad(:)=tau/scale_factor
162
163    CALL getin("niterdivgrad",niterdivgrad)
164
165
166    cgraddiv=-1
167    cdivgrad=-1
168    cgradrot=-1
169
170!$OMP BARRIER
171!$OMP MASTER
172    DO ind=1,ndomain
173      CALL swap_dimensions(ind)
174      CALL swap_geometry(ind)
175      u=f_u(ind)
176
177! set random seed to get reproductibility when using a different number of process
178      CALL RANDOM_SEED(size=M)
179      CALL RANDOM_SEED(put=(/(i,i=1,M)/))
180
181      ! This cannot be ported on GPU due to compiler limitations
182      DO j=jj_begin,jj_end
183        DO i=ii_begin,ii_end
184          ij=(j-1)*iim+i   
185          CALL RANDOM_NUMBER(r)
186          u(ij+u_right)=r-0.5
187          CALL RANDOM_NUMBER(r)
188          u(ij+u_lup)=r-0.5
189          CALL RANDOM_NUMBER(r)
190          u(ij+u_ldown)=r-0.5
191        ENDDO
192      ENDDO       
193    ENDDO
194!$OMP END MASTER 
195!$OMP BARRIER
196
197    DO it=1,20
198     
199      dumax=0
200      DO iter=1,nitergdiv
201        CALL update_device_field(f_u)
202        CALL transfert_request(f_u,req_e1_vect)
203         
204        DO ind=1,ndomain
205          IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE
206          CALL swap_dimensions(ind)
207          CALL swap_geometry(ind)
208          u=f_u(ind)
209          du=f_du(ind)
210          CALL compute_gradiv_inplace(u,Ai,ne,le,de,1,1)
211          ! This should be ported on GPU but we are running into compiler issues...
212          !$acc update host(u(:)) wait
213          du=u
214        ENDDO
215      ENDDO
216
217      CALL update_device_field(f_du)
218      CALL transfert_request(f_du,req_e1_vect)
219      CALL update_host_field(f_du)
220
221      DO ind=1,ndomain
222        IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE
223        CALL swap_dimensions(ind)
224        CALL swap_geometry(ind)
225        du=f_du(ind)
226         
227        ! Not ported on GPU because all the other kernels cannot be ported
228        DO j=jj_begin,jj_end
229          DO i=ii_begin,ii_end
230            ij=(j-1)*iim+i   
231            if (le(ij+u_right)>1e-100) dumax=MAX(dumax,ABS(du(ij+u_right)))
232            if (le(ij+u_lup)>1e-100)   dumax=MAX(dumax,ABS(du(ij+u_lup)))
233            if (le(ij+u_ldown)>1e-100) dumax=MAX(dumax,ABS(du(ij+u_ldown)))
234          ENDDO
235        ENDDO
236      ENDDO
237
238      IF (using_mpi) THEN
239        CALL reduce_max_omp(dumax,dumax1)
240!$OMP MASTER       
241        CALL MPI_ALLREDUCE(dumax1,dumax,1,MPI_REAL8,MPI_MAX,comm_icosa,ierr)
242!$OMP END MASTER     
243        CALL bcast_omp(dumax) 
244      ELSE
245        CALL allreduce_max_omp(dumax,dumax1)
246        dumax=dumax1
247      ENDIF 
248                       
249      DO ind=1,ndomain
250        IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE
251        CALL swap_dimensions(ind)
252        CALL swap_geometry(ind)
253        u=f_u(ind)
254        du=f_du(ind)
255        ! This should be ported on GPU but we are running into compiler issues...
256        u=du/dumax
257      ENDDO
258      IF (is_master) PRINT *,"gradiv : it :",it ,": dumax",dumax
259
260    ENDDO 
261    IF (is_master) PRINT *,"gradiv : dumax",dumax
262    IF (is_master) PRINT *, 'mean T-cell edge size (km)', 1.45*radius/iim_glo/1000., &
263                              'effective T-cell half-edge size (km)', dumax**(-.5/nitergdiv)/1000
264    IF (is_master)  PRINT *, 'Max. time step assuming c=340 m/s and Courant number=3 (ARK2.3) :', &
265                               3./340.*dumax**(-.5/nitergdiv)
266
267    cgraddiv=dumax**(-1./nitergdiv)
268    IF (is_master) PRINT *,"cgraddiv : ",cgraddiv
269
270!$OMP BARRIER
271!$OMP MASTER
272    DO ind=1,ndomain
273      CALL swap_dimensions(ind)
274      CALL swap_geometry(ind)
275      u=f_u(ind)
276
277! set random seed to get reproductibility when using a different number of process
278      CALL RANDOM_SEED(size=M)
279      CALL RANDOM_SEED(put=(/(i,i=1,M)/))
280 
281      ! This cannot be ported on GPU due to compiler limitations
282       DO j=jj_begin,jj_end
283        DO i=ii_begin,ii_end
284          ij=(j-1)*iim+i   
285          CALL RANDOM_NUMBER(r)
286          u(ij+u_right)=r-0.5
287          CALL RANDOM_NUMBER(r)
288          u(ij+u_lup)=r-0.5
289          CALL RANDOM_NUMBER(r)
290          u(ij+u_ldown)=r-0.5
291        ENDDO
292      ENDDO       
293    ENDDO
294!$OMP END MASTER 
295!$OMP BARRIER
296
297    DO it=1,20
298 
299      dumax=0
300      DO iter=1,nitergrot
301        CALL update_device_field(f_u)
302        CALL transfert_request(f_u,req_e1_vect)
303        DO ind=1,ndomain
304          IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE
305          CALL swap_dimensions(ind)
306          CALL swap_geometry(ind)
307          u=f_u(ind)
308          du=f_du(ind)
309          CALL compute_gradrot_inplace(u,Av,ne,le,de,1,1)
310          ! This should be ported on GPU but we are running into compiler issues...
311          !$acc update host(u(:)) wait
312          du=u
313        ENDDO
314      ENDDO
315     
316      CALL update_device_field(f_du)
317      CALL transfert_request(f_du,req_e1_vect)
318      CALL update_host_field(f_du)
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        du=f_du(ind)
325       
326        ! Not ported on GPU because all the other kernels cannot be ported
327        DO j=jj_begin,jj_end
328          DO i=ii_begin,ii_end
329            ij=(j-1)*iim+i   
330            if (le(ij+u_right)>1e-100) dumax=MAX(dumax,ABS(du(ij+u_right)))
331            if (le(ij+u_lup)>1e-100) dumax=MAX(dumax,ABS(du(ij+u_lup)))
332            if (le(ij+u_ldown)>1e-100) dumax=MAX(dumax,ABS(du(ij+u_ldown)))
333          ENDDO
334        ENDDO
335      ENDDO
336
337      IF (using_mpi) THEN
338        CALL reduce_max_omp(dumax,dumax1)
339!$OMP MASTER       
340        CALL MPI_ALLREDUCE(dumax1,dumax,1,MPI_REAL8,MPI_MAX,comm_icosa,ierr)
341!$OMP END MASTER     
342        CALL bcast_omp(dumax) 
343      ELSE
344        CALL allreduce_max_omp(dumax,dumax1)
345        dumax=dumax1
346      ENDIF 
347
348     
349      DO ind=1,ndomain
350        IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE
351        CALL swap_dimensions(ind)
352        CALL swap_geometry(ind)
353        u=f_u(ind)
354        du=f_du(ind)
355        ! This should be ported on GPU but we are running into compiler issues...
356        u=du/dumax
357      ENDDO
358     
359      IF (is_master) PRINT *,"gradrot : it :",it ,": dumax",dumax
360
361    ENDDO 
362    IF (is_master) PRINT *,"gradrot : dumax",dumax
363 
364    cgradrot=dumax**(-1./nitergrot) 
365    IF (is_master) PRINT *,"cgradrot : ",cgradrot
366   
367
368
369!$OMP BARRIER
370!$OMP MASTER
371    DO ind=1,ndomain
372      CALL swap_dimensions(ind)
373      CALL swap_geometry(ind)
374      theta=f_theta(ind)
375
376! set random seed to get reproductibility when using a different number of process
377      CALL RANDOM_SEED(size=M)
378      CALL RANDOM_SEED(put=(/(i,i=1,M)/))
379
380      ! This cannot be ported on GPU due to compiler limitations
381      DO j=jj_begin,jj_end
382        DO i=ii_begin,ii_end
383          ij=(j-1)*iim+i   
384          CALL RANDOM_NUMBER(r)
385          theta(ij)=r-0.5
386        ENDDO
387      ENDDO 
388    ENDDO
389!$OMP END MASTER
390!$OMP BARRIER
391
392    DO it=1,20
393
394      dthetamax=0
395      DO iter=1,niterdivgrad
396        CALL update_device_field(f_theta)
397        CALL transfert_request(f_theta,req_i1)
398        DO ind=1,ndomain
399          IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE
400          CALL swap_dimensions(ind)
401          CALL swap_geometry(ind)
402          theta=f_theta(ind)
403          dtheta=f_dtheta(ind)
404          CALL compute_divgrad_inplace(theta,Ai,ne,le,de,1,1)
405          ! This should be ported on GPU but we are running into compiler issues...
406          !$acc update host(theta(:)) wait
407          dtheta=theta
408        ENDDO
409      ENDDO
410
411      CALL update_device_field(f_dtheta)
412      CALL transfert_request(f_dtheta,req_i1)
413      CALL update_host_field(f_dtheta)
414     
415      DO ind=1,ndomain
416        IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE
417        CALL swap_dimensions(ind)
418        CALL swap_geometry(ind)
419        dtheta=f_dtheta(ind)
420       
421        ! Not ported on GPU because all the other kernels cannot be ported
422        DO j=jj_begin,jj_end
423          DO i=ii_begin,ii_end
424            ij=(j-1)*iim+i   
425            dthetamax=MAX(dthetamax,ABS(dtheta(ij)))
426          ENDDO
427        ENDDO
428      ENDDO
429
430      IF (using_mpi) THEN
431        CALL reduce_max_omp(dthetamax ,dthetamax1)
432!$OMP MASTER       
433        CALL MPI_ALLREDUCE(dthetamax1,dthetamax,1,MPI_REAL8,MPI_MAX,comm_icosa,ierr)
434!$OMP END MASTER     
435        CALL bcast_omp(dthetamax) 
436      ELSE
437        CALL allreduce_max_omp(dthetamax,dthetamax1)
438        dumax=dumax1
439      ENDIF 
440     
441      IF (is_master) PRINT *,"divgrad : it :",it ,": dthetamax",dthetamax
442
443      DO ind=1,ndomain
444        IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE
445        CALL swap_dimensions(ind)
446        CALL swap_geometry(ind)
447        theta=f_theta(ind)
448        dtheta=f_dtheta(ind)
449        ! This should be ported on GPU but we are running into compiler issues...
450        theta=dtheta/dthetamax
451      ENDDO
452    ENDDO 
453
454    IF (is_master) PRINT *,"divgrad : divgrad",dthetamax
455 
456    cdivgrad=dthetamax**(-1./niterdivgrad) 
457    IF (is_master) PRINT *,"cdivgrad : ",cdivgrad
458
459    ! Added CMIP5 dissipation vertical profile (SF, 19/09/18)
460    vert_prof_dissip = 1
461    CALL getin('vert_prof_dissip',vert_prof_dissip)
462    IF (vert_prof_dissip == 1) then
463       dissip_zref = 30.
464       dissip_deltaz = 10.
465       dissip_factz = 4.
466       CALL getin('dissip_factz',dissip_factz )
467       CALL getin('dissip_deltaz',dissip_deltaz )
468       CALL getin('dissip_zref',dissip_zref )
469       DO l=1,llm
470          pseudoz=8.*log(preff/presnivs(l))
471          zvert=1+ &
472               (tanh((pseudoz-dissip_zref)/dissip_deltaz)+1.)/2. &
473               *(dissip_factz-1.)
474          tau_graddiv(l) = tau_graddiv(l)/zvert
475          tau_gradrot(l) = tau_gradrot(l)/zvert
476          tau_divgrad(l) = tau_divgrad(l)/zvert
477       ENDDO
478    ELSE     
479       fact=2
480       DO l=1,llm
481          IF(ap_bp_present) THEN ! height-dependent dissipation
482             zz= 1. - preff/presnivs(l)
483          ELSE
484             zz = 0.
485          END IF
486          zvert=fact-(fact-1)/(1+zz*zz)
487          tau_graddiv(l) = tau_graddiv(l)/zvert
488          tau_gradrot(l) = tau_gradrot(l)/zvert
489          tau_divgrad(l) = tau_divgrad(l)/zvert
490       ENDDO
491    ENDIF
492
493    mintau=tau_graddiv(1)
494    DO l=1,llm
495      mintau=MIN(mintau,tau_graddiv(l))
496      mintau=MIN(mintau,tau_gradrot(l))
497      mintau=MIN(mintau,tau_divgrad(l))
498    ENDDO
499
500    IF(rayleigh_friction_type>0) mintau=MIN(mintau,rayleigh_tau)
501
502    IF(mintau>0) THEN
503       IF (itau_dissip==0) THEN
504         IF (is_master) PRINT *,"init_dissip: Automatic computation of itau_dissip..."
505         itau_dissip=INT(mintau/dt)
506       ENDIF
507    ELSE
508       IF (is_master) PRINT *,"init_dissip: No dissipation time set, setting itau_dissip to 1000000000"
509       itau_dissip=100000000
510    END IF
511    itau_dissip=MAX(1,itau_dissip)
512    dtdissip=itau_dissip*dt
513    IF (is_master) THEN
514      PRINT *,"init_dissip: rayleigh_tau",rayleigh_tau, "mintau ",mintau
515      PRINT *,"init_dissip: itau_dissip",itau_dissip," dtdissip ",dtdissip
516    ENDIF
517
518    !$acc update device(tau_graddiv(:),tau_gradrot(:),tau_divgrad(:)) async
519
520  END SUBROUTINE init_dissip 
521 
522 
523  SUBROUTINE dissip(f_ps,f_mass,f_phis,f_geopot,f_theta_rhodz,f_ue, f_dtheta_rhodz,f_due)
524  USE icosa
525  USE theta2theta_rhodz_mod
526  USE pression_mod
527  USE exner_mod
528  USE geopotential_mod
529  USE trace
530  USE time_mod
531  USE omp_para
532  USE abort_mod
533  IMPLICIT NONE
534    TYPE(t_field),POINTER :: f_ps(:), f_mass(:), f_phis(:), f_geopot(:)
535    TYPE(t_field),POINTER :: f_theta_rhodz(:), f_dtheta_rhodz(:)
536    TYPE(t_field),POINTER :: f_ue(:), f_due(:)
537
538    REAL(rstd),POINTER         :: due(:,:)
539    REAL(rstd),POINTER         :: phi(:,:), ue(:,:)
540    REAL(rstd),POINTER         :: due_diss1(:,:)
541    REAL(rstd),POINTER         :: due_diss2(:,:)
542    REAL(rstd),POINTER         :: dtheta_rhodz(:,:,:)
543    REAL(rstd),POINTER         :: dtheta_rhodz_diss(:,:)
544
545    INTEGER :: ind, shear
546    INTEGER :: l,ij,nn
547
548!$OMP BARRIER
549   
550    CALL trace_start("dissip")
551    CALL gradiv(f_ue,f_due_diss_gradiv)
552    CALL gradrot(f_ue,f_due_diss_gradrot)
553
554    CALL divgrad_theta_rhodz(f_mass,f_theta_rhodz,f_dtheta_rhodz_diss)
555   
556    DO ind=1,ndomain
557      IF (.NOT. assigned_domain(ind)) CYCLE
558      CALL swap_dimensions(ind)
559      CALL swap_geometry(ind)
560      due=f_due(ind) 
561      due_diss1=f_due_diss_gradiv(ind)
562      due_diss2=f_due_diss_gradrot(ind)
563      dtheta_rhodz=f_dtheta_rhodz(ind)
564      dtheta_rhodz_diss=f_dtheta_rhodz_diss(ind)
565
566      !$acc parallel loop collapse(2) present(due(:,:), dtheta_rhodz(:,:,:), due_diss1(:,:), due_diss2(:,:), dtheta_rhodz_diss(:,:), tau_graddiv(:), tau_gradrot(:), tau_divgrad(:)) async
567      DO l=ll_begin,ll_end
568!DIR$ SIMD
569        DO ij=ij_begin,ij_end
570
571            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 
572            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
573            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
574
575            dtheta_rhodz(ij,l,1) = -0.5*dtheta_rhodz_diss(ij,l)/tau_divgrad(l)*itau_dissip
576        ENDDO
577      ENDDO
578      !$acc end parallel loop
579
580!      dtheta_rhodz=0
581!      due=0
582
583      IF(rayleigh_friction_type>0) THEN
584        CALL abort_acc("dissip/rayleigh_friction_type>0")
585       IF(rayleigh_friction_type<99) THEN
586         phi=f_geopot(ind)
587         ue=f_ue(ind)
588         DO l=ll_begin,ll_end
589            DO ij=ij_begin,ij_end
590               CALL relax(t_right, u_right)
591               CALL relax(t_lup,   u_lup)
592               CALL relax(t_ldown, u_ldown)
593            ENDDO
594         END DO
595       ELSE
596         ue=f_ue(ind)
597            !$acc parallel loop present(ue(:,:), due(:,:), lat_e(:))
598            DO ij=ij_begin,ij_end
599              nn = ij+u_right
600              IF (ABS(lat_e(nn)) .gt. rayleigh_limlat) THEN
601                !print*, "latitude", lat_e(nn)*180./3.14159
602                due(nn,ll_begin:ll_begin+1) = due(nn,ll_begin:ll_begin+1) - (ue(nn,ll_begin:ll_begin+1)/rayleigh_tau)
603              ENDIF
604              nn = ij+u_lup
605              IF (ABS(lat_e(nn)) .gt. rayleigh_limlat) THEN
606                due(nn,ll_begin:ll_begin+1) = due(nn,ll_begin:ll_begin+1) - (ue(nn,ll_begin:ll_begin+1)/rayleigh_tau)
607              ENDIF
608              nn = ij+u_ldown
609              IF (ABS(lat_e(nn)) .gt. rayleigh_limlat) THEN
610                due(nn,ll_begin:ll_begin+1) = due(nn,ll_begin:ll_begin+1) - (ue(nn,ll_begin:ll_begin+1)/rayleigh_tau)
611              ENDIF
612            ENDDO
613            !$acc end parallel loop
614       ENDIF
615      END IF
616   END DO
617   CALL trace_end("dissip")
618
619   !CALL write_dissip_tendencies
620!$OMP BARRIER
621
622    CONTAINS
623
624      SUBROUTINE relax(shift_t, shift_u)
625        USE dcmip_initial_conditions_test_1_2_3
626        REAL(rstd) :: z, ulon,ulat, & ! input to test2_schaer_mountain
627             p,hyam,hybm,w,t,phis,ps,rho,q, &   ! unused input/output to test2_schaer_mountain
628             fz, u3d(3), uref
629        REAL(rstd), PARAMETER :: zh=2e4,ztop=3e4  ! DCMIP values
630        LOGICAL :: hybrid_eta
631        INTEGER :: shift_u, shift_t, zcoords, nn
632        z = (phi(ij,l)+phi(ij+shift_t,l))/(2.*g)
633        IF(z>zh) THEN  ! relax only in the sponge layer z>zh
634
635           nn = ij+shift_u
636           zcoords = 1 ; hybrid_eta = .FALSE. ! use z instead of p or hyam/hybm
637           CALL test2_schaer_mountain(lon_e(nn),lat_e(nn),p,z,zcoords,hybrid_eta, &
638                hyam,hybm,shear,ulon,ulat,w,t,phis,ps,rho,q)
639!           u3d = ulon*elon_e(nn,:) + ulat*elat_e(nn,:)
640           u3d = ulon*elon_e(nn,:) ! ulat=0
641           uref = sum(u3d*ep_e(nn,:))
642
643           fz = sin((pi/2)*(z-zh)/(ztop-zh))
644           fz = fz*fz/rayleigh_tau
645           due(nn,l) = due(nn,l) - itau_dissip*fz*(ue(nn,l)-uref)
646!           fz = 1./rayleigh_tau
647!           due(nn,l) = due(nn,l) - fz*ue(nn,l)
648         END IF
649       END SUBROUTINE relax
650     
651       SUBROUTINE write_dissip_tendencies
652         USE observable_mod, ONLY : f_buf_ulon, f_buf_ulat
653         USE wind_mod
654         USE output_field_mod
655
656         CALL transfert_request(f_due_diss_gradiv,req_e1_vect)
657         CALL un2ulonlat(f_due_diss_gradiv, f_buf_ulon, f_buf_ulat, (1./(tau_graddiv(1))))
658         CALL output_field("dulon_diss1",f_buf_ulon)
659         CALL output_field("dulat_diss1",f_buf_ulat)
660!
661         CALL transfert_request(f_due_diss_gradrot,req_e1_vect)
662         CALL un2ulonlat(f_due_diss_gradrot, f_buf_ulon, f_buf_ulat, (1./(tau_graddiv(1))))
663         CALL output_field("dulon_diss2",f_buf_ulon)
664         CALL output_field("dulat_diss2",f_buf_ulat)
665       END SUBROUTINE write_dissip_tendencies
666
667  END SUBROUTINE dissip
668
669  SUBROUTINE gradiv(f_ue,f_due)
670  USE icosa
671  USE trace
672  USE omp_para
673  IMPLICIT NONE
674    TYPE(t_field),POINTER :: f_ue(:)
675    TYPE(t_field),POINTER :: f_due(:)
676    REAL(rstd),POINTER    :: ue(:,:)
677    REAL(rstd),POINTER    :: due(:,:)
678    INTEGER :: ind
679    INTEGER :: it,l,ij
680       
681    CALL trace_start("gradiv")
682
683    DO ind=1,ndomain
684      IF (.NOT. assigned_domain(ind)) CYCLE
685      CALL swap_dimensions(ind)
686      CALL swap_geometry(ind)
687      ue=f_ue(ind)
688      due=f_due(ind) 
689      !$acc parallel loop present(ue(:,:), due(:,:)) async
690      DO  l = ll_begin, ll_end
691        !$acc loop
692!DIR$ SIMD
693        DO ij=ij_begin,ij_end
694             due(ij+u_right,l)=ue(ij+u_right,l)
695             due(ij+u_lup,l)=ue(ij+u_lup,l)
696             due(ij+u_ldown,l)=ue(ij+u_ldown,l)
697        ENDDO
698      ENDDO
699    ENDDO
700
701    DO it=1,nitergdiv
702      CALL send_message(f_due,req_due_gradiv)
703      CALL wait_message(req_due_gradiv)
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_gradrot)
752      CALL wait_message(req_due_gradrot)
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.