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

Last change on this file since 354 was 354, checked in by dubos, 9 years ago

Moved output of dyn fields out of caldyn_gcm
File size: 2.2 KB
Line 
1MODULE geopotential_mod
2  IMPLICIT NONE
3  PRIVATE
4
5  PUBLIC :: compute_geopotential
6CONTAINS
7
8  SUBROUTINE geopotential(f_phis,f_pks,f_pk,f_theta,f_phi) ! ORPHAN
9  USE icosa
10    TYPE(t_field), POINTER :: f_phis(:)
11    TYPE(t_field), POINTER :: f_pks(:)
12    TYPE(t_field), POINTER :: f_pk(:)
13    TYPE(t_field), POINTER :: f_theta(:)
14    TYPE(t_field), POINTER :: f_phi(:)
15 
16    REAL(rstd), POINTER :: phis(:)
17    REAL(rstd), POINTER :: pks(:)
18    REAL(rstd), POINTER :: pk(:,:)
19    REAL(rstd), POINTER :: theta(:,:)
20    REAL(rstd), POINTER :: phi(:,:)
21    INTEGER :: ind
22
23    DO ind=1,ndomain
24      IF (.NOT. assigned_domain(ind)) CYCLE
25      CALL swap_dimensions(ind)
26      CALL swap_geometry(ind)
27      phis=f_phis(ind)
28      pks=f_pks(ind)
29      pk=f_pk(ind)
30      theta=f_theta(ind)
31      phi=f_phi(ind)
32      CALL compute_geopotential(phis,pks,pk,theta,phi,0)
33    ENDDO
34 
35  END SUBROUTINE geopotential
36 
37  SUBROUTINE compute_geopotential(phis,pks,pk,theta,phi,offset)
38  USE icosa
39  USE omp_para
40  IMPLICIT NONE
41    REAL(rstd),INTENT(IN) :: phis(iim*jjm)
42    REAL(rstd),INTENT(IN) :: pks(iim*jjm)
43    REAL(rstd),INTENT(IN) :: pk(iim*jjm,llm)
44    REAL(rstd),INTENT(IN) :: theta(iim*jjm,llm)
45    REAL(rstd),INTENT(OUT) :: phi(iim*jjm,llm)
46    INTEGER,INTENT(IN) :: offset
47    INTEGER :: i,j,ij,l
48
49!!! Compute geopotential
50  ! LMD-Z : dh=pi.dtheta-vdp  => (1/rho) dp = pi.dtheta - d(h=pi.theta) = -dtheta.gradpi
51  ! g.dz =  -(1/rho)dp =  theta.gradpi
52  ! direct (variational)
53  ! dz = -(1/rho.g)dp, rho = RT/p = (R/C_p).theta.pi/p
54   
55  ! for first layer
56
57! flush pks, pk thetha, phis
58!$OMP BARRIER
59   IF(is_omp_level_master) THEN
60     DO j=jj_begin-offset,jj_end+offset
61       DO i=ii_begin-offset,ii_end+offset
62         ij=(j-1)*iim+i
63         phi( ij,1 ) = phis( ij ) + theta(ij,1) * ( pks(ij) - pk(ij,1) )
64       ENDDO
65     ENDDO
66         
67    ! for other layers
68     DO l = 2, llm
69       DO j=jj_begin-offset,jj_end+offset
70         DO i=ii_begin-offset,ii_end+offset
71           ij=(j-1)*iim+i
72           phi(ij,l) = phi(ij,l-1) + 0.5 * ( theta(ij,l)  + theta(ij,l-1) )  & 
73                                         * (  pk(ij,l-1) -  pk(ij,l)    )
74         ENDDO
75       ENDDO
76     ENDDO
77   ENDIF
78! flush phi
79!$OMP BARRIER
80   
81  END SUBROUTINE compute_geopotential
82
83END MODULE geopotential_mod
Note: See TracBrowser for help on using the repository browser.