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

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

Various optimisations

  • dissipation is not called every timestep (similar way than LMDZ)
  • transfert size of halos has been reduced : need to synchronise redondant data between tiles at itau_sync timestep

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