Morrison-Gettelman microphysics

CMakeLists.txt

@@ -129,15 +129,31 @@ elseif (${CMAKE_Fortran_COMPILER_ID} STREQUAL "Intel")
                                 ./physics/cu_gf_sh.F90
                                 PROPERTIES COMPILE_FLAGS "${CMAKE_Fortran_FLAGS} -r8 -ftz")
     # Replace -xHost or -xCORE-AVX2 with -xCORE-AVX-I for certain files (following FV3/gfsphysics/makefile)
+    # for bit-for-bit reproducibility with non-CCPP builds. These may go in the future once the CCPP solution
+    # is fully accepted.
     set(CMAKE_Fortran_FLAGS_LOPT ${CMAKE_Fortran_FLAGS})
-    string(REPLACE "xHost" "xCORE-AVX-I"
+    string(REPLACE "-xHOST" "-xCORE-AVX-I"
            CMAKE_Fortran_FLAGS_LOPT
            "${CMAKE_Fortran_FLAGS_LOPT}")
-    string(REPLACE "xCORE-AVX2" "xCORE-AVX-I"
+    string(REPLACE "-xCORE-AVX2" "-xCORE-AVX-I"
            CMAKE_Fortran_FLAGS_LOPT
            "${CMAKE_Fortran_FLAGS_LOPT}")
     SET_SOURCE_FILES_PROPERTIES(./physics/radiation_aerosols.f
                                 PROPERTIES COMPILE_FLAGS "${CMAKE_Fortran_FLAGS_LOPT} -r8 -ftz")
+    # Force consistent results of math calculations for MG microphysics;
+    # in Debug/Bitforbit) mode; without this flag, the results of the
+    # intrinsic gamma function are different for the non-CCPP and CCPP
+    # version (on Theia with Intel 18). Note this is only required with
+    # dynamic CCPP builds (hybrid, standalone), not with static CCPP builds.
+    if (${CMAKE_BUILD_TYPE} MATCHES "Debug")
+      SET_SOURCE_FILES_PROPERTIES(./physics/micro_mg2_0.F90
+                                  ./physics/micro_mg3_0.F90
+                                  PROPERTIES COMPILE_FLAGS "${CMAKE_Fortran_FLAGS} -fimf-arch-consistency=true")
+    elseif (${CMAKE_BUILD_TYPE} MATCHES "Bitforbit")
+      SET_SOURCE_FILES_PROPERTIES(./physics/micro_mg2_0.F90
+                                  ./physics/micro_mg3_0.F90
+                                  PROPERTIES COMPILE_FLAGS "${CMAKE_Fortran_FLAGS} -fimf-arch-consistency=true")
+    endif (${CMAKE_BUILD_TYPE} MATCHES "Debug")
   else (PROJECT STREQUAL "CCPP-FV3")
     SET_SOURCE_FILES_PROPERTIES(./physics/module_bfmicrophysics.f ./physics/rascnvv2.f ./physics/sflx.f ./physics/sfc_diff.f ./physics/sfc_diag.f PROPERTIES COMPILE_FLAGS -r8)
     SET_SOURCE_FILES_PROPERTIES(./physics/module_nst_model.f90 ./physics/calpreciptype.f90 PROPERTIES COMPILE_FLAGS "-r8 -free")

physics/GFS_debug.F90

@@ -6,7 +6,11 @@ module GFS_diagtoscreen
 
       public GFS_diagtoscreen_init, GFS_diagtoscreen_run, GFS_diagtoscreen_finalize
 
-#define PRINT_SUM
+      public print_my_stuff, chksum_int, chksum_real
+
+#define PRINT_CHKSUM
+!#define PRINT_SUM
+
       interface print_var
         module procedure print_logic_0d
         module procedure print_int_0d
@@ -18,10 +22,10 @@ module GFS_diagtoscreen
       end interface
 
       integer, parameter :: ISTART = 1
-      integer, parameter :: IEND = 11
+      integer, parameter :: IEND = 9999999
 
       integer, parameter :: KSTART = 1
-      integer, parameter :: KEND = 11
+      integer, parameter :: KEND = 9999999
 
       contains
 
@@ -86,7 +90,7 @@ subroutine GFS_diagtoscreen_run (Model, Statein, Stateout, Sfcprop, Coupling,
          integer,                    intent(out) :: errflg
 
          !--- local variables
-         integer :: impi, iomp, ierr
+         integer :: impi, iomp, ierr, n
          integer :: mpirank, mpisize, mpicomm
          integer :: omprank, ompsize
 
@@ -115,7 +119,7 @@ subroutine GFS_diagtoscreen_run (Model, Statein, Stateout, Sfcprop, Coupling,
 !$OMP BARRIER
 #endif
 #ifdef MPI
-         call MPI_BARRIER(mpicomm,ierr)
+!         call MPI_BARRIER(mpicomm,ierr)
 #endif
 
          do impi=0,mpisize-1
@@ -230,6 +234,9 @@ subroutine GFS_diagtoscreen_run (Model, Statein, Stateout, Sfcprop, Coupling,
                      call print_var(mpirank,omprank, Tbd%blkno, 'Tbd%phy_fctd'        , Tbd%phy_fctd)
                      call print_var(mpirank,omprank, Tbd%blkno, 'Tbd%phy_f2d'         , Tbd%phy_f2d)
                      call print_var(mpirank,omprank, Tbd%blkno, 'Tbd%phy_f3d'         , Tbd%phy_f3d)
+                     do n=1,size(Tbd%phy_f3d(1,1,:))
+                         call print_var(mpirank,omprank, Tbd%blkno, 'Tbd%phy_f3d_n'   , Tbd%phy_f3d(:,:,n))
+                     end do
                      call print_var(mpirank,omprank, Tbd%blkno, 'Tbd%blkno'           , Tbd%blkno)
                      ! Diag (incomplete)
                      call print_var(mpirank,omprank, Tbd%blkno, 'Diag%topfsw%upfxc',    Diag%topfsw%upfxc)
@@ -410,15 +417,15 @@ subroutine GFS_diagtoscreen_run (Model, Statein, Stateout, Sfcprop, Coupling,
 #endif
              end do
 #ifdef MPI
-             call MPI_BARRIER(mpicomm,ierr)
+!             call MPI_BARRIER(mpicomm,ierr)
 #endif
          end do
 
 #ifdef OPENMP
 !$OMP BARRIER
 #endif
 #ifdef MPI
-         call MPI_BARRIER(mpicomm,ierr)
+!         call MPI_BARRIER(mpicomm,ierr)
 #endif
 
       end subroutine GFS_diagtoscreen_run
@@ -431,7 +438,7 @@ subroutine print_logic_0d(mpirank,omprank,blkno,name,var)
           character(len=*), intent(in) :: name
           logical, intent(in) :: var
 
-          write(0,'(2a,3i4,1x,l)') 'XXX: ', trim(name), mpirank, omprank, blkno, var
+          write(0,'(2a,3i6,1x,l)') 'XXX: ', trim(name), mpirank, omprank, blkno, var
 
       end subroutine print_logic_0d
 
@@ -443,7 +450,7 @@ subroutine print_int_0d(mpirank,omprank,blkno,name,var)
           character(len=*), intent(in) :: name
           integer, intent(in) :: var
 
-          write(0,'(2a,3i4,i15)') 'XXX: ', trim(name), mpirank, omprank, blkno, var
+          write(0,'(2a,3i6,i15)') 'XXX: ', trim(name), mpirank, omprank, blkno, var
 
       end subroutine print_int_0d
 
@@ -460,10 +467,12 @@ subroutine print_int_1d(mpirank,omprank,blkno,name,var)
           integer :: i
 
 #ifdef PRINT_SUM
-          write(0,'(2a,3i4,3i15)') 'XXX: ', trim(name), mpirank, omprank, blkno, sum(var), minval(var), maxval(var)
+          write(0,'(2a,3i6,3i15)') 'XXX: ', trim(name), mpirank, omprank, blkno, sum(var), minval(var), maxval(var)
+#elif defined(PRINT_CHKSUM)
+          write(0,'(2a,3i6,i20,2i15)') 'XXX: ', trim(name), mpirank, omprank, blkno, chksum_int(size(var),var), minval(var), maxval(var)
 #else
           do i=ISTART,min(IEND,size(var(:)))
-              write(0,'(2a,3i4,i6,i15)') 'XXX: ', trim(name), mpirank, omprank, blkno, i, var(i)
+              write(0,'(2a,3i6,i6,i15)') 'XXX: ', trim(name), mpirank, omprank, blkno, i, var(i)
           end do
 #endif
 
@@ -479,7 +488,7 @@ subroutine print_real_0d(mpirank,omprank,blkno,name,var)
           character(len=*), intent(in) :: name
           real(kind_phys), intent(in) :: var
 
-          write(0,'(2a,3i4,e35.25)') 'XXX: ', trim(name), mpirank, omprank, blkno, var
+          write(0,'(2a,3i6,e35.25)') 'XXX: ', trim(name), mpirank, omprank, blkno, var
 
       end subroutine print_real_0d
 
@@ -496,10 +505,12 @@ subroutine print_real_1d(mpirank,omprank,blkno,name,var)
           integer :: i
 
 #ifdef PRINT_SUM
-          write(0,'(2a,3i4,3e35.25)') 'XXX: ', trim(name), mpirank, omprank, blkno, sum(var), minval(var), maxval(var)
+          write(0,'(2a,3i6,3e35.25)') 'XXX: ', trim(name), mpirank, omprank, blkno, sum(var), minval(var), maxval(var)
+#elif defined(PRINT_CHKSUM)
+          write(0,'(2a,3i6,i20,2e35.25)') 'XXX: ', trim(name), mpirank, omprank, blkno, chksum_real(size(var),var), minval(var), maxval(var)
 #else
           do i=ISTART,min(IEND,size(var(:)))
-              write(0,'(2a,3i4,i6,e35.25)') 'XXX: ', trim(name), mpirank, omprank, blkno, i, var(i)
+              write(0,'(2a,3i6,i6,e35.25)') 'XXX: ', trim(name), mpirank, omprank, blkno, i, var(i)
           end do
 #endif
 
@@ -518,11 +529,13 @@ subroutine print_real_2d(mpirank,omprank,blkno,name,var)
           integer :: k, i
 
 #ifdef PRINT_SUM
-          write(0,'(2a,3i4,3e35.25)') 'XXX: ', trim(name), mpirank, omprank, blkno, sum(var), minval(var), maxval(var)
+          write(0,'(2a,3i6,3e35.25)') 'XXX: ', trim(name), mpirank, omprank, blkno, sum(var), minval(var), maxval(var)
+#elif defined(PRINT_CHKSUM)
+          write(0,'(2a,3i6,i20,2e35.25)') 'XXX: ', trim(name), mpirank, omprank, blkno, chksum_real(size(var),reshape(var,(/size(var)/))), minval(var), maxval(var)
 #else
           do i=ISTART,min(IEND,size(var(:,1)))
               do k=KSTART,min(KEND,size(var(1,:)))
-                  write(0,'(2a,3i4,2i6,e35.25)') 'XXX: ', trim(name), mpirank, omprank, blkno, i, k, var(i,k)
+                  write(0,'(2a,3i6,2i6,e35.25)') 'XXX: ', trim(name), mpirank, omprank, blkno, i, k, var(i,k)
               end do
           end do
 #endif
@@ -542,19 +555,97 @@ subroutine print_real_3d(mpirank,omprank,blkno,name,var)
           integer :: k, i, l
 
 #ifdef PRINT_SUM
-          write(0,'(2a,3i4,3e35.25)') 'XXX: ', trim(name), mpirank, omprank, blkno, sum(var), minval(var), maxval(var)
+          write(0,'(2a,3i6,3e35.25)') 'XXX: ', trim(name), mpirank, omprank, blkno, sum(var), minval(var), maxval(var)
+#elif defined(PRINT_CHKSUM)
+          write(0,'(2a,3i6,i20,2e35.25)') 'XXX: ', trim(name), mpirank, omprank, blkno, chksum_real(size(var),reshape(var,(/size(var)/))), minval(var), maxval(var)
 #else
           do i=ISTART,min(IEND,size(var(:,1,1)))
               do k=KSTART,min(KEND,size(var(1,:,1)))
                   do l=1,size(var(1,1,:))
-                      write(0,'(2a,3i4,3i6,e35.25)') 'XXX: ', trim(name), mpirank, omprank, blkno, i, k, l, var(i,k,l)
+                      write(0,'(2a,3i6,3i6,e35.25)') 'XXX: ', trim(name), mpirank, omprank, blkno, i, k, l, var(i,k,l)
                   end do
               end do
           end do
 #endif
 
       end subroutine print_real_3d
 
+      function chksum_int(N, var) result(hash)
+          implicit none
+          integer, intent(in) :: N
+          integer, dimension(1:N), intent(in) :: var
+          integer*8, dimension(1:N) :: int_var
+          integer*8 :: a, b, i, hash
+          integer*8, parameter :: mod_adler=65521
+
+          a=1
+          b=0
+          i=1
+          hash = 0
+          int_var = TRANSFER(var, a, N)
+
+          do i= 1, N
+              a = MOD(a + int_var(i), mod_adler)
+              b = MOD(b+a, mod_adler)
+          end do
+
+          hash = ior(b * 65536, a)
+
+     end function chksum_int
+
+      function chksum_real(N, var) result(hash)
+          use machine,               only: kind_phys
+          implicit none
+          integer, intent(in) :: N
+          real(kind_phys), dimension(1:N), intent(in) :: var
+          integer*8, dimension(1:N) :: int_var
+          integer*8 :: a, b, i, hash
+          integer*8, parameter :: mod_adler=65521
+
+          a=1
+          b=0
+          i=1
+          hash = 0
+          int_var = TRANSFER(var, a, N)
+
+          do i= 1, N
+              a = MOD(a + int_var(i), mod_adler)
+              b = MOD(b+a, mod_adler)
+          end do
+
+          hash = ior(b * 65536, a)
+
+     end function chksum_real
+
+     function print_my_stuff(mpitoprint,omptoprint) result(flag)
+#ifdef MPI
+         use mpi
+#endif
+#ifdef OPENMP
+         use omp_lib
+#endif
+         implicit none
+         integer, intent(in) :: mpitoprint, omptoprint
+         logical :: flag
+         integer :: ompthread, mpirank, ierr
+#ifdef MPI
+         call MPI_COMM_RANK(MPI_COMM_WORLD, mpirank, ierr)
+#else
+         mpirank = 0
+#endif
+#ifdef OPENMP
+         ompthread = OMP_GET_THREAD_NUM()
+#else
+         ompthread = 0
+#endif
+
+         if (mpitoprint==mpirank .and. omptoprint==ompthread) then
+             flag = .true.
+         else
+             flag = .false.
+         end if
+     end function print_my_stuff
+
     end module GFS_diagtoscreen