source: codes/icosagcm/trunk/src/dissip_gcm.f90 @ 123

Last change on this file since 123 was 109, checked in by dubos, 12 years ago

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