source: codes/icosagcm/trunk/src/diagnostics/geopotential_mod.f90 @ 899

Last change on this file since 899 was 899, checked in by adurocher, 5 years ago

trunk : Fixed GCC warnings

Fixed iso c bindings
fixed warnings with -Wall -Wno-aliasing -Wno-unused -Wno-unused-dummy-argument -Wno-maybe-uninitialized -Wno-tabs warnings
Removed all unused variables (-Wunused-variable)
vector%dot_product is now dot_product_3d to avoid compilation warning "dot_product shadows intrinsic" with GCC
File size: 2.8 KB
Line 
1MODULE geopotential_mod
2  IMPLICIT NONE
3  PRIVATE
4 
5  PUBLIC :: compute_geopotential
6CONTAINS
7 
8  SUBROUTINE thetarhodz2geopot(f_ps,f_phis,f_theta_rhodz, f_p,f_theta,f_phi) ! FIXME : never called, dry only
9    USE icosa
10    USE omp_para
11    USE pression_mod
12    USE theta2theta_rhodz_mod
13    TYPE(t_field), POINTER :: f_ps(:), f_phis(:), f_theta_rhodz(:), &  ! IN                                       
14         f_p(:), f_theta(:), f_phi(:)               ! OUT                                                         
15    REAL(rstd),POINTER :: p(:,:), theta(:,:,:), theta_rhodz(:,:,:), &
16         phi(:,:), phis(:), ps(:)
17    INTEGER :: ind
18   
19    DO ind=1,ndomain
20       IF (.NOT. assigned_domain(ind)) CYCLE
21       CALL swap_dimensions(ind)
22       CALL swap_geometry(ind)
23       ps = f_ps(ind)
24       p  = f_p(ind)
25!$OMP BARRIER                                                                                                     
26       CALL compute_pression(ps,p,0)
27!$OMP BARRIER                                                                                                     
28       theta_rhodz = f_theta_rhodz(ind)
29       theta = f_theta(ind)
30       CALL compute_theta_rhodz2theta(ps, theta_rhodz,theta,0)
31       phis = f_phis(ind)
32       phi = f_phi(ind)
33       CALL compute_geopotential(phis,p,theta,phi,0)
34    END DO
35
36  END SUBROUTINE thetarhodz2geopot
37
38
39  SUBROUTINE compute_geopotential(phis,p,theta,phi,offset)
40    USE icosa
41    USE omp_para
42    REAL(rstd),INTENT(IN) :: p(iim*jjm,llm+1)
43    REAL(rstd),INTENT(IN) :: phis(iim*jjm)
44    REAL(rstd),INTENT(IN) :: theta(iim*jjm,llm,nqdyn)
45    REAL(rstd),INTENT(OUT) :: phi(iim*jjm,llm)
46    INTEGER,INTENT(IN) :: offset
47    INTEGER :: i,j,ij,l
48    REAL(rstd) :: mg_ij, p_ij, exner_ij
49
50!!! Compute geopotential
51  ! LMD-Z : dh=pi.dtheta-vdp  => (1/rho) dp = pi.dtheta - d(h=pi.theta) = -dtheta.gradpi
52  ! g.dz =  -(1/rho)dp =  theta.gradpi
53  ! direct (variational)
54  ! dz = -(1/rho.g)dp, rho = RT/p = (R/C_p).theta.pi/p
55   
56  ! for first layer
57
58! flush pks, pk thetha, phis
59!$OMP BARRIER
60    IF(is_omp_level_master) THEN
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,1 ) = phis( ij )
65          END DO
66       END DO
67         
68    ! for other layers
69     DO l = 1, llm-1
70       DO j=jj_begin-offset,jj_end+offset
71         DO i=ii_begin-offset,ii_end+offset
72           ij=(j-1)*iim+i
73           mg_ij = p(ij,l+1)-p(ij,l)
74           ! v=RT/p=(kappa.cpp).Theta.(p/preff)**kappa /p
75           p_ij=.5*(p(ij,l+1)+p(ij,l))
76           exner_ij = (p_ij/preff)**kappa ! NB : no cpp factor
77           phi(ij,l+1) = phi(ij,l) + (kappa*cpp)*mg_ij*theta(ij,l,1)*exner_ij/p_ij
78         ENDDO
79       ENDDO
80     ENDDO
81   ENDIF
82! flush phi
83!$OMP BARRIER
84   
85  END SUBROUTINE compute_geopotential
86
87END MODULE geopotential_mod
Note: See TracBrowser for help on using the repository browser.