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

Last change on this file since 141 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
RevLine 
[15]1MODULE dissip_gcm_mod
[19]2  USE icosa
[15]3
[26]4  PRIVATE
5
6  TYPE(t_field),POINTER,SAVE :: f_due_diss1(:)
7  TYPE(t_field),POINTER,SAVE :: f_due_diss2(:)
8
[109]9  TYPE(t_field),POINTER,SAVE :: f_theta(:), f_phi(:), f_pk(:), f_pks(:), f_p(:)
[26]10  TYPE(t_field),POINTER,SAVE :: f_dtheta_diss(:)
11  TYPE(t_field),POINTER,SAVE :: f_dtheta_rhodz_diss(:)
12 
[15]13 
[26]14  INTEGER,SAVE :: nitergdiv=1
15  INTEGER,SAVE :: nitergrot=1
16  INTEGER,SAVE :: niterdivgrad=1
[15]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
[109]25
26  INTEGER, SAVE :: rayleigh_friction_type, rayleigh_shear
27  REAL, save    :: rayleigh_tau
[15]28 
[26]29  INTEGER,SAVE :: idissip
30  REAL,SAVE    :: dtdissip
[15]31 
[26]32  PUBLIC init_dissip, dissip
[15]33CONTAINS
34
35  SUBROUTINE allocate_dissip
[19]36  USE icosa
[15]37  IMPLICIT NONE 
[26]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)
[109]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)
[26]48   
[15]49    ALLOCATE(tau_graddiv(llm))
50    ALLOCATE(tau_gradrot(llm))   
51    ALLOCATE(tau_divgrad(llm))
52  END SUBROUTINE allocate_dissip
53 
[98]54  SUBROUTINE init_dissip
[19]55  USE icosa
[15]56  USE disvert_mod
[26]57  USE mpi_mod
58  USE mpipara
[15]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(:)
[59]70  REAL(rstd)            :: dumax,dumax1
71  REAL(rstd)            :: dthetamax,dthetamax1
[15]72  REAL(rstd)            :: r
73  REAL(rstd)            :: tau
74  REAL(rstd)            :: zz, zvert, fact
75  INTEGER               :: l
[109]76  CHARACTER(len=255)    :: rayleigh_friction_key
[15]77 
78           
[26]79  INTEGER :: i,j,n,ind,it,iter
[15]80
[109]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
[131]86      IF (is_mpi_root) PRINT *, 'No Rayleigh friction'
[109]87   CASE('dcmip2_schaer_noshear')
88      rayleigh_friction_type=1
89      rayleigh_shear=0
[131]90      IF (is_mpi_root) PRINT *, 'Rayleigh friction : Schaer-like mountain without shear DCMIP2.1'
[109]91   CASE('dcmip2_schaer_shear')
92      rayleigh_shear=1
93      rayleigh_friction_type=2
[131]94      IF (is_mpi_root) PRINT *, 'Rayleigh friction : Schaer-like mountain with shear DCMIP2.2'
[109]95   CASE DEFAULT
[131]96      IF (is_mpi_root) PRINT *, 'Bad selector : rayleigh_friction_type =', TRIM(rayleigh_friction_key), ' in dissip_gcm.f90/init_dissip'
[109]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
[131]105         IF (is_mpi_root) PRINT *, 'Forbidden : negative value for rayleigh_friction_tau =',rayleigh_tau
[109]106         STOP
107      END IF
108   END IF
109
[15]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)
[32]118    tau_graddiv(:)=tau/scale_factor
[26]119
120    CALL getin("nitergdiv",nitergdiv)
[15]121   
122    tau_gradrot(:)=5000
[26]123    CALL getin("tau_gradrot",tau_gradrot)
[32]124    tau_gradrot(:)=tau/scale_factor
[26]125
126    CALL getin("nitergrot",nitergrot)
[15]127   
128
129    tau_divgrad(:)=5000
130    CALL getin("tau_divgrad",tau)
[32]131    tau_divgrad(:)=tau/scale_factor
[26]132
133    CALL getin("niterdivgrad",niterdivgrad)
[15]134 
[26]135!   CALL create_request(field_u,req_dissip)
[15]136
[26]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
[15]142
[26]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   
[15]152
[26]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   
[15]157
[26]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   
[15]167
[26]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
[15]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
[29]203    DO it=1,20
[26]204     
[15]205      dumax=0
[26]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
[15]216      ENDDO
[26]217
218      CALL transfert_request(f_du,req_e1)
[15]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)
[26]225         
[15]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
[59]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                       
[15]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
[131]248      IF (is_mpi_root) PRINT *,"gradiv : it :",it ,": dumax",dumax
[15]249
250    ENDDO 
[131]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)
[129]256
[15]257    cgraddiv=dumax**(-1./nitergdiv)
[131]258    IF (is_mpi_root) PRINT *,"cgraddiv : ",cgraddiv
[15]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
[29]280    DO it=1,20
[15]281 
282      dumax=0
[26]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
[15]293      ENDDO
294
[26]295      CALL transfert_request(f_du,req_e1)
[15]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
[26]312
[59]313      IF (using_mpi) THEN
314        CALL MPI_ALLREDUCE(dumax,dumax1,1,MPI_REAL8,MPI_MAX,comm_icosa,ierr)
315        dumax=dumax1
316      ENDIF
[15]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     
[131]326      IF (is_mpi_root) PRINT *,"gradrot : it :",it ,": dumax",dumax
[15]327
328    ENDDO 
[131]329    IF (is_mpi_root) PRINT *,"gradrot : dumax",dumax
[15]330 
[26]331    cgradrot=dumax**(-1./nitergrot) 
[131]332    IF (is_mpi_root) PRINT *,"cgradrot : ",cgradrot
[15]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
[29]351    DO it=1,20
[15]352
353      dthetamax=0
[26]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
[15]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
[59]382      IF (using_mpi) THEN
383        CALL MPI_ALLREDUCE(dthetamax,dthetamax1,1,MPI_REAL8,MPI_MAX,comm_icosa,ierr)
384        dthetamax=dthetamax1
385      ENDIF
[131]386      IF (is_mpi_root) PRINT *,"divgrad : it :",it ,": dthetamax",dthetamax
[15]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)
[131]398    IF (is_mpi_root) PRINT *,"divgrad : divgrad",dthetamax
[15]399 
[26]400    cdivgrad=dthetamax**(-1./niterdivgrad) 
[131]401    IF (is_mpi_root) PRINT *,"cdivgrad : ",cdivgrad
[15]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 
[109]426  SUBROUTINE dissip(f_ue,f_due,f_ps,f_phis,f_theta_rhodz,f_dtheta_rhodz)
[19]427  USE icosa
[15]428  USE theta2theta_rhodz_mod
[109]429  USE pression_mod
430  USE exner_mod
431  USE geopotential_mod
[15]432  IMPLICIT NONE
433    TYPE(t_field),POINTER :: f_ue(:)
434    TYPE(t_field),POINTER :: f_due(:)
[109]435    TYPE(t_field),POINTER :: f_ps(:), f_phis(:)
[15]436    TYPE(t_field),POINTER :: f_theta_rhodz(:)
437    TYPE(t_field),POINTER :: f_dtheta_rhodz(:)
[26]438
[15]439    REAL(rstd),POINTER         :: due(:,:)
[109]440    REAL(rstd),POINTER         :: phi(:,:), ue(:,:)
[26]441    REAL(rstd),POINTER         :: due_diss1(:,:)
442    REAL(rstd),POINTER         :: due_diss2(:,:)
[15]443    REAL(rstd),POINTER         :: dtheta_rhodz(:,:)
[26]444    REAL(rstd),POINTER         :: dtheta_rhodz_diss(:,:)
[15]445
[109]446    INTEGER :: ind, shear
[15]447    INTEGER :: l,i,j,n
448   
[26]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   
[109]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
[15]461    DO ind=1,ndomain
462      CALL swap_dimensions(ind)
463      CALL swap_geometry(ind)
464      due=f_due(ind) 
[26]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
[15]474
[129]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) )
[26]478
[129]479            dtheta_rhodz(n,l) = -0.5*dtheta_rhodz_diss(n,l)/tau_divgrad(l)
[26]480          ENDDO
481        ENDDO
482      ENDDO
[109]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
[15]529     
530  END SUBROUTINE dissip
531
[26]532  SUBROUTINE gradiv(f_ue,f_due)
[19]533  USE icosa
[15]534  IMPLICIT NONE
[26]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
[15]549
[26]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
[15]561
[26]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(:,:)
[15]572    INTEGER :: ind
[26]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
[15]582
[26]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)
[15]592      ENDDO
593
[26]594    ENDDO
[15]595
[26]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
[15]615
[26]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
[15]626
[26]627    ENDDO
628
629  END SUBROUTINE divgrad
630   
[15]631 
[26]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
[15]711 
[26]712 
713  SUBROUTINE compute_gradiv(ue,gradivu_e,ll)
[19]714  USE icosa
[15]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
[26]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
[15]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   
[26]780  END SUBROUTINE compute_gradiv
[15]781 
[26]782  SUBROUTINE compute_divgrad(theta,divgrad_i,ll)
[19]783  USE icosa
[15]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
[26]836          divgrad_i(n,l)=-divgrad_i(n,l)*cdivgrad
[15]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
[26]847  END SUBROUTINE compute_divgrad
[15]848
849   
[26]850  SUBROUTINE compute_gradrot(ue,gradrot_e,ll)
[19]851  USE icosa
[15]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   
[26]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
[15]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   
[26]920  END SUBROUTINE compute_gradrot
[15]921
922
923END MODULE dissip_gcm_mod
924       
Note: See TracBrowser for help on using the repository browser.