Changeset 201 for codes/icosagcm/trunk

Timestamp:

07/08/14 15:03:32 (10 years ago)

Author:

dubos

Message:

Simplified etat0 : isothermal, jablonowsky06

Location:

codes/icosagcm/trunk/src

Files:

• etat0.f90 (modified) (5 diffs)
• etat0_jablonowsky06.f90 (modified) (3 diffs)
• vector.f90 (modified) (1 diff)

codes/icosagcm/trunk/src/etat0.f90

@@ -13 +13 @@
   
   SUBROUTINE etat0(f_ps,f_mass,f_phis,f_theta_rhodz,f_u, f_q)
-    USE icosa
     USE mpipara, ONLY : is_mpi_root
     USE disvert_mod
@@ -54 +53 @@
     CASE ('isothermal')
        CALL getin_etat0_isothermal
-       CALL etat0_collocated(init_mass, f_phis,f_ps,f_mass,f_theta_rhodz,f_u, f_q)
+       CALL etat0_collocated(f_phis,f_ps,f_mass,f_theta_rhodz,f_u, f_q)
     CASE ('williamson91.6')
        init_mass=.FALSE.
        CALL etat0_williamson_new(f_phis,f_mass,f_theta_rhodz,f_u, f_q)
     CASE ('jablonowsky06')
-       CALL etat0_jablonowsky06(f_ps,f_phis,f_theta_rhodz,f_u, f_q)
+!       CALL etat0_jablonowsky06(f_ps,f_phis,f_theta_rhodz,f_u, f_q)
+       CALL etat0_collocated(f_phis,f_ps,f_mass,f_theta_rhodz,f_u, f_q)
     CASE ('academic')
        CALL etat0_academic(f_ps,f_phis,f_theta_rhodz,f_u, f_q)
@@ -106 +106 @@
   END SUBROUTINE etat0
 
-  SUBROUTINE etat0_collocated(init_mass, f_ps,f_mass,f_phis,f_theta_rhodz,f_u, f_q)
+  SUBROUTINE etat0_collocated(f_phis,f_ps,f_mass,f_theta_rhodz,f_u, f_q)
+    USE mpipara
     IMPLICIT NONE
-    LOGICAL, INTENT(IN) :: init_mass
     TYPE(t_field),POINTER :: f_ps(:)
     TYPE(t_field),POINTER :: f_mass(:)
@@ -134 +134 @@
       u=f_u(ind)
       q=f_q(ind)
-      CALL compute_etat0_collocated(init_mass, ps,mass, phis, theta_rhodz, u, q)
+      CALL compute_etat0_collocated(ps,mass, phis, theta_rhodz, u, q)
+
     ENDDO
   END SUBROUTINE etat0_collocated
 
-  SUBROUTINE compute_etat0_collocated(init_mass, ps,mass, phis, theta_rhodz, u, q)
-    USE icosa
+  SUBROUTINE compute_etat0_collocated(ps,mass, phis, theta_rhodz, u, q)
     USE disvert_mod
     USE theta2theta_rhodz_mod
     USE wind_mod
+    USE etat0_jablonowsky06_mod, ONLY : compute_jablonowsky06 => compute_etat0_new
     IMPLICIT NONE
-    LOGICAL, INTENT(IN) :: init_mass
     REAL(rstd),INTENT(INOUT) :: ps(iim*jjm)
     REAL(rstd),INTENT(INOUT) :: mass(iim*jjm,llm)
@@ -182 +182 @@
     END DO
 
-    CALL compute_etat0_ngrid(iim*jjm,lon_i,lat_i, phis, ps, mass, temp_i, ulon_i, ulat_i, q)
-    IF(init_mass) THEN
-       CALL compute_rhodz(.TRUE., ps, mass)
-    END IF
-    CALL compute_temperature2theta_rhodz(ps,temp_i,theta_rhodz,0) ! FIXME - works with shallow-water ?
-
-    CALL compute_etat0_ngrid(3*iim*jjm,lon_e,lat_e, phis_e,ps_e,mass_e, temp_e, ulon_e, ulat_e, q_e)
-    CALL compute_wind_from_lonlat_compound(ulon_e, ulat_e, u)
-
-  END SUBROUTINE compute_etat0_collocated
-
-  SUBROUTINE compute_etat0_ngrid(ngrid,lon,lat, phis, ps, mass, temp, ulon, ulat, q)
-    IMPLICIT NONE  
-    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(INOUT) :: ps(ngrid)
-    REAL(rstd),INTENT(INOUT) :: mass(ngrid,llm)
-    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)
-
     SELECT CASE (TRIM(etat0_type))
     CASE ('isothermal')
-       CALL compute_etat0_isothermal(ngrid, phis, ps, temp, ulon, ulat, q)
+       CALL compute_etat0_isothermal(iim*jjm, phis, ps, temp_i, ulon_i, ulat_i, q)
+       CALL compute_etat0_isothermal(3*iim*jjm, phis_e, ps_e, temp_e, ulon_e, ulat_e, q_e)
+    CASE('jablonowsky06')
+       CALL compute_jablonowsky06(iim*jjm,lon_i,lat_i, phis, ps, temp_i, ulon_i, ulat_i)
+       CALL compute_jablonowsky06(3*iim*jjm,lon_e,lat_e, phis_e, ps_e, temp_e, ulon_e, ulat_e)
     END SELECT
 
-  END SUBROUTINE compute_etat0_ngrid
+    CALL compute_temperature2theta_rhodz(ps,temp_i,theta_rhodz,0)       
+    CALL compute_wind_perp_from_lonlat_compound(ulon_e, ulat_e, u)
+
+  END SUBROUTINE compute_etat0_collocated
 
 !----------------------------- Resting isothermal state --------------------------------
 
   SUBROUTINE getin_etat0_isothermal
-    etat0_temp=300;
-    CALL getin("etat0_temp",etat0_temp)
+    etat0_temp=300
+    CALL getin("etat0_isothermal_temp",etat0_temp)
   END SUBROUTINE getin_etat0_isothermal
 

codes/icosagcm/trunk/src/etat0_jablonowsky06.f90

@@ -6 +6 @@
   REAL(rstd),PARAMETER :: ps0=1e5
   REAL(rstd),PARAMETER :: u0=35
-!  REAL(rstd),PARAMETER :: u0=0
   REAL(rstd),PARAMETER :: T0=288
   REAL(rstd),PARAMETER :: DeltaT=4.8e5
@@ -12 +11 @@
   REAL(rstd),PARAMETER :: Gamma=0.005
   REAL(rstd),PARAMETER :: up0=1
-  PUBLIC  test_etat0_jablonowsky06, etat0, compute_etat0_jablonowsky06
+  PUBLIC  test_etat0_jablonowsky06, etat0, compute_etat0_jablonowsky06, compute_etat0_new
 CONTAINS
   
@@ -248 +247 @@
 
   END SUBROUTINE compute_etat0_jablonowsky06
+
+  SUBROUTINE compute_etat0_new(ngrid,lat,lon, phis, ps, temp, ulon, ulat)
+    USE disvert_mod
+    IMPLICIT NONE  
+    INTEGER, INTENT(IN) :: ngrid
+    REAL(rstd),INTENT(IN) :: lat(ngrid)
+    REAL(rstd),INTENT(IN) :: lon(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)
+    
+    INTEGER :: l,ij
+    REAL(rstd) :: eta(llm)
+    REAL(rstd) :: etav(llm)
+    REAL(rstd) :: etas, etavs, Tave, phis_ave, r2
+    
+    DO l=1,llm
+       eta(l)  = 0.5 *( ap(l)/ps0+bp(l) + ap(l+1)/ps0+bp(l+1) )
+       etav(l) = (eta(l)-eta0)*Pi/2
+    ENDDO
+    etas  = ap(1)+bp(1)
+    etavs = (etas-eta0)*Pi/2
+        
+    phis_ave=T0*g/Gamma*(1-etas**(Rd*Gamma/g))
+    DO ij=1,ngrid
+       ps(ij)=ps0
+       phis(ij) = phis_ave + u0*cos(etavs)**1.5*( (-2*sin(lat(ij))**6 * (cos(lat(ij))**2+1./3) + 10./63 )*u0*cos(etavs)**1.5  &
+            +(8./5*cos(lat(ij))**3 * (sin(lat(ij))**2 + 2./3) - Pi/4)*radius*Omega )
+    ENDDO
+
+    DO l=1,llm
+       Tave=T0*eta(l)**(Rd*Gamma/g)
+       IF (etat>eta(l)) Tave=Tave+DeltaT*(etat-eta(l))**5
+       DO ij=1,ngrid
+          r2 = arc(Pi/9,2*Pi/9, lon(ij),lat(ij))**2
+          temp(ij,l) = Tave + 0.75*(eta(l)*Pi*u0/Rd)*sin(etav(l))*cos(etav(l))**0.5  &
+               * ( (-2*sin(lat(ij))**6*(cos(lat(ij))**2+1./3)+10./63)*2*u0*cos(etav(l))**1.5 &
+               + (8./5*cos(lat(ij))**3*(sin(lat(ij))**2+2./3)-Pi/4)*radius*Omega)
+          ulon(ij,l) = u0*cos(etav(l))**1.5*sin(2*lat(ij))**2 + up0*exp(-r2/0.01)
+          ulat(ij,l) = 0
+       ENDDO
+    ENDDO
+    
+    
+  END SUBROUTINE compute_etat0_new
   
 END MODULE etat0_jablonowsky06_mod

codes/icosagcm/trunk/src/vector.f90

@@ -47 +47 @@
      
     END SUBROUTINE cross_product2
-    
+
+    FUNCTION arc(lon1,lat1, lon2,lat2)
+      REAL(rstd) :: lon1,lat1, lon2,lat2, arc2, &
+           x1,x2,x3,y1,y2,y3,z1,z2,z3
+      x1=cos(lat1)*cos(lon1)
+      x2=cos(lat2)*cos(lon2)
+      y1=cos(lat1)*sin(lon1)
+      y2=cos(lat2)*sin(lon2)
+      z1=sin(lat1)
+      z2=sin(lat2)
+      x3=y1*z2-y2*z1
+      y3=z1*x2-z2*x1
+      z3=x1*y2-x2*y1
+      arc=ASIN(SQRT(x3*x3+y3*y3+z3*z3))
+    END FUNCTION arc
+
 END MODULE vector