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

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

Bug fix in dissipation.
=> better conservation of the mass

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