AWICM3 changes to coupling, restart file writing and DMOC

cmake/FindOASIS.cmake

@@ -7,5 +7,5 @@ find_library(MCT_Fortran_LIBRARIES mct HINTS ${TOPLEVEL_DIR}/../oasis/build/lib/
 find_path(MPEU_Fortran_INCLUDE_DIRECTORIES m_mpout.mod HINTS ${TOPLEVEL_DIR}/../oasis/build/lib/psmile/mct/mpeu)
 find_library(MPEU_Fortran_LIBRARIES mpeu HINTS ${TOPLEVEL_DIR}/../oasis/build/lib/psmile/mct/mpeu)
 
-find_path(SCRIP_Fortran_INCLUDE_DIRECTORIES remap_bicubic.mod HINTS ${TOPLEVEL_DIR}/../oasis/build/lib/psmile/scrip)
+find_path(SCRIP_Fortran_INCLUDE_DIRECTORIES remap_bicubic_reduced.mod HINTS ${TOPLEVEL_DIR}/../oasis/build/lib/psmile/scrip)
 find_library(SCRIP_Fortran_LIBRARIES scrip HINTS ${TOPLEVEL_DIR}/../oasis/build/lib/psmile/scrip)

src/cpl_driver.F90

@@ -24,10 +24,11 @@ module cpl_driver
 
 #if defined (__oifs)
   integer, parameter         :: nsend = 5
+  integer, parameter         :: nrecv = 12
 #else
   integer, parameter         :: nsend = 4
-#endif
   integer, parameter         :: nrecv = 12
+#endif
 
   integer, dimension(nsend)  :: send_id
   integer, dimension(nrecv)  :: recv_id
@@ -222,7 +223,7 @@ subroutine cpl_oasis3mct_define_unstr(mesh)
 
     real(kind=WP), allocatable :: all_x_coords(:, :)     ! longitude coordinates
     real(kind=WP), allocatable :: all_y_coords(:, :)     ! latitude  coordinates
-    real(kind=WP), allocatable :: all_elem_area(:,:)    
+    real(kind=WP), allocatable :: all_area(:,:)    
 
 #include "associate_mesh.h"
 
@@ -304,11 +305,11 @@ subroutine cpl_oasis3mct_define_unstr(mesh)
     if (mype .eq. localroot) then
       ALLOCATE(all_x_coords(number_of_all_points, 1))
       ALLOCATE(all_y_coords(number_of_all_points, 1))
-      ALLOCATE(all_elem_area(number_of_all_points, 1))
+      ALLOCATE(all_area(number_of_all_points, 1))
     else 
       ALLOCATE(all_x_coords(1, 1))
       ALLOCATE(all_y_coords(1, 1))
-      ALLOCATE(all_elem_area(1, 1))
+      ALLOCATE(all_area(1, 1))
     endif
 
     displs_from_all_pes(1) = 0
@@ -331,7 +332,7 @@ subroutine cpl_oasis3mct_define_unstr(mesh)
     if (mype .eq. 0) then 
       print *, 'FESOM before 3rd GatherV'
     endif
-    CALL MPI_GATHERV(elem_area, my_number_of_points, MPI_DOUBLE_PRECISION, all_elem_area,  &
+    CALL MPI_GATHERV(area(1,:), my_number_of_points, MPI_DOUBLE_PRECISION, all_area,  &
                     counts_from_all_pes, displs_from_all_pes, MPI_DOUBLE_PRECISION, localroot, MPI_COMM_FESOM, ierror)
 
     if (mype .eq. 0) then 
@@ -358,7 +359,7 @@ subroutine cpl_oasis3mct_define_unstr(mesh)
           DEALLOCATE(unstr_mask)
 
           print *, 'FESOM before write area'
-          CALL oasis_write_area(grid_name, number_of_all_points, 1, all_elem_area)
+          CALL oasis_write_area(grid_name, number_of_all_points, 1, all_area)
 
        end if
       print *, 'FESOM before terminate_grids_writing'
@@ -377,12 +378,17 @@ subroutine cpl_oasis3mct_define_unstr(mesh)
 ! ... Define symbolic names for the transient fields send by the ocean
 !     These must be identical to the names specified in the SMIOC file.
 !
+#if defined (__oifs)
     cpl_send( 1)='sst_feom' ! 1. sea surface temperature [K]       ->
+    cpl_send( 2)='sie_feom' ! 2. sea ice extent [%-100]            ->
+    cpl_send( 3)='snt_feom' ! 3. snow thickness [m]                ->
+    cpl_send( 4)='ist_feom' ! 4. sea ice surface temperature [K]   ->
+    cpl_send( 5)='sia_feom' ! 5. sea ice albedo [%-100]            ->
+#else
+    cpl_send( 1)='sst_feom' ! 1. sea surface temperature [°C]      ->
     cpl_send( 2)='sit_feom' ! 2. sea ice thickness [m]             ->
     cpl_send( 3)='sie_feom' ! 3. sea ice extent [%-100]            ->
     cpl_send( 4)='snt_feom' ! 4. snow thickness [m]                ->
-#if defined (__oifs)
-    cpl_send( 5)='sia_feom' ! 5. sea ice albedo [%-100]            ->
 #endif
 
 
@@ -391,6 +397,7 @@ subroutine cpl_oasis3mct_define_unstr(mesh)
 ! ...  Define symbolic names for transient fields received by the ocean.
 !      These must be identical to the names specified in the SMIOC file.
 !
+#if defined (__oifs)
     cpl_recv(1)  = 'taux_oce'
     cpl_recv(2)  = 'tauy_oce'
     cpl_recv(3)  = 'taux_ico'
@@ -403,6 +410,20 @@ subroutine cpl_oasis3mct_define_unstr(mesh)
     cpl_recv(10) = 'heat_ico'
     cpl_recv(11) = 'heat_swo'    
     cpl_recv(12) = 'hydr_oce'
+#else
+    cpl_recv(1)  = 'taux_oce'
+    cpl_recv(2)  = 'tauy_oce'
+    cpl_recv(3)  = 'taux_ico'
+    cpl_recv(4)  = 'tauy_ico'    
+    cpl_recv(5)  = 'prec_oce'
+    cpl_recv(6)  = 'snow_oce'    
+    cpl_recv(7)  = 'evap_oce'
+    cpl_recv(8)  = 'subl_oce'
+    cpl_recv(9)  = 'heat_oce'
+    cpl_recv(10) = 'heat_ico'
+    cpl_recv(11) = 'heat_swo'    
+    cpl_recv(12) = 'hydr_oce'
+#endif
 
     if (mype .eq. 0) then 
        print *, 'FESOM after declaring the transient variables'

src/gen_forcing_couple.F90

@@ -58,25 +58,35 @@ subroutine update_atm_forcing(istep, mesh)
      do i=1,nsend
          exchange  =0.
          if (i.eq.1) then
+#if defined (__oifs) 
+            ! AWI-CM3 outgoing state vectors
             do n=1,myDim_nod2D+eDim_nod2D
-#if defined (__oifs)
             exchange(n)=tr_arr(1, n, 1)+tmelt	                    ! sea surface temperature [K]
+            end do
+            elseif (i.eq.2) then
+            exchange(:) = a_ice(:)                                  ! ice concentation [%]
+            elseif (i.eq.3) then
+            exchange(:) = m_snow(:)                                 ! snow thickness
+            elseif (i.eq.4) then
+            exchange(:) = ice_temp(:)                               ! ice surface temperature
+            elseif (i.eq.5) then
+            exchange(:) = ice_alb(:)                                ! ice albedo
+            else	    
+            print *, 'not installed yet or error in cpl_oasis3mct_send', mype
 #else
-            exchange(n)=tr_arr(1, n, 1)		                    ! sea surface temperature [°C]
-#endif
+            ! AWI-CM2 outgoing state vectors
+            do n=1,myDim_nod2D+eDim_nod2D
+            exchange(n)=tr_arr(1, n, 1)                             ! sea surface temperature [°C]
             end do
             elseif (i.eq.2) then
             exchange(:) = m_ice(:)                                  ! ice thickness [m]
             elseif (i.eq.3) then
             exchange(:) = a_ice(:)                                  ! ice concentation [%]
             elseif (i.eq.4) then
             exchange(:) = m_snow(:)                                 ! snow thickness
-#if defined (__oifs)
-            elseif (i.eq.5) then
-            exchange(:) = ice_alb(:)                                ! ice albedo
-#endif
             else	    
             print *, 'not installed yet or error in cpl_oasis3mct_send', mype
+#endif
          endif
          call cpl_oasis3mct_send(i, exchange, action)
       enddo
@@ -170,11 +180,11 @@ subroutine update_atm_forcing(istep, mesh)
 	     call force_flux_consv(shortwave, mask, i, 0,action, mesh)
          elseif (i.eq.12) then
              if (action) then
-	     runoff(:)                   =  exchange(:)        ! runoff + calving
+	     runoff(:)            =  exchange(:)        ! runoff + calving
     	     mask=1.
 	     call force_flux_consv(runoff, mask, i, 0,action, mesh)
              end if
-	  end if  	  
+	 end if  	  
 #ifdef VERBOSE
 	  if (mype==0) then
 		write(*,*) 'FESOM RECV: flux ', i, ', max val: ', maxval(exchange)

src/gen_forcing_init.F90

@@ -54,14 +54,12 @@ subroutine forcing_array_setup(mesh)
   allocate(runoff(n2), evaporation(n2),ice_sublimation(n2))
 
 #if defined (__oasis)
-  allocate(sublimation(n2), evap_no_ifrac(n2))
   allocate(tmp_sublimation(n2),tmp_evap_no_ifrac(n2), tmp_shortwave(n2))
+  allocate(sublimation(n2),evap_no_ifrac(n2))
   allocate(atm_net_fluxes_north(nrecv), atm_net_fluxes_south(nrecv))
   allocate(oce_net_fluxes_north(nrecv), oce_net_fluxes_south(nrecv))
   allocate(flux_correction_north(nrecv), flux_correction_south(nrecv))
   allocate(flux_correction_total(nrecv))
-  sublimation=0.0_WP
-  evap_no_ifrac=0.0_WP
   tmp_sublimation = 0.0_WP
   tmp_evap_no_ifrac = 0.0_WP
   tmp_shortwave = 0.0_WP
@@ -72,6 +70,8 @@ subroutine forcing_array_setup(mesh)
   flux_correction_north=0.0_WP
   flux_correction_south=0.0_WP
   flux_correction_total=0.0_WP  
+  evap_no_ifrac=0.0_WP
+  sublimation=0.0_WP
 #endif