Changeset 179

Timestamp:

12/13/13 18:33:56 (11 years ago)

Author:

mtort

Message:

Prepared structure for non-traditional shallow-water equations

Location:

codes/icosagcm/trunk/src

Files:

• caldyn_gcm.f90 (modified) (12 diffs)
• etat0_williamson.f90 (modified) (1 diff)
• wind.f90 (modified) (2 diffs)

codes/icosagcm/trunk/src/caldyn_gcm.f90

@@ -16 +16 @@
   TYPE(t_field),POINTER :: f_geopot(:)
   TYPE(t_field),POINTER :: f_wwuu(:)
+  TYPE(t_field),POINTER :: f_planetvel(:)
    
   INTEGER :: caldyn_conserv
@@ -32 +33 @@
     IMPLICIT NONE
     CHARACTER(len=255) :: def
+    INTEGER            :: ind
+    REAL(rstd),POINTER :: planetvel(:)
   
     def='enstrophy'
@@ -48 +51 @@
 
     CALL allocate_caldyn
-  
+
+    DO ind=1,ndomain
+       CALL swap_dimensions(ind)
+       CALL swap_geometry(ind)
+       planetvel=f_planetvel(ind)
+       CALL compute_planetvel(planetvel)
+    END DO
+
   END SUBROUTINE init_caldyn
- 
+
   SUBROUTINE allocate_caldyn
   USE icosa
@@ -71 +81 @@
     CALL allocate_field(f_geopot,field_t,type_real,llm+1,name='geopot')  ! geopotential
     CALL allocate_field(f_wwuu,  field_u,type_real,llm+1,name='wwuu')
+    CALL allocate_field(f_planetvel,  field_u,type_real, name='planetvel') ! planetary velocity at r=a
 
   END SUBROUTINE allocate_caldyn
@@ -210 +221 @@
           pk = f_pk(ind)
           geopot = f_geopot(ind)  
-          CALL compute_geopot(ps,mass,theta, pk,geopot)
+          CALL compute_geopot(ps,mass,theta,u, pk,geopot)
           hflux=f_hflux(ind)
           convm = f_dmass(ind)
@@ -237 +248 @@
           pk = f_pk(ind)
           geopot = f_geopot(ind)  
-          CALL compute_geopot(ps,mass,theta, pk,geopot)
+          CALL compute_geopot(ps,mass,theta,u, pk,geopot)
           hflux=f_hflux(ind)
           convm = f_dmass(ind)
@@ -280 +291 @@
     
 END SUBROUTINE caldyn
+
+SUBROUTINE compute_planetvel(planetvel)
+  USE wind_mod
+  REAL(rstd),INTENT(OUT)  :: planetvel(iim*3*jjm)
+  REAL(rstd) :: ulon(iim*3*jjm)
+  REAL(rstd) :: ulat(iim*3*jjm)
+  REAL(rstd) :: lon,lat
+  INTEGER :: ij
+  DO ij=ij_begin_ext,ij_end_ext
+     CALL xyz2lonlat(xyz_e(ij+u_right,:),lon,lat)
+     ulon(ij+u_right)=a*omega*cos(lat)
+     ulat(ij+u_right)=0
+
+     CALL xyz2lonlat(xyz_e(ij+u_lup,:),lon,lat)
+     ulon(ij+u_lup)=a*omega*cos(lat)
+     ulat(ij+u_lup)=0
+
+     CALL xyz2lonlat(xyz_e(ij+u_ldown,:),lon,lat)
+     ulon(ij+u_ldown)=a*omega*cos(lat)
+     ulat(ij+u_ldown)=0
+  END DO
+  CALL compute_wind2D_perp_from_lonlat_compound(ulon, ulat, planetvel)
+END SUBROUTINE compute_planetvel
     
 SUBROUTINE compute_pvort(ps,u,theta_rhodz, rhodz,theta,qu,qv)
@@ -334 +368 @@
   DO l = ll_begin,ll_end
 !DIR$ SIMD
-     DO ij=ij_begin_ext,ij_end_ext         
+     DO ij=ij_begin_ext,ij_end_ext
           etav= 1./Av(ij+z_up)*(  ne_rup        * u(ij+u_rup,l)        * de(ij+u_rup)         &
                                 + ne_left * u(ij+t_rup+u_left,l) * de(ij+t_rup+u_left)  &
@@ -370 +404 @@
   END SUBROUTINE compute_pvort
   
-  SUBROUTINE compute_geopot(ps,rhodz,theta,pk,geopot) 
+  SUBROUTINE compute_geopot(ps,rhodz,theta,u, pk,geopot) 
   USE icosa
   USE disvert_mod
@@ -380 +414 @@
     REAL(rstd),INTENT(IN)    :: rhodz(iim*jjm,llm)
     REAL(rstd),INTENT(IN)    :: theta(iim*jjm,llm)    ! potential temperature
+    REAL(rstd),INTENT(IN)    :: u(3*iim*jjm,llm)      ! prognostic velocity
     REAL(rstd),INTENT(INOUT) :: pk(iim*jjm,llm)       ! Exner function
     REAL(rstd),INTENT(INOUT) :: geopot(iim*jjm,llm+1) ! geopotential
     
+    REAL(rstd) :: urel(3*iim*jjm,llm)  ! relative velocity
+    REAL(rstd) :: cor_NT(iim*jjm,llm)  ! NT coriolis force u.(du/dPhi)
     INTEGER :: i,j,ij,l
     REAL(rstd) :: p_ik, exner_ik
@@ -474 +511 @@
     REAL(rstd),INTENT(OUT) :: du(iim*3*jjm,llm)
 
+    REAL(rstd) :: urel(3*iim*jjm,llm)  ! relative velocity
     REAL(rstd) :: Ftheta(3*iim*jjm,llm) ! theta flux
     REAL(rstd) :: berni(iim*jjm,llm)  ! Bernoulli function
@@ -629 +667 @@
           DO ij=ij_begin,ij_end         
                 
-                berni(ij,l) = pk(ij,l) + &
+                berni(ij,l) = pk(ij,l) &
                      + 1/(4*Ai(ij))*(le(ij+u_right)*de(ij+u_right)*u(ij+u_right,l)**2 +    &
                                      le(ij+u_rup)*de(ij+u_rup)*u(ij+u_rup,l)**2 +          &

codes/icosagcm/trunk/src/etat0_williamson.f90

@@ -12 +12 @@
   
   
-  SUBROUTINE etat0_williamson(f_h,f_u)
+  SUBROUTINE etat0_williamson(f_h,f_u) ! Deprecated / to be removed
   USE icosa
   IMPLICIT NONE

codes/icosagcm/trunk/src/wind.f90

@@ -184 +184 @@
   END SUBROUTINE compute_wind_from_lonlat_compound
  
+  SUBROUTINE compute_wind2D_from_lonlat_compound(ulon, ulat, u)
+  USE icosa  
+  
+  IMPLICIT NONE
+  REAL(rstd) :: u(3*iim*jjm,3)
+  REAL(rstd) :: ulon(3*iim*jjm)
+  REAL(rstd) :: ulat(3*iim*jjm)
+  
+  INTEGER :: i,j,ij
+  
+  DO j=jj_begin-1,jj_end+1
+     DO i=ii_begin-1,ii_end+1
+        ij=(j-1)*iim+i
+        u(ij+u_right,:)=ulon(ij+u_right)*elon_e(ij+u_right,:)+ ulat(ij+u_right)*elat_e(ij+u_right,:)
+        u(ij+u_lup,:)=ulon(ij+u_lup)*elon_e(ij+u_lup,:) + ulat(ij+u_lup)*elat_e(ij+u_lup,:)
+        u(ij+u_ldown,:)=ulon(ij+u_ldown)*elon_e(ij+u_ldown,:) + ulat(ij+u_ldown)*elat_e(ij+u_ldown,:)
+     ENDDO
+  ENDDO
+    
+  END SUBROUTINE compute_wind2D_from_lonlat_compound
+ 
   SUBROUTINE compute_wind_perp_from_lonlat_compound(ulon, ulat, up)
   USE icosa  
@@ -210 +231 @@
   END SUBROUTINE compute_wind_perp_from_lonlat_compound
    
+  SUBROUTINE compute_wind2D_perp_from_lonlat_compound(ulon, ulat, up)
+  USE icosa  
+    
+  IMPLICIT NONE
+  REAL(rstd) :: up(3*iim*jjm)
+  REAL(rstd) :: ulon(3*iim*jjm)
+  REAL(rstd) :: ulat(3*iim*jjm)
+  REAL(rstd) :: u(3*iim*jjm,3)
+
+  INTEGER :: i,j,ij 
+  
+  CALL compute_wind2D_from_lonlat_compound(ulon, ulat, u)
+  DO j=jj_begin-1,jj_end+1
+     DO i=ii_begin-1,ii_end+1
+        ij=(j-1)*iim+i
+        up(ij+u_right)=sum(u(ij+u_right,:)*ep_e(ij+u_right,:))
+        up(ij+u_lup)=sum(u(ij+u_lup,:)*ep_e(ij+u_lup,:))
+        up(ij+u_ldown)=sum(u(ij+u_ldown,:)*ep_e(ij+u_ldown,:))
+     ENDDO
+  ENDDO
+   
+  END SUBROUTINE compute_wind2D_perp_from_lonlat_compound
+   
  SUBROUTINE compute_wind_centered_lonlat_compound(uc, ulon, ulat)
   USE icosa