Changeset 668

Timestamp:

01/27/18 00:10:07 (6 years ago)

Author:

dubos

Message:

trunk : compute and output physics tendencies on the dynamics side

Location:

codes/icosagcm/trunk/src

Files:

• diagnostics/observable.f90 (modified) (1 diff)
• diagnostics/wind.F90 (modified) (4 diffs)
• icosa_init.f90 (modified) (1 diff)
• physics/physics.f90 (modified) (12 diffs)

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

@@ -15 +15 @@
   TYPE(t_field),POINTER, SAVE :: f_theta(:)
 
-  PUBLIC init_observable, write_output_fields_basic, f_theta, f_buf_i
+  PUBLIC init_observable, write_output_fields_basic, & 
+       f_theta, f_buf_i, f_buf_ulon, f_buf_ulat
 
 CONTAINS

codes/icosagcm/trunk/src/diagnostics/wind.F90

@@ -12 +12 @@
 CONTAINS
 
-  SUBROUTINE un2ulonlat(f_u, f_ulon, f_ulat)
+  SUBROUTINE un2ulonlat(f_u, f_ulon, f_ulat, scale_)
     TYPE(t_field), POINTER :: f_u(:) ! IN  : normal velocity components on edges
     TYPE(t_field), POINTER :: f_ulon(:), f_ulat(:) ! OUT : velocity reconstructed at hexagons    
-    REAL(rstd),POINTER :: u(:,:),  ulon(:,:), ulat(:,:)
+    REAL(rstd),POINTER     :: u(:,:),  ulon(:,:), ulat(:,:)
+    REAL(rstd), OPTIONAL   :: scale_ 
+    REAL(rstd)             :: scale
     INTEGER :: ind
-
+    scale = MERGE(scale_, 1., PRESENT(scale_))
     DO ind=1,ndomain
        IF (.NOT. assigned_domain(ind)) CYCLE
@@ -25 +27 @@
        ulon=f_ulon(ind)
        ulat=f_ulat(ind)
-       CALL compute_un2ulonlat(u,ulon, ulat)
+       CALL compute_un2ulonlat(u,ulon, ulat, scale)
     END DO
 
@@ -49 +51 @@
   END SUBROUTINE ulonlat2un
  
-  SUBROUTINE compute_wind_centered(ue,ucenter)
+  SUBROUTINE compute_wind_centered(ue,ucenter,scale_)
   REAL(rstd) :: ue(3*iim*jjm,llm)
   REAL(rstd) :: ucenter(iim*jjm,llm,3)
+  REAL(rstd), INTENT(IN), OPTIONAL :: scale_
   INTEGER :: ij,l
-  REAL(rstd), PARAMETER :: scale=1.
-  REAL(rstd) :: fac, ue_le, cx,cy,cz, ux,uy,uz
+  REAL(rstd) :: scale,fac, ue_le, cx,cy,cz, ux,uy,uz
+  scale = MERGE(scale_, 1., PRESENT(scale_))
 #include "../kernels/wind_centered.k90"
  END SUBROUTINE compute_wind_centered
@@ -328 +331 @@
 
 
- SUBROUTINE compute_un2ulonlat(un, ulon, ulat)
+ SUBROUTINE compute_un2ulonlat(un, ulon, ulat, scale)
   REAL(rstd),INTENT(IN)  :: un(3*iim*jjm,llm)
   REAL(rstd),INTENT(OUT) :: ulon(iim*jjm,llm)
   REAL(rstd),INTENT(OUT) :: ulat(iim*jjm,llm)
-
   REAL(rstd)             :: uc(iim*jjm,llm,3)
-    
-  CALL compute_wind_centered(un,uc)  
+  REAL(rstd)             :: scale
+  CALL compute_wind_centered(un,uc,scale)
   CALL compute_wind_centered_lonlat_compound(uc, ulon, ulat)
  

codes/icosagcm/trunk/src/icosa_init.f90

@@ -58 +58 @@
 
     CALL init_diagflux
+    CALL zero_du_phys
     CALL timeloop
     CALL switch_omp_no_distrib_level

codes/icosagcm/trunk/src/physics/physics.f90

@@ -1 +1 @@
 MODULE physics_mod
-
+  USE icosa
   USE field_mod
-
+  USE physics_interface_mod
+  USE omp_para
+  IMPLICIT NONE
   PRIVATE
 
@@ -14 +16 @@
   TYPE(t_field),POINTER,SAVE :: f_p(:), f_pk(:)
   TYPE(t_field),POINTER,SAVE :: f_temp(:)
+  TYPE(t_field),POINTER,SAVE :: f_du_phys(:)
 
   CHARACTER(LEN=255),SAVE :: physics_type
 !$OMP THREADPRIVATE(physics_type)
 
-  PUBLIC :: physics, init_physics
+  PUBLIC :: physics, init_physics, zero_du_phys
 
 CONTAINS
@@ -25 +28 @@
     USE mpipara
     USE etat0_mod
-    USE icosa
-    USE physics_interface_mod
     USE physics_dcmip_mod, ONLY : init_physics_dcmip=>init_physics
     USE physics_dcmip2016_mod, ONLY : init_physics_dcmip2016=>init_physics
@@ -32 +33 @@
     USE physics_lmdz_generic_mod, ONLY : init_physics_lmdz_generic=>init_physics
     USE physics_external_mod, ONLY : init_physics_external=>init_physics
-    IMPLICIT NONE
 
     physics_inout%dt_phys = dt*itau_physics
@@ -84 +84 @@
     END SELECT
 
+    CALL allocate_field(f_du_phys,field_u,type_real,llm, name='du_phys')
+
     IF(is_mpi_root) PRINT *, 'phys_type = ',phys_type
   END SUBROUTINE init_physics
 
+  SUBROUTINE zero_du_phys()
+    REAL(rstd), DIMENSION(:,:), POINTER :: du
+    INTEGER :: ind
+    DO ind=1,ndomain
+       IF (.NOT. assigned_domain(ind)) CYCLE
+       CALL swap_dimensions(ind)
+       CALL swap_geometry(ind)
+       du=f_du_phys(ind)
+       du(:,ll_begin:ll_end) = 0.
+    END DO
+  END SUBROUTINE zero_du_phys
+
+  SUBROUTINE add_du_phys(coef, f_u)
+    REAL(rstd), INTENT(IN) :: coef  ! -1 before physics, +1 after physics
+    TYPE(t_field),POINTER :: f_u(:) ! velocity field before/after call to physics
+    REAL(rstd), DIMENSION(:,:), POINTER :: u, du
+    INTEGER :: ind
+    DO ind=1,ndomain
+       IF (.NOT. assigned_domain(ind)) CYCLE
+       CALL swap_dimensions(ind)
+       CALL swap_geometry(ind)
+       du=f_du_phys(ind)
+       u=f_u(ind)
+       du(:,ll_begin:ll_end) = du(:,ll_begin:ll_end) + coef*u(:,ll_begin:ll_end)
+    END DO
+  END SUBROUTINE add_du_phys
+
   SUBROUTINE physics(it,f_phis, f_ps, f_theta_rhodz, f_ue, f_wflux, f_q)
-    USE icosa
-    USE physics_interface_mod
     USE physics_lmdz_generic_mod, ONLY : physics_lmdz_generic => physics
     USE physics_external_mod, ONLY : physics_external => physics
@@ -96 +123 @@
     USE etat0_heldsz_mod
     USE etat0_venus_mod, ONLY : phys_venus => physics
-    IMPLICIT NONE
     INTEGER, INTENT(IN)   :: it
     TYPE(t_field),POINTER :: f_phis(:)
@@ -109 +135 @@
     TYPE(t_physics_inout) :: args
 
-    IF(MOD(it,itau_physics)==0) THEN
-    
+    IF(MOD(it,itau_physics)==0 .AND. phys_type/=phys_none) THEN
+
+       ! as a result of the the two calls to add_du_phys,
+       ! du_phys increases by u(after physics) - u (before physics)
+       CALL add_du_phys(-1., f_ue)
+
        SELECT CASE(phys_type)
-       CASE (phys_none)
-          ! No physics, do nothing
        CASE(phys_HS94)
           CALL held_suarez(f_ps,f_theta_rhodz,f_ue) 
@@ -129 +157 @@
        CALL transfert_request(f_ue,req_e0_vect)
        CALL transfert_request(f_q,req_i0)
+
+       CALL add_du_phys(1., f_ue)
     END IF
 
     IF (mod(it,itau_out)==0 ) THEN
+       CALL write_physics_tendencies
+       CALL zero_du_phys
        SELECT CASE(phys_type)
        CASE (phys_DCMIP)
@@ -142 +174 @@
   END SUBROUTINE physics
 
+  SUBROUTINE write_physics_tendencies
+    USE observable_mod, ONLY : f_buf_ulon, f_buf_ulat
+    USE wind_mod
+    USE output_field_mod
+    CALL transfert_request(f_du_phys,req_e1_vect)
+    CALL un2ulonlat(f_du_phys, f_buf_ulon, f_buf_ulat, (1./(dt*itau_out)))
+    CALL output_field("dulon_phys",f_buf_ulon)
+    CALL output_field("dulat_phys",f_buf_ulat)
+  END SUBROUTINE write_physics_tendencies
+    
   SUBROUTINE physics_column(it, f_phis, f_ps, f_theta_rhodz, f_ue, f_q)
-    USE icosa
-    USE physics_interface_mod
     USE physics_dcmip_mod, ONLY : full_physics_dcmip => full_physics
     USE physics_dcmip2016_mod, ONLY : full_physics_dcmip2016 => full_physics
     USE theta2theta_rhodz_mod
     USE mpipara
-    USE omp_para
     USE checksum_mod
-    IMPLICIT NONE
     TYPE(t_field),POINTER :: f_phis(:)
     TYPE(t_field),POINTER :: f_ps(:)
@@ -229 +267 @@
 
   SUBROUTINE pack_physics(info, phis, ps, temp, ue, q, p, pk, ulon, ulat )
-    USE icosa
     USE wind_mod
     USE pression_mod
     USE theta2theta_rhodz_mod
-    USE physics_interface_mod
     USE exner_mod
-    USE omp_para
-    IMPLICIT NONE
     TYPE(t_pack_info) :: info
     REAL(rstd) :: phis(iim*jjm)
@@ -272 +306 @@
 
   SUBROUTINE unpack_physics(info, ps,temp, q, dulon, dulat)
-    USE icosa
-    USE physics_interface_mod
     USE theta2theta_rhodz_mod
-    USE omp_para
-    IMPLICIT NONE
     TYPE(t_pack_info) :: info
     REAL(rstd) :: ps(iim*jjm)
@@ -303 +333 @@
 
   SUBROUTINE compute_update_velocity(dulon, dulat, ue)
-    USE icosa
-    USE physics_interface_mod
     USE wind_mod
-    USE omp_para
-    IMPLICIT NONE
     REAL(rstd) :: dulon(iim*jjm,llm)
     REAL(rstd) :: dulat(iim*jjm,llm)