Changeset 592

Timestamp:

10/17/17 23:11:14 (7 years ago)

Author:

dubos

Message:

devel : finalize diagnostics of tracer fluxes

Location:

codes/icosagcm/devel/src

Files:

• diagnostics/diagflux.f90 (modified) (5 diffs)
• diagnostics/observable.f90 (modified) (1 diff)
• time/timeloop_gcm.f90 (modified) (3 diffs)
• transport/advect_tracer.f90 (modified) (7 diffs)

codes/icosagcm/devel/src/diagnostics/diagflux.f90

@@ -4 +4 @@
   SAVE
   
-  TYPE(t_field),POINTER :: f_qfluxt(:), f_qfluxt_lon(:), f_qfluxt_lat(:) ! time-integrated flux of scalars and its reconstruction at cell centers
+  TYPE(t_field),POINTER :: &
+       f_masst(:), f_qmasst(:), & ! time-integrated mass, tracer mass, 
+       f_massfluxt(:), f_qfluxt(:), & ! mass flux and tracer flux
+       f_qfluxt_lon(:), f_qfluxt_lat(:) ! scalar flux reconstructed  cell centers
   LOGICAL :: diagflux_on
   !$OMP THREADPRIVATE(diagflux_on)
@@ -15 +18 @@
     CALL getin("diagflux", diagflux_on)
     IF(diagflux_on) THEN
+       CALL allocate_field(f_masst,  field_t,type_real,llm, name="masst")
+       CALL allocate_field(f_qmasst,  field_t,type_real,llm,nqtot, name="qmasst")
+       CALL allocate_field(f_massfluxt, field_u,type_real,llm, name="massfluxt")
        CALL allocate_field(f_qfluxt,  field_u,type_real,llm,nqtot, name="qfluxt")
        CALL allocate_field(f_qfluxt_lon,  field_t,type_real,llm,nqtot, name="qfluxt_lon")
@@ -20 +26 @@
        CALL zero_qfluxt
     ELSE
+       CALL allocate_field(f_masst,  field_t,type_real,0, name="masst")
+       CALL allocate_field(f_qmasst,  field_t,type_real,llm,0, name="qmasst")
+       CALL allocate_field(f_massfluxt, field_u,type_real,0, name="massfluxt")
        CALL allocate_field(f_qfluxt,  field_u,type_real,llm,0, name="qfluxt")
        CALL allocate_field(f_qfluxt_lon,  field_t,type_real,llm,0, name="qfluxt_lon")
@@ -30 +39 @@
     USE omp_para
     INTEGER :: ind
-    REAL(rstd), POINTER :: qfluxt(:,:,:)
+    REAL(rstd), POINTER :: buf2(:,:),buf3(:,:,:)
     DO ind=1,ndomain
        IF (.NOT. assigned_domain(ind)) CYCLE
        CALL swap_dimensions(ind)
-       qfluxt=f_qfluxt(ind)
-       qfluxt(:,ll_begin:ll_end,:)=0.
+       buf2=f_masst(ind)
+       buf2(:,ll_begin:ll_end)=0.
+       buf2=f_massfluxt(ind)
+       buf2(:,ll_begin:ll_end)=0.
+       buf3=f_qmasst(ind)
+       buf3(:,ll_begin:ll_end,:)=0.
+       buf3=f_qfluxt(ind)
+       buf3(:,ll_begin:ll_end,:)=0.
     END DO
   END SUBROUTINE zero_qfluxt
 
-  SUBROUTINE flux_centered_lonlat(scale, f_flux, f_flux_lon, f_flux_lat)
+  SUBROUTINE flux_centered_lonlat(scale, f_massflux, f_flux, f_massflux_lon, f_massflux_lat, f_flux_lon, f_flux_lat)
     REAL(rstd), INTENT(IN) :: scale
-    TYPE(t_field),POINTER :: f_flux(:), f_flux_lon(:), f_flux_lat(:)
-    REAL(rstd), POINTER :: flux(:,:,:), flux_lon(:,:,:), flux_lat(:,:,:)
+    TYPE(t_field),POINTER :: f_flux(:), f_flux_lon(:), f_flux_lat(:), &
+         f_massflux(:), f_massflux_lon(:), f_massflux_lat(:)
+    REAL(rstd), POINTER :: flux(:,:,:), flux_lon(:,:,:), flux_lat(:,:,:), &
+         massflux(:,:), massflux_lon(:,:), massflux_lat(:,:)
     INTEGER :: ind, itrac
     DO ind=1,ndomain
@@ -54 +71 @@
           CALL compute_flux_centered_lonlat(scale, flux(:,:,itrac), flux_lon(:,:,itrac), flux_lat(:,:,itrac))
        END DO
+       massflux=f_massflux(ind)
+       massflux_lon=f_massflux_lon(ind)
+       massflux_lat=f_massflux_lat(ind)
+       CALL compute_flux_centered_lonlat(scale, massflux, massflux_lon, massflux_lat)
     END DO
   END SUBROUTINE flux_centered_lonlat

codes/icosagcm/devel/src/diagnostics/observable.f90

@@ -163 +163 @@
     IF(.NOT. init) THEN
        IF(diagflux_on) THEN
-          CALL flux_centered_lonlat(1./(itau_out*dt) , f_qfluxt, f_qfluxt_lon, f_qfluxt_lat)
+          CALL output_field("mass_t", f_masst)
+          CALL output_field("qmass_t", f_qmasst)
+          CALL flux_centered_lonlat(1./(itau_out*dt) , f_massfluxt, f_qfluxt, &
+               f_buf_ulon, f_buf_ulat, f_qfluxt_lon, f_qfluxt_lat)
+          CALL output_field("massflux_lon",f_buf_ulon)
+          CALL output_field("massflux_lat",f_buf_ulat)
           CALL output_field("qflux_lon",f_qfluxt_lon)
           CALL output_field("qflux_lat",f_qfluxt_lat)

codes/icosagcm/devel/src/time/timeloop_gcm.f90

@@ -182 +182 @@
     REAL(rstd),POINTER :: rhodz(:,:), mass(:,:), ps(:)
 
-    INTEGER :: ind
-    INTEGER :: it,i,j,l,n,  stage
-    LOGICAL :: fluxt_zero(ndomain) ! set to .TRUE. to start accumulating fluxes in time
+    REAL(rstd) :: adv_over_out ! ratio itau_adv/itau_out
+    INTEGER :: ind, it,i,j,l,n,  stage
+    LOGICAL :: fluxt_zero(ndomain) ! set to .TRUE. to start accumulating mass fluxes in time
     LOGICAL, PARAMETER :: check_rhodz=.FALSE.
     INTEGER :: start_clock, stop_clock, rate_clock
@@ -211 +211 @@
 
     IF(positive_theta) CALL copy_theta_to_q(f_theta_rhodz,f_rhodz,f_q)
+    IF(diagflux_on) THEN
+       adv_over_out = itau_adv*(1./itau_out)
+    ELSE
+       adv_over_out = 0.
+    END IF
 
     CALL check_conserve(f_ps,f_dps,f_u,f_theta_rhodz,f_phis,itau0)  
@@ -299 +304 @@
        
        IF(MOD(it,itau_adv)==0) THEN
-          CALL advect_tracer(diagflux_on, f_hfluxt,f_wfluxt,f_u, f_q,f_rhodz,f_qfluxt)  ! update q and rhodz after RK step
+          CALL advect_tracer(f_hfluxt,f_wfluxt,f_u, f_q,f_rhodz, & ! update q and rhodz after RK step
+               adv_over_out, f_masst,f_qmasst,f_massfluxt, f_qfluxt)  ! accumulate mass and fluxes if diagflux_on
           fluxt_zero=.TRUE.
           ! FIXME : check that rhodz is consistent with ps

codes/icosagcm/devel/src/transport/advect_tracer.f90

@@ -54 +54 @@
   END SUBROUTINE init_advect_tracer
 
-  SUBROUTINE advect_tracer(diagflux_on, f_hfluxt, f_wfluxt,f_u, f_q,f_rhodz,f_qfluxt)
-    USE mpipara
+  SUBROUTINE advect_tracer(f_hfluxt, f_wfluxt,f_u, f_q,f_rhodz,&
+       frac, f_masst,f_qmasst,f_massfluxt,f_qfluxt)
+    USE omp_para
     USE trace
     USE write_field_mod
     USE tracer_mod
-    LOGICAL, INTENT(IN) :: diagflux_on
     TYPE(t_field),POINTER :: f_hfluxt(:)   ! time-integrated horizontal mass flux
     TYPE(t_field),POINTER :: f_wfluxt(:)   ! time-integrated vertical mass flux
@@ -65 +65 @@
     TYPE(t_field),POINTER :: f_q(:)        ! tracer
     TYPE(t_field),POINTER :: f_rhodz(:)    ! mass field at beginning of macro time step
+    REAL(rstd), INTENT(in):: frac          ! ratio itau_adv/itau_out or 0. if not diagflux_on
+    TYPE(t_field),POINTER :: f_masst(:)    ! time-integrated mass
+    TYPE(t_field),POINTER :: f_qmasst(:)   ! time-integrated tracer mass
+    TYPE(t_field),POINTER :: f_massfluxt(:)! time-integrated horizontal mass flux
     TYPE(t_field),POINTER :: f_qfluxt(:)   ! time-integrated horizontal tracer flux
 
     REAL(rstd),POINTER :: q(:,:,:), normal(:,:), tangent(:,:), sqrt_leng(:), gradq3d(:,:,:), cc(:,:,:)
-    REAL(rstd),POINTER :: hfluxt(:,:), wfluxt(:,:), qfluxt(:,:,:)
+    REAL(rstd),POINTER :: hfluxt(:,:), wfluxt(:,:), masst(:,:), qmasst(:,:,:), massfluxt(:,:), qfluxt(:,:,:)
     REAL(rstd),POINTER :: rhodz(:,:), u(:,:) 
 ! temporary shared variable for vlz
@@ -76 +80 @@
     REAL(rstd),POINTER ::  wq(:,:)           ! time-integrated flux of q
     
-     INTEGER :: ind,k, nq_last
+    INTEGER :: ind,k, nq_last
     LOGICAL,SAVE :: first=.TRUE.
 !$OMP THREADPRIVATE(first)
@@ -140 +144 @@
     CALL send_message(f_cc,req_cc)
 
-
     ! horizontal transport - split in two to place transfer of gradq3d
     DO k = 1, nqtot
@@ -153 +156 @@
           sqrt_leng=f_sqrt_leng(ind)
           CALL compute_gradq3d(q(:,:,k),sqrt_leng,gradq3d,xyz_i,xyz_v)
-
         END DO
 
@@ -159 +161 @@
         CALL wait_message(req_cc)
         CALL wait_message(req_gradq3d)
-
 
         DO ind=1,ndomain
@@ -171 +172 @@
           qfluxt  = f_qfluxt(ind) 
           gradq3d = f_gradq3d(ind)
-          CALL compute_advect_horiz(k==nq_last,diagflux_on, hfluxt,cc,gradq3d, rhodz,q(:,:,k),qfluxt(:,:,k))
+
+          IF(frac>0.) THEN ! accumulate mass, mass flux and tracer mass
+             qmasst  = f_qmasst(ind) 
+             qmasst(:,ll_begin:ll_end,k) = qmasst(:,ll_begin:ll_end,k) + &
+                  frac*rhodz(:,ll_begin:ll_end)*q(:,ll_begin:ll_end,k)
+             IF(k==nq_last) THEN
+                masst  = f_masst(ind)
+                massfluxt  = f_massfluxt(ind)
+                masst(:,ll_begin:ll_end) = masst(:,ll_begin:ll_end)+frac*rhodz(:,ll_begin:ll_end)
+                massfluxt(:,ll_begin:ll_end) = massfluxt(:,ll_begin:ll_end)+hfluxt(:,ll_begin:ll_end)
+             END IF
+          END IF
+          CALL compute_advect_horiz(k==nq_last,frac>0., hfluxt,cc,gradq3d, rhodz, q(:,:,k), qfluxt(:,:,k))
         END DO
+
       ENDIF
     END DO