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

Last change on this file since 78 was 50, checked in by dubos, 12 years ago

Temporary fix to problem with Exner function
Now also outputs velocity lat-lon, geopotential, pressure (at interfaces), Exner pressure, potential temp
Tested : test case DCMIP-3-1 with nbp=20
File size: 3.3 KB
Line 
1MODULE exner_mod
2
3CONTAINS
4
5  SUBROUTINE exner(f_ps,f_p,f_pks,f_pk)
6  USE icosa
7  IMPLICIT NONE
8    TYPE(t_field), POINTER :: f_ps(:)
9    TYPE(t_field), POINTER :: f_p(:)
10    TYPE(t_field), POINTER :: f_pks(:)
11    TYPE(t_field), POINTER :: f_pk(:)
12 
13    REAL(rstd), POINTER :: ps(:)
14    REAL(rstd), POINTER :: p(:,:)
15    REAL(rstd), POINTER :: pks(:)
16    REAL(rstd), POINTER :: pk(:,:)
17    INTEGER :: ind
18
19    DO ind=1,ndomain
20      CALL swap_dimensions(ind)
21      CALL swap_geometry(ind)
22      ps=f_ps(ind)
23      p=f_p(ind)
24      pks=f_pks(ind)
25      pk=f_pk(ind)
26      CALL compute_exner(ps, p, pks, pk, 0)
27    ENDDO
28 
29  END SUBROUTINE exner
30 
31  SUBROUTINE compute_exner(ps,p,pks,pk,offset)
32  USE icosa
33  USE disvert_mod
34  USE pression_mod
35  IMPLICIT NONE
36    REAL(rstd),INTENT(IN) :: ps(iim*jjm)
37    REAL(rstd),INTENT(IN) :: p(iim*jjm,llm+1)
38    REAL(rstd),INTENT(OUT) :: pks(iim*jjm)
39    REAL(rstd),INTENT(OUT) :: pk(iim*jjm,llm)
40    INTEGER,INTENT(IN) :: offset
41    INTEGER :: i,j,ij,l
42    REAL(rstd) :: alpha(iim*jjm,llm),beta(iim*jjm,llm)
43    REAL(rstd) :: delta 
44   
45    IF(.FALSE.) THEN          ! LMD-Z style calculation of Exner pressure
46       !! Compute Alpha and Beta
47       
48       ! for llm layer
49       !$OMP DO
50       DO j=jj_begin-offset,jj_end+offset
51          DO i=ii_begin-offset,ii_end+offset
52             ij=(j-1)*iim+i
53             alpha(ij,llm) = 0.
54             beta (ij,llm) = 1./ (1+ 2*kappa)
55          ENDDO
56       ENDDO
57       
58       ! for other layer   
59       DO l = llm-1 , 2 , -1
60          !$OMP DO
61          DO j=jj_begin-offset,jj_end+offset
62             DO i=ii_begin-offset,ii_end+offset
63                ij=(j-1)*iim+i
64                delta = p(ij,l)* (1+2*kappa) + p(ij,l+1)* ( beta(ij,l+1)- (1+2*kappa) )
65                alpha(ij,l)  = - p(ij,l+1) / delta * alpha(ij,l+1)
66                beta (ij,l)  =   p(ij,l  ) / delta   
67             ENDDO
68          ENDDO
69       ENDDO
70       
71       !! Compute pk
72       
73       ! for first layer
74       !$OMP DO
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          !$OMP DO
88          DO j=jj_begin-offset,jj_end+offset
89             DO i=ii_begin-offset,ii_end+offset
90                ij=(j-1)*iim+i
91                pk(ij,l) = alpha(ij,l) + beta(ij,l) * pk(ij,l-1)
92             ENDDO
93          ENDDO
94       ENDDO
95       
96    ELSE ! Simple calculation of Exner pressure based on centered average
97       ! surface : pks
98       !$OMP DO
99       DO j=jj_begin-offset,jj_end+offset
100          DO i=ii_begin-offset,ii_end+offset
101             ij=(j-1)*iim+i
102             pks(ij) = cpp * ( ps(ij)/preff ) ** kappa
103          ENDDO
104       ENDDO
105
106       ! 3D : pk
107       DO l = 1, llm
108          !$OMP DO
109          DO j=jj_begin-offset,jj_end+offset
110             DO i=ii_begin-offset,ii_end+offset
111                ij=(j-1)*iim+i
112                pk(ij,l) = cpp * ((.5/preff)*(p(ij,l)+p(ij,l+1))) ** kappa
113             ENDDO
114          ENDDO
115       ENDDO
116    END IF
117   
118  END SUBROUTINE compute_exner
119 
120END MODULE exner_mod
Note: See TracBrowser for help on using the repository browser.