source: codes/icosagcm/trunk/src/dissip/dissip_gcm.f90 @ 781

Last change on this file since 781 was 781, checked in by sfromang, 6 years ago

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