Changeset 914 for codes/icosagcm/devel

Timestamp:

06/17/19 18:05:22 (5 years ago)

Author:

dubos

Message:

devel : moved diagnostics of temperature and velocity into separate modules

Location:

codes/icosagcm/devel/src/diagnostics

Files:

• compute_temperature.F90 (added)
• compute_velocity.F90 (added)
• observable.f90 (modified) (4 diffs)

codes/icosagcm/devel/src/diagnostics/observable.f90

@@ -44 +44 @@
     USE earth_const
     USE compute_pression_mod, ONLY : pression_mid
+    USE compute_temperature_mod
+    USE compute_velocity_mod
     USE vertical_interp_mod
     USE theta2theta_rhodz_mod
@@ -97 +99 @@
 
     CALL pression_mid(f_ps, f_pmid)
-    CALL diagnose_temperature(f_pmid, f_q, f_buf_i) ! f_buf_i : IN = theta, out = T
+    CALL temperature(f_pmid, f_q, f_buf_i) ! f_buf_i : IN = theta, out = T
 
     IF(init) THEN
@@ -125 +127 @@
     IF(grid_type == grid_unst) RETURN
 
-    CALL progonostic_vel_to_horiz(f_geopot, f_ps, f_mass, f_u, f_W, f_buf_uh, f_buf_i)
+    CALL velocity(f_geopot, f_ps, f_mass, f_u, f_W, f_buf_Fel, f_buf_uh, f_buf_i)
     CALL transfert_request(f_buf_uh,req_e1_vect) 
     CALL un2ulonlat(f_buf_uh, f_buf_ulon, f_buf_ulat)
@@ -210 +212 @@
   END SUBROUTINE output_energyflux
   
- !------------------- Conversion from prognostic to observable variables ------------------
-
-  SUBROUTINE progonostic_vel_to_horiz(f_geopot, f_ps, f_rhodz, f_u, f_W, f_uh, f_uz)
-    USE disvert_mod, ONLY : caldyn_eta, eta_mass
-    USE compute_diagnostics_mod, ONLY : compute_rhodz
-    TYPE(t_field), POINTER :: f_geopot(:), f_ps(:), f_rhodz(:), &
-         f_u(:), f_W(:), f_uz(:), &  ! IN
-         f_uh(:)                         ! OUT
-    REAL(rstd),POINTER :: geopot(:,:), ps(:), rhodz(:,:), u(:,:), W(:,:), uh(:,:), uz(:,:), F_el(:,:)
-    INTEGER :: ind
-    
-    DO ind=1,ndomain
-       IF (.NOT. assigned_domain(ind)) CYCLE
-       CALL swap_dimensions(ind)
-       CALL swap_geometry(ind)
-       geopot = f_geopot(ind)
-       rhodz = f_rhodz(ind)
-       u = f_u(ind)
-       W = f_W(ind)
-       uh  = f_uh(ind)
-       F_el  = f_buf_Fel(ind)
-       IF(caldyn_eta==eta_mass) THEN
-          ps=f_ps(ind)
-          CALL compute_rhodz(.TRUE., ps, rhodz)
-       END IF
-       uz = f_uz(ind)
-       !$OMP BARRIER
-       CALL compute_prognostic_vel_to_horiz(geopot,rhodz,u,W, F_el, uh,uz)
-       !$OMP BARRIER
-    END DO
-  END SUBROUTINE progonostic_vel_to_horiz
-  
-  SUBROUTINE compute_prognostic_vel_to_horiz(Phi, rhodz, u, W, F_el, uh, uz)
-    USE omp_para
-    REAL(rstd), INTENT(IN) :: Phi(iim*jjm,llm+1)
-    REAL(rstd), INTENT(IN) :: rhodz(iim*jjm,llm)
-    REAL(rstd), INTENT(IN) :: u(3*iim*jjm,llm)
-    REAL(rstd), INTENT(IN) :: W(iim*jjm,llm+1)
-    REAL(rstd), INTENT(OUT) :: uh(3*iim*jjm,llm)
-    REAL(rstd), INTENT(OUT) :: uz(iim*jjm,llm)
-    INTEGER :: ij,l
-    REAL(rstd) :: F_el(3*iim*jjm,llm+1)
-    REAL(rstd) :: uu_right, uu_lup, uu_ldown, W_el, DePhil
-    ! NB : u and uh are not in DEC form, they are normal components    
-    ! => we must divide by de
-    IF(hydrostatic) THEN
-       uh(:,:)=u(:,:)
-       uz(:,:)=0.
-    ELSE
-       DO l=ll_begin, ll_endp1 ! compute on l levels (interfaces)
-          DO ij=ij_begin_ext, ij_end_ext
-             ! Compute on edge 'right'
-             W_el   = .5*( W(ij,l)+W(ij+t_right,l) )
-             DePhil = ne_right*(Phi(ij+t_right,l)-Phi(ij,l))
-             F_el(ij+u_right,l)   = DePhil*W_el/de(ij+u_right)
-             ! Compute on edge 'lup'
-             W_el   = .5*( W(ij,l)+W(ij+t_lup,l) )
-             DePhil = ne_lup*(Phi(ij+t_lup,l)-Phi(ij,l))
-             F_el(ij+u_lup,l)   = DePhil*W_el/de(ij+u_lup)
-             ! Compute on edge 'ldown'
-             W_el   = .5*( W(ij,l)+W(ij+t_ldown,l) )
-             DePhil = ne_ldown*(Phi(ij+t_ldown,l)-Phi(ij,l))
-             F_el(ij+u_ldown,l) = DePhil*W_el/de(ij+u_ldown)
-          END DO
-       END DO
-
-       ! We need a barrier here because we compute F_el above and do a vertical average below
-       !$OMP BARRIER
-
-       DO l=ll_begin, ll_end ! compute on k levels (full levels)
-          DO ij=ij_begin_ext, ij_end_ext
-             ! w = vertical momentum = g^-2*dPhi/dt = uz/g
-             uz(ij,l) = (.5*g)*(W(ij,l)+W(ij,l+1))/rhodz(ij,l)
-             ! uh = u-w.grad(Phi) = u - uz.grad(z)
-             uh(ij+u_right,l) = u(ij+u_right,l) - (F_el(ij+u_right,l)+F_el(ij+u_right,l+1)) / (rhodz(ij,l)+rhodz(ij+t_right,l))
-             uh(ij+u_lup,l)   = u(ij+u_lup,l)   - (F_el(ij+u_lup,l)+F_el(ij+u_lup,l+1))     / (rhodz(ij,l)+rhodz(ij+t_lup,l))
-             uh(ij+u_ldown,l) = u(ij+u_ldown,l) - (F_el(ij+u_ldown,l)+F_el(ij+u_ldown,l+1)) / (rhodz(ij,l)+rhodz(ij+t_ldown,l))
-          END DO
-       END DO
-
-    END IF
-  END SUBROUTINE compute_prognostic_vel_to_horiz
-
-  SUBROUTINE diagnose_temperature(f_pmid,f_q,f_temp)
-    USE icosa
-    IMPLICIT NONE
-    TYPE(t_field), POINTER :: f_pmid(:)         ! IN
-    TYPE(t_field), POINTER :: f_q(:)            ! IN
-    TYPE(t_field), POINTER :: f_temp(:)         ! INOUT
-    
-    REAL(rstd), POINTER :: pmid(:,:)
-    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)
-       pmid=f_pmid(ind)
-       q=f_q(ind)
-       temp=f_temp(ind)
-       CALL compute_diagnose_temp(pmid,q,temp)
-    END DO
-  END SUBROUTINE diagnose_temperature
-  
-  SUBROUTINE compute_diagnose_temp(pmid,q,temp)
-    USE omp_para
-    REAL(rstd),INTENT(IN)    :: pmid(iim*jjm,llm)
-    REAL(rstd),INTENT(IN)    :: q(iim*jjm,llm,nqtot)
-    REAL(rstd),INTENT(INOUT) :: temp(iim*jjm,llm)
-
-    REAL(rstd) :: Rd, p_ik, theta_ik, temp_ik, qv, chi, Rmix
-    INTEGER :: ij,l
-
-    Rd = kappa*cpp
-    DO l=ll_begin,ll_end
-       DO ij=ij_begin,ij_end
-          p_ik = pmid(ij,l)
-          theta_ik = temp(ij,l)
-          qv = q(ij,l,1) ! water vapor 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(:)
-    TYPE(t_field), POINTER :: f_q(:)
-    TYPE(t_field), POINTER :: f_T(:)
-    
-    REAL(rstd),POINTER :: Tv(:,:), q(:,:,:), T(:,:)
-    INTEGER :: ind
-    
-    DO ind=1,ndomain
-       IF (.NOT. assigned_domain(ind)) CYCLE
-       CALL swap_dimensions(ind)
-       CALL swap_geometry(ind)
-       Tv=f_Tv(ind)
-       T=f_T(ind)
-       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
-
   SUBROUTINE divide_by_mass(iq, f_mass, f_theta_rhodz, f_theta)
     INTEGER, INTENT(IN) :: iq