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

Last change on this file since 169 was 167, checked in by dubos, 11 years ago

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