Changeset 434 for codes/icosagcm/trunk/src/observable.f90

Timestamp:

06/15/16 22:35:17 (8 years ago)

Author:

dubos

Message:

Dry/moist output cleanup

File:

• codes/icosagcm/trunk/src/observable.f90 (modified) (3 diffs)

codes/icosagcm/trunk/src/observable.f90

@@ -93 +93 @@
     END IF
 
-    CALL theta_rhodz2temperature(f_ps,f_theta_rhodz,f_buf_i) ; 
-    CALL Tv2T(f_buf_i,f_q,f_buf1_i)
+!    CALL theta_rhodz2temperature(f_ps,f_theta_rhodz,f_buf_i)
+!    CALL Tv2T(f_buf_i,f_q,f_buf1_i)
+    CALL diagnose_temperature(f_ps, f_theta_rhodz, f_q, f_buf_i)
+
     IF(init) THEN
        CALL output_field("temp_init",f_buf_i)
@@ -234 +236 @@
   END SUBROUTINE compute_prognostic_vel_to_horiz
 
+  SUBROUTINE diagnose_temperature(f_ps,f_theta_rhodz,f_q,f_temp)
+    USE icosa
+    USE pression_mod
+    IMPLICIT NONE
+    TYPE(t_field), POINTER :: f_ps(:)           ! IN
+    TYPE(t_field), POINTER :: f_theta_rhodz(:)  ! IN
+    TYPE(t_field), POINTER :: f_q(:)            ! IN
+    TYPE(t_field), POINTER :: f_temp(:)         ! OUT
+    
+    REAL(rstd), POINTER :: ps(:)
+    REAL(rstd), POINTER :: theta_rhodz(:,:,:)
+    REAL(rstd), POINTER :: q(:,:,:)
+    REAL(rstd), POINTER :: temp(:,:)
+    INTEGER :: ind
+    
+    DO ind=1,ndomain
+       IF (.NOT. assigned_domain(ind)) CYCLE
+       CALL swap_dimensions(ind)
+       CALL swap_geometry(ind)
+       ps=f_ps(ind)
+       theta_rhodz=f_theta_rhodz(ind)
+       q=f_q(ind)
+       temp=f_temp(ind)
+       CALL compute_diagnose_temp(ps,theta_rhodz,q,temp)
+    END DO
+  END SUBROUTINE diagnose_temperature
+  
+  SUBROUTINE compute_diagnose_temp(ps,theta_rhodz,q,temp)
+    USE omp_para
+    USE pression_mod
+    REAL(rstd),INTENT(IN)  :: ps(iim*jjm)
+    REAL(rstd),INTENT(IN)  :: theta_rhodz(iim*jjm,llm,nqdyn)
+    REAL(rstd),INTENT(IN)  :: q(iim*jjm,llm,nqtot)
+    REAL(rstd),INTENT(OUT) :: temp(iim*jjm,llm)
+
+    REAL(rstd)  :: p(iim*jjm,llm+1)
+    REAL(rstd) :: Rd, p_ik, theta_ik, temp_ik, qv, chi, Rmix
+    INTEGER :: ij,l
+
+    Rd = kappa*cpp
+    CALL compute_pression(ps,p,0)
+    DO l=ll_begin,ll_end
+       DO ij=ij_begin,ij_end
+          p_ik = .5*(p(ij,l)+p(ij,l+1))
+          theta_ik = g*theta_rhodz(ij,l,1)/(p(ij,l)-p(ij,l+1))
+          qv = q(ij,l,1) ! water vaper mixing ratio = mv/md
+          SELECT CASE(caldyn_thermo)
+          CASE(thermo_theta)
+             temp_ik = theta_ik*((p_ik/preff)**kappa)
+          CASE(thermo_entropy)
+             temp_ik = Treff*exp((theta_ik + Rd*log(p_ik/preff))/cpp)
+          CASE(thermo_moist)
+             Rmix = Rd+qv*Rv
+             chi = ( theta_ik + Rmix*log(p_ik/preff) ) / (cpp + qv*cppv) ! log(T/Treff)
+             temp_ik = Treff*exp(chi)
+          END SELECT
+          IF(physics_thermo==thermo_fake_moist) temp_ik=temp_ik/(1+0.608*qv)
+          temp(ij,l)=temp_ik
+       END DO
+    END DO
+  END SUBROUTINE compute_diagnose_temp
+
   SUBROUTINE Tv2T(f_Tv, f_q, f_T)
     TYPE(t_field), POINTER :: f_TV(:)
@@ -247 +311 @@
        CALL swap_geometry(ind)
        Tv=f_Tv(ind)
-       q=f_q(ind)
        T=f_T(ind)
-       T=Tv/(1+0.608*q(:,:,1))
-    END DO
-    
+       SELECT CASE(physics_thermo)
+       CASE(thermo_dry)
+          T=Tv
+       CASE(thermo_fake_moist)
+          q=f_q(ind)
+          T=Tv/(1+0.608*q(:,:,1))
+       END SELECT
+    END DO
   END SUBROUTINE Tv2T