Changeset 78

Timestamp:

08/02/12 23:14:56 (12 years ago)

Author:

ymipsl

Message:

update dcmip4 testcase

YM

Location:

codes/icosagcm/trunk/src

Files:

• etat0.f90 (modified) (2 diffs)
• etat0_dcmip41.f90 (modified) (6 diffs)

codes/icosagcm/trunk/src/etat0.f90

@@ -10 +10 @@
     USE etat0_dcmip2_mod, ONLY : etat0_dcmip2=>etat0
     USE etat0_dcmip3_mod, ONLY : etat0_dcmip3=>etat0  
-    USE etat0_dcmip41_mod, ONLY : etat0_dcmip41=>etat0  
+    USE etat0_dcmip4_mod, ONLY : etat0_dcmip4=>etat0  
     IMPLICIT NONE
     TYPE(t_field),POINTER :: f_ps(:)
@@ -33 +33 @@
     CASE ('dcmip3')
        CALL etat0_dcmip3(f_ps,f_phis,f_theta_rhodz,f_u, f_q)
-     CASE ('dcmip41')
-       CALL etat0_dcmip41(f_ps,f_phis,f_theta_rhodz,f_u, f_q)
+     CASE ('dcmip4')
+       CALL etat0_dcmip4(f_ps,f_phis,f_theta_rhodz,f_u, f_q)
    CASE DEFAULT
        PRINT*, 'Bad selector for variable etat0 <',etat0_type, &

codes/icosagcm/trunk/src/etat0_dcmip41.f90

@@ -1 @@
-MODULE etat0_dcmip41_mod
+MODULE etat0_dcmip4_mod
   USE icosa
   PRIVATE
@@ -11 +11 @@
   REAL(rstd),PARAMETER :: Gamma=0.005
   REAL(rstd),PARAMETER :: up0=1
+  REAL(rstd) :: latw=2*Pi/9
+  REAL(rstd) :: pw=34000
+  REAL(rstd) :: q0=0.021
   REAL(rstd) :: lonc
   REAL(rstd) :: latc
@@ -42 +45 @@
       u=f_u(ind)
       q=f_q(ind)
-      CALL compute_etat0_dcmip41(ps, phis, theta_rhodz, u, q)
+      CALL compute_etat0_dcmip4(ps, phis, theta_rhodz, u, q)
     ENDDO
 
   END SUBROUTINE etat0
   
-  SUBROUTINE compute_etat0_dcmip41(ps, phis, theta_rhodz, u, q)
+  SUBROUTINE compute_etat0_dcmip4(ps, phis, theta_rhodz, u, q)
   USE icosa
   USE disvert_mod
@@ -78 +81 @@
   REAL(rstd) :: lonx,latx
   REAL(rstd) :: dthetaodeta_ave, dthetaodeta, dthetaodlat, duodeta, K, r
-  lonc=Pi/9
-  latc=2*Pi/9 
+  INTEGER :: testcase
+  
+    lonc=Pi/9
+    latc=2*Pi/9 
+  
+    testcase=1
+    CALL getin("dcmip4_testcase",testcase)
+  
   
     DO l=1,llm
@@ -131 +140 @@
        Tave=T0*eta(l)**(Rd*Gamma/g)
        IF (etat>eta(l)) Tave=Tave+DeltaT*(etat-eta(l))**5
-       DO j=jj_begin,jj_end
-         DO i=ii_begin,ii_end
+       DO j=jj_begin-1,jj_end+1
+         DO i=ii_begin-1,ii_end+1
            ij=(j-1)*iim+i
            CALL xyz2lonlat(xyz_i(ij,:)/radius,lon,lat)
@@ -161 +170 @@
     CALL compute_theta2theta_rhodz(ps,theta,theta_rhodz,0)
 
-    q(:,:,1)=theta(:,:)
-    
-
-    DO l=1,llm
-       dthetaodeta_ave = T0 *( Rd*Gamma/g - kappa)* eta(l)**(Rd*Gamma/g-kappa-1)
-       IF (etat>eta(l)) dthetaodeta_ave = dthetaodeta_ave - DeltaT * ( 5*(etat-eta(l))**4 * eta(l)**(-kappa)  &
+
+    IF (testcase==1) THEN
+    
+      q(:,:,1)=theta(:,:)
+
+      DO l=1,llm
+         dthetaodeta_ave = T0 *( Rd*Gamma/g - kappa)* eta(l)**(Rd*Gamma/g-kappa-1)
+         IF (etat>eta(l)) dthetaodeta_ave = dthetaodeta_ave - DeltaT * ( 5*(etat-eta(l))**4 * eta(l)**(-kappa)  &
                                                           + kappa * (etat-eta(l))**5 * eta(l)**(-kappa-1))
-       DO j=jj_begin,jj_end
-         DO i=ii_begin,ii_end
-           ij=(j-1)*iim+i
-           CALL xyz2lonlat(xyz_i(ij,:)/radius,lon,lat)
-           dthetaodeta=dthetaodeta_ave + 3./4. * Pi * u0/Rd*(1-kappa)*eta(l)**(-kappa)*sin(etav(l))*cos(etav(l))**0.5 * Y(ij,l)  & 
-                                       + 3/8. * Pi**2*u0/Rd * eta(l)**(1-kappa) * cos(etav(l))**1.5 * Y(ij,l)                    & 
-                                  - 3./16. * Pi**2 * u0 /Rd * eta(l)**(1-kappa) * sin(etav(l))**2 * cos(etav(l))**(-0.5) *Y(ij,l)&
-                                  - 9./8.  * Pi**2 * u0 /Rd * eta(l)**(1-kappa) * sin(etav(l))**2 * cos(etav(l))   &
+         DO j=jj_begin,jj_end
+           DO i=ii_begin,ii_end
+             ij=(j-1)*iim+i
+             CALL xyz2lonlat(xyz_i(ij,:)/radius,lon,lat)
+             dthetaodeta=dthetaodeta_ave + 3./4. * Pi * u0/Rd*(1-kappa)*eta(l)**(-kappa)*sin(etav(l))*cos(etav(l))**0.5 * Y(ij,l) & 
+                                        + 3/8. * Pi**2*u0/Rd * eta(l)**(1-kappa) * cos(etav(l))**1.5 * Y(ij,l)                    & 
+                                   - 3./16. * Pi**2 * u0 /Rd * eta(l)**(1-kappa) * sin(etav(l))**2 * cos(etav(l))**(-0.5) *Y(ij,l)&
+                                   - 9./8.  * Pi**2 * u0 /Rd * eta(l)**(1-kappa) * sin(etav(l))**2 * cos(etav(l))   &
                                                                                   * (-2*sin(lat)**6*(cos(lat)**2+1./3.)+10./63.) 
-           dthetaodlat=3./4.*Pi*u0/Rd*eta(l)**(1-kappa)*sin(etav(l))*cos(etav(l))**0.5                                   &
-                       *( 2*u0*cos(etav(l))**1.5 * ( -12 * cos(lat)*sin(lat)**5*(cos(lat)**2+1./3.)+4*cos(lat)*sin(lat)**7) &
-                      + radius*omega*(-24./5. * sin(lat) * cos(lat)**2 * (sin(lat)**2 + 2./3.) + 16./5. * cos(lat)**4 * sin(lat)))
+             dthetaodlat=3./4.*Pi*u0/Rd*eta(l)**(1-kappa)*sin(etav(l))*cos(etav(l))**0.5                                   &
+                         *( 2*u0*cos(etav(l))**1.5 * ( -12 * cos(lat)*sin(lat)**5*(cos(lat)**2+1./3.)+4*cos(lat)*sin(lat)**7) &
+                        + radius*omega*(-24./5. * sin(lat) * cos(lat)**2 * (sin(lat)**2 + 2./3.) + 16./5. * cos(lat)**4 * sin(lat)))
                           
-           duodeta=-u0 * sin(2*lat)**2 * 3./4.*Pi * cos(etav(l))**0.5 * sin(etav(l))
+             duodeta=-u0 * sin(2*lat)**2 * 3./4.*Pi * cos(etav(l))**0.5 * sin(etav(l))
            
-           K=sin(latc)*sin(lat)+cos(latc)*cos(lat)*cos(lon-lonc)
-           r=radius*acos(K)
-           vort = -4*u0/radius*cos(etav(l))**1.5 * sin(lat) * cos(lat) * (2.-5.*sin(lat)**2)                  & 
-                  + up0/radius*exp(-(r/(0.1*radius))**2) * (tan(lat)-2*(radius/(0.1*radius))**2 * acos(K) * (sin(latc)*cos(lat) &
+             K=sin(latc)*sin(lat)+cos(latc)*cos(lat)*cos(lon-lonc)
+             r=radius*acos(K)
+             vort = -4*u0/radius*cos(etav(l))**1.5 * sin(lat) * cos(lat) * (2.-5.*sin(lat)**2)                  & 
+                    + up0/radius*exp(-(r/(0.1*radius))**2) * (tan(lat)-2*(radius/(0.1*radius))**2 * acos(K) * (sin(latc)*cos(lat) &
                                                             -cos(latc)*sin(lat)*cos(lon-lonc))/(sqrt(1-K**2)))
-           q(ij,l,2)=ABS(g/preff*(-1./radius*duodeta*dthetaodlat-(2*sin(lat)*omega+vort)*dthetaodeta))
+             q(ij,l,2)=ABS(g/preff*(-1./radius*duodeta*dthetaodlat-(2*sin(lat)*omega+vort)*dthetaodeta))
+           ENDDO
          ENDDO
-       ENDDO
-     ENDDO       
+       ENDDO       
+
+    ELSE IF (testcase==2) THEN
+    
+      DO l=1,llm
+         DO j=jj_begin,jj_end
+           DO i=ii_begin,ii_end
+             ij=(j-1)*iim+i
+             CALL xyz2lonlat(xyz_i(ij,:)/radius,lon,lat)
+             q(ij,l,1)=q0*exp(-(lat/latw)**4)*exp(-((eta(l)-1)*preff/pw)**2)
+           ENDDO
+         ENDDO
+       ENDDO
+    
+    ENDIF       
 
         
-  END SUBROUTINE compute_etat0_dcmip41
-  
-END MODULE etat0_dcmip41_mod
+  END SUBROUTINE compute_etat0_dcmip4
+  
+END MODULE etat0_dcmip4_mod