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

Last change on this file since 915 was 856, checked in by dubos, 5 years ago

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