Changeset 428 for codes/icosagcm/trunk/src/geopotential_mod.f90

Timestamp:

06/14/16 21:54:26 (8 years ago)

Author:

dubos

Message:

theta-related cleanup

File:

• codes/icosagcm/trunk/src/geopotential_mod.f90 (modified) (2 diffs)

codes/icosagcm/trunk/src/geopotential_mod.f90

@@ -2 +2 @@
   IMPLICIT NONE
   PRIVATE
-
+  
   PUBLIC :: compute_geopotential
 CONTAINS
+  
+  SUBROUTINE thetarhodz2geopot(f_ps,f_phis,f_theta_rhodz, f_p,f_theta,f_phi) ! FIXME : never called, dry only
+    USE icosa
+    USE omp_para
+    USE pression_mod
+    USE theta2theta_rhodz_mod
+    TYPE(t_field), POINTER :: f_ps(:), f_phis(:), f_theta_rhodz(:), &  ! IN                                       
+         f_p(:), f_theta(:), f_phi(:)               ! OUT                                                         
+    REAL(rstd),POINTER :: pk(:,:), p(:,:), theta(:,:,:), theta_rhodz(:,:,:), &
+         phi(:,:), phis(:), ps(:)
+    INTEGER :: ind
+    
+    DO ind=1,ndomain
+       IF (.NOT. assigned_domain(ind)) CYCLE
+       CALL swap_dimensions(ind)
+       CALL swap_geometry(ind)
+       ps = f_ps(ind)
+       p  = f_p(ind)
+!$OMP BARRIER                                                                                                     
+       CALL compute_pression(ps,p,0)
+!$OMP BARRIER                                                                                                     
+       theta_rhodz = f_theta_rhodz(ind)
+       theta = f_theta(ind)
+       CALL compute_theta_rhodz2theta(ps, theta_rhodz,theta,0)
+       phis = f_phis(ind)
+       phi = f_phi(ind)
+       CALL compute_geopotential(phis,p,theta,phi,0)
+    END DO
 
-  SUBROUTINE geopotential(f_phis,f_pks,f_pk,f_theta,f_phi) ! ORPHAN
-  USE icosa
-    TYPE(t_field), POINTER :: f_phis(:)
-    TYPE(t_field), POINTER :: f_pks(:)
-    TYPE(t_field), POINTER :: f_pk(:)
-    TYPE(t_field), POINTER :: f_theta(:)
-    TYPE(t_field), POINTER :: f_phi(:)
-  
-    REAL(rstd), POINTER :: phis(:)
-    REAL(rstd), POINTER :: pks(:)
-    REAL(rstd), POINTER :: pk(:,:)
-    REAL(rstd), POINTER :: theta(:,:)
-    REAL(rstd), POINTER :: phi(:,:)
-    INTEGER :: ind
+  END SUBROUTINE thetarhodz2geopot
 
-    DO ind=1,ndomain
-      IF (.NOT. assigned_domain(ind)) CYCLE
-      CALL swap_dimensions(ind)
-      CALL swap_geometry(ind)
-      phis=f_phis(ind)
-      pks=f_pks(ind)
-      pk=f_pk(ind)
-      theta=f_theta(ind)
-      phi=f_phi(ind)
-      CALL compute_geopotential(phis,pks,pk,theta,phi,0)
-    ENDDO
-  
-  END SUBROUTINE geopotential
-  
-  SUBROUTINE compute_geopotential(phis,pks,pk,theta,phi,offset)
-  USE icosa
-  USE omp_para
-  IMPLICIT NONE
+
+  SUBROUTINE compute_geopotential(phis,p,theta,phi,offset)
+    USE icosa
+    USE omp_para
+    REAL(rstd),INTENT(IN) :: p(iim*jjm,llm+1)
     REAL(rstd),INTENT(IN) :: phis(iim*jjm)
-    REAL(rstd),INTENT(IN) :: pks(iim*jjm)
-    REAL(rstd),INTENT(IN) :: pk(iim*jjm,llm)
-    REAL(rstd),INTENT(IN) :: theta(iim*jjm,llm)
+    REAL(rstd),INTENT(IN) :: theta(iim*jjm,llm,nqdyn)
     REAL(rstd),INTENT(OUT) :: phi(iim*jjm,llm)
     INTEGER,INTENT(IN) :: offset
     INTEGER :: i,j,ij,l
+    REAL(rstd) :: mg_ij, p_ij, exner_ij
 
 !!! Compute geopotential
@@ -57 +58 @@
 ! flush pks, pk thetha, phis
 !$OMP BARRIER
-   IF(is_omp_level_master) THEN
-     DO j=jj_begin-offset,jj_end+offset
-       DO i=ii_begin-offset,ii_end+offset
-         ij=(j-1)*iim+i
-         phi( ij,1 ) = phis( ij ) + theta(ij,1) * ( pks(ij) - pk(ij,1) )
-       ENDDO
-     ENDDO
+    IF(is_omp_level_master) THEN
+       DO j=jj_begin-offset,jj_end+offset
+          DO i=ii_begin-offset,ii_end+offset
+             ij=(j-1)*iim+i
+             phi( ij,1 ) = phis( ij )
+          END DO
+       END DO
           
     ! for other layers
-     DO l = 2, llm
+     DO l = 1, llm-1
        DO j=jj_begin-offset,jj_end+offset
          DO i=ii_begin-offset,ii_end+offset
            ij=(j-1)*iim+i
-           phi(ij,l) = phi(ij,l-1) + 0.5 * ( theta(ij,l)  + theta(ij,l-1) )  & 
-                                         * (  pk(ij,l-1) -  pk(ij,l)    )
+           mg_ij = p(ij,l+1)-p(ij,l)
+           ! v=RT/p=(kappa.cpp).Theta.(p/preff)**kappa /p
+           p_ij=.5*(p(ij,l+1)+p(ij,l))
+           exner_ij = (p_ij/preff)**kappa ! NB : no cpp factor
+           phi(ij,l+1) = phi(ij,l) + (kappa*cpp)*mg_ij*theta(ij,l,1)*exner_ij/p_ij
          ENDDO
        ENDDO