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

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

Add new openMP parallelism based on distribution of domains on threads. There is no more limitation of number of threads by MPI process.

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.