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

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

devel : added safeguard for Courant criterion for horizontal dissipation
File size: 32.0 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       CALL dissip_timescale
129    END IF
130
131  END SUBROUTINE init_dissip
132
133  SUBROUTINE dissip_constants
134    ! The SAVE attribute is required here, otherwise allocate_field will not work with OpenMP
135    TYPE(t_field),POINTER, SAVE  :: f_u(:)
136    TYPE(t_field),POINTER, SAVE  :: f_du(:)
137    TYPE(t_field),POINTER, SAVE  :: f_theta(:)
138    TYPE(t_field),POINTER, SAVE  :: f_dtheta(:)
139    REAL(rstd),POINTER :: u(:)
140    REAL(rstd),POINTER :: du(:)
141    REAL(rstd),POINTER :: theta(:)
142    REAL(rstd),POINTER :: dtheta(:)
143    REAL(rstd) :: dumax, dthetamax
144    INTEGER :: it, iter, ind
145    CALL allocate_field(f_u,field_u,type_real)
146    CALL allocate_field(f_du,field_u,type_real)
147    CALL allocate_field(f_theta,field_t,type_real)
148    CALL allocate_field(f_dtheta,field_t,type_real)
149
150    PRINT *, '=========', SHAPE(f_u)
151
152    !--------------------------- Compute cgraddiv ----------------------------
153    cgraddiv=-1.
154    CALL random_vel(f_u)   
155    DO it=1,20
156       DO iter=1,nitergdiv
157          CALL transfert_request(f_u,req_e1_vect)
158          DO ind=1,ndomain
159             IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE
160             CALL swap_dimensions(ind)
161             CALL swap_geometry(ind)
162             u=f_u(ind)
163             du=f_du(ind)
164             CALL compute_gradiv(u,du,1,1)
165          ENDDO
166       ENDDO
167       CALL transfert_request(f_du,req_e1_vect)       
168       dumax = max_vel(f_du)       
169       CALL rescale_field(1./dumax, f_du, f_u)
170       IF (is_master) PRINT *,"gradiv : it :",it ,": dumax",dumax
171    ENDDO
172
173    cgraddiv=dumax**(-1./nitergdiv)
174    IF (is_master) THEN
175       PRINT *,"gradiv : dumax",dumax
176       PRINT *, 'mean T-cell edge size (km)', 1.45*radius/iim_glo/1000., &
177         'effective T-cell half-edge size (km)', dumax**(-.5/nitergdiv)/1000
178       PRINT *, 'Max. time step assuming c=340 m/s and Courant number=3 (ARK2.3) :', &
179         3./340.*dumax**(-.5/nitergdiv)
180       PRINT *,"cgraddiv : ",cgraddiv
181    END IF
182   
183    !----------------- Compute cgradrot --------------------
184    cgradrot=-1.
185    CALL random_vel(f_u)
186    DO it=1,20
187       DO iter=1,nitergrot
188          CALL transfert_request(f_u,req_e1_vect)
189          DO ind=1,ndomain
190             IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE
191             CALL swap_dimensions(ind)
192             CALL swap_geometry(ind)
193             u=f_u(ind)
194             du=f_du(ind)
195             CALL compute_gradrot(u,du,1,1)
196             u=du
197          ENDDO
198       ENDDO
199       CALL transfert_request(f_du,req_e1_vect)
200       dumax = max_vel(f_du)
201       CALL rescale_field(1./dumax, f_du, f_u)
202       IF (is_master) PRINT *,"gradrot : it :",it ,": dumax",dumax
203    ENDDO
204
205    cgradrot=dumax**(-1./nitergrot) 
206    IF (is_master) PRINT *,"gradrot : dumax",dumax 
207    IF (is_master) PRINT *,"cgradrot : ",cgradrot
208   
209    !----------------- Compute cdivgrad --------------------
210
211    cdivgrad=-1.
212    CALL random_scalar(f_theta)   
213    DO it=1,20       
214       DO iter=1,niterdivgrad
215          CALL transfert_request(f_theta,req_i1)
216          DO ind=1,ndomain
217             IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE
218             CALL swap_dimensions(ind)
219             CALL swap_geometry(ind)
220             theta=f_theta(ind)
221             dtheta=f_dtheta(ind)
222             CALL compute_divgrad(theta,dtheta,1,1)
223          ENDDO
224       ENDDO       
225       CALL transfert_request(f_dtheta,req_i1)
226       dthetamax = max_scalar(f_dtheta)     
227       IF (is_master) PRINT *,"divgrad : it :",it ,": dthetamax",dthetamax
228       CALL rescale_field(1./dthetamax, f_dtheta, f_theta)
229    END DO
230   
231    cdivgrad=dthetamax**(-1./niterdivgrad) 
232    IF (is_master) PRINT *,"divgrad : divgrad",dthetamax
233    IF (is_master) PRINT *,"cdivgrad : ",cdivgrad
234
235  END SUBROUTINE dissip_constants
236
237  SUBROUTINE dissip_profile
238    USE disvert_mod
239    ! parameters used by the various profiles
240    ! IF planet_type == earth
241    !    IF dissip_vert_prof == 0
242    !      none
243    !    IF dissip_vert_prof == 1
244    !      dissip_zref, dissip_deltaz, dissip_factz
245    ! IF planet_type == other
246    !    IF dissip_vert_prof == 0
247    !      dissip_fac_mid
248    !      + dissip_deltaz, dissip_hdelta, dissip_fac_up, dissip_pupstart IF ok_strato     
249    !    IF dissip_vert_prof == 1
250    !      fac_mid, fac_up, startalt, delta => middle (hardcoded), scaleheight
251    !     
252
253    REAL(rstd), PARAMETER :: fact=2., &
254         fac_mid=3.,   & ! coefficient for lower/middle atmosphere
255         fac_up=30.,   & ! coefficient for upper atmosphere
256         startalt=70., & ! altitude (in km) for the transition from middle to upper atm.
257         delta=30.,    & ! Size (in km) of the transition region between middle/upper atm.
258         middle=startalt+delta/2.
259    REAL(rstd) :: dissip_zref, dissip_deltaz, dissip_factz, &
260         dissip_hdelta, dissip_fac_up, dissip_fac_mid, dissip_pupstart, &
261         scaleheight, &
262         zz, pseudoz, pup, &
263         sigma(llm), zvert(llm)
264    CHARACTER(LEN=255) :: planet_type ! earth, other, other_strato ; other_strato = other + ok_strato
265    LOGICAL :: ok_strato
266    INTEGER    :: l, dissip_vert_prof
267
268    ! select vertical profile of horizontal dissipation coefficients, see also [781]
269
270    planet_type='earth'
271    CALL getin('dissip_planet_type', planet_type)
272    SELECT CASE(planet_type)
273    CASE('earth','other')
274       ok_strato=.FALSE.
275    CASE('other_strato')
276       planet_type='other'
277       ok_strato=.TRUE.
278    CASE DEFAULT
279       STOP "Invalid value of dissip_planet_type, valid values are <earth>, <other>, <other_strato>"
280    END SELECT
281
282    dissip_vert_prof = 0
283    CALL getin('dissip_vert_prof',dissip_vert_prof)
284
285    SELECT CASE(dissip_vert_prof)
286    CASE(0)
287       IF(TRIM(planet_type)=='other') THEN
288          ! Default values given below are for a Venus-like planet
289          CALL getin('dissip_fac_mid',dissip_fac_mid )
290          dissip_fac_mid=2.
291          IF(ok_strato) THEN
292             dissip_fac_up=50.
293             ! deltaz et hdelta in km
294             dissip_deltaz=30.
295             dissip_hdelta=5.
296             ! pupstart in Pa
297             dissip_pupstart=1.e4         
298             CALL getin('dissip_deltaz',dissip_deltaz )
299             CALL getin('dissip_hdelta',dissip_hdelta )
300             CALL getin('dissip_fac_up',dissip_fac_up )
301             CALL getin('dissip_pupstart',dissip_pupstart)
302          END IF
303       END IF
304    CASE(1)
305       IF(TRIM(planet_type)=='earth') THEN
306          dissip_zref = 30.
307          dissip_deltaz = 10.
308          dissip_factz = 4.
309          CALL getin('dissip_zref',dissip_zref )
310          CALL getin('dissip_deltaz',dissip_deltaz )
311          CALL getin('dissip_factz',dissip_factz )
312       ELSE
313          ! fac_mid, fac_up, startalt, delta => middle are hardcoded
314          CALL getin('dissip_scaleheight', scaleheight)
315       END IF
316    CASE DEFAULT
317       STOP 'Invalid value of dissip_vert_prof : valid values are 0,1'
318    END SELECT
319
320    IF(ap_bp_present) THEN
321       sigma(:) = preff/presnivs(:)
322    ELSE
323       sigma(:) = 1.
324    END IF
325
326    SELECT CASE(TRIM(planet_type))
327    CASE('earth')
328       DO l=1,llm
329          IF(dissip_vert_prof == 1) THEN
330             pseudoz=8.*LOG(sigma(l))
331             zvert(l)=1+ &
332                  (TANH((pseudoz-dissip_zref)/dissip_deltaz)+1.)/2. &
333                  *(dissip_factz-1.)
334          ELSE
335             zz = 1.-sigma(l)
336             zvert(l)= fact -( fact-1.)/( 1.+zz*zz )
337          END IF
338       END DO
339    CASE('other')
340       SELECT CASE(dissip_vert_prof)
341       CASE(0)
342          ! First step: going from 1 to dissip_fac_mid (in gcm.def)
343          !============
344          DO l=1,llm
345             zz       = 1. - sigma(l)
346             zvert(l) = dissip_fac_mid -( dissip_fac_mid-1.)/( 1.+zz*zz )
347          END DO
348          !
349          ! Second step if ok_strato:  from dissip_fac_mid to dissip_fac_up (in gcm.def)
350          !==========================
351          ! Utilisation de la fonction tangente hyperbolique pour augmenter
352          ! arbitrairement la dissipation et donc la stabilite du modele a
353          ! partir d'une certaine altitude.
354          !
355          !   Le facteur multiplicatif de basse atmosphere etant deja pris
356          !   en compte, il faut diviser le facteur multiplicatif de haute
357          !   atmosphere par celui-ci.
358
359          IF(ok_strato) THEN
360             Pup = dissip_pupstart*exp(-0.5*dissip_deltaz/dissip_hdelta)
361             DO l=1,llm
362                zvert(l)= zvert(l)*(1.0+( (dissip_fac_up/dissip_fac_mid-1) &
363                     *(1.-(0.5*(1+tanh(-6./dissip_deltaz*              &
364                     (-dissip_hdelta*log(presnivs(l)/Pup))  ))))  ))
365             END DO
366          END IF
367       CASE(1)
368          DO l=1,llm
369             zz      = 1. - sigma(l)
370             zvert(l)= fac_mid -( fac_mid-1.)/( 1.+zz*zz )     
371             zvert(l)= zvert(l)*(1.0+((fac_up/fac_mid-1)*    &
372                  (1.-(0.5*(1+tanh(-6./                 &
373                  delta*(scaleheight*(-log(presnivs(l)/preff))-middle))))) &
374                  ))
375          END DO
376       END SELECT
377    END SELECT
378
379    IF(is_master) PRINT *, 'vertical profile of horizontal dissipation : ', zvert(:)
380
381    DO l=1,llm
382       tau_graddiv(l) = tau_graddiv(l)/zvert(l)
383       tau_gradrot(l) = tau_gradrot(l)/zvert(l)
384       tau_divgrad(l) = tau_divgrad(l)/zvert(l)
385    END DO
386  END SUBROUTINE dissip_profile
387 
388  SUBROUTINE dissip_timescale   
389    INTEGER :: l
390    REAL(rstd) :: mintau
391    mintau=tau_graddiv(1)
392    DO l=1,llm
393       mintau=MIN(mintau,tau_graddiv(l))
394       mintau=MIN(mintau,tau_gradrot(l))
395       mintau=MIN(mintau,tau_divgrad(l))
396    END DO
397   
398    IF(rayleigh_friction_type>0) mintau=MIN(mintau,rayleigh_tau)
399   
400    IF(mintau>0) THEN
401       IF (itau_dissip==0) THEN
402          IF (is_master) PRINT *,"init_dissip: Automatic computation of itau_dissip..."
403          itau_dissip=INT(mintau/dt)
404       ENDIF
405    ELSE
406       IF (is_master) PRINT *,"init_dissip: No dissipation time set, setting itau_dissip to 1000000000"
407       itau_dissip=100000000
408    END IF
409    itau_dissip=MAX(1,itau_dissip)
410    dtdissip=itau_dissip*dt
411
412    IF (is_master) THEN
413       PRINT *,"init_dissip: rayleigh_tau",rayleigh_tau, "mintau ",mintau
414       PRINT *,"init_dissip: itau_dissip",itau_dissip," dtdissip ",dtdissip
415    ENDIF
416   
417    IF (dtdissip>2.*mintau) THEN
418       IF(is_master) PRINT *, 'The CFL condition for dissipation dtdissip<2*mintau is violated : dtdissip, mintau ', dtdissip, mintau
419       STOP
420    END IF
421
422  END SUBROUTINE dissip_timescale
423
424  SUBROUTINE dissip(f_ps,f_mass,f_phis,f_geopot,f_theta_rhodz,f_ue, f_dtheta_rhodz,f_due)
425    TYPE(t_field),POINTER :: f_ps(:), f_mass(:), f_phis(:), f_geopot(:)
426    TYPE(t_field),POINTER :: f_theta_rhodz(:), f_dtheta_rhodz(:)
427    TYPE(t_field),POINTER :: f_ue(:), f_due(:)
428
429    REAL(rstd),POINTER         :: due(:,:)
430    REAL(rstd),POINTER         :: phi(:,:), ue(:,:)
431    REAL(rstd),POINTER         :: due_diss1(:,:)
432    REAL(rstd),POINTER         :: due_diss2(:,:)
433    REAL(rstd),POINTER         :: dtheta_rhodz(:,:,:)
434    REAL(rstd),POINTER         :: dtheta_rhodz_diss(:,:)
435
436    INTEGER :: ind, l,ij,nn
437
438!$OMP BARRIER
439   
440    CALL trace_start("dissip")
441    CALL gradiv(f_ue,f_due_diss1)
442    CALL gradrot(f_ue,f_due_diss2)
443    CALL divgrad_theta_rhodz(f_mass,f_theta_rhodz,f_dtheta_rhodz_diss)
444   
445    DO ind=1,ndomain
446      IF (.NOT. assigned_domain(ind)) CYCLE
447      CALL swap_dimensions(ind)
448      CALL swap_geometry(ind)
449      due=f_due(ind) 
450      due_diss1=f_due_diss1(ind)
451      due_diss2=f_due_diss2(ind)
452      dtheta_rhodz=f_dtheta_rhodz(ind)
453      dtheta_rhodz_diss=f_dtheta_rhodz_diss(ind)
454
455      DO l=ll_begin,ll_end
456!DIR$ SIMD
457        DO ij=ij_begin,ij_end
458            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 
459            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
460            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
461            dtheta_rhodz(ij,l,1) = -0.5*dtheta_rhodz_diss(ij,l)/tau_divgrad(l)*itau_dissip
462        ENDDO
463      ENDDO
464     
465      IF(rayleigh_friction_type>0) THEN
466         IF(rayleigh_friction_type<99) THEN
467            phi=f_geopot(ind)
468            ue=f_ue(ind)
469            DO l=ll_begin,ll_end
470               DO ij=ij_begin,ij_end
471                  CALL relax(t_right, u_right)
472                  CALL relax(t_lup,   u_lup)
473                  CALL relax(t_ldown, u_ldown)
474               ENDDO
475            END DO
476         ELSE
477            ue=f_ue(ind)
478            DO ij=ij_begin,ij_end
479               nn = ij+u_right
480               IF (ABS(lat_e(nn)) .gt. rayleigh_limlat) THEN
481                  !print*, "latitude", lat_e(nn)*180./3.14159
482                  due(nn,ll_begin:ll_begin+1) = due(nn,ll_begin:ll_begin+1) - (ue(nn,ll_begin:ll_begin+1)/rayleigh_tau)
483               ENDIF
484               nn = ij+u_lup
485               IF (ABS(lat_e(nn)) .gt. rayleigh_limlat) THEN
486                  due(nn,ll_begin:ll_begin+1) = due(nn,ll_begin:ll_begin+1) - (ue(nn,ll_begin:ll_begin+1)/rayleigh_tau)
487               ENDIF
488               nn = ij+u_ldown
489               IF (ABS(lat_e(nn)) .gt. rayleigh_limlat) THEN
490                  due(nn,ll_begin:ll_begin+1) = due(nn,ll_begin:ll_begin+1) - (ue(nn,ll_begin:ll_begin+1)/rayleigh_tau)
491               ENDIF
492            ENDDO
493         ENDIF
494      END IF
495   END DO
496   
497   CALL trace_end("dissip")
498   
499   CALL write_dissip_tendencies
500!$OMP BARRIER
501
502 CONTAINS
503   
504   SUBROUTINE relax(shift_t, shift_u)
505     USE dcmip_initial_conditions_test_1_2_3
506     REAL(rstd) :: z, ulon,ulat, lon,lat, & ! input to test2_schaer_mountain
507          p,hyam,hybm,w,t,phis,ps,rho,q, &   ! unused input/output to test2_schaer_mountain
508          fz, u3d(3), uref
509     REAL(rstd), PARAMETER :: zh=2e4,ztop=3e4  ! DCMIP values
510     LOGICAL :: hybrid_eta
511     INTEGER :: shift_u, shift_t, zcoords, nn
512     z = (phi(ij,l)+phi(ij+shift_t,l))/(2.*g)
513     IF(z>zh) THEN  ! relax only in the sponge layer z>zh
514        nn = ij+shift_u
515        zcoords = 1 ; hybrid_eta = .FALSE. ! use z instead of p or hyam/hybm
516        CALL test2_schaer_mountain(lon_e(nn),lat_e(nn),p,z,zcoords,hybrid_eta, &
517             hyam,hybm,rayleigh_shear,ulon,ulat,w,t,phis,ps,rho,q)
518        u3d = ulon*elon_e(nn,:) ! ulat=0
519        uref = sum(u3d*ep_e(nn,:))
520       
521        fz = sin((pi/2)*(z-zh)/(ztop-zh))
522        fz = fz*fz/rayleigh_tau
523        due(nn,l) = due(nn,l) - itau_dissip*fz*(ue(nn,l)-uref)
524     END IF
525   END SUBROUTINE relax
526   
527   SUBROUTINE write_dissip_tendencies
528     USE observable_mod, ONLY : f_buf_ulon, f_buf_ulat
529     USE wind_mod
530     USE output_field_mod
531     
532     CALL transfert_request(f_due_diss1,req_e1_vect)
533     CALL un2ulonlat(f_due_diss1, f_buf_ulon, f_buf_ulat, (1./(tau_graddiv(1))))
534     CALL output_field("dulon_diss1",f_buf_ulon)
535     CALL output_field("dulat_diss1",f_buf_ulat)
536     !
537     CALL transfert_request(f_due_diss2,req_e1_vect)
538     CALL un2ulonlat(f_due_diss2, f_buf_ulon, f_buf_ulat, (1./(tau_graddiv(1))))
539     CALL output_field("dulon_diss2",f_buf_ulon)
540     CALL output_field("dulat_diss2",f_buf_ulat)
541   END SUBROUTINE write_dissip_tendencies
542   
543 END SUBROUTINE dissip
544
545 
546 SUBROUTINE gradiv(f_ue,f_due)
547   TYPE(t_field), POINTER :: f_ue(:)
548   TYPE(t_field), POINTER :: f_due(:)
549   REAL(rstd),POINTER :: ue(:,:)
550   REAL(rstd),POINTER :: due(:,:)
551   INTEGER :: ind
552   INTEGER :: it,l,ij
553   
554   CALL trace_start("gradiv")
555   
556   DO ind=1,ndomain
557      IF (.NOT. assigned_domain(ind)) CYCLE
558      CALL swap_dimensions(ind)
559      CALL swap_geometry(ind)
560      ue=f_ue(ind)
561      due=f_due(ind) 
562      DO  l = ll_begin, ll_end
563         !DIR$ SIMD
564         DO ij=ij_begin,ij_end
565            due(ij+u_right,l)=ue(ij+u_right,l)
566            due(ij+u_lup,l)=ue(ij+u_lup,l)
567            due(ij+u_ldown,l)=ue(ij+u_ldown,l)
568         ENDDO
569      ENDDO
570   ENDDO
571   
572   DO it=1,nitergdiv     
573      CALL send_message(f_due,req_due)
574      CALL wait_message(req_due)
575      DO ind=1,ndomain
576         IF (.NOT. assigned_domain(ind)) CYCLE
577         CALL swap_dimensions(ind)
578         CALL swap_geometry(ind)
579         due=f_due(ind) 
580         CALL compute_gradiv(due,due,ll_begin,ll_end)
581      ENDDO
582   ENDDO
583   
584   CALL trace_end("gradiv")
585 END SUBROUTINE gradiv
586 
587 
588 SUBROUTINE gradrot(f_ue,f_due)
589   TYPE(t_field), POINTER :: f_ue(:)
590    TYPE(t_field), POINTER :: f_due(:)
591    REAL(rstd),POINTER :: ue(:,:)
592    REAL(rstd),POINTER :: due(:,:)
593    INTEGER :: ind, it,l,ij
594    CALL trace_start("gradrot")
595
596    DO ind=1,ndomain
597      IF (.NOT. assigned_domain(ind)) CYCLE
598      CALL swap_dimensions(ind)
599      CALL swap_geometry(ind)
600      ue=f_ue(ind)
601      due=f_due(ind) 
602      DO  l = ll_begin, ll_end
603!DIR$ SIMD
604        DO ij=ij_begin,ij_end
605             due(ij+u_right,l)=ue(ij+u_right,l)
606             due(ij+u_lup,l)=ue(ij+u_lup,l)
607             due(ij+u_ldown,l)=ue(ij+u_ldown,l)
608        ENDDO
609      ENDDO
610    ENDDO
611
612    DO it=1,nitergrot
613       CALL send_message(f_due,req_due)
614       CALL wait_message(req_due)
615       DO ind=1,ndomain
616          IF (.NOT. assigned_domain(ind)) CYCLE
617          CALL swap_dimensions(ind)
618          CALL swap_geometry(ind)
619          due=f_due(ind) 
620          CALL compute_gradrot(due,due,ll_begin,ll_end)
621       ENDDO
622    ENDDO
623   
624    CALL trace_end("gradrot")
625  END SUBROUTINE gradrot
626 
627  SUBROUTINE divgrad(f_theta,f_dtheta)
628    TYPE(t_field), POINTER :: f_theta(:)
629    TYPE(t_field), POINTER :: f_dtheta(:)
630    REAL(rstd),POINTER :: theta(:,:)
631    REAL(rstd),POINTER :: dtheta(:,:)
632    INTEGER :: ind, it
633
634    CALL trace_start("divgrad")
635       
636    DO ind=1,ndomain
637      IF (.NOT. assigned_domain(ind)) CYCLE
638      CALL swap_dimensions(ind)
639      CALL swap_geometry(ind)
640      theta=f_theta(ind)
641      dtheta=f_dtheta(ind) 
642      dtheta=theta
643    ENDDO
644
645    DO it=1,niterdivgrad       
646      CALL transfert_request(f_dtheta,req_i1)       
647      DO ind=1,ndomain
648        IF (.NOT. assigned_domain(ind)) CYCLE
649        CALL swap_dimensions(ind)
650        CALL swap_geometry(ind)
651        dtheta=f_dtheta(ind) 
652        CALL compute_divgrad(dtheta,dtheta,ll_begin,ll_end)
653      ENDDO
654    ENDDO
655
656    CALL trace_end("divgrad")
657  END SUBROUTINE divgrad
658   
659  SUBROUTINE divgrad_theta_rhodz(f_mass,f_theta_rhodz,f_dtheta_rhodz)
660    TYPE(t_field), POINTER :: f_mass(:)
661    TYPE(t_field), POINTER :: f_theta_rhodz(:)
662    TYPE(t_field), POINTER :: f_dtheta_rhodz(:)
663   
664    REAL(rstd),POINTER :: mass(:,:)
665    REAL(rstd),POINTER :: theta_rhodz(:,:,:)
666    REAL(rstd),POINTER :: dtheta_rhodz(:,:)
667
668    INTEGER :: ind
669    INTEGER :: it,l,ij
670
671    CALL trace_start("divgrad")
672
673    DO ind=1,ndomain
674      IF (.NOT. assigned_domain(ind)) CYCLE
675      CALL swap_dimensions(ind)
676      CALL swap_geometry(ind)
677      mass=f_mass(ind)
678      theta_rhodz=f_theta_rhodz(ind)
679      dtheta_rhodz=f_dtheta_rhodz(ind) 
680      DO  l = ll_begin, ll_end
681!DIR$ SIMD
682        DO ij=ij_begin,ij_end
683            dtheta_rhodz(ij,l) = theta_rhodz(ij,l,1) / mass(ij,l)
684        ENDDO
685      ENDDO
686    ENDDO
687
688    DO it=1,niterdivgrad
689       
690      CALL send_message(f_dtheta_rhodz,req_dtheta)
691      CALL wait_message(req_dtheta)
692      DO ind=1,ndomain
693        IF (.NOT. assigned_domain(ind)) CYCLE
694        CALL swap_dimensions(ind)
695        CALL swap_geometry(ind)
696        dtheta_rhodz=f_dtheta_rhodz(ind) 
697        CALL compute_divgrad(dtheta_rhodz,dtheta_rhodz,ll_begin,ll_end)
698      ENDDO
699
700    ENDDO
701
702    DO ind=1,ndomain
703      IF (.NOT. assigned_domain(ind)) CYCLE
704      CALL swap_dimensions(ind)
705      CALL swap_geometry(ind)
706      dtheta_rhodz=f_dtheta_rhodz(ind) 
707      mass=f_mass(ind)
708   
709      DO  l = ll_begin, ll_end
710!DIR$ SIMD
711        DO ij=ij_begin,ij_end
712            dtheta_rhodz(ij,l) = dtheta_rhodz(ij,l) * mass(ij,l)
713        ENDDO
714      ENDDO
715    ENDDO
716
717    CALL trace_end("divgrad")
718  END SUBROUTINE divgrad_theta_rhodz 
719 
720  SUBROUTINE compute_gradiv(ue,gradivu_e,llb,lle)
721    INTEGER,INTENT(IN)     :: llb
722    INTEGER,INTENT(IN)     :: lle
723    REAL(rstd),INTENT(IN)  :: ue(iim*3*jjm,llm)
724    REAL(rstd),INTENT(OUT) :: gradivu_e(iim*3*jjm,llm)
725    REAL(rstd) :: divu_i(iim*jjm,llb:lle)
726   
727    INTEGER :: ij,l
728
729    DO l=llb,lle
730!DIR$ SIMD
731      DO ij=ij_begin,ij_end
732          divu_i(ij,l)=1./Ai(ij)*(ne(ij,right)*ue(ij+u_right,l)*le(ij+u_right)  +  &
733                             ne(ij,rup)*ue(ij+u_rup,l)*le(ij+u_rup)        +  & 
734                             ne(ij,lup)*ue(ij+u_lup,l)*le(ij+u_lup)        +  & 
735                             ne(ij,left)*ue(ij+u_left,l)*le(ij+u_left)     +  & 
736                             ne(ij,ldown)*ue(ij+u_ldown,l)*le(ij+u_ldown)  +  & 
737                             ne(ij,rdown)*ue(ij+u_rdown,l)*le(ij+u_rdown))
738      ENDDO
739    ENDDO
740   
741    DO l=llb,lle
742!DIR$ SIMD
743      DO ij=ij_begin,ij_end 
744          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) )       
745          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))   
746          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) )
747      ENDDO
748    ENDDO
749
750    DO l=llb,lle
751!DIR$ SIMD
752      DO ij=ij_begin,ij_end
753          gradivu_e(ij+u_right,l)=-gradivu_e(ij+u_right,l)*cgraddiv
754          gradivu_e(ij+u_lup,l)=-gradivu_e(ij+u_lup,l)*cgraddiv
755          gradivu_e(ij+u_ldown,l)=-gradivu_e(ij+u_ldown,l)*cgraddiv
756      ENDDO
757    ENDDO
758
759   
760  END SUBROUTINE compute_gradiv
761 
762  SUBROUTINE compute_divgrad(theta,divgrad_i,llb,lle)
763    INTEGER,INTENT(IN)     :: llb
764    INTEGER,INTENT(IN)     :: lle
765    REAL(rstd),INTENT(IN)  :: theta(iim*jjm,llm)
766    REAL(rstd),INTENT(OUT) :: divgrad_i(iim*jjm,llm)
767    REAL(rstd)  :: grad_e(3*iim*jjm,llb:lle)
768    INTEGER :: ij,l       
769    DO l=llb,lle
770!DIR$ SIMD
771      DO ij=ij_begin_ext,ij_end_ext
772          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) )
773          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 ))
774          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) )
775      ENDDO
776    ENDDO   
777   
778    DO l=llb,lle
779!DIR$ SIMD
780      DO ij=ij_begin,ij_end
781          divgrad_i(ij,l)=1./Ai(ij)*(ne(ij,right)*grad_e(ij+u_right,l)*le(ij+u_right)  +  &
782                             ne(ij,rup)*grad_e(ij+u_rup,l)*le(ij+u_rup)              +  & 
783                             ne(ij,lup)*grad_e(ij+u_lup,l)*le(ij+u_lup)              +  & 
784                             ne(ij,left)*grad_e(ij+u_left,l)*le(ij+u_left)           +  & 
785                             ne(ij,ldown)*grad_e(ij+u_ldown,l)*le(ij+u_ldown)        +  & 
786                             ne(ij,rdown)*grad_e(ij+u_rdown,l)*le(ij+u_rdown))
787      ENDDO
788    ENDDO
789   
790    DO l=llb,lle
791      DO ij=ij_begin,ij_end
792          divgrad_i(ij,l)=-divgrad_i(ij,l)*cdivgrad
793      ENDDO
794    ENDDO
795
796  END SUBROUTINE compute_divgrad
797   
798  SUBROUTINE compute_gradrot(ue,gradrot_e,llb,lle)
799    INTEGER,INTENT(IN)     :: llb
800    INTEGER,INTENT(IN)     :: lle
801    REAL(rstd),INTENT(IN)  :: ue(iim*3*jjm,llm)
802    REAL(rstd),INTENT(OUT) :: gradrot_e(iim*3*jjm,llm)
803    REAL(rstd) :: rot_v(2*iim*jjm,llb:lle)
804
805    INTEGER :: ij,l
806     
807    DO l=llb,lle
808!DIR$ SIMD
809      DO ij=ij_begin_ext,ij_end_ext
810       
811          rot_v(ij+z_up,l)= 1./Av(ij+z_up)*(  ne(ij,rup)*ue(ij+u_rup,l)*de(ij+u_rup)                     &
812                                + ne(ij+t_rup,left)*ue(ij+t_rup+u_left,l)*de(ij+t_rup+u_left)          &
813                                - ne(ij,lup)*ue(ij+u_lup,l)*de(ij+u_lup) ) 
814                             
815          rot_v(ij+z_down,l) = 1./Av(ij+z_down)*( ne(ij,ldown)*ue(ij+u_ldown,l)*de(ij+u_ldown)                 &
816                                     + ne(ij+t_ldown,right)*ue(ij+t_ldown+u_right,l)*de(ij+t_ldown+u_right)  &
817                                     - ne(ij,rdown)*ue(ij+u_rdown,l)*de(ij+u_rdown) )
818 
819      ENDDO
820    ENDDO                             
821 
822    DO l=llb,lle
823!DIR$ SIMD
824      DO ij=ij_begin,ij_end
825          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))         
826          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))
827          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))       
828      ENDDO
829    ENDDO
830
831    DO l=llb,lle
832!DIR$ SIMD
833      DO ij=ij_begin,ij_end   
834          gradrot_e(ij+u_right,l)=-gradrot_e(ij+u_right,l)*cgradrot       
835          gradrot_e(ij+u_lup,l)=-gradrot_e(ij+u_lup,l)*cgradrot
836          gradrot_e(ij+u_ldown,l)=-gradrot_e(ij+u_ldown,l)*cgradrot
837      ENDDO
838    ENDDO 
839   
840  END SUBROUTINE compute_gradrot
841
842!----------------------- Utility routines ------------------
843
844  SUBROUTINE global_max(dumax)
845    USE mpi_mod
846    USE mpipara
847    REAL(rstd) :: dumax, dumax1
848    IF (using_mpi) THEN
849       CALL reduce_max_omp(dumax,dumax1)
850       !$OMP MASTER       
851       CALL MPI_ALLREDUCE(dumax1,dumax,1,MPI_REAL8,MPI_MAX,comm_icosa,ierr)
852       !$OMP END MASTER     
853       CALL bcast_omp(dumax) 
854    ELSE
855       CALL allreduce_max_omp(dumax,dumax1)
856       dumax=dumax1
857    ENDIF
858  END SUBROUTINE global_max   
859
860  FUNCTION max_vel(f_du)
861    TYPE(t_field)       :: f_du(:)
862    REAL(rstd),POINTER  :: du(:)
863    INTEGER :: ind, i,j,ij
864    REAL(rstd) :: max_vel, dumax
865
866    dumax=0.
867    DO ind=1,ndomain
868       IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE
869       CALL swap_dimensions(ind)
870       CALL swap_geometry(ind)
871       du=f_du(ind)
872       
873       DO j=jj_begin,jj_end
874          DO i=ii_begin,ii_end
875             ij=(j-1)*iim+i   
876             if (le(ij+u_right)>1e-100) dumax=MAX(dumax,ABS(du(ij+u_right)))
877             if (le(ij+u_lup)>1e-100) dumax=MAX(dumax,ABS(du(ij+u_lup)))
878             if (le(ij+u_ldown)>1e-100) dumax=MAX(dumax,ABS(du(ij+u_ldown)))
879          ENDDO
880       ENDDO
881    ENDDO
882    CALL global_max(dumax)
883    max_vel=dumax
884  END FUNCTION max_vel
885
886  FUNCTION max_scalar(f_dtheta)
887    TYPE(t_field)       :: f_dtheta(:)
888    REAL(rstd),POINTER  :: dtheta(:)
889    INTEGER :: ind, i,j,ij
890    REAL(rstd) :: max_scalar, dthetamax
891    DO ind=1,ndomain
892       IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE
893       CALL swap_dimensions(ind)
894       dtheta=f_dtheta(ind)
895       DO j=jj_begin,jj_end
896          DO i=ii_begin,ii_end
897             ij=(j-1)*iim+i
898             dthetamax=MAX(dthetamax,ABS(dtheta(ij)))
899          ENDDO
900       ENDDO
901    ENDDO
902    CALL global_max(dthetamax)
903    max_scalar=dthetamax
904  END FUNCTION max_scalar
905
906  SUBROUTINE random_vel(f_u)
907    TYPE(t_field)      :: f_u(:)
908    REAL(rstd),POINTER :: u(:)
909    INTEGER :: M, ind, i,j,ij
910    REAL(rstd) :: r
911
912    !$OMP BARRIER
913    !$OMP MASTER
914    DO ind=1,ndomain
915       CALL swap_dimensions(ind)
916       CALL swap_geometry(ind)
917       u=f_u(ind)
918       
919       ! set random seed to get reproductibility when using a different number of process
920       CALL RANDOM_SEED(size=M)
921       CALL RANDOM_SEED(put=(/(i,i=1,M)/))
922       
923       DO j=jj_begin,jj_end
924          DO i=ii_begin,ii_end
925             ij=(j-1)*iim+i   
926             CALL RANDOM_NUMBER(r)
927             u(ij+u_right)=r-0.5
928             CALL RANDOM_NUMBER(r)
929             u(ij+u_lup)=r-0.5
930             CALL RANDOM_NUMBER(r)
931             u(ij+u_ldown)=r-0.5
932          ENDDO
933       ENDDO
934    ENDDO
935    !$OMP END MASTER 
936    !$OMP BARRIER
937   
938  END SUBROUTINE random_vel
939 
940  SUBROUTINE random_scalar(f_theta)
941    TYPE(t_field)      :: f_theta(:)
942    REAL(rstd),POINTER :: theta(:)
943    INTEGER :: M, ind, i,j,ij
944    REAL(rstd) :: r
945    !$OMP BARRIER
946    !$OMP MASTER
947    DO ind=1,ndomain
948       CALL swap_dimensions(ind)
949       CALL swap_geometry(ind)
950       theta=f_theta(ind)
951       ! set random seed to get reproductibility when using a different number of process
952       CALL RANDOM_SEED(size=M)
953       CALL RANDOM_SEED(put=(/(i,i=1,M)/))
954       
955       DO j=jj_begin,jj_end
956          DO i=ii_begin,ii_end
957             ij=(j-1)*iim+i   
958             CALL RANDOM_NUMBER(r)
959             theta(ij)=r-0.5
960          ENDDO
961       ENDDO
962    ENDDO
963    !$OMP END MASTER
964    !$OMP BARRIER
965  END SUBROUTINE random_scalar
966 
967  SUBROUTINE rescale_field(scal, f_du, f_u)
968    TYPE(t_field)      :: f_du(:), f_u(:)
969    REAL(rstd),POINTER :: du(:), u(:)
970    REAL(rstd) :: scal
971    INTEGER :: ind
972    DO ind=1,ndomain
973       IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE
974       CALL swap_dimensions(ind)
975       u=f_u(ind)
976       du=f_du(ind)
977       u=scal*du
978    ENDDO
979  END SUBROUTINE rescale_field
980
981END MODULE dissip_gcm_mod
Note: See TracBrowser for help on using the repository browser.