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

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