Changeset 199 for codes/icosagcm/trunk

Timestamp:

07/07/14 15:16:06 (10 years ago)

Author:

dubos

Message:

Isothermal state of rest - untested

File:

• codes/icosagcm/trunk/src/etat0.f90 (modified) (3 diffs)

codes/icosagcm/trunk/src/etat0.f90

@@ -1 +1 @@
 MODULE etat0_mod
+  USE icosa
+  PRIVATE
+
     CHARACTER(len=255),SAVE :: etat0_type
 !$OMP THREADPRIVATE(etat0_type)
+
+    REAL(rstd) :: etat0_temp
+
+    PUBLIC :: etat0, etat0_type
 
 CONTAINS
@@ -45 +52 @@
     
     SELECT CASE (TRIM(etat0_type))
+    CASE ('isothermal')
+       CALL getin_etat0_isothermal
+       CALL etat0_collocated(init_mass, 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)
-       init_mass=.FALSE.
     CASE ('jablonowsky06')
        CALL etat0_jablonowsky06(f_ps,f_phis,f_theta_rhodz,f_u, f_q)
@@ -95 +105 @@
 
   END SUBROUTINE etat0
-         
+
+  SUBROUTINE etat0_collocated(init_mass, f_ps,f_mass,f_phis,f_theta_rhodz,f_u, f_q)
+    IMPLICIT NONE
+    LOGICAL, INTENT(IN) :: init_mass
+    TYPE(t_field),POINTER :: f_ps(:)
+    TYPE(t_field),POINTER :: f_mass(:)
+    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 :: mass(:,:)
+    REAL(rstd),POINTER :: phis(:)
+    REAL(rstd),POINTER :: theta_rhodz(:,:)
+    REAL(rstd),POINTER :: u(:,:)
+    REAL(rstd),POINTER :: q(:,:,:)
+    INTEGER :: ind
+
+    DO ind=1,ndomain
+      IF (.NOT. assigned_domain(ind)) CYCLE
+      CALL swap_dimensions(ind)
+      CALL swap_geometry(ind)
+      ps=f_ps(ind)
+      mass=f_mass(ind)
+      phis=f_phis(ind)
+      theta_rhodz=f_theta_rhodz(ind)
+      u=f_u(ind)
+      q=f_q(ind)
+      CALL compute_etat0_collocated(init_mass, 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
+    USE disvert_mod
+    USE theta2theta_rhodz_mod
+    USE wind_mod
+    IMPLICIT NONE
+    LOGICAL, INTENT(IN) :: init_mass
+    REAL(rstd),INTENT(INOUT) :: ps(iim*jjm)
+    REAL(rstd),INTENT(INOUT) :: mass(iim*jjm,llm)
+    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)
+
+    REAL(rstd) :: lon_i(iim*jjm)
+    REAL(rstd) :: lat_i(iim*jjm)
+    REAL(rstd) :: temp_i(iim*jjm,llm)
+    REAL(rstd) :: ulon_i(iim*jjm,llm)
+    REAL(rstd) :: ulat_i(iim*jjm,llm)
+
+    REAL(rstd) :: lon_e(3*iim*jjm)
+    REAL(rstd) :: lat_e(3*iim*jjm)
+    REAL(rstd) :: ps_e(3*iim*jjm)
+    REAL(rstd) :: mass_e(3*iim*jjm,llm)
+    REAL(rstd) :: phis_e(3*iim*jjm)
+    REAL(rstd) :: temp_e(3*iim*jjm,llm)
+    REAL(rstd) :: ulon_e(3*iim*jjm,llm)
+    REAL(rstd) :: ulat_e(3*iim*jjm,llm)
+    REAL(rstd) :: q_e(3*iim*jjm,llm,nqtot)
+
+    INTEGER :: l,i,j,ij
+
+    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
+             CALL xyz2lonlat(xyz_i(ij,:)/radius,lon_i(ij),lat_i(ij))
+             CALL xyz2lonlat(xyz_e(ij+u_right,:)/radius,lon_e(ij+u_right),lat_e(ij+u_right))
+             CALL xyz2lonlat(xyz_e(ij+u_lup,:)/radius,lon_e(ij+u_lup),lat_e(ij+u_lup))
+             CALL xyz2lonlat(xyz_e(ij+u_ldown,:)/radius,lon_e(ij+u_ldown),lat_e(ij+u_ldown))
+          END DO
+       END DO
+    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)
+    END SELECT
+
+  END SUBROUTINE compute_etat0_ngrid
+
+!----------------------------- Resting isothermal state --------------------------------
+
+  SUBROUTINE getin_etat0_isothermal
+    etat0_temp=300;
+    CALL getin("etat0_temp",etat0_temp)
+  END SUBROUTINE getin_etat0_isothermal
+
+  SUBROUTINE compute_etat0_isothermal(ngrid, phis, ps, temp, ulon, ulat, q)
+    IMPLICIT NONE  
+    INTEGER, INTENT(IN)    :: 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)
+    phis(:)=0
+    ps(:)=preff
+    temp(:,:)=etat0_temp
+    ulon(:,:)=0
+    ulat(:,:)=0
+    q(:,:,:)=0
+  END SUBROUTINE compute_etat0_isothermal
+
 END MODULE etat0_mod
-