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

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

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