Changeset 374 for codes/icosagcm

Timestamp:

04/08/16 09:41:16 (8 years ago)

Author:

dubos

Message:

More NH diagnostics - bugfix needed

Location:

codes/icosagcm/trunk/src

Files:

• observable.f90 (modified) (5 diffs)
• timeloop_gcm.f90 (modified) (2 diffs)

codes/icosagcm/trunk/src/observable.f90

@@ -4 +4 @@
   PRIVATE
 
-  TYPE(t_field),POINTER, SAVE :: f_buf_i(:), f_buf_ulon(:), f_buf_ulat(:), f_buf_u3d(:)
+  TYPE(t_field),POINTER, SAVE :: f_buf_i(:), &
+       f_buf_uh(:), & ! horizontal velocity, different from prognostic velocity if NH
+       f_buf_ulon(:), f_buf_ulat(:), &
+       f_buf_u3d(:) ! unused, remove ?
   TYPE(t_field),POINTER, SAVE :: f_buf1_i(:), f_buf2_i(:)
   TYPE(t_field),POINTER, SAVE :: f_buf_v(:), f_buf_s(:), f_buf_p(:)
@@ -21 +24 @@
     CALL allocate_field(f_buf_ulon,field_t,type_real,llm, name="buf_ulon")
     CALL allocate_field(f_buf_ulat,field_t,type_real,llm, name="buf_ulat")
-    CALL allocate_field(f_buf_v,   field_z,type_real,llm)
-    CALL allocate_field(f_buf_s,   field_t,type_real)
+    CALL allocate_field(f_buf_uh,  field_u,type_real,llm, name="buf_uh")
+    CALL allocate_field(f_buf_v,   field_z,type_real,llm, name="buf_v")
+    CALL allocate_field(f_buf_s,   field_t,type_real, name="buf_s")
 
     CALL allocate_field(f_theta, field_t,type_real,llm,  name='theta')   ! potential temperature
   END SUBROUTINE init_observable
   
-  SUBROUTINE write_output_fields_basic(f_ps, f_u, f_q)
+  SUBROUTINE write_output_fields_basic(f_ps, f_mass, f_geopot, f_u, f_W, f_q)
     USE wind_mod
     USE output_field_mod
     USE omp_para
-    TYPE(t_field),POINTER :: f_ps(:), f_u(:), f_q(:)
+    TYPE(t_field),POINTER :: f_ps(:), f_mass(:), f_geopot(:), f_u(:), f_W(:), f_q(:)
 !    IF (is_master) PRINT *,'CALL write_output_fields_basic'
-    CALL un2ulonlat(f_u, f_buf_ulon, f_buf_ulat)
+    CALL progonostic_vel_to_horiz(f_geopot, f_ps, f_mass, f_u, f_W, f_buf_uh, f_buf_i)
+    CALL transfert_request(f_buf_uh,req_e1_vect) 
+    CALL output_field("uz",f_buf_i)
+    CALL un2ulonlat(f_buf_uh, f_buf_ulon, f_buf_ulat)
     CALL output_field("ulon",f_buf_ulon)
     CALL output_field("ulat",f_buf_ulat)
@@ -39 +46 @@
     CALL output_field("Ai",geom%Ai)
     !       CALL output_field("dps",f_dps)
-    !CALL output_field("mass",f_mass)
+    CALL output_field("mass",f_mass)
+    CALL output_field("geopot",f_geopot)
     !       CALL output_field("dmass",f_dmass)
     !       CALL output_field("vort",f_qv)
@@ -50 +58 @@
   SUBROUTINE write_output_fields(f_ps, f_phis, f_dps, f_u, f_theta_rhodz, f_q, &
        f_buf_i, f_buf_v, f_buf_i3, f_buf1_i, f_buf2_i, f_buf_s, f_buf_p)
-    USE icosa
     USE vorticity_mod
     USE theta2theta_rhodz_mod
@@ -122 +129 @@
   
   END SUBROUTINE write_output_fields
-  
+
+!------------------- Conversion from prognostic to observable variables ------------------
+
+  SUBROUTINE progonostic_vel_to_horiz(f_geopot, f_ps, f_rhodz, f_u, f_W, f_uh, f_uz)
+    USE disvert_mod
+    TYPE(t_field), POINTER :: f_geopot(:), f_ps(:), f_rhodz(:), &
+         f_u(:), f_W(:), f_uz(:), &  ! IN
+         f_uh(:)                         ! OUT
+    REAL(rstd),POINTER :: geopot(:,:), ps(:), rhodz(:,:), u(:,:), W(:,:), uh(:,:), uz(:,:)
+    INTEGER :: ind
+    
+    DO ind=1,ndomain
+       IF (.NOT. assigned_domain(ind)) CYCLE
+       CALL swap_dimensions(ind)
+       CALL swap_geometry(ind)
+       geopot = f_geopot(ind)
+       rhodz = f_rhodz(ind)
+       u = f_u(ind)
+       W = f_W(ind)
+       uh  = f_uh(ind)
+       IF(caldyn_eta==eta_mass) THEN
+          ps=f_ps(ind)
+          CALL compute_rhodz(.TRUE., ps, rhodz)
+       END IF
+       uz = f_uz(ind)
+       CALL compute_prognostic_vel_to_horiz(geopot,rhodz,u,W,uh,uz)
+    END DO
+  END SUBROUTINE progonostic_vel_to_horiz
+  
+  SUBROUTINE compute_prognostic_vel_to_horiz(Phi, rhodz, u, W, uh, uz)
+    USE omp_para
+    REAL(rstd), INTENT(IN) :: Phi(iim*jjm,llm+1)
+    REAL(rstd), INTENT(IN) :: rhodz(iim*jjm,llm)
+    REAL(rstd), INTENT(IN) :: u(3*iim*jjm,llm)
+    REAL(rstd), INTENT(IN) :: W(iim*jjm,llm+1)
+    REAL(rstd), INTENT(OUT) :: uh(3*iim*jjm,llm)
+    REAL(rstd), INTENT(OUT) :: uz(iim*jjm,llm)
+    INTEGER :: ij,l
+    REAL(rstd) :: F_el(3*iim*jjm,llm+1)
+    REAL(rstd) :: uu_right, uu_lup, uu_ldown, W_el, DePhil
+    IF(hydrostatic) THEN
+       uh(:,:)=u(:,:)
+       uz(:,:)=0.
+    ELSE
+       DO l=ll_begin, ll_endp1 ! compute on l levels (interfaces)
+          DO ij=ij_begin_ext, ij_end_ext
+             ! Compute on edge 'right'
+             W_el   = .5*( W(ij,l)+W(ij+t_right,l) )
+             DePhil = ne_right*(Phi(ij+t_right,l)-Phi(ij,l))
+             F_el(ij+u_right,l)   = DePhil*W_el
+             ! Compute on edge 'lup'
+             W_el   = .5*( W(ij,l)+W(ij+t_lup,l) )
+             DePhil = ne_lup*(Phi(ij+t_lup,l)-Phi(ij,l))
+             F_el(ij+u_lup,l)   = DePhil*W_el
+             ! Compute on edge 'ldown'
+             W_el   = .5*( W(ij,l)+W(ij+t_ldown,l) )
+             DePhil = ne_ldown*(Phi(ij+t_ldown,l)-Phi(ij,l))
+             F_el(ij+u_ldown,l) = DePhil*W_el
+          END DO
+       END DO
+
+       DO l=ll_begin, ll_end ! compute on k levels (full levels)
+          DO ij=ij_begin_ext, ij_end_ext
+             ! w = vertical momentum = g^-2*dPhi/dt = g*uz
+             uz(ij,l) = (.5*g)*(W(ij,l)+W(ij,l+1))/rhodz(ij,l)
+             ! uh = u-w.grad(Phi) = u - uz.grad(z)
+             uh(ij+u_right,l) = u(ij+u_right,l) - (F_el(ij+u_right,l)+F_el(ij+u_right,l+1)) / (rhodz(ij,l)+rhodz(ij+t_right,l))
+             uh(ij+u_lup,l)   = u(ij+u_lup,l)   - (F_el(ij+u_lup,l)+F_el(ij+u_lup,l+1))     / (rhodz(ij,l)+rhodz(ij+t_lup,l))
+             uh(ij+u_ldown,l) = u(ij+u_ldown,l) - (F_el(ij+u_ldown,l)+F_el(ij+u_ldown,l+1)) / (rhodz(ij,l)+rhodz(ij+t_ldown,l))
+          END DO
+       END DO
+
+    END IF
+  END SUBROUTINE compute_prognostic_vel_to_horiz
+
   SUBROUTINE thetarhodz2geopot(f_ps,f_phis,f_theta_rhodz, f_pks,f_p,f_theta,f_pk,f_phi) 
     USE field_mod

codes/icosagcm/trunk/src/timeloop_gcm.f90

@@ -34 +34 @@
     IF (.NOT. enable_io) itau_out=HUGE(itau_out)
     
-    def='ARK2.3'
+    def='RK2.5'
     CALL getin('time_scheme',def)
     SELECT CASE (TRIM(def))
@@ -231 +231 @@
        IF (mod(it,itau_out)==0 ) THEN
           CALL transfert_request(f_u,req_e1_vect)
-          CALL write_output_fields_basic(f_ps, f_u, f_q)
+          CALL write_output_fields_basic(f_ps, f_mass, f_geopot, f_u, f_W, f_q)
        ENDIF