Changeset 366 for codes/icosagcm/trunk/src/etat0.f90

Timestamp:

10/30/15 15:41:06 (9 years ago)

Author:

dubos

Message:

Progress towards NH

File:

• codes/icosagcm/trunk/src/etat0.f90 (modified) (16 diffs)

codes/icosagcm/trunk/src/etat0.f90

@@ -13 +13 @@
 CONTAINS
   
-  SUBROUTINE etat0(f_ps,f_mass,f_phis,f_theta_rhodz,f_u, f_q)
+  SUBROUTINE etat0(f_ps,f_mass,f_phis,f_theta_rhodz,f_u, f_geopot,f_w, f_q)
     USE mpipara, ONLY : is_mpi_root
     USE disvert_mod
@@ -35 +35 @@
     TYPE(t_field),POINTER :: f_theta_rhodz(:)
     TYPE(t_field),POINTER :: f_u(:)
+    TYPE(t_field),POINTER :: f_geopot(:)
+    TYPE(t_field),POINTER :: f_w(:)
     TYPE(t_field),POINTER :: f_q(:)
     
     REAL(rstd),POINTER :: ps(:), mass(:,:)
-    LOGICAL :: init_mass, collocated
+    LOGICAL :: autoinit_mass, autoinit_geopot, collocated
     INTEGER :: ind,i,j,ij,l
 
@@ -45 +47 @@
     ! the initial distribution of mass is taken to be the same
     ! as what the mass coordinate would dictate
-    ! however if etat0_XXX defines mass then the flag init_mass must be set to .FALSE.
+    ! however if etat0_XXX defines mass then the flag autoinit_mass must be set to .FALSE.
     ! otherwise mass will be overwritten 
-    init_mass = (caldyn_eta == eta_lag)
+    autoinit_mass = (caldyn_eta == eta_lag)
 
     etat0_type='jablonowsky06'
@@ -70 +72 @@
         CALL getin_etat0_dcmip5
     CASE ('williamson91.6')
-       init_mass=.FALSE.
+       autoinit_mass=.FALSE.
        CALL getin_etat0_williamson
     CASE DEFAULT
        collocated=.FALSE.
+       autoinit_mass = .FALSE.
     END SELECT
 
@@ -90 +93 @@
    CASE DEFAULT
       IF(collocated) THEN
-         CALL etat0_collocated(f_phis,f_ps,f_mass,f_theta_rhodz,f_u, f_q)
+         CALL etat0_collocated(f_phis,f_ps,f_mass,f_theta_rhodz,f_u, f_geopot,f_W, f_q)
       ELSE
          PRINT*, 'Bad selector for variable etat0 <',etat0_type, &
@@ -98 +101 @@
     END SELECT
 
-    IF(init_mass) THEN ! initialize mass distribution using ps
 !       !$OMP BARRIER
+    IF(autoinit_mass) THEN
        DO ind=1,ndomain
           IF (.NOT. assigned_domain(ind)) CYCLE
@@ -105 +108 @@
           CALL swap_geometry(ind)
           mass=f_mass(ind); ps=f_ps(ind)
-          CALL compute_rhodz(.TRUE., ps, mass)
+          CALL compute_rhodz(.TRUE., ps, mass) ! initialize mass distribution using ps
        END DO
     END IF
-
+ 
   END SUBROUTINE etat0
 
-  SUBROUTINE etat0_collocated(f_phis,f_ps,f_mass,f_theta_rhodz,f_u, f_q)
+  SUBROUTINE etat0_collocated(f_phis,f_ps,f_mass,f_theta_rhodz,f_u, f_geopot,f_W, f_q)
     USE theta2theta_rhodz_mod
     IMPLICIT NONE
@@ -119 +122 @@
     TYPE(t_field),POINTER :: f_theta_rhodz(:)
     TYPE(t_field),POINTER :: f_u(:)
+    TYPE(t_field),POINTER :: f_geopot(:)
+    TYPE(t_field),POINTER :: f_W(:)
     TYPE(t_field),POINTER :: f_q(:)
   
@@ -128 +133 @@
     REAL(rstd),POINTER :: temp(:,:)
     REAL(rstd),POINTER :: u(:,:)
+    REAL(rstd),POINTER :: geopot(:,:)
+    REAL(rstd),POINTER :: W(:,:)
     REAL(rstd),POINTER :: q(:,:,:)
     INTEGER :: ind
@@ -143 +150 @@
       temp=f_temp(ind)
       u=f_u(ind)
+      geopot=f_geopot(ind)
+      w=f_w(ind)
       q=f_q(ind)
 
-      IF( TRIM(etat0_type)=='williamson91.6' ) THEN 
-       CALL compute_etat0_collocated(ps,mass, phis, theta_rhodz, u, q)
+      IF( TRIM(etat0_type)=='williamson91.6' ) THEN
+       CALL compute_etat0_collocated(ps,mass, phis, theta_rhodz, u, geopot, W, q)
       ELSE
-       CALL compute_etat0_collocated(ps,mass, phis, temp, u, q)
+       CALL compute_etat0_collocated(ps,mass, phis, temp, u, geopot, W, q)
       ENDIF
     ENDDO
@@ -158 +167 @@
   END SUBROUTINE etat0_collocated
 
-  SUBROUTINE compute_etat0_collocated(ps,mass, phis, temp_i, u, q)
+  SUBROUTINE compute_etat0_collocated(ps,mass,phis,temp_i,u, geopot,W, q)
     USE wind_mod
+    USE disvert_mod
     USE etat0_jablonowsky06_mod, ONLY : compute_jablonowsky06 => compute_etat0
     USE etat0_dcmip1_mod, ONLY : compute_dcmip1 => compute_etat0
@@ -174 +184 @@
     REAL(rstd),INTENT(OUT) :: temp_i(iim*jjm,llm)
     REAL(rstd),INTENT(OUT) :: u(3*iim*jjm,llm)
+    REAL(rstd),INTENT(OUT) :: W(iim*jjm,llm+1)
+    REAL(rstd),INTENT(OUT) :: geopot(iim*jjm,llm+1)
     REAL(rstd),INTENT(OUT) :: q(iim*jjm,llm,nqtot)
 
@@ -183 +195 @@
     REAL(rstd) :: phis_e(3*iim*jjm)
     REAL(rstd) :: temp_e(3*iim*jjm,llm)
+    REAL(rstd) :: geopot_e(3*iim*jjm,llm+1)
     REAL(rstd) :: ulon_e(3*iim*jjm,llm)
     REAL(rstd) :: ulat_e(3*iim*jjm,llm)
@@ -188 +201 @@
 
     INTEGER :: l,i,j,ij
+    REAL :: p_ik, v_ik, mass_ik
+    LOGICAL :: autoinit_mass, autoinit_NH
+
+    ! For NH geopotential and vertical momentum must be initialized.
+    ! Unless autoinit_mass is set to .FALSE. , they will be initialized
+    ! to hydrostatic geopotential and zero
+    autoinit_mass = .TRUE.
+    autoinit_NH = .NOT. hydrostatic
+    w(:,:) = 0
 
     !$OMP BARRIER
@@ -208 +230 @@
        CALL compute_dcmip2(3*iim*jjm,lon_e,lat_e, phis_e, ps_e, temp_e, ulon_e, ulat_e)      
     CASE('dcmip3')
-       CALL compute_dcmip3(iim*jjm,lon_i,lat_i, phis, ps, temp_i, ulon_i, ulat_i, q)
-       CALL compute_dcmip3(3*iim*jjm,lon_e,lat_e, phis_e, ps_e, temp_e, ulon_e, ulat_e, q_e)
+       CALL compute_dcmip3(iim*jjm,lon_i,lat_i, phis, ps, temp_i, ulon_i, ulat_i, geopot, q)
+       CALL compute_dcmip3(3*iim*jjm,lon_e,lat_e, phis_e, ps_e, temp_e, ulon_e, ulat_e, geopot_e, q_e)
+       autoinit_NH = .FALSE. ! compute_dcmip3 initializes geopot
     CASE('dcmip4')
        CALL compute_dcmip4(iim*jjm,lon_i,lat_i, phis, ps, temp_i, ulon_i, ulat_i, q)
@@ -219 +242 @@
        CALL compute_w91_6(iim*jjm,lon_i,lat_i, phis, mass(:,1), temp_i(:,1), ulon_i(:,1), ulat_i(:,1))
        CALL compute_w91_6(3*iim*jjm,lon_e,lat_e, phis_e, mass_e(:,1), temp_e(:,1), ulon_e(:,1), ulat_e(:,1))
+       autoinit_mass = .FALSE. ! do not overwrite mass
     END SELECT
+
+    IF(autoinit_mass) CALL compute_rhodz(.TRUE., ps, mass) ! initialize mass distribution using ps
+    IF(autoinit_NH) THEN
+       geopot(:,1) = phis(:) ! surface geopotential
+       DO l = 1, llm
+          DO ij=1,iim*jjm
+             ! hybrid pressure coordinate
+             p_ik = ptop + mass_ak(l) + mass_bk(l)*ps(ij)
+             mass_ik = (mass_dak(l) + mass_dbk(l)*ps(ij))/g
+             ! v=R.T/p, R=kappa*cpp
+             v_ik = kappa*cpp*temp_i(ij,l)/p_ik
+             geopot(ij,l+1) = geopot(ij,l) + mass_ik*v_ik*g
+          END DO
+       END DO
+    END IF
 
     !$OMP BARRIER