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

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

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