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

Timestamp:

06/08/16 01:51:21 (8 years ago)

Author:

dubos

Message:

Introduced entropy as prognostic variable - tested with JW06

File:

• codes/icosagcm/trunk/src/etat0.f90 (modified) (1 diff)

codes/icosagcm/trunk/src/etat0.f90

@@ -172 +172 @@
          CALL compute_etat0_collocated(ps,mass, phis, temp, u, geopot, W, q)
       ENDIF
+
+      IF( TRIM(etat0_type)/='williamson91.6' ) CALL compute_temperature2entropy(ps,temp,q,theta_rhodz, 1)
+    
     ENDDO
     
-    IF( TRIM(etat0_type)/='williamson91.6' ) CALL temperature2theta_rhodz(f_ps,f_temp,f_theta_rhodz)
-    
     CALL deallocate_field(f_temp)
     
   END SUBROUTINE etat0_collocated
+
+  SUBROUTINE compute_temperature2entropy(ps,temp,q,theta_rhodz,offset)
+    USE icosa
+    USE pression_mod
+    USE exner_mod
+    USE omp_para
+    IMPLICIT NONE
+    REAL(rstd),INTENT(IN)  :: ps(iim*jjm)
+    REAL(rstd),INTENT(IN)  :: temp(iim*jjm,llm)
+    REAL(rstd),INTENT(IN)  :: q(iim*jjm,llm,nqtot)
+    REAL(rstd),INTENT(OUT) :: theta_rhodz(iim*jjm,llm)
+    INTEGER,INTENT(IN) :: offset
+
+    REAL(rstd) :: p(iim*jjm,llm+1)
+    REAL(rstd) :: cppd,Rd, mass, p_ij, q_ij,r_ij, chi,nu, entropy, theta
+    INTEGER :: i,j,ij,l
+
+    cppd=cpp
+    Rd=kappa*cppd
+
+    CALL compute_pression(ps,p,offset)
+    ! flush p
+    !$OMP BARRIER
+    DO    l    = ll_begin, ll_end
+       DO j=jj_begin-offset,jj_end+offset
+          DO i=ii_begin-offset,ii_end+offset
+             ij=(j-1)*iim+i
+             mass = (p(ij,l)-p(ij,l+1))/g ! dry+moist mass
+             p_ij = .5*(p(ij,l)+p(ij,l+1))  ! pressure at full level
+             SELECT CASE(caldyn_thermo)
+             CASE(thermo_theta)
+                theta = temp(ij,l)*(p_ij/preff)**(-kappa) 
+                theta_rhodz(ij,l) = mass * theta
+             CASE(thermo_entropy)
+                nu = log(p_ij/preff)
+                chi = log(temp(ij,l)/Treff)
+                entropy = cppd*chi-Rd*nu
+                theta_rhodz(ij,l) = mass * entropy
+!             CASE(thermo_moist)
+!                q_ij=q(ij,l,1)
+!                r_ij=1.-q_ij
+!                mass=mass*(1-q_ij) ! dry mass
+!                nu = log(p_ij/preff)
+!                chi = log(temp(ij,l)/Treff)
+!                entropy = r_ij*(cppd*chi-Rd*nu) + q_ij*(cppv*chi-Rv*nu)
+!                theta_rhodz(ij,l) = mass * entropy                
+                CASE DEFAULT
+                   STOP
+             END SELECT
+          ENDDO
+       ENDDO
+    ENDDO
+    !$OMP BARRIER  
+  END SUBROUTINE compute_temperature2entropy
 
   SUBROUTINE compute_etat0_collocated(ps,mass,phis,temp_i,u, geopot,W, q)