source: codes/icosagcm/trunk/src/geopotential_mod.f90 @ 121

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

Options for calculation of Exner pressure
File size: 2.0 KB
Line 
1MODULE geopotential_mod
2
3CONTAINS
4
5  SUBROUTINE geopotential(f_phis,f_pks,f_pk,f_theta,f_phi)
6  USE icosa
7  IMPLICIT NONE
8    TYPE(t_field), POINTER :: f_phis(:)
9    TYPE(t_field), POINTER :: f_pks(:)
10    TYPE(t_field), POINTER :: f_pk(:)
11    TYPE(t_field), POINTER :: f_theta(:)
12    TYPE(t_field), POINTER :: f_phi(:)
13 
14    REAL(rstd), POINTER :: phis(:)
15    REAL(rstd), POINTER :: pks(:)
16    REAL(rstd), POINTER :: pk(:,:)
17    REAL(rstd), POINTER :: theta(:,:)
18    REAL(rstd), POINTER :: phi(:,:)
19    INTEGER :: ind
20
21    DO ind=1,ndomain
22      CALL swap_dimensions(ind)
23      CALL swap_geometry(ind)
24      phis=f_phis(ind)
25      pks=f_pks(ind)
26      pk=f_pk(ind)
27      theta=f_theta(ind)
28      phi=f_phi(ind)
29      CALL compute_geopotential(phis,pks,pk,theta,phi,0)
30    ENDDO
31 
32  END SUBROUTINE geopotential
33 
34  SUBROUTINE compute_geopotential(phis,pks,pk,theta,phi,offset)
35  USE icosa
36  IMPLICIT NONE
37    REAL(rstd),INTENT(IN) :: phis(iim*jjm)
38    REAL(rstd),INTENT(IN) :: pks(iim*jjm)
39    REAL(rstd),INTENT(IN) :: pk(iim*jjm,llm)
40    REAL(rstd),INTENT(IN) :: theta(iim*jjm,llm)
41    REAL(rstd),INTENT(OUT) :: phi(iim*jjm,llm)
42    INTEGER,INTENT(IN) :: offset
43    INTEGER :: i,j,ij,l
44
45!!! Compute geopotential
46  ! LMD-Z : dh=pi.dtheta-vdp  => (1/rho) dp = pi.dtheta - d(h=pi.theta) = -dtheta.gradpi
47  ! g.dz =  -(1/rho)dp =  theta.gradpi
48  ! direct (variational)
49  ! dz = -(1/rho.g)dp, rho = RT/p = (R/C_p).theta.pi/p
50   
51  ! for first 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       phi( ij,1 ) = phis( ij ) + theta(ij,1) * ( pks(ij) - pk(ij,1) )
56     ENDDO
57   ENDDO
58         
59  ! for other layers
60   DO l = 2, llm
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         phi(ij,l) = phi(ij,l-1) + 0.5 * ( theta(ij,l)  + theta(ij,l-1) )  & 
65                                       * (  pk(ij,l-1) -  pk(ij,l)    )
66       ENDDO
67     ENDDO
68   ENDDO
69   
70  END SUBROUTINE compute_geopotential
71
72END MODULE geopotential_mod
Note: See TracBrowser for help on using the repository browser.