source: codes/icosagcm/branches/SATURN_DYNAMICO/ICOSAGCM/src/geopotential_mod.f90 @ 316

Last change on this file since 316 was 221, checked in by ymipsl, 10 years ago

Creating temporary dynamico/lmdz/saturn branche

YM
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      IF (.NOT. assigned_domain(ind)) CYCLE
23      CALL swap_dimensions(ind)
24      CALL swap_geometry(ind)
25      phis=f_phis(ind)
26      pks=f_pks(ind)
27      pk=f_pk(ind)
28      theta=f_theta(ind)
29      phi=f_phi(ind)
30      CALL compute_geopotential(phis,pks,pk,theta,phi,0)
31    ENDDO
32 
33  END SUBROUTINE geopotential
34 
35  SUBROUTINE compute_geopotential(phis,pks,pk,theta,phi,offset)
36  USE icosa
37  IMPLICIT NONE
38    REAL(rstd),INTENT(IN) :: phis(iim*jjm)
39    REAL(rstd),INTENT(IN) :: pks(iim*jjm)
40    REAL(rstd),INTENT(IN) :: pk(iim*jjm,llm)
41    REAL(rstd),INTENT(IN) :: theta(iim*jjm,llm)
42    REAL(rstd),INTENT(OUT) :: phi(iim*jjm,llm)
43    INTEGER,INTENT(IN) :: offset
44    INTEGER :: i,j,ij,l
45
46!!! Compute geopotential
47  ! LMD-Z : dh=pi.dtheta-vdp  => (1/rho) dp = pi.dtheta - d(h=pi.theta) = -dtheta.gradpi
48  ! g.dz =  -(1/rho)dp =  theta.gradpi
49  ! direct (variational)
50  ! dz = -(1/rho.g)dp, rho = RT/p = (R/C_p).theta.pi/p
51   
52  ! for first layer
53   DO j=jj_begin-offset,jj_end+offset
54     DO i=ii_begin-offset,ii_end+offset
55       ij=(j-1)*iim+i
56       phi( ij,1 ) = phis( ij ) + theta(ij,1) * ( pks(ij) - pk(ij,1) )
57     ENDDO
58   ENDDO
59         
60  ! for other layers
61   DO l = 2, llm
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         phi(ij,l) = phi(ij,l-1) + 0.5 * ( theta(ij,l)  + theta(ij,l-1) )  & 
66                                       * (  pk(ij,l-1) -  pk(ij,l)    )
67       ENDDO
68     ENDDO
69   ENDDO
70   
71  END SUBROUTINE compute_geopotential
72
73END MODULE geopotential_mod
Note: See TracBrowser for help on using the repository browser.