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

Last change on this file since 899 was 899, checked in by adurocher, 5 years ago

trunk : Fixed GCC warnings

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