source: codes/icosagcm/trunk/src/exner.f90 @ 155

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

Implementation of mixte parallelism MPI/OpenMP into src directory

YM
File size: 3.3 KB
Line 
1MODULE exner_mod
2
3  INTEGER :: caldyn_exner
4  INTEGER, PARAMETER :: lmdz=3, direct=4
5
6CONTAINS
7 
8  SUBROUTINE exner(f_ps,f_p,f_pks,f_pk)
9  USE icosa
10  IMPLICIT NONE
11    TYPE(t_field), POINTER :: f_ps(:)
12    TYPE(t_field), POINTER :: f_p(:)
13    TYPE(t_field), POINTER :: f_pks(:)
14    TYPE(t_field), POINTER :: f_pk(:)
15 
16    REAL(rstd), POINTER :: ps(:)
17    REAL(rstd), POINTER :: p(:,:)
18    REAL(rstd), POINTER :: pks(:)
19    REAL(rstd), POINTER :: pk(:,:)
20    INTEGER :: ind
21
22    DO ind=1,ndomain
23      CALL swap_dimensions(ind)
24      CALL swap_geometry(ind)
25      ps=f_ps(ind)
26      p=f_p(ind)
27      pks=f_pks(ind)
28      pk=f_pk(ind)
29      CALL compute_exner(ps, p, pks, pk, 0)
30    ENDDO
31 
32  END SUBROUTINE exner
33 
34  SUBROUTINE compute_exner(ps,p,pks,pk,offset)
35  USE icosa
36  USE disvert_mod
37  USE pression_mod
38  IMPLICIT NONE
39    REAL(rstd),INTENT(IN) :: ps(iim*jjm)
40    REAL(rstd),INTENT(IN) :: p(iim*jjm,llm+1)
41    REAL(rstd),INTENT(OUT) :: pks(iim*jjm)
42    REAL(rstd),INTENT(OUT) :: pk(iim*jjm,llm)
43    INTEGER,INTENT(IN) :: offset
44    INTEGER :: i,j,ij,l
45    REAL(rstd) :: alpha(iim*jjm,llm),beta(iim*jjm,llm)
46    REAL(rstd) :: delta 
47   
48    IF(caldyn_exner == lmdz) THEN          ! LMD-Z style calculation of Exner pressure
49       !! Compute Alpha and Beta
50       
51       ! for llm layer
52       DO j=jj_begin-offset,jj_end+offset
53          DO i=ii_begin-offset,ii_end+offset
54             ij=(j-1)*iim+i
55             alpha(ij,llm) = 0.
56             beta (ij,llm) = 1./ (1+ 2*kappa)
57          ENDDO
58       ENDDO
59       
60       ! for other layer   
61       DO l = llm-1 , 2 , -1
62          DO j=jj_begin-offset,jj_end+offset
63             DO i=ii_begin-offset,ii_end+offset
64                ij=(j-1)*iim+i
65                delta = p(ij,l)* (1+2*kappa) + p(ij,l+1)* ( beta(ij,l+1)- (1+2*kappa) )
66                alpha(ij,l)  = - p(ij,l+1) / delta * alpha(ij,l+1)
67                beta (ij,l)  =   p(ij,l  ) / delta   
68             ENDDO
69          ENDDO
70       ENDDO
71       
72       !! Compute pk
73       
74       ! for first layer
75       DO j=jj_begin-offset,jj_end+offset
76          DO i=ii_begin-offset,ii_end+offset
77             ij=(j-1)*iim+i
78             pks(ij) = cpp * ( ps(ij)/preff ) ** kappa
79             pk(ij,1) = ( p(ij,1)*pks(ij) - 0.5*alpha(ij,2)*p(ij,2) )     /    &
80                  (  p(ij,1)* (1.+kappa) + 0.5*( beta(ij,2)-(1.+2*kappa) )* p(ij,2) )
81          ENDDO
82       ENDDO
83       
84       ! for other layers
85       
86       DO l = 2, llm
87          DO j=jj_begin-offset,jj_end+offset
88             DO i=ii_begin-offset,ii_end+offset
89                ij=(j-1)*iim+i
90                pk(ij,l) = alpha(ij,l) + beta(ij,l) * pk(ij,l-1)
91             ENDDO
92          ENDDO
93       ENDDO
94       
95    ELSE ! Simple calculation of Exner pressure based on centered average
96       ! surface : pks
97       DO j=jj_begin-offset,jj_end+offset
98          DO i=ii_begin-offset,ii_end+offset
99             ij=(j-1)*iim+i
100             pks(ij) = cpp * ( ps(ij)/preff ) ** kappa
101          ENDDO
102       ENDDO
103
104       ! 3D : pk
105       DO l = 1, llm
106          DO j=jj_begin-offset,jj_end+offset
107             DO i=ii_begin-offset,ii_end+offset
108                ij=(j-1)*iim+i
109                pk(ij,l) = cpp * ((.5/preff)*(p(ij,l)+p(ij,l+1))) ** kappa
110             ENDDO
111          ENDDO
112       ENDDO
113    END IF
114   
115  END SUBROUTINE compute_exner
116 
117END MODULE exner_mod
Note: See TracBrowser for help on using the repository browser.