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

Last change on this file since 352 was 295, checked in by ymipsl, 10 years ago

Merging OpenMP parallisme mode : by subdomain and on vertical level.
This feature is actually experimental but may be retro-compatible with the last method based only on subdomain

YM
File size: 2.2 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  USE omp_para
38  IMPLICIT NONE
39    REAL(rstd),INTENT(IN) :: phis(iim*jjm)
40    REAL(rstd),INTENT(IN) :: pks(iim*jjm)
41    REAL(rstd),INTENT(IN) :: pk(iim*jjm,llm)
42    REAL(rstd),INTENT(IN) :: theta(iim*jjm,llm)
43    REAL(rstd),INTENT(OUT) :: phi(iim*jjm,llm)
44    INTEGER,INTENT(IN) :: offset
45    INTEGER :: i,j,ij,l
46
47!!! Compute geopotential
48  ! LMD-Z : dh=pi.dtheta-vdp  => (1/rho) dp = pi.dtheta - d(h=pi.theta) = -dtheta.gradpi
49  ! g.dz =  -(1/rho)dp =  theta.gradpi
50  ! direct (variational)
51  ! dz = -(1/rho.g)dp, rho = RT/p = (R/C_p).theta.pi/p
52   
53  ! for first layer
54
55! flush pks, pk thetha, phis
56!$OMP BARRIER
57   IF(is_omp_level_master) THEN
58     DO j=jj_begin-offset,jj_end+offset
59       DO i=ii_begin-offset,ii_end+offset
60         ij=(j-1)*iim+i
61         phi( ij,1 ) = phis( ij ) + theta(ij,1) * ( pks(ij) - pk(ij,1) )
62       ENDDO
63     ENDDO
64         
65    ! for other layers
66     DO l = 2, llm
67       DO j=jj_begin-offset,jj_end+offset
68         DO i=ii_begin-offset,ii_end+offset
69           ij=(j-1)*iim+i
70           phi(ij,l) = phi(ij,l-1) + 0.5 * ( theta(ij,l)  + theta(ij,l-1) )  & 
71                                         * (  pk(ij,l-1) -  pk(ij,l)    )
72         ENDDO
73       ENDDO
74     ENDDO
75   ENDIF
76! flush phi
77!$OMP BARRIER
78   
79  END SUBROUTINE compute_geopotential
80
81END MODULE geopotential_mod
Note: See TracBrowser for help on using the repository browser.