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

Last change on this file since 123 was 121, checked in by dubos, 11 years ago

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