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

Last change on this file since 849 was 845, checked in by dubos, 5 years ago

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