Index: codes/icosagcm/trunk/src/etat0_dcmip4.f90
===================================================================
--- codes/icosagcm/trunk/src/etat0_dcmip4.f90	(revision 286)
+++ codes/icosagcm/trunk/src/etat0_dcmip4.f90	(revision 346)
@@ -1,5 +1,7 @@
 MODULE etat0_dcmip4_mod
   USE icosa
+  IMPLICIT NONE
   PRIVATE
+  
   REAL(rstd),PARAMETER :: eta0=0.252
   REAL(rstd),PARAMETER :: etat=0.2
@@ -11,220 +13,109 @@
   REAL(rstd),PARAMETER :: Gamma=0.005
   REAL(rstd),PARAMETER :: up0=1
-  REAL(rstd) :: latw=2*Pi/9
-!$OMP THREADPRIVATE(latw)
-  REAL(rstd) :: pw=34000
-!$OMP THREADPRIVATE(pw)
-  REAL(rstd) :: q0=0.021
-!$OMP THREADPRIVATE(q0)
-  REAL(rstd) :: lonc
-!$OMP THREADPRIVATE(lonc)
-  REAL(rstd) :: latc
-!$OMP THREADPRIVATE(latc)
+  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)  
+  
+  PUBLIC getin_etat0, compute_etat0
 
-  INTEGER,SAVE :: testcase
-!$OMP THREADPRIVATE(testcase)  
+CONTAINS
 
-  PUBLIC  etat0
-CONTAINS
-  
-   
-    
-  SUBROUTINE etat0(f_ps,f_phis,f_theta_rhodz,f_u, f_q)
-  USE icosa
-  IMPLICIT NONE
-    TYPE(t_field),POINTER :: f_ps(:)
-    TYPE(t_field),POINTER :: f_phis(:)
-    TYPE(t_field),POINTER :: f_theta_rhodz(:)
-    TYPE(t_field),POINTER :: f_u(:)
-    TYPE(t_field),POINTER :: f_q(:)
-  
-    REAL(rstd),POINTER :: ps(:)
-    REAL(rstd),POINTER :: phis(:)
-    REAL(rstd),POINTER :: theta_rhodz(:,:)
-    REAL(rstd),POINTER :: u(:,:)
-    REAL(rstd),POINTER :: q(:,:,:)
-    INTEGER :: ind
-
+  SUBROUTINE getin_etat0
+    USE mpipara, ONLY : is_mpi_root
+    IF(nqtot<2) THEN
+       IF (is_mpi_root)  THEN
+          PRINT *, "nqtot must be at least 2 for test case DCMIP4"
+       END IF
+       STOP
+    END IF
     testcase=1
     CALL getin("dcmip4_testcase",testcase)
+  END SUBROUTINE getin_etat0
+
+  SUBROUTINE compute_etat0(ngrid,lon,lat, phis,ps,temp,ulon,ulat,q)
+    USE icosa
+    USE disvert_mod
+    INTEGER, INTENT(IN) :: ngrid
+    REAL(rstd),INTENT(IN) :: lon(ngrid)
+    REAL(rstd),INTENT(IN) :: lat(ngrid)
+    REAL(rstd),INTENT(OUT) :: phis(ngrid)
+    REAL(rstd),INTENT(OUT) :: ps(ngrid)
+    REAL(rstd),INTENT(OUT) :: temp(ngrid,llm)
+    REAL(rstd),INTENT(OUT) :: ulon(ngrid,llm)
+    REAL(rstd),INTENT(OUT) :: ulat(ngrid,llm)
+    REAL(rstd),INTENT(OUT) :: q(ngrid,llm,nqtot)
     
-    DO ind=1,ndomain
-      IF (.NOT. assigned_domain(ind)) CYCLE
-      CALL swap_dimensions(ind)
-      CALL swap_geometry(ind)
-      ps=f_ps(ind)
-      phis=f_phis(ind)
-      theta_rhodz=f_theta_rhodz(ind)
-      u=f_u(ind)
-      q=f_q(ind)
-      CALL compute_etat0_dcmip4(ps, phis, theta_rhodz, u, q)
-    ENDDO
-
-  END SUBROUTINE etat0
-  
-  SUBROUTINE compute_etat0_dcmip4(ps, phis, theta_rhodz, u, q)
-  USE icosa
-  USE disvert_mod
-  USE pression_mod
-  USE exner_mod
-  USE geopotential_mod
-  USE theta2theta_rhodz_mod
-  IMPLICIT NONE  
-  REAL(rstd),INTENT(OUT) :: ps(iim*jjm)
-  REAL(rstd),INTENT(OUT) :: phis(iim*jjm)
-  REAL(rstd),INTENT(OUT) :: theta_rhodz(iim*jjm,llm)
-  REAL(rstd),INTENT(OUT) :: u(3*iim*jjm,llm)
-  REAL(rstd),INTENT(OUT) :: q(iim*jjm,llm,nqtot)
-  
-  INTEGER :: i,j,l,ij
-  REAL(rstd) :: theta(iim*jjm,llm)
-  REAL(rstd) :: Y(iim*jjm,llm)
-  REAL(rstd) :: vort
-  REAL(rstd) :: eta(llm)
-  REAL(rstd) :: etav(llm)
-  REAL(rstd) :: etas, etavs
-  REAL(rstd) :: lon,lat
-  REAL(rstd) :: ulon(3)
-  REAL(rstd) :: ep(3), norm_ep
-  REAL(rstd) :: Tave, T
-  REAL(rstd) :: phis_ave
-  REAL(rstd) :: V0(3)
-  REAL(rstd) :: r2
-  REAL(rstd) :: utot
-  REAL(rstd) :: lonx,latx
-  REAL(rstd) :: dthetaodeta_ave, dthetaodeta, dthetaodlat, duodeta, K, r
-  
-    lonc=Pi/9
-    latc=2*Pi/9 
-  
-    DO l=1,llm
-      eta(l)= 0.5 *( ap(l)/preff+bp(l) + ap(l+1)/preff+bp(l+1) )
-      etav(l)=(eta(l)-eta0)*Pi/2
-    ENDDO
+    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
+    
     etas=ap(1)/preff+bp(1)
     etavs=(etas-eta0)*Pi/2
-
-    DO j=jj_begin,jj_end
-      DO i=ii_begin,ii_end
-        ij=(j-1)*iim+i
-        ps(ij)=ps0
-      ENDDO
+    phis_ave=T0*g/Gamma*(1-etas**(Rd*Gamma/g))
+    
+    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
     ENDDO
     
+    DO l=1,llm
+       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
+          
+          IF (testcase==1) THEN
+             q(ij,l,1)=T*etal**(-kappa)
+             IF(nqtot>2) q(ij,l,3)=1.
+             dthetaodeta=dthetaodeta_ave + 3./4. * Pi * u0/Rd*(1-kappa)*etal**(-kappa)*sin(etavl)*cos(etavl)**0.5 * Y & 
+                  + 3/8. * Pi**2*u0/Rd * etal**(1-kappa) * cos(etavl)**1.5 * Y                    & 
+                  - 3./16. * Pi**2 * u0 /Rd * etal**(1-kappa) * sin(etavl)**2 * cos(etavl)**(-0.5) *Y &
+                  - 9./8.  * Pi**2 * u0 /Rd * etal**(1-kappa) * sin(etavl)**2 * cos(etavl)   &
+                  * (-2*sinlat**6*(coslat**2+1./3.)+10./63.) 
+             dthetaodlat=3./4.*Pi*u0/Rd*etal**(1-kappa)*sin(etavl)*cos(etavl)**0.5                                   &
+                  *( 2*u0*cos(etavl)**1.5 * ( -12 * coslat*sinlat**5*(coslat**2+1./3.)+4*coslat*sinlat**7) &
+                  + radius*omega*(-24./5. * sinlat * coslat**2 * (sinlat**2 + 2./3.) + 16./5. * coslat**4 * sinlat))
+             
+             duodeta=-u0 * sin(2*lat(ij))**2 * 3./4.*Pi * cos(etavl)**0.5 * sin(etavl)
+             
+             vort = -4*u0/radius*cos(etavl)**1.5 * sinlat * coslat * (2.-5.*sinlat**2)                  & 
+                  + up0/radius*exp(-(r/(0.1*radius))**2) * (tan(lat(ij))-2*(radius/(0.1*radius))**2 * acos(K) * (sin(latc)*coslat &
+                  -cos(latc)*sinlat*cos(lon(ij)-lonc))/(sqrt(1-K**2)))
+             q(ij,l,2)=ABS(g/preff*(-1./radius*duodeta*dthetaodlat-(2*sinlat*omega+vort)*dthetaodeta))
+             IF(nqtot>3) q(ij,l,4)=cos(lon(ij))*coslat
+          ELSE IF (testcase==2) THEN
+             q(ij,l,1)=q0*exp(-(lat(ij)/latw)**4)*exp(-((etal-1)*preff/pw)**2)
+          END IF
+       END DO
+    END DO
     
-    CALL lonlat2xyz(lonc,latc,V0)
-
-    u(:,:)=1e10      
-    DO l=1,llm
-      DO j=jj_begin-1,jj_end+1
-        DO i=ii_begin-1,ii_end+1
-          ij=(j-1)*iim+i
-          
-          lon=lon_e(ij+u_right) ; lat=lat_e(ij+u_right)
-          K=sin(latc)*sin(lat)+cos(latc)*cos(lat)*cos(lon-lonc)
-          r=radius*acos(K)
-          utot=u0*cos(etav(l))**1.5*sin(2*lat)**2 + up0*exp(-(r/(0.1*radius))**2)
-          u(ij+u_right,l) = utot * sum(elon_e(ij+u_right,:) * ep_e(ij+u_right,:))
-
-          lon=lon_e(ij+u_lup) ; lat=lat_e(ij+u_lup)
-          K=sin(latc)*sin(lat)+cos(latc)*cos(lat)*cos(lon-lonc)
-          r=radius*acos(K)
-          utot=u0*cos(etav(l))**1.5*sin(2*lat)**2 + up0*exp(-(r/(0.1*radius))**2)
-          u(ij+u_lup,l) = utot * sum(elon_e(ij+u_lup,:) * ep_e(ij+u_lup,:))
-
-          lon=lon_e(ij+u_ldown) ; lat=lat_e(ij+u_ldown)
-          K=sin(latc)*sin(lat)+cos(latc)*cos(lat)*cos(lon-lonc)
-          r=radius*acos(K)
-          utot=u0*cos(etav(l))**1.5*sin(2*lat)**2 + up0*exp(-(r/(0.1*radius))**2)
-          u(ij+u_ldown,l) = utot * sum(elon_e(ij+u_ldown,:) * ep_e(ij+u_ldown,:))
-
-       ENDDO
-     ENDDO
-    ENDDO 
-     
-     
-     DO l=1,llm
-       Tave=T0*eta(l)**(Rd*Gamma/g)
-       IF (etat>eta(l)) Tave=Tave+DeltaT*(etat-eta(l))**5
-       DO j=jj_begin-1,jj_end+1
-         DO i=ii_begin-1,ii_end+1
-           ij=(j-1)*iim+i
-           lat=lat_i(ij)
-            Y(ij,l)=((-2*sin(lat)**6*(cos(lat)**2+1./3)+10./63)*2*u0*cos(etav(l))**1.5     &
-                                + (8./5*cos(lat)**3*(sin(lat)**2+2./3)-Pi/4)*radius*Omega)
-            T=Tave+ 0.75*(eta(l)*Pi*u0/Rd)*sin(etav(l))*cos(etav(l))**0.5 * Y(ij,l)
-            
-            theta(ij,l)=T*eta(l)**(-kappa)
-
-          ENDDO
-       ENDDO
-     ENDDO
-     
-     
-     phis_ave=T0*g/Gamma*(1-etas**(Rd*Gamma/g))
-     DO j=jj_begin,jj_end
-       DO i=ii_begin,ii_end
-         ij=(j-1)*iim+i
-         lat=lat_i(ij)
-         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
-
-    CALL compute_theta2theta_rhodz(ps,theta,theta_rhodz,0)
-
-
-    IF (testcase==1) THEN
-    
-      q(:,:,1)=theta(:,:)
-      IF(nqtot>2) q(:,:,3)=1.
-
-      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
-             lon=lon_i(ij) ; lat=lat_i(ij)
-             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)))
-                          
-             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) &
-                                                            -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))
-             IF(nqtot>3) q(ij,l,4)=cos(lon)*cos(lat)
-           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
-             q(ij,l,1)=q0*exp(-(lat_i(ij)/latw)**4)*exp(-((eta(l)-1)*preff/pw)**2)
-           ENDDO
-         ENDDO
-       ENDDO
-    
-    ENDIF       
-
-        
-  END SUBROUTINE compute_etat0_dcmip4
+  END SUBROUTINE compute_etat0
   
 END MODULE etat0_dcmip4_mod
