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

Last change on this file since 140 was 131, checked in by ymipsl, 11 years ago

Some operations must be only done by the mpi master task.

YM
File size: 25.4 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      IF (is_mpi_root) PRINT *, 'No Rayleigh friction'
87   CASE('dcmip2_schaer_noshear')
88      rayleigh_friction_type=1
89      rayleigh_shear=0
90      IF (is_mpi_root) 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      IF (is_mpi_root) PRINT *, 'Rayleigh friction : Schaer-like mountain with shear DCMIP2.2'
95   CASE DEFAULT
96      IF (is_mpi_root) 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         IF (is_mpi_root) 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      IF (is_mpi_root) PRINT *,"gradiv : it :",it ,": dumax",dumax
249
250    ENDDO 
251    IF (is_mpi_root) PRINT *,"gradiv : dumax",dumax
252    IF (is_mpi_root) PRINT *, 'mean T-cell edge size (km)', 1.45*radius/iim_glo/1000., &
253                              'effective T-cell half-edge size (km)', dumax**(-.5/nitergdiv)/1000
254    IF (is_mpi_root)  PRINT *, 'Max. time step assuming c=340 m/s and Courant number=2.8 :', &
255                               2.8/340.*dumax**(-.5/nitergdiv)
256
257    cgraddiv=dumax**(-1./nitergdiv)
258    IF (is_mpi_root) PRINT *,"cgraddiv : ",cgraddiv
259
260    DO ind=1,ndomain
261      CALL swap_dimensions(ind)
262      CALL swap_geometry(ind)
263      u=f_u(ind)
264
265     DO j=jj_begin,jj_end
266        DO i=ii_begin,ii_end
267          n=(j-1)*iim+i
268          CALL RANDOM_NUMBER(r)
269          u(n+u_right)=r-0.5
270          CALL RANDOM_NUMBER(r)
271          u(n+u_lup)=r-0.5
272          CALL RANDOM_NUMBER(r)
273          u(n+u_ldown)=r-0.5
274       ENDDO
275      ENDDO
276       
277    ENDDO
278
279
280    DO it=1,20
281 
282      dumax=0
283      DO iter=1,nitergrot
284        CALL transfert_request(f_u,req_e1)
285        DO ind=1,ndomain
286          CALL swap_dimensions(ind)
287          CALL swap_geometry(ind)
288          u=f_u(ind)
289          du=f_du(ind)
290          CALL compute_gradrot(u,du,1)
291          u=du
292        ENDDO
293      ENDDO
294
295      CALL transfert_request(f_du,req_e1)
296     
297      DO ind=1,ndomain
298        CALL swap_dimensions(ind)
299        CALL swap_geometry(ind)
300        u=f_u(ind)
301        du=f_du(ind)
302       
303        DO j=jj_begin,jj_end
304          DO i=ii_begin,ii_end
305            n=(j-1)*iim+i
306            if (le(n+u_right)>1e-100) dumax=MAX(dumax,ABS(du(n+u_right)))
307            if (le(n+u_lup)>1e-100) dumax=MAX(dumax,ABS(du(n+u_lup)))
308            if (le(n+u_ldown)>1e-100) dumax=MAX(dumax,ABS(du(n+u_ldown)))
309          ENDDO
310        ENDDO
311      ENDDO
312
313      IF (using_mpi) THEN
314        CALL MPI_ALLREDUCE(dumax,dumax1,1,MPI_REAL8,MPI_MAX,comm_icosa,ierr)
315        dumax=dumax1
316      ENDIF
317     
318      DO ind=1,ndomain
319        CALL swap_dimensions(ind)
320        CALL swap_geometry(ind)
321        u=f_u(ind)
322        du=f_du(ind)
323        u=du/dumax
324      ENDDO
325     
326      IF (is_mpi_root) PRINT *,"gradrot : it :",it ,": dumax",dumax
327
328    ENDDO 
329    IF (is_mpi_root) PRINT *,"gradrot : dumax",dumax
330 
331    cgradrot=dumax**(-1./nitergrot) 
332    IF (is_mpi_root) PRINT *,"cgradrot : ",cgradrot
333   
334
335
336    DO ind=1,ndomain
337      CALL swap_dimensions(ind)
338      CALL swap_geometry(ind)
339      theta=f_theta(ind)
340 
341       DO j=jj_begin,jj_end
342          DO i=ii_begin,ii_end
343            n=(j-1)*iim+i
344            CALL RANDOM_NUMBER(r)
345            theta(n)=r-0.5
346         ENDDO
347        ENDDO
348       
349    ENDDO
350
351    DO it=1,20
352
353      dthetamax=0
354      DO iter=1,niterdivgrad
355        CALL transfert_request(f_theta,req_i1)
356        DO ind=1,ndomain
357          CALL swap_dimensions(ind)
358          CALL swap_geometry(ind)
359          theta=f_theta(ind)
360          dtheta=f_dtheta(ind)
361          CALL compute_divgrad(theta,dtheta,1)
362          theta=dtheta
363        ENDDO
364      ENDDO
365!      CALL writefield("divgrad",f_dtheta)
366
367      CALL transfert_request(f_dtheta,req_i1)
368     
369      DO ind=1,ndomain
370        CALL swap_dimensions(ind)
371        CALL swap_geometry(ind)
372        theta=f_theta(ind)
373        dtheta=f_dtheta(ind)
374       
375        DO j=jj_begin,jj_end
376          DO i=ii_begin,ii_end
377            n=(j-1)*iim+i
378            dthetamax=MAX(dthetamax,ABS(dtheta(n)))
379          ENDDO
380        ENDDO
381      ENDDO
382      IF (using_mpi) THEN
383        CALL MPI_ALLREDUCE(dthetamax,dthetamax1,1,MPI_REAL8,MPI_MAX,comm_icosa,ierr)
384        dthetamax=dthetamax1
385      ENDIF
386      IF (is_mpi_root) PRINT *,"divgrad : it :",it ,": dthetamax",dthetamax
387
388      DO ind=1,ndomain
389        CALL swap_dimensions(ind)
390        CALL swap_geometry(ind)
391        theta=f_theta(ind)
392        dtheta=f_dtheta(ind)
393        theta=dtheta/dthetamax
394      ENDDO
395    ENDDO 
396
397!    CALL writefield("divgrad",f_dtheta)
398    IF (is_mpi_root) PRINT *,"divgrad : divgrad",dthetamax
399 
400    cdivgrad=dthetamax**(-1./niterdivgrad) 
401    IF (is_mpi_root) PRINT *,"cdivgrad : ",cdivgrad
402   
403
404
405
406
407
408     
409    fact=2
410    DO l=1,llm
411       zz= 1. - preff/presnivs(l)
412       zvert=fact-(fact-1)/(1+zz*zz)
413       tau_graddiv(l) = tau_graddiv(l)/zvert
414       tau_gradrot(l) = tau_gradrot(l)/zvert
415       tau_divgrad(l) = tau_divgrad(l)/zvert
416    ENDDO
417       
418!    idissip=INT(MIN(tetagdiv,tetagrot)/(2.*dt))
419!    idissip=MAX(1,idissip)
420!    dtdissip=idissip*dt
421!    PRINT *,"idissip",idissip," dtdissip ",dtdissip
422
423  END SUBROUTINE init_dissip 
424 
425 
426  SUBROUTINE dissip(f_ue,f_due,f_ps,f_phis,f_theta_rhodz,f_dtheta_rhodz)
427  USE icosa
428  USE theta2theta_rhodz_mod
429  USE pression_mod
430  USE exner_mod
431  USE geopotential_mod
432  IMPLICIT NONE
433    TYPE(t_field),POINTER :: f_ue(:)
434    TYPE(t_field),POINTER :: f_due(:)
435    TYPE(t_field),POINTER :: f_ps(:), f_phis(:)
436    TYPE(t_field),POINTER :: f_theta_rhodz(:)
437    TYPE(t_field),POINTER :: f_dtheta_rhodz(:)
438
439    REAL(rstd),POINTER         :: due(:,:)
440    REAL(rstd),POINTER         :: phi(:,:), ue(:,:)
441    REAL(rstd),POINTER         :: due_diss1(:,:)
442    REAL(rstd),POINTER         :: due_diss2(:,:)
443    REAL(rstd),POINTER         :: dtheta_rhodz(:,:)
444    REAL(rstd),POINTER         :: dtheta_rhodz_diss(:,:)
445
446    INTEGER :: ind, shear
447    INTEGER :: l,i,j,n
448   
449    CALL gradiv(f_ue,f_due_diss1)
450    CALL gradrot(f_ue,f_due_diss2)
451    CALL theta_rhodz2theta(f_ps,f_theta_rhodz,f_theta)
452    CALL divgrad(f_theta,f_dtheta_diss)
453    CALL theta2theta_rhodz(f_ps,f_dtheta_diss,f_dtheta_rhodz_diss)
454   
455    IF(rayleigh_friction_type>0) THEN
456       CALL pression(f_ps, f_p)
457       CALL exner(f_ps, f_p, f_pks, f_pk)
458       CALL geopotential(f_phis,f_pks,f_pk,f_theta,f_phi)
459    END IF
460
461    DO ind=1,ndomain
462      CALL swap_dimensions(ind)
463      CALL swap_geometry(ind)
464      due=f_due(ind) 
465      due_diss1=f_due_diss1(ind)
466      due_diss2=f_due_diss2(ind)
467      dtheta_rhodz=f_dtheta_rhodz(ind)
468      dtheta_rhodz_diss=f_dtheta_rhodz_diss(ind)
469           
470      DO l=1,llm
471        DO j=jj_begin,jj_end
472          DO i=ii_begin,ii_end
473            n=(j-1)*iim+i
474
475            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) )
476            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) )
477            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) )
478
479            dtheta_rhodz(n,l) = -0.5*dtheta_rhodz_diss(n,l)/tau_divgrad(l)
480          ENDDO
481        ENDDO
482      ENDDO
483
484      IF(rayleigh_friction_type>0) THEN
485         phi=f_phi(ind)
486         ue=f_ue(ind)
487         DO l=1,llm
488            DO j=jj_begin,jj_end
489               DO i=ii_begin,ii_end
490                  n=(j-1)*iim+i
491                  CALL relax(t_right, u_right)
492                  CALL relax(t_lup,   u_lup)
493                  CALL relax(t_ldown, u_ldown)
494               ENDDO
495            ENDDO
496         END DO
497      END IF
498   END DO
499
500    CONTAINS
501      SUBROUTINE relax(shift_t, shift_u)
502        USE dcmip_initial_conditions_test_1_2_3
503        REAL(rstd) :: z, ulon,ulat, lon,lat, & ! input to test2_schaer_mountain
504             p,hyam,hybm,w,t,phis,ps,rho,q, &   ! unused input/output to test2_schaer_mountain
505             fz, u3d(3), uref
506        REAL(rstd), PARAMETER :: zh=2e4,ztop=3e4  ! DCMIP values
507        LOGICAL :: hybrid_eta
508        INTEGER :: shift_u, shift_t, zcoords, nn
509        z = (phi(n,l)+phi(n+shift_t,l))/(2.*g)
510        IF(z>zh) THEN  ! relax only in the sponge layer z>zh
511!           PRINT *, 'z>zh : z,zh,l',z,zh,l
512!           STOP
513           nn = n+shift_u
514           CALL xyz2lonlat(xyz_e(nn,:),lon,lat)
515           zcoords = 1 ; hybrid_eta = .FALSE. ! use z instead of p or hyam/hybm
516           CALL test2_schaer_mountain(lon,lat,p,z,zcoords,hybrid_eta, &
517                hyam,hybm,shear,ulon,ulat,w,t,phis,ps,rho,q)
518!           u3d = ulon*elon_e(nn,:) + ulat*elat_e(nn,:)
519           u3d = ulon*elon_e(nn,:) ! ulat=0
520           uref = sum(u3d*ep_e(nn,:))
521
522           fz = sin((pi/2)*(z-zh)/(ztop-zh))
523           fz = fz*fz/rayleigh_tau
524!           fz = 1/rayleigh_tau
525           due(nn,l) = due(nn,l) - fz*(ue(nn,l)-uref)
526!           due(nn,l) = due(nn,l) - fz*ue(nn,l)
527         END IF
528       END SUBROUTINE relax
529     
530  END SUBROUTINE dissip
531
532  SUBROUTINE gradiv(f_ue,f_due)
533  USE icosa
534  IMPLICIT NONE
535    TYPE(t_field),POINTER :: f_ue(:)
536    TYPE(t_field),POINTER :: f_due(:)
537    REAL(rstd),POINTER    :: ue(:,:)
538    REAL(rstd),POINTER    :: due(:,:)
539    INTEGER :: ind
540    INTEGER :: it
541       
542    DO ind=1,ndomain
543      CALL swap_dimensions(ind)
544      CALL swap_geometry(ind)
545      ue=f_ue(ind)
546      due=f_due(ind) 
547      due=ue
548    ENDDO
549
550    DO it=1,nitergdiv
551       
552      CALL transfert_request(f_due,req_e1)
553       
554      DO ind=1,ndomain
555        CALL swap_dimensions(ind)
556        CALL swap_geometry(ind)
557        due=f_due(ind) 
558        CALL compute_gradiv(due,due,llm)
559      ENDDO
560    ENDDO
561
562  END SUBROUTINE gradiv
563 
564
565  SUBROUTINE gradrot(f_ue,f_due)
566  USE icosa
567  IMPLICIT NONE
568    TYPE(t_field),POINTER :: f_ue(:)
569    TYPE(t_field),POINTER :: f_due(:)
570    REAL(rstd),POINTER    :: ue(:,:)
571    REAL(rstd),POINTER    :: due(:,:)
572    INTEGER :: ind
573    INTEGER :: it
574       
575    DO ind=1,ndomain
576      CALL swap_dimensions(ind)
577      CALL swap_geometry(ind)
578      ue=f_ue(ind)
579      due=f_due(ind) 
580      due=ue
581    ENDDO
582
583    DO it=1,nitergrot
584       
585      CALL transfert_request(f_due,req_e1)
586       
587      DO ind=1,ndomain
588        CALL swap_dimensions(ind)
589        CALL swap_geometry(ind)
590        due=f_due(ind) 
591        CALL compute_gradrot(due,due,llm)
592      ENDDO
593
594    ENDDO
595
596  END SUBROUTINE gradrot
597 
598  SUBROUTINE divgrad(f_theta,f_dtheta)
599  USE icosa
600  IMPLICIT NONE
601    TYPE(t_field),POINTER :: f_theta(:)
602    TYPE(t_field),POINTER :: f_dtheta(:)
603    REAL(rstd),POINTER    :: theta(:,:)
604    REAL(rstd),POINTER    :: dtheta(:,:)
605    INTEGER :: ind
606    INTEGER :: it
607       
608    DO ind=1,ndomain
609      CALL swap_dimensions(ind)
610      CALL swap_geometry(ind)
611      theta=f_theta(ind)
612      dtheta=f_dtheta(ind) 
613      dtheta=theta
614    ENDDO
615
616    DO it=1,niterdivgrad
617       
618      CALL transfert_request(f_dtheta,req_i1)
619       
620      DO ind=1,ndomain
621        CALL swap_dimensions(ind)
622        CALL swap_geometry(ind)
623        dtheta=f_dtheta(ind) 
624        CALL compute_divgrad(dtheta,dtheta,llm)
625      ENDDO
626
627    ENDDO
628
629  END SUBROUTINE divgrad
630   
631 
632
633!  SUBROUTINE compute_dissip(ue,due,ps,theta_rhodz,dtheta_rhodz)
634!  USE icosa
635!  USE theta2theta_rhodz_mod
636!  IMPLICIT NONE
637!    REAL(rstd)         :: ue(3*iim*jjm,llm)
638!    REAL(rstd)         :: due(3*iim*jjm,llm)
639!    REAL(rstd)         :: ps(iim*jjm)
640!    REAL(rstd)         :: theta_rhodz(iim*jjm,llm)
641!    REAL(rstd)         :: dtheta_rhodz(iim*jjm,llm)
642!
643!    REAL(rstd),SAVE,ALLOCATABLE :: theta(:,:)
644!    REAL(rstd),SAVE,ALLOCATABLE :: du_dissip(:,:)
645!    REAL(rstd),SAVE,ALLOCATABLE :: dtheta_dissip(:,:)
646!    REAL(rstd),SAVE,ALLOCATABLE :: dtheta_rhodz_dissip(:,:)
647!
648!    INTEGER :: ind
649!    INTEGER :: l,i,j,n
650!
651!!$OMP BARRIER
652!!$OMP MASTER
653!      ALLOCATE(theta(iim*jjm,llm))
654!      ALLOCATE(du_dissip(3*iim*jjm,llm))
655!      ALLOCATE(dtheta_dissip(iim*jjm,llm))
656!      ALLOCATE(dtheta_rhodz_dissip(iim*jjm,llm))
657!!$OMP END MASTER
658!!$OMP BARRIER
659!   
660!      CALL gradiv(ue,du_dissip,llm)
661!      DO l=1,llm
662!!$OMP DO
663!        DO j=jj_begin,jj_end
664!          DO i=ii_begin,ii_end
665!            n=(j-1)*iim+i
666!            due(n+u_right,l)=due(n+u_right,l)+du_dissip(n+u_right,l)/tau_graddiv(l)*0.5
667!            due(n+u_lup,l)=due(n+u_lup,l)+du_dissip(n+u_lup,l)/tau_graddiv(l)*0.5
668!            due(n+u_ldown,l)=due(n+u_ldown,l)+du_dissip(n+u_ldown,l)/tau_graddiv(l)*0.5
669!          ENDDO
670!        ENDDO
671!      ENDDO
672!     
673!      CALL gradrot(ue,du_dissip,llm)
674!
675!      DO l=1,llm
676!!$OMP DO
677!        DO j=jj_begin,jj_end
678!          DO i=ii_begin,ii_end
679!            n=(j-1)*iim+i
680!            due(n+u_right,l)=due(n+u_right,l)+du_dissip(n+u_right,l)/tau_gradrot(l)*0.5
681!            due(n+u_lup,l)=due(n+u_lup,l)+du_dissip(n+u_lup,l)/tau_gradrot(l)*0.5
682!            due(n+u_ldown,l)=due(n+u_ldown,l)+du_dissip(n+u_ldown,l)/tau_gradrot(l)*0.5
683!          ENDDO
684!        ENDDO
685!      ENDDO
686!     
687!      CALL compute_theta_rhodz2theta(ps,theta_rhodz,theta,1)
688!      CALL divgrad(theta,dtheta_dissip,llm)
689!      CALL compute_theta2theta_rhodz(ps,dtheta_dissip,dtheta_rhodz_dissip,0)
690!
691!      DO l=1,llm
692!!$OMP DO
693!        DO j=jj_begin,jj_end
694!          DO i=ii_begin,ii_end
695!            n=(j-1)*iim+i
696!            dtheta_rhodz(n,l)=dtheta_rhodz(n,l)+dtheta_rhodz_dissip(n,l)/tau_divgrad(l)*0.5
697!          ENDDO
698!        ENDDO
699!      ENDDO
700!
701!!$OMP BARRIER
702!!$OMP MASTER
703!      DEALLOCATE(theta)
704!      DEALLOCATE(du_dissip)
705!      DEALLOCATE(dtheta_dissip)
706!      DEALLOCATE(dtheta_rhodz_dissip)
707!!$OMP END MASTER
708!!$OMP BARRIER     
709!
710!  END SUBROUTINE compute_dissip
711 
712 
713  SUBROUTINE compute_gradiv(ue,gradivu_e,ll)
714  USE icosa
715  IMPLICIT NONE
716    INTEGER,INTENT(IN)     :: ll
717    REAL(rstd),INTENT(IN)  :: ue(iim*3*jjm,ll)
718    REAL(rstd),INTENT(OUT) :: gradivu_e(iim*3*jjm,ll)
719    REAL(rstd),SAVE,ALLOCATABLE :: divu_i(:,:)
720   
721    INTEGER :: i,j,n,l
722
723!$OMP BARRIER
724!$OMP MASTER
725      ALLOCATE(divu_i(iim*jjm,ll))
726!$OMP END MASTER
727!$OMP BARRIER
728
729    DO l=1,ll
730!$OMP DO
731      DO j=jj_begin,jj_end
732        DO i=ii_begin,ii_end
733          n=(j-1)*iim+i
734          divu_i(n,l)=1./Ai(n)*(ne(n,right)*ue(n+u_right,l)*le(n+u_right)  +  &
735                             ne(n,rup)*ue(n+u_rup,l)*le(n+u_rup)        +  & 
736                             ne(n,lup)*ue(n+u_lup,l)*le(n+u_lup)        +  & 
737                             ne(n,left)*ue(n+u_left,l)*le(n+u_left)     +  & 
738                             ne(n,ldown)*ue(n+u_ldown,l)*le(n+u_ldown)  +  & 
739                             ne(n,rdown)*ue(n+u_rdown,l)*le(n+u_rdown))
740        ENDDO
741      ENDDO
742    ENDDO
743   
744    DO l=1,ll
745!$OMP DO
746      DO j=jj_begin,jj_end
747        DO i=ii_begin,ii_end
748 
749          n=(j-1)*iim+i
750 
751          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) )       
752
753          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))       
754   
755          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) )
756
757        ENDDO
758      ENDDO
759    ENDDO
760
761    DO l=1,ll
762!$OMP DO
763      DO j=jj_begin,jj_end
764        DO i=ii_begin,ii_end
765          n=(j-1)*iim+i
766          gradivu_e(n+u_right,l)=-gradivu_e(n+u_right,l)*cgraddiv
767          gradivu_e(n+u_lup,l)=-gradivu_e(n+u_lup,l)*cgraddiv
768          gradivu_e(n+u_ldown,l)=-gradivu_e(n+u_ldown,l)*cgraddiv
769        ENDDO
770      ENDDO
771    ENDDO
772
773!$OMP BARRIER
774!$OMP MASTER
775      DEALLOCATE(divu_i)
776!$OMP END MASTER
777!$OMP BARRIER
778
779   
780  END SUBROUTINE compute_gradiv
781 
782  SUBROUTINE compute_divgrad(theta,divgrad_i,ll)
783  USE icosa
784  IMPLICIT NONE
785    INTEGER,INTENT(IN)     :: ll
786    REAL(rstd),INTENT(IN)  :: theta(iim*jjm,ll)
787    REAL(rstd),INTENT(OUT) :: divgrad_i(iim*jjm,ll)
788    REAL(rstd),SAVE,ALLOCATABLE :: grad_e(:,:)
789
790    INTEGER :: i,j,n,l
791
792!$OMP BARRIER
793!$OMP MASTER
794      ALLOCATE(grad_e(3*iim*jjm,ll))
795!$OMP END MASTER
796!$OMP BARRIER
797       
798    DO l=1,ll
799!$OMP DO
800      DO j=jj_begin-1,jj_end+1
801        DO i=ii_begin-1,ii_end+1
802   
803          n=(j-1)*iim+i
804 
805          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) )       
806
807          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 ))       
808   
809          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) )
810
811        ENDDO
812      ENDDO
813    ENDDO
814   
815   
816    DO l=1,ll
817!$OMP DO
818      DO j=jj_begin,jj_end
819        DO i=ii_begin,ii_end
820          n=(j-1)*iim+i
821          divgrad_i(n,l)=1./Ai(n)*(ne(n,right)*grad_e(n+u_right,l)*le(n+u_right)  +  &
822                             ne(n,rup)*grad_e(n+u_rup,l)*le(n+u_rup)              +  & 
823                             ne(n,lup)*grad_e(n+u_lup,l)*le(n+u_lup)              +  & 
824                             ne(n,left)*grad_e(n+u_left,l)*le(n+u_left)           +  & 
825                             ne(n,ldown)*grad_e(n+u_ldown,l)*le(n+u_ldown)        +  & 
826                             ne(n,rdown)*grad_e(n+u_rdown,l)*le(n+u_rdown))
827        ENDDO
828      ENDDO
829    ENDDO
830   
831    DO l=1,ll
832!$OMP DO
833      DO j=jj_begin,jj_end
834        DO i=ii_begin,ii_end
835          n=(j-1)*iim+i
836          divgrad_i(n,l)=-divgrad_i(n,l)*cdivgrad
837        ENDDO
838      ENDDO
839    ENDDO
840
841!$OMP BARRIER
842!$OMP MASTER
843      DEALLOCATE(grad_e)
844!$OMP END MASTER
845!$OMP BARRIER
846
847  END SUBROUTINE compute_divgrad
848
849   
850  SUBROUTINE compute_gradrot(ue,gradrot_e,ll)
851  USE icosa
852  IMPLICIT NONE
853    INTEGER,INTENT(IN)     :: ll
854    REAL(rstd),INTENT(IN)  :: ue(iim*3*jjm,ll)
855    REAL(rstd),INTENT(OUT) :: gradrot_e(iim*3*jjm,ll)
856    REAL(rstd),SAVE,ALLOCATABLE  :: rot_v(:,:)
857
858    INTEGER :: i,j,n,l
859
860!$OMP BARRIER
861!$OMP MASTER
862      ALLOCATE(rot_v(2*iim*jjm,ll))
863!$OMP END MASTER
864!$OMP BARRIER
865     
866    DO l=1,ll
867!$OMP DO
868      DO j=jj_begin-1,jj_end+1
869        DO i=ii_begin-1,ii_end+1
870          n=(j-1)*iim+i
871       
872          rot_v(n+z_up,l)= 1./Av(n+z_up)*(  ne(n,rup)*ue(n+u_rup,l)*de(n+u_rup)                     &
873                                + ne(n+t_rup,left)*ue(n+t_rup+u_left,l)*de(n+t_rup+u_left)          &
874                                - ne(n,lup)*ue(n+u_lup,l)*de(n+u_lup) ) 
875                             
876          rot_v(n+z_down,l) = 1./Av(n+z_down)*( ne(n,ldown)*ue(n+u_ldown,l)*de(n+u_ldown)                 &
877                                     + ne(n+t_ldown,right)*ue(n+t_ldown+u_right,l)*de(n+t_ldown+u_right)  &
878                                     - ne(n,rdown)*ue(n+u_rdown,l)*de(n+u_rdown) )
879 
880        ENDDO
881      ENDDO
882    ENDDO                             
883 
884    DO l=1,ll
885!$OMP DO
886      DO j=jj_begin,jj_end
887        DO i=ii_begin,ii_end
888          n=(j-1)*iim+i
889 
890          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)) 
891         
892          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)) 
893       
894          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))
895       
896        ENDDO
897      ENDDO
898    ENDDO
899
900    DO l=1,ll
901!$OMP DO
902      DO j=jj_begin,jj_end
903        DO i=ii_begin,ii_end
904          n=(j-1)*iim+i
905   
906          gradrot_e(n+u_right,l)=-gradrot_e(n+u_right,l)*cgradrot       
907          gradrot_e(n+u_lup,l)=-gradrot_e(n+u_lup,l)*cgradrot
908          gradrot_e(n+u_ldown,l)=-gradrot_e(n+u_ldown,l)*cgradrot
909       
910        ENDDO
911      ENDDO
912    ENDDO 
913   
914!$OMP BARRIER
915!$OMP MASTER
916      DEALLOCATE(rot_v)
917!$OMP END MASTER
918!$OMP BARRIER
919   
920  END SUBROUTINE compute_gradrot
921
922
923END MODULE dissip_gcm_mod
924       
Note: See TracBrowser for help on using the repository browser.