Changeset 714 for codes/icosagcm/devel/src/physics/physics.f90

Timestamp:

08/03/18 16:53:37 (6 years ago)

Author:

dubos

Message:

devel : backported from trunk commits r607,r648,r649,r667,r668,r669,r706

File:

• codes/icosagcm/devel/src/physics/physics.f90 (modified) (12 diffs)

codes/icosagcm/devel/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)