Changeset 397

Timestamp:

06/06/16 20:41:52 (8 years ago)

Author:

ymipsl

Message:

Prepare DCMIP2016 output by XIOS2

YM

Location:

codes/icosagcm/trunk/src

Files:

• observable.f90 (modified) (5 diffs)
• physics_dcmip2016.f90 (modified) (4 diffs)
• pression.f90 (modified) (2 diffs)
• timeloop_gcm.f90 (modified) (1 diff)

codes/icosagcm/trunk/src/observable.f90

@@ -10 +10 @@
   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(:)
+  TYPE(t_field),POINTER, SAVE :: f_pmid(:)
 
 ! temporary shared variable for caldyn
@@ -15 +16 @@
 
   PUBLIC init_observable, write_output_fields_basic, f_theta
-
+  LOGICAL,SAVE :: first_output=.TRUE.
+!$OMP THREADPRIVATE(first_output)
+  
 CONTAINS
   
   SUBROUTINE init_observable
     CALL allocate_field(f_buf_i,   field_t,type_real,llm,name="buffer_i")
+    CALL allocate_field(f_buf1_i,   field_t,type_real,llm,name="buffer1_i")
+    CALL allocate_field(f_buf2_i,   field_t,type_real,llm,name="buffer2_i")
     CALL allocate_field(f_buf_p,   field_t,type_real,llm+1) 
     CALL allocate_field(f_buf_u3d, field_t,type_real,3,llm)  ! 3D vel at cell centers
@@ -29 +34 @@
 
     CALL allocate_field(f_theta, field_t,type_real,llm,  name='theta')   ! potential temperature
+    CALL allocate_field(f_pmid,  field_t,type_real,llm,  name='pmid')   ! mid layer pressure
   END SUBROUTINE init_observable
   
-  SUBROUTINE write_output_fields_basic(f_ps, f_mass, f_geopot, f_u, f_W, f_q)
+  SUBROUTINE write_output_fields_basic(f_ps, f_mass, f_geopot, f_theta_rhodz, f_u, f_W, f_q)
     USE wind_mod
     USE output_field_mod
     USE omp_para
-    TYPE(t_field),POINTER :: f_ps(:), f_mass(:), f_geopot(:), f_u(:), f_W(:), f_q(:)
+    USE time_mod
+    USE xios
+    USE disvert_mod
+    USE earth_const
+    USE pression_mod
+    USE vertical_interp_mod
+    USE theta2theta_rhodz_mod
+    USE wind_mod
+    USE omega_mod
+    
+    TYPE(t_field),POINTER :: f_ps(:), f_mass(:), f_geopot(:), f_theta_rhodz(:), f_u(:), f_W(:), f_q(:)
 !    IF (is_master) PRINT *,'CALL write_output_fields_basic'
+    REAL :: scalar(1)
+    REAL :: mid_ap(llm)
+    REAL :: mid_bp(llm)
+    INTEGER :: l
+
+    IF (first_output) THEN
+      scalar(1)=dt
+      CALL xios_send_field("timestep", scalar)
+      scalar(1)=preff
+      CALL xios_send_field("preff", scalar)
+      CALL xios_send_field("ap",ap)
+      CALL xios_send_field("bp",bp)
+      DO l=1,llm
+        mid_ap(l)=(ap(l)+ap(l+1))/2
+        mid_bp(l)=(bp(l)+bp(l+1))/2
+      ENDDO
+      CALL xios_send_field("mid_ap",mid_ap)
+      CALL xios_send_field("mid_bp",mid_bp)
+        
+      first_output=.FALSE.
+    ENDIF
+    
     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 vertical_interp(f_ps,f_buf_i,f_buf_s,85000.)
+    CALL output_field("w850",f_buf_s)
+    CALL vertical_interp(f_ps,f_buf_i,f_buf_s,50000.)
+    CALL output_field("w500",f_buf_s)
+    
+    
     CALL un2ulonlat(f_buf_uh, f_buf_ulon, f_buf_ulat)
     CALL output_field("ulon",f_buf_ulon)
@@ -45 +89 @@
     CALL output_field("ps",f_ps)
     CALL output_field("Ai",geom%Ai)
+
     !       CALL output_field("dps",f_dps)
     CALL output_field("mass",f_mass)
@@ -50 +95 @@
     !       CALL output_field("dmass",f_dmass)
     !       CALL output_field("vort",f_qv)
-    CALL output_field("theta",f_theta)
+    
+    
+    CALL theta_rhodz2temperature(f_ps,f_theta_rhodz,f_buf_i) ; 
+    CALL Tv2T(f_buf_i,f_q,f_buf1_i) 
+    CALL output_field("temp",f_buf_i)
+    CALL vertical_interp(f_ps,f_buf_i,f_buf_s,85000.)
+    CALL output_field("t850",f_buf_s)
+    CALL vertical_interp(f_ps,f_buf_i,f_buf_s,50000.)
+    CALL output_field("t500",f_buf_s)
+    CALL vertical_interp(f_ps,f_buf_i,f_buf_s,preff)
+    CALL output_field("SST",f_buf_s)
+
+
+    CALL extract_slice(f_theta_rhodz, f_buf_i,1)
+    CALL output_field("theta",f_buf_i)
+           
     !       CALL output_field("exner",f_pk)
     !       CALL output_field("pv",f_qv)
     CALL output_field("q",f_q)
+    CALL pression_mid(f_ps, f_pmid)
+    CALL output_field("p",f_pmid)
+
+    CALL vertical_interp(f_ps,f_buf_ulon,f_buf_s,85000.)
+    CALL output_field("u850",f_buf_s)
+    CALL vertical_interp(f_ps,f_buf_ulon,f_buf_s,50000.)
+    CALL output_field("u500",f_buf_s)
+
+    CALL vertical_interp(f_ps,f_buf_ulat,f_buf_s,85000.)
+    CALL output_field("v850",f_buf_s)
+    CALL vertical_interp(f_ps,f_buf_ulat,f_buf_s,50000.)
+    CALL output_field("v500",f_buf_s)
+    
+    CALL w_omega(f_ps, f_u, f_buf_i)
+    CALL output_field("omega",f_buf_i)
+    CALL vertical_interp(f_ps,f_buf_i,f_buf_s,85000.)
+    CALL output_field("omega850",f_buf_s)
+    CALL vertical_interp(f_ps,f_buf_i,f_buf_s,50000.)
+    CALL output_field("omega500",f_buf_s)
+
+
+    
   END SUBROUTINE write_output_fields_basic
   

codes/icosagcm/trunk/src/physics_dcmip2016.f90

@@ -120 +120 @@
     REAL(rstd) :: qfi(ngrid,llm,5)
 
-    REAL(rstd)  :: rho(llm), z(llm), theta(llm)
+    REAL(rstd)  :: rho(llm), z(llm), theta(llm), qv(llm),qc(llm),qr(llm)
     REAL(rstd)  :: lastz
     REAL(rstd)  :: dcl1,dcl2
-    INTEGER :: l,ll,ij
+     INTEGER :: l,ll,ij
     REAL(rstd) :: dt_phys, inv_dt
-
+    INTEGER :: simple_physic_testcase
+    
     ! prepare input fields and mirror vertical index      
     ps(:) = p(:,1) ! surface pressure
@@ -147 +148 @@
     ENDDO
     
+    
     IF (testcase==moist_baroclinic .OR. testcase==cyclone  ) THEN
-      CALL simple_physics(ngrid, llm, dt_phys, lat, tfi, qfi(:,:,1) , ufi, vfi, pmid, pint, pdel, 1./pdel, ps, precl, 1, .FALSE., .FALSE.)
+      IF (testcase==moist_baroclinic) simple_physic_testcase=1
+      IF (testcase==cyclone) simple_physic_testcase=0
+      CALL simple_physics(ngrid, llm, dt_phys, lat, tfi, qfi(:,:,1) , ufi, vfi, pmid, pint, pdel, 1./pdel, ps, precl, &
+                          simple_physic_testcase, .TRUE., .FALSE.)
     ENDIF
 
@@ -164 +169 @@
           ENDDO
 
-          CALL KESSLER(theta(:), qfi(ij,llm:1:-1,1), qfi(ij,llm:1:-1,2), qfi(ij,llm:1:-1,3), rho(:),  &
-                       pk(ij,:), dt_phys, z(:), llm, precl(ij)) 
+          qv(:)=qfi(ij,llm:1:-1,1)
+          qc(:)=qfi(ij,llm:1:-1,2)
+          qr(:)=qfi(ij,llm:1:-1,3)
+          
+!          CALL KESSLER(theta(:), qv, qc, qr, rho(:),  &
+!                       pk(ij,:), dt_phys, z(:), llm, precl(ij)) 
+          
           
           DO l=1,llm
@@ -172 +182 @@
            Tfi(ij,ll) = theta(l)  * ( pk(ij,l) / cpp)
           ENDDO
+
+          qfi(ij,:,1)=qv(llm:1:-1)
+          qfi(ij,:,2)=qc(llm:1:-1)
+          qfi(ij,:,3)=qr(llm:1:-1)
+
        ENDDO
     ENDIF

codes/icosagcm/trunk/src/pression.f90

@@ -25 +25 @@
   
   END SUBROUTINE pression
+
+  SUBROUTINE pression_mid(f_ps,f_pmid)
+  USE icosa
+  IMPLICIT NONE
+    TYPE(t_field), POINTER :: f_ps(:)
+    TYPE(t_field), POINTER :: f_pmid(:)
   
+    REAL(rstd), POINTER :: ps(:)
+    REAL(rstd), POINTER :: pmid(:,:)
+    INTEGER :: ind
+
+!$OMP BARRIER
+    DO ind=1,ndomain
+      IF (.NOT. assigned_domain(ind)) CYCLE
+      CALL swap_dimensions(ind)
+      CALL swap_geometry(ind)
+      ps=f_ps(ind)
+      pmid=f_pmid(ind)
+      CALL compute_pression_mid(ps, pmid,0)
+    ENDDO
+!$OMP BARRIER
+  
+  END SUBROUTINE pression_mid
+
   SUBROUTINE compute_pression(ps,p,offset)
   USE icosa
@@ -49 +72 @@
 
   END SUBROUTINE compute_pression
+  
+  SUBROUTINE compute_pression_mid(ps,pmid,offset)
+  USE icosa
+  USE disvert_mod
+  USE omp_para
+  IMPLICIT NONE
+    REAL(rstd),INTENT(IN) :: ps(iim*jjm)
+    REAL(rstd),INTENT(OUT) :: pmid(iim*jjm,llm)
+    INTEGER,INTENT(IN) :: offset
+    INTEGER :: i,j,ij,l
+
+    IF(ap_bp_present) THEN
+      DO    l    = ll_begin, ll_end
+        DO j=jj_begin-offset,jj_end+offset
+          DO i=ii_begin-offset,ii_end+offset
+            ij=(j-1)*iim+i
+            pmid(ij,l) = 0.5*(ap(l)+ap(l+1) + (bp(l)+bp(l+1)) * ps(ij))
+          ENDDO
+        ENDDO
+      ENDDO
+    END IF
+
+  END SUBROUTINE compute_pression_mid
 
 END MODULE pression_mod

codes/icosagcm/trunk/src/timeloop_gcm.f90

@@ -251 +251 @@
        IF (mod(it,itau_out)==0 ) THEN
           CALL transfert_request(f_u,req_e1_vect)
-          CALL write_output_fields_basic(f_ps, f_mass, f_geopot, f_u, f_W, f_q)
+          CALL write_output_fields_basic(f_ps, f_mass, f_geopot, f_theta_rhodz, f_u, f_W, f_q)
        ENDIF