Changeset 434

Timestamp:

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

Author:

dubos

Message:

Dry/moist output cleanup

Location:

codes/icosagcm/trunk/src

Files:

• observable.f90 (modified) (3 diffs)
• physics_dcmip2016.f90 (modified) (3 diffs)
• theta_rhodz.f90 (modified) (1 diff)

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
 

codes/icosagcm/trunk/src/physics_dcmip2016.f90

@@ -34 +34 @@
   SUBROUTINE init_physics
     USE physics_interface_mod
+    USE omp_para
     IMPLICIT NONE
     INTEGER :: ngrid
     CHARACTER(LEN=255) :: testcase_str 
     
+    testcase_str='undefined'
     CALL getin("physics_dcmip2016",testcase_str)
    
@@ -53 +55 @@
       CASE DEFAULT
          PRINT*, 'Bad selector for dcmip2016 test case <', testcase_str, &
-              '> options are <dry_baroclinic>, <moist_baroclinic>, <cyclone>, <supercell>'
+              '> options are <dry_baroclinic>, <moist_baroclinic_full>, <moist_baroclinic_kessler>, <cyclone>, <supercell>'
          STOP
     END SELECT
@@ -59 +61 @@
     PBL=.FALSE.
     CALL getin("physics_dcmip2016_PBL",PBL)
+    IF(is_master) THEN
+       IF(PBL) THEN
+          PRINT *, 'PBL=.TRUE., Bryan PBL activated'
+       ELSE
+          PRINT *, 'PBL=.FALSE., Reed & Jablonowski PBL activated'
+       END IF
+    END IF
 
     ngrid = physics_inout%ngrid

codes/icosagcm/trunk/src/theta_rhodz.f90

@@ -48 +48 @@
   USE icosa
   USE pression_mod
-  USE exner_mod
   IMPLICIT NONE
     TYPE(t_field), POINTER :: f_ps(:)