Changeset 408

Timestamp:

06/09/16 02:15:40 (8 years ago)

Author:

ymipsl

Message:

Modif of dcmip16 baroclinic wave etat0

YM

File:

• codes/icosagcm/trunk/src/etat0_dcmip2016_baroclinic_wave.f90 (modified) (5 diffs)

codes/icosagcm/trunk/src/etat0_dcmip2016_baroclinic_wave.f90

@@ -3 +3 @@
   IMPLICIT NONE
   PRIVATE
-  
-  REAL(rstd),PARAMETER :: eta0=0.252
-  REAL(rstd),PARAMETER :: etat=0.2
-  REAL(rstd),PARAMETER :: ps0=1e5
-  REAL(rstd),PARAMETER :: u0=35
-  REAL(rstd),PARAMETER :: T0=288
-  REAL(rstd),PARAMETER :: DeltaT=4.8e5
-  REAL(rstd),PARAMETER :: Rd=287
-  REAL(rstd),PARAMETER :: Gamma=0.005
-  REAL(rstd),PARAMETER :: up0=1
-  REAL(rstd),PARAMETER :: lonc=Pi/9, latc=2*Pi/9, latw=2*Pi/9
-  REAL(rstd),PARAMETER :: pw=34000
-  REAL(rstd),PARAMETER :: q0=0.021
-  
+ 
   INTEGER,SAVE :: testcase
   !$OMP THREADPRIVATE(testcase)  
   
+  INTEGER :: perturbation
+  !$OMP THREADPRIVATE(perturbation)  
+
   PUBLIC getin_etat0, compute_etat0
 
@@ -27 +17 @@
     USE mpipara, ONLY : is_mpi_root
     USE tracer_mod
+    IMPLICIT NONE
+    LOGICAL :: is_moist
+    CHARACTER(LEN=255) :: str_perturbation
+    
     IF(nqtot<5) THEN
        IF (is_mpi_root)  THEN
@@ -33 +27 @@
        STOP
     END IF
-!    CALL set_advection_scheme(1,advect_none)
-!    CALL set_advection_scheme(2,advect_none)
-!    CALL set_advection_scheme(3,advect_none)
     
+    str_perturbation="exponential"
+    CALL getin("dcmip2016_baroclinic_wave_perturbation",str_perturbation)
+    IF (TRIM(str_perturbation)=="exponential") THEN
+      perturbation=0
+    ELSE IF (TRIM(str_perturbation)=="stream") THEN
+      perturbation=1
+    ENDIF
+    
+  
   END SUBROUTINE getin_etat0
 
@@ -43 +43 @@
     USE disvert_mod
     USE omp_para
-    USE terminator
-    USE tracer_mod
+    USE dcmip2016_baroclinic_wave_mod, ONLY : baroclinic_wave_test
+    USE earth_const
+    USE terminator, ONLY: initial_value_Terminator
+    IMPLICIT NONE
     INTEGER, INTENT(IN) :: ngrid
     REAL(rstd),INTENT(IN) :: lon(ngrid)
@@ -55 +57 @@
     REAL(rstd),INTENT(OUT) :: q(ngrid,llm,nqtot)
     
-    INTEGER :: l,ij
-    REAL(rstd) :: etal, etavl, etas, etavs, sinlat, coslat, &
-         Y, Tave, T, phis_ave, vort, r2, utot, &
-         dthetaodeta_ave, dthetaodeta, dthetaodlat, duodeta, K, r
+    INTEGER :: deep=0
+    INTEGER :: zcoords
+    REAL :: p,z
+    REAL :: rho, thetav
+    INTEGER :: ij,l
+    INTEGER :: moist
     
-    etas=ap(1)/preff+bp(1)
-    etavs=(etas-eta0)*Pi/2
-    phis_ave=T0*g/Gamma*(1-etas**(Rd*Gamma/g))
-    
+    moist=0
+    IF (physics_thermo==thermo_moist .OR. physics_thermo==thermo_fake_moist) moist=1
+   
     DO ij=1,ngrid
-       sinlat=SIN(lat(ij))
-       coslat=COS(lat(ij))
-       phis(ij)=phis_ave+u0*cos(etavs)**1.5*( (-2*sinlat**6 * (coslat**2+1./3) + 10./63 )*u0*cos(etavs)**1.5  &
-            +(8./5*coslat**3 * (sinlat**2 + 2./3) - Pi/4)*radius*Omega )
-       ps(ij)=ps0
+       z=0.
+       zcoords=1
+       CALL baroclinic_wave_test(deep,moist,perturbation,scale_factor,lon(ij),lat(ij),p,z,zcoords,ulon(ij,1),ulat(ij,1), &
+                                 temp(ij,1),thetav,phis(ij),ps(ij),rho,q(ij,1,1))
+       
+       zcoords=0
+       DO l=ll_begin,ll_end
+         p=0.5*(ap(l)+ap(l+1) + (bp(l)+bp(l+1)) * ps(ij))
+         CALL baroclinic_wave_test(deep,moist,perturbation,scale_factor,lon(ij),lat(ij),p,z,zcoords,ulon(ij,l),ulat(ij,l), &
+                                   temp(ij,l),thetav,phis(ij),ps(ij),rho,q(ij,l,1))
+                                   
+        IF (physics_thermo==thermo_fake_moist) temp(ij,l)=Temp(ij,l)*(1+0.608*q(ij,l,1))
+        q(ij,l,2)=0.
+        q(ij,l,3)=0.
+        CALL initial_value_Terminator(lat(ij),lon(ij),q(ij,l,4),q(ij,l,5))
+       END DO
     ENDDO
-    
-    DO l=ll_begin,ll_end
-       etal = 0.5 *( ap(l)/preff+bp(l) + ap(l+1)/preff+bp(l+1) )
-       etavl=(etal-eta0)*Pi/2
-       
-       Tave=T0*etal**(Rd*Gamma/g)
-       dthetaodeta_ave = T0 *( Rd*Gamma/g - kappa)* etal**(Rd*Gamma/g-kappa-1)
-       IF (etat>etal) THEN
-          Tave=Tave+DeltaT*(etat-etal)**5
-          dthetaodeta_ave = dthetaodeta_ave - DeltaT * ( 5*(etat-etal)**4 * etal**(-kappa)  &
-               + kappa * (etat-etal)**5 * etal**(-kappa-1))
-       END IF
-       
-       DO ij=1,ngrid
-          sinlat=SIN(lat(ij))
-          coslat=COS(lat(ij))
-          
-          K=sin(latc)*sinlat+cos(latc)*coslat*cos(lon(ij)-lonc)
-          r=radius*acos(K)
-          utot=u0*cos(etavl)**1.5*sin(2*lat(ij))**2 + up0*exp(-(r/(0.1*radius))**2)
-          ulon(ij,l) = utot
-          ulat(ij,l) = 0.
-          Y = ((-2*sinlat**6*(coslat**2+1./3)+10./63)*2*u0*cos(etavl)**1.5     &
-               + (8./5*coslat**3*(sinlat**2+2./3)-Pi/4)*radius*Omega)
-          T = Tave + 0.75*(etal*Pi*u0/Rd)*sin(etavl)*cos(etavl)**0.5 * Y
-          temp(ij,l)=T
-          
-          q(ij,l,1)=q0*exp(-(lat(ij)/latw)**4)*exp(-((etal-1)*preff/pw)**2)
-          q(ij,l,2)=0.
-          q(ij,l,3)=0.
-          CALL initial_value_Terminator(lat(ij),lon(ij),q(ij,l,4),q(ij,l,5))
-       END DO
-    END DO
-    
+
   END SUBROUTINE compute_etat0
-  
+
 END MODULE etat0_dcmip2016_baroclinic_wave_mod