Changeset 67

Timestamp:

08/01/12 06:58:27 (12 years ago)

Author:

ymipsl

Message:

bug fix for dcmip4.1 testcase

YM

Location:

codes/icosagcm/trunk/src

Files:

• etat0.f90 (modified) (1 diff)
• etat0_dcmip41.f90 (modified) (4 diffs)

codes/icosagcm/trunk/src/etat0.f90

@@ -37 +37 @@
    CASE DEFAULT
        PRINT*, 'Bad selector for variable etat0 <',etat0_type, &
-            '> options are <jablonowsky06>, <academic>, <dcmip[1-3]> '
+            '> options are <jablonowsky06>, <academic>, <dcmip[1-4]> '
        STOP
     END SELECT

codes/icosagcm/trunk/src/etat0_dcmip41.f90

@@ -82 +82 @@
   
     DO l=1,llm
-      eta(l)= 0.5 *( ap(l)/ps0+bp(l) + ap(l+1)/ps0+bp(l+1) )
+      eta(l)= 0.5 *( ap(l)/preff+bp(l) + ap(l+1)/preff+bp(l+1) )
       etav(l)=(eta(l)-eta0)*Pi/2
     ENDDO
-    etas=ap(1)+bp(1)
+    etas=ap(1)/preff+bp(1)
     etavs=(etas-eta0)*Pi/2
 
@@ -154 +154 @@
          phis(ij)=phis_ave+u0*cos(etavs)**1.5*( (-2*sin(lat)**6 * (cos(lat)**2+1./3) + 10./63 )*u0*cos(etavs)**1.5  &
                                            +(8./5*cos(lat)**3 * (sin(lat)**2 + 2./3) - Pi/4)*radius*Omega )
+!         phis(ij)=phis_ave+u0*cos(etavs)**1.5
+
        ENDDO
      ENDDO
@@ -170 +172 @@
            ij=(j-1)*iim+i
            CALL xyz2lonlat(xyz_i(ij,:)/radius,lon,lat)
-           dthetaodeta=dthetaodeta_ave + 0.25 * Pi * u0/Rd*(1-kappa)*eta(l)**(-kappa)*sin(etav(l))*cos(etav(l))**0.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))   &
+           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                                   &
@@ -182 +185 @@
            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*(r/(0.1*radius))**2 * acos(K) * (sin(latc)*cos(lat) &
+           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)))
-           duodeta= -u0*sin(2*lat)*3*Pi/4*cos(etav(l))**0.5*sin(etav(l))
-           
-           q(ij,l,2)=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