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

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

Set constant sign for wind way :
ne(ij,right)==ne_right=1
ne(ij,rup)==ne_rup=-1
ne(ij,lup)==ne_lup=1
ne(ij,left)==ne_left=-1
ne(ij,ldown)==ne_ldown=1
ne(ij,rdown)==ne_rdown=-1

Modified transfert function to be compliant for this convention.

YM
File size: 25.7 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_vect)
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_vect)
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_vect)
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_vect)
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  USE trace
433  IMPLICIT NONE
434    TYPE(t_field),POINTER :: f_ue(:)
435    TYPE(t_field),POINTER :: f_due(:)
436    TYPE(t_field),POINTER :: f_ps(:), f_phis(:)
437    TYPE(t_field),POINTER :: f_theta_rhodz(:)
438    TYPE(t_field),POINTER :: f_dtheta_rhodz(:)
439
440    REAL(rstd),POINTER         :: due(:,:)
441    REAL(rstd),POINTER         :: phi(:,:), ue(:,:)
442    REAL(rstd),POINTER         :: due_diss1(:,:)
443    REAL(rstd),POINTER         :: due_diss2(:,:)
444    REAL(rstd),POINTER         :: dtheta_rhodz(:,:)
445    REAL(rstd),POINTER         :: dtheta_rhodz_diss(:,:)
446
447    INTEGER :: ind, shear
448    INTEGER :: l,i,j,n
449   
450    CALL trace_start("dissip")
451    CALL gradiv(f_ue,f_due_diss1)
452    CALL gradrot(f_ue,f_due_diss2)
453    CALL theta_rhodz2theta(f_ps,f_theta_rhodz,f_theta)
454    CALL divgrad(f_theta,f_dtheta_diss)
455    CALL theta2theta_rhodz(f_ps,f_dtheta_diss,f_dtheta_rhodz_diss)
456   
457    IF(rayleigh_friction_type>0) THEN
458       CALL pression(f_ps, f_p)
459       CALL exner(f_ps, f_p, f_pks, f_pk)
460       CALL geopotential(f_phis,f_pks,f_pk,f_theta,f_phi)
461    END IF
462
463    DO ind=1,ndomain
464      CALL swap_dimensions(ind)
465      CALL swap_geometry(ind)
466      due=f_due(ind) 
467      due_diss1=f_due_diss1(ind)
468      due_diss2=f_due_diss2(ind)
469      dtheta_rhodz=f_dtheta_rhodz(ind)
470      dtheta_rhodz_diss=f_dtheta_rhodz_diss(ind)
471           
472      DO l=1,llm
473        DO j=jj_begin,jj_end
474          DO i=ii_begin,ii_end
475            n=(j-1)*iim+i
476
477            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) )
478            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) )
479            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) )
480
481            dtheta_rhodz(n,l) = -0.5*dtheta_rhodz_diss(n,l)/tau_divgrad(l)
482          ENDDO
483        ENDDO
484      ENDDO
485
486      IF(rayleigh_friction_type>0) THEN
487         phi=f_phi(ind)
488         ue=f_ue(ind)
489         DO l=1,llm
490            DO j=jj_begin,jj_end
491               DO i=ii_begin,ii_end
492                  n=(j-1)*iim+i
493                  CALL relax(t_right, u_right)
494                  CALL relax(t_lup,   u_lup)
495                  CALL relax(t_ldown, u_ldown)
496               ENDDO
497            ENDDO
498         END DO
499      END IF
500   END DO
501
502   CALL trace_end("dissip")
503
504    CONTAINS
505      SUBROUTINE relax(shift_t, shift_u)
506        USE dcmip_initial_conditions_test_1_2_3
507        REAL(rstd) :: z, ulon,ulat, lon,lat, & ! input to test2_schaer_mountain
508             p,hyam,hybm,w,t,phis,ps,rho,q, &   ! unused input/output to test2_schaer_mountain
509             fz, u3d(3), uref
510        REAL(rstd), PARAMETER :: zh=2e4,ztop=3e4  ! DCMIP values
511        LOGICAL :: hybrid_eta
512        INTEGER :: shift_u, shift_t, zcoords, nn
513        z = (phi(n,l)+phi(n+shift_t,l))/(2.*g)
514        IF(z>zh) THEN  ! relax only in the sponge layer z>zh
515!           PRINT *, 'z>zh : z,zh,l',z,zh,l
516!           STOP
517           nn = n+shift_u
518           CALL xyz2lonlat(xyz_e(nn,:),lon,lat)
519           zcoords = 1 ; hybrid_eta = .FALSE. ! use z instead of p or hyam/hybm
520           CALL test2_schaer_mountain(lon,lat,p,z,zcoords,hybrid_eta, &
521                hyam,hybm,shear,ulon,ulat,w,t,phis,ps,rho,q)
522!           u3d = ulon*elon_e(nn,:) + ulat*elat_e(nn,:)
523           u3d = ulon*elon_e(nn,:) ! ulat=0
524           uref = sum(u3d*ep_e(nn,:))
525
526           fz = sin((pi/2)*(z-zh)/(ztop-zh))
527           fz = fz*fz/rayleigh_tau
528!           fz = 1/rayleigh_tau
529           due(nn,l) = due(nn,l) - fz*(ue(nn,l)-uref)
530!           due(nn,l) = due(nn,l) - fz*ue(nn,l)
531         END IF
532       END SUBROUTINE relax
533     
534  END SUBROUTINE dissip
535
536  SUBROUTINE gradiv(f_ue,f_due)
537  USE icosa
538  USE trace
539  IMPLICIT NONE
540    TYPE(t_field),POINTER :: f_ue(:)
541    TYPE(t_field),POINTER :: f_due(:)
542    REAL(rstd),POINTER    :: ue(:,:)
543    REAL(rstd),POINTER    :: due(:,:)
544    INTEGER :: ind
545    INTEGER :: it
546       
547    CALL trace_start("gradiv")
548
549    DO ind=1,ndomain
550      CALL swap_dimensions(ind)
551      CALL swap_geometry(ind)
552      ue=f_ue(ind)
553      due=f_due(ind) 
554      due=ue
555    ENDDO
556
557    DO it=1,nitergdiv
558       
559      CALL transfert_request(f_due,req_e1_vect)
560       
561      DO ind=1,ndomain
562        CALL swap_dimensions(ind)
563        CALL swap_geometry(ind)
564        due=f_due(ind) 
565        CALL compute_gradiv(due,due,llm)
566      ENDDO
567    ENDDO
568
569   CALL trace_end("gradiv")
570
571  END SUBROUTINE gradiv
572 
573
574  SUBROUTINE gradrot(f_ue,f_due)
575  USE icosa
576  USE trace
577  IMPLICIT NONE
578    TYPE(t_field),POINTER :: f_ue(:)
579    TYPE(t_field),POINTER :: f_due(:)
580    REAL(rstd),POINTER    :: ue(:,:)
581    REAL(rstd),POINTER    :: due(:,:)
582    INTEGER :: ind
583    INTEGER :: it
584       
585    CALL trace_start("gradrot")
586
587    DO ind=1,ndomain
588      CALL swap_dimensions(ind)
589      CALL swap_geometry(ind)
590      ue=f_ue(ind)
591      due=f_due(ind) 
592      due=ue
593    ENDDO
594
595    DO it=1,nitergrot
596       
597      CALL transfert_request(f_due,req_e1_vect)
598       
599      DO ind=1,ndomain
600        CALL swap_dimensions(ind)
601        CALL swap_geometry(ind)
602        due=f_due(ind) 
603        CALL compute_gradrot(due,due,llm)
604      ENDDO
605
606    ENDDO
607
608    CALL trace_end("gradrot")
609
610  END SUBROUTINE gradrot
611 
612  SUBROUTINE divgrad(f_theta,f_dtheta)
613  USE icosa
614  USE trace
615  IMPLICIT NONE
616    TYPE(t_field),POINTER :: f_theta(:)
617    TYPE(t_field),POINTER :: f_dtheta(:)
618    REAL(rstd),POINTER    :: theta(:,:)
619    REAL(rstd),POINTER    :: dtheta(:,:)
620    INTEGER :: ind
621    INTEGER :: it
622
623    CALL trace_start("divgrad")
624       
625    DO ind=1,ndomain
626      CALL swap_dimensions(ind)
627      CALL swap_geometry(ind)
628      theta=f_theta(ind)
629      dtheta=f_dtheta(ind) 
630      dtheta=theta
631    ENDDO
632
633    DO it=1,niterdivgrad
634       
635      CALL transfert_request(f_dtheta,req_i1)
636       
637      DO ind=1,ndomain
638        CALL swap_dimensions(ind)
639        CALL swap_geometry(ind)
640        dtheta=f_dtheta(ind) 
641        CALL compute_divgrad(dtheta,dtheta,llm)
642      ENDDO
643
644    ENDDO
645
646    CALL trace_end("divgrad")
647
648  END SUBROUTINE divgrad
649   
650 
651
652!  SUBROUTINE compute_dissip(ue,due,ps,theta_rhodz,dtheta_rhodz)
653!  USE icosa
654!  USE theta2theta_rhodz_mod
655!  IMPLICIT NONE
656!    REAL(rstd)         :: ue(3*iim*jjm,llm)
657!    REAL(rstd)         :: due(3*iim*jjm,llm)
658!    REAL(rstd)         :: ps(iim*jjm)
659!    REAL(rstd)         :: theta_rhodz(iim*jjm,llm)
660!    REAL(rstd)         :: dtheta_rhodz(iim*jjm,llm)
661!
662!    REAL(rstd),SAVE,ALLOCATABLE :: theta(:,:)
663!    REAL(rstd),SAVE,ALLOCATABLE :: du_dissip(:,:)
664!    REAL(rstd),SAVE,ALLOCATABLE :: dtheta_dissip(:,:)
665!    REAL(rstd),SAVE,ALLOCATABLE :: dtheta_rhodz_dissip(:,:)
666!
667!    INTEGER :: ind
668!    INTEGER :: l,i,j,n
669!
670!!$OMP BARRIER
671!!$OMP MASTER
672!      ALLOCATE(theta(iim*jjm,llm))
673!      ALLOCATE(du_dissip(3*iim*jjm,llm))
674!      ALLOCATE(dtheta_dissip(iim*jjm,llm))
675!      ALLOCATE(dtheta_rhodz_dissip(iim*jjm,llm))
676!!$OMP END MASTER
677!!$OMP BARRIER
678!   
679!      CALL gradiv(ue,du_dissip,llm)
680!      DO l=1,llm
681!!$OMP DO
682!        DO j=jj_begin,jj_end
683!          DO i=ii_begin,ii_end
684!            n=(j-1)*iim+i
685!            due(n+u_right,l)=due(n+u_right,l)+du_dissip(n+u_right,l)/tau_graddiv(l)*0.5
686!            due(n+u_lup,l)=due(n+u_lup,l)+du_dissip(n+u_lup,l)/tau_graddiv(l)*0.5
687!            due(n+u_ldown,l)=due(n+u_ldown,l)+du_dissip(n+u_ldown,l)/tau_graddiv(l)*0.5
688!          ENDDO
689!        ENDDO
690!      ENDDO
691!     
692!      CALL gradrot(ue,du_dissip,llm)
693!
694!      DO l=1,llm
695!!$OMP DO
696!        DO j=jj_begin,jj_end
697!          DO i=ii_begin,ii_end
698!            n=(j-1)*iim+i
699!            due(n+u_right,l)=due(n+u_right,l)+du_dissip(n+u_right,l)/tau_gradrot(l)*0.5
700!            due(n+u_lup,l)=due(n+u_lup,l)+du_dissip(n+u_lup,l)/tau_gradrot(l)*0.5
701!            due(n+u_ldown,l)=due(n+u_ldown,l)+du_dissip(n+u_ldown,l)/tau_gradrot(l)*0.5
702!          ENDDO
703!        ENDDO
704!      ENDDO
705!     
706!      CALL compute_theta_rhodz2theta(ps,theta_rhodz,theta,1)
707!      CALL divgrad(theta,dtheta_dissip,llm)
708!      CALL compute_theta2theta_rhodz(ps,dtheta_dissip,dtheta_rhodz_dissip,0)
709!
710!      DO l=1,llm
711!!$OMP DO
712!        DO j=jj_begin,jj_end
713!          DO i=ii_begin,ii_end
714!            n=(j-1)*iim+i
715!            dtheta_rhodz(n,l)=dtheta_rhodz(n,l)+dtheta_rhodz_dissip(n,l)/tau_divgrad(l)*0.5
716!          ENDDO
717!        ENDDO
718!      ENDDO
719!
720!!$OMP BARRIER
721!!$OMP MASTER
722!      DEALLOCATE(theta)
723!      DEALLOCATE(du_dissip)
724!      DEALLOCATE(dtheta_dissip)
725!      DEALLOCATE(dtheta_rhodz_dissip)
726!!$OMP END MASTER
727!!$OMP BARRIER     
728!
729!  END SUBROUTINE compute_dissip
730 
731 
732  SUBROUTINE compute_gradiv(ue,gradivu_e,ll)
733  USE icosa
734  IMPLICIT NONE
735    INTEGER,INTENT(IN)     :: ll
736    REAL(rstd),INTENT(IN)  :: ue(iim*3*jjm,ll)
737    REAL(rstd),INTENT(OUT) :: gradivu_e(iim*3*jjm,ll)
738    REAL(rstd),SAVE,ALLOCATABLE :: divu_i(:,:)
739   
740    INTEGER :: i,j,n,l
741
742!$OMP BARRIER
743!$OMP MASTER
744      ALLOCATE(divu_i(iim*jjm,ll))
745!$OMP END MASTER
746!$OMP BARRIER
747
748    DO l=1,ll
749!$OMP DO
750      DO j=jj_begin,jj_end
751        DO i=ii_begin,ii_end
752          n=(j-1)*iim+i
753          divu_i(n,l)=1./Ai(n)*(ne(n,right)*ue(n+u_right,l)*le(n+u_right)  +  &
754                             ne(n,rup)*ue(n+u_rup,l)*le(n+u_rup)        +  & 
755                             ne(n,lup)*ue(n+u_lup,l)*le(n+u_lup)        +  & 
756                             ne(n,left)*ue(n+u_left,l)*le(n+u_left)     +  & 
757                             ne(n,ldown)*ue(n+u_ldown,l)*le(n+u_ldown)  +  & 
758                             ne(n,rdown)*ue(n+u_rdown,l)*le(n+u_rdown))
759        ENDDO
760      ENDDO
761    ENDDO
762   
763    DO l=1,ll
764!$OMP DO
765      DO j=jj_begin,jj_end
766        DO i=ii_begin,ii_end
767 
768          n=(j-1)*iim+i
769 
770          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) )       
771
772          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))       
773   
774          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) )
775
776        ENDDO
777      ENDDO
778    ENDDO
779
780    DO l=1,ll
781!$OMP DO
782      DO j=jj_begin,jj_end
783        DO i=ii_begin,ii_end
784          n=(j-1)*iim+i
785          gradivu_e(n+u_right,l)=-gradivu_e(n+u_right,l)*cgraddiv
786          gradivu_e(n+u_lup,l)=-gradivu_e(n+u_lup,l)*cgraddiv
787          gradivu_e(n+u_ldown,l)=-gradivu_e(n+u_ldown,l)*cgraddiv
788        ENDDO
789      ENDDO
790    ENDDO
791
792!$OMP BARRIER
793!$OMP MASTER
794      DEALLOCATE(divu_i)
795!$OMP END MASTER
796!$OMP BARRIER
797
798   
799  END SUBROUTINE compute_gradiv
800 
801  SUBROUTINE compute_divgrad(theta,divgrad_i,ll)
802  USE icosa
803  IMPLICIT NONE
804    INTEGER,INTENT(IN)     :: ll
805    REAL(rstd),INTENT(IN)  :: theta(iim*jjm,ll)
806    REAL(rstd),INTENT(OUT) :: divgrad_i(iim*jjm,ll)
807    REAL(rstd),SAVE,ALLOCATABLE :: grad_e(:,:)
808
809    INTEGER :: i,j,n,l
810
811!$OMP BARRIER
812!$OMP MASTER
813      ALLOCATE(grad_e(3*iim*jjm,ll))
814!$OMP END MASTER
815!$OMP BARRIER
816       
817    DO l=1,ll
818!$OMP DO
819      DO j=jj_begin-1,jj_end+1
820        DO i=ii_begin-1,ii_end+1
821   
822          n=(j-1)*iim+i
823 
824          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) )       
825
826          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 ))       
827   
828          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) )
829
830        ENDDO
831      ENDDO
832    ENDDO
833   
834   
835    DO l=1,ll
836!$OMP DO
837      DO j=jj_begin,jj_end
838        DO i=ii_begin,ii_end
839          n=(j-1)*iim+i
840          divgrad_i(n,l)=1./Ai(n)*(ne(n,right)*grad_e(n+u_right,l)*le(n+u_right)  +  &
841                             ne(n,rup)*grad_e(n+u_rup,l)*le(n+u_rup)              +  & 
842                             ne(n,lup)*grad_e(n+u_lup,l)*le(n+u_lup)              +  & 
843                             ne(n,left)*grad_e(n+u_left,l)*le(n+u_left)           +  & 
844                             ne(n,ldown)*grad_e(n+u_ldown,l)*le(n+u_ldown)        +  & 
845                             ne(n,rdown)*grad_e(n+u_rdown,l)*le(n+u_rdown))
846        ENDDO
847      ENDDO
848    ENDDO
849   
850    DO l=1,ll
851!$OMP DO
852      DO j=jj_begin,jj_end
853        DO i=ii_begin,ii_end
854          n=(j-1)*iim+i
855          divgrad_i(n,l)=-divgrad_i(n,l)*cdivgrad
856        ENDDO
857      ENDDO
858    ENDDO
859
860!$OMP BARRIER
861!$OMP MASTER
862      DEALLOCATE(grad_e)
863!$OMP END MASTER
864!$OMP BARRIER
865
866  END SUBROUTINE compute_divgrad
867
868   
869  SUBROUTINE compute_gradrot(ue,gradrot_e,ll)
870  USE icosa
871  IMPLICIT NONE
872    INTEGER,INTENT(IN)     :: ll
873    REAL(rstd),INTENT(IN)  :: ue(iim*3*jjm,ll)
874    REAL(rstd),INTENT(OUT) :: gradrot_e(iim*3*jjm,ll)
875    REAL(rstd),SAVE,ALLOCATABLE  :: rot_v(:,:)
876
877    INTEGER :: i,j,n,l
878
879!$OMP BARRIER
880!$OMP MASTER
881      ALLOCATE(rot_v(2*iim*jjm,ll))
882!$OMP END MASTER
883!$OMP BARRIER
884     
885    DO l=1,ll
886!$OMP DO
887      DO j=jj_begin-1,jj_end+1
888        DO i=ii_begin-1,ii_end+1
889          n=(j-1)*iim+i
890       
891          rot_v(n+z_up,l)= 1./Av(n+z_up)*(  ne(n,rup)*ue(n+u_rup,l)*de(n+u_rup)                     &
892                                + ne(n+t_rup,left)*ue(n+t_rup+u_left,l)*de(n+t_rup+u_left)          &
893                                - ne(n,lup)*ue(n+u_lup,l)*de(n+u_lup) ) 
894                             
895          rot_v(n+z_down,l) = 1./Av(n+z_down)*( ne(n,ldown)*ue(n+u_ldown,l)*de(n+u_ldown)                 &
896                                     + ne(n+t_ldown,right)*ue(n+t_ldown+u_right,l)*de(n+t_ldown+u_right)  &
897                                     - ne(n,rdown)*ue(n+u_rdown,l)*de(n+u_rdown) )
898 
899        ENDDO
900      ENDDO
901    ENDDO                             
902 
903    DO l=1,ll
904!$OMP DO
905      DO j=jj_begin,jj_end
906        DO i=ii_begin,ii_end
907          n=(j-1)*iim+i
908 
909          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)) 
910         
911          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)) 
912       
913          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))
914       
915        ENDDO
916      ENDDO
917    ENDDO
918
919    DO l=1,ll
920!$OMP DO
921      DO j=jj_begin,jj_end
922        DO i=ii_begin,ii_end
923          n=(j-1)*iim+i
924   
925          gradrot_e(n+u_right,l)=-gradrot_e(n+u_right,l)*cgradrot       
926          gradrot_e(n+u_lup,l)=-gradrot_e(n+u_lup,l)*cgradrot
927          gradrot_e(n+u_ldown,l)=-gradrot_e(n+u_ldown,l)*cgradrot
928       
929        ENDDO
930      ENDDO
931    ENDDO 
932   
933!$OMP BARRIER
934!$OMP MASTER
935      DEALLOCATE(rot_v)
936!$OMP END MASTER
937!$OMP BARRIER
938   
939  END SUBROUTINE compute_gradrot
940
941
942END MODULE dissip_gcm_mod
943       
Note: See TracBrowser for help on using the repository browser.