Correct confusion of gamma's with gamma's_tra

.github/workflows/fesom2.1.yml

@@ -1,7 +1,7 @@
 
 name: FESOM2 main test
 
-# Controls when the action will run. Triggers the workflow on push or pull request.
+# Controls when the action will run. Triggers the workflow on push or pull request. 
 
 on: [push, pull_request]
 

.github/workflows/fesom2_icepack.yml

@@ -1,7 +1,7 @@
 
 name: FESOM2_icepack
 
-# Controls when the action will run. Triggers the workflow on push or pull request.
+# Controls when the action will run. Triggers the workflow on push or pull request. 
 
 on: [push, pull_request]
 

config/namelist.forcing.ncep2

@@ -0,0 +1,56 @@
+! This is the namelist file for forcing
+
+&forcing_exchange_coeff
+Ce_atm_oce=0.00175 ! exchange coeff. of latent heat over open water
+Ch_atm_oce=0.00175 ! exchange coeff. of sensible heat over open water
+Cd_atm_oce=0.001  ! drag coefficient between atmosphere and water
+Ce_atm_ice=0.00175 ! exchange coeff. of latent heat over ice
+Ch_atm_ice=0.00175 ! exchange coeff. of sensible heat over ice
+Cd_atm_ice=0.0012  ! drag coefficient between atmosphere and ice 
+Swind     =0.0     ! parameterization for coupled current feedback
+/
+
+&forcing_bulk
+AOMIP_drag_coeff=.false.
+ncar_bulk_formulae=.true.
+ncar_bulk_z_wind=10.0 ! height at which wind forcing is located (CORE, JRA-do: 10m, JRA, NCEP:2m)
+ncar_bulk_z_tair=2.0 ! height at which temp forcing is located (CORE, JRA-do: 10m, JRA, NCEP:2m)
+ncar_bulk_z_shum=2.0 ! height at which humi forcing is located (CORE, JRA-do: 10m, JRA, NCEP:2m)
+/
+
+&land_ice
+use_landice_water=.false.
+landice_start_mon=5
+landice_end_mon=10
+/
+
+&nam_sbc
+   nm_xwind_file = '/work/ollie/clidyn/forcing/NCEP2/uwnd.10m.gauss.'   ! name of file with winds, if nm_sbc=2
+   nm_ywind_file = '/work/ollie/clidyn/forcing/NCEP2/vwnd.10m.gauss.'   ! name of file with winds, if nm_sbc=2
+   nm_humi_file  = '/work/ollie/clidyn/forcing/NCEP2/shum.2m.gauss.'    ! name of file with 2m specific humidity
+   nm_qsr_file   = '/work/ollie/clidyn/forcing/NCEP2/dswrf.sfc.gauss.'  ! name of file with solar heat
+   nm_qlw_file   = '/work/ollie/clidyn/forcing/NCEP2/dlwrf.sfc.gauss.'  ! name of file with Long wave
+   nm_tair_file  = '/work/ollie/clidyn/forcing/NCEP2/air.2m.gauss.'     ! name of file with 2m air temperature
+   nm_prec_file  = '/work/ollie/clidyn/forcing/NCEP2/prate.sfc.gauss.'  ! name of file with rain fall
+   nm_snow_file  = ''                                                   ! name of file with snow fall
+   nm_mslp_file  = ''                                                   ! air_pressure_at_sea_level
+   nm_xwind_var  = 'uwnd'   	! name of variable in file with wind
+   nm_ywind_var  = 'vwnd'   	! name of variable in file with wind
+   nm_humi_var   = 'shum'   	! name of variable in file with humidity
+   nm_qsr_var    = 'dswrf'  	! name of variable in file with solar heat
+   nm_qlw_var    = 'dlwrf'  	! name of variable in file with Long wave
+   nm_tair_var   = 'air'    	! name of variable in file with 2m air temperature
+   nm_prec_var   = 'prate'  	! name of variable in file with total precipitation
+   nm_snow_var   = ''       	! name of variable in file with total precipitation
+   nm_mslp_var   = ''       	! name of variable in file with air_pressure_at_sea_level
+   nm_nc_iyear   = 1800
+   nm_nc_imm     = 1            ! initial month of time axis in netCDF
+   nm_nc_idd     = 1            ! initial day of time axis in netCDF
+   nm_nc_freq    = 24           ! data points per day (i.e. 86400 if the time axis is in seconds)
+   nm_nc_tmid    = 0            ! 1 if the time stamps are given at the mid points of the netcdf file, 0 otherwise (i.e. 1 in CORE1, CORE2; 0 in JRA55)
+   l_xwind=.true., l_ywind=.true., l_humi=.true., l_qsr=.true., l_qlw=.true., l_tair=.true., l_prec=.true., l_mslp=.false., l_cloud=.false., l_snow=.false.
+   nm_runoff_file     ='/work/ollie/clidyn/forcing/JRA55-do-v1.4.0/CORE2_runoff.nc'
+   runoff_data_source ='CORE2'  !Dai09, CORE2, JRA55
+   nm_sss_data_file   ='/work/ollie/clidyn/forcing/JRA55-do-v1.4.0/PHC2_salx.nc'
+   sss_data_source    ='CORE2'
+/

env/bsc/shell

@@ -1,5 +1,5 @@
 
-module load intel/2017.4 impi/2017.4 mkl/2017.4 bsc/1.0 netcdf/4.2
+module load intel/2020.1 impi/2018.4 mkl/2020.1 bsc/1.0 netcdf/4.2
 
 export FC=mpiifort CC=mpiicc CXX=mpiicpc
 

setups/test_pi/setup.yml

@@ -61,12 +61,12 @@ namelist.io:
                 prec: 8
 
 fcheck:
-    a_ice: 0.26911274194532003
-    salt: 23.9440531023692
-    temp: 1.7017743034836539
-    sst: 8.532529081624512
-    u: -0.0014065854610620704
-    v: 0.00014195144238082126
+    a_ice: 0.2691276598479261
+    salt: 23.944024679303666
+    temp: 1.701768750033021
+    sst: 8.531528640978305
+    u: -0.0014072137861434184
+    v: 0.00014184602459601167
 
 
 

setups/test_pi_icepack/setup.yml

@@ -73,13 +73,13 @@ namelist.io:
                 prec: 8
 
 fcheck:
-    a_ice: 0.3059942958760058
-    salt: 23.866224273520945
-    temp: 1.7172059436119271
-    sst: 8.725966058658427
-    u: -0.0014448488412238854
-    v: 0.00018596541127645607
-    aicen: 0.061198859175201174
+    a_ice: 0.30599570824298994
+    salt: 23.866195774787034
+    temp: 1.717206693389919
+    sst: 8.725991935766256
+    u: -0.0014448974204450153
+    v: 0.00018600030457097512
+    aicen: 0.06119914164859799
 
 
 

src/CMakeLists.txt

@@ -1,5 +1,4 @@
 cmake_minimum_required(VERSION 3.4)
-set(CMAKE_OSX_DEPLOYMENT_TARGET "10.9")
 
 project(fesom C Fortran)
 
@@ -68,15 +67,15 @@ if(${VERBOSE})
 endif()
 # CMAKE_Fortran_COMPILER_ID will also work if a wrapper is being used (e.g. mpif90 wraps ifort -> compiler id is Intel)
 if(${CMAKE_Fortran_COMPILER_ID} STREQUAL  Intel )
-	target_compile_options(${PROJECT_NAME} PRIVATE -r8 -i4 -fp-model precise -no-prec-div -no-prec-sqrt -fast-transcendentals -xHost -ip -init=zero -no-wrap-margin)
-#	target_compile_options(${PROJECT_NAME} PRIVATE -r8 -i4 -fp-model precise -no-prec-div -no-prec-sqrt -fast-transcendentals -xHost -ip -g -traceback -check all,noarg_temp_created,bounds,uninit ) #-ftrapuv ) #-init=zero)
+	target_compile_options(${PROJECT_NAME} PRIVATE -r8 -i4 -fp-model precise -no-prec-div -no-prec-sqrt -fimf-use-svml -xHost -ip -init=zero -no-wrap-margin)
+#	target_compile_options(${PROJECT_NAME} PRIVATE -r8 -i4 -fp-model precise -no-prec-div -no-prec-sqrt -fimf-use-svml -xHost -ip -g -traceback -check all,noarg_temp_created,bounds,uninit ) #-ftrapuv ) #-init=zero)
 elseif(${CMAKE_Fortran_COMPILER_ID} STREQUAL  GNU )
 	target_compile_options(${PROJECT_NAME} PRIVATE -O3 -finit-local-zero  -finline-functions -march=native -fimplicit-none  -fdefault-real-8 -ffree-line-length-none)
 	if(CMAKE_Fortran_COMPILER_VERSION VERSION_GREATER_EQUAL 10 )
 		target_compile_options(${PROJECT_NAME} PRIVATE -fallow-argument-mismatch) # gfortran v10 is strict about erroneous API calls: "Rank mismatch between actual argument at (1) and actual argument at (2) (scalar and rank-1)"
 	endif()
 elseif(${CMAKE_Fortran_COMPILER_ID} STREQUAL Cray )
-       target_compile_options(${PROJECT_NAME} PRIVATE -c -emf -hbyteswapio -hflex_mp=conservative -hfp1 -hadd_paren -Ounroll0 -hipa0 -r am -s real64)
+       target_compile_options(${PROJECT_NAME} PRIVATE -c -emf -hbyteswapio -hflex_mp=conservative -hfp1 -hadd_paren -Ounroll0 -hipa0 -r am -s real64 -hnoomp)
 endif()
 target_include_directories(${PROJECT_NAME} PRIVATE ${NETCDF_Fortran_INCLUDE_DIRECTORIES} ${OASIS_Fortran_INCLUDE_DIRECTORIES})
 target_include_directories(${PROJECT_NAME} PRIVATE ${MCT_Fortran_INCLUDE_DIRECTORIES} ${MPEU_Fortran_INCLUDE_DIRECTORIES})

src/gen_forcing_couple.F90

@@ -141,113 +141,115 @@ subroutine update_atm_forcing(istep, mesh)
          call cpl_oasis3mct_recv (i,exchange,action)
 	 !if (.not. action) cycle
 	 !Do not apply a correction at first time step!
-	 if (i==1 .and. action .and. istep/=1) call net_rec_from_atm(action)
-         if (i.eq.1) then
-     	     if (.not. action) cycle
-             stress_atmoce_x(:) =  exchange(:)                    ! taux_oce
-	     do_rotate_oce_wind=.true.
-         elseif (i.eq.2) then
-     	     if (.not. action) cycle
-             stress_atmoce_y(:) =  exchange(:)                    ! tauy_oce
-	     do_rotate_oce_wind=.true.
-         elseif (i.eq.3) then
-     	     if (.not. action) cycle	
-             stress_atmice_x(:) =  exchange(:)                    ! taux_ice
-	     do_rotate_ice_wind=.true.
-         elseif (i.eq.4) then
-     	     if (.not. action) cycle	
-             stress_atmice_y(:) =  exchange(:)                    ! tauy_ice
-	     do_rotate_ice_wind=.true.	     
-         elseif (i.eq.5) then
-             if (action) then 
-	        prec_rain(:)    =  exchange(:)	                  ! tot_prec
-		mask=1.
-		call force_flux_consv(prec_rain, mask, i, 0,action, mesh)
-	     end if
-         elseif (i.eq.6) then 
-	     if (action) then
-	        prec_snow(:)    =  exchange(:)                    ! snowfall
-		mask=1.
-		call force_flux_consv(prec_snow, mask,i,1,action, mesh) ! Northern hemisphere
-		call force_flux_consv(prec_snow, mask,i,2,action, mesh) ! Southern Hemisphere
-             end if
+    if (i==1 .and. action .and. istep/=1) call net_rec_from_atm(action)
+        if (i.eq.1) then
+            if (.not. action) cycle
+            stress_atmoce_x(:) =  exchange(:)                    ! taux_oce
+            do_rotate_oce_wind=.true.
+        elseif (i.eq.2) then
+            if (.not. action) cycle
+            stress_atmoce_y(:) =  exchange(:)                    ! tauy_oce
+            do_rotate_oce_wind=.true.
+        elseif (i.eq.3) then
+            if (.not. action) cycle	
+            stress_atmice_x(:) =  exchange(:)                    ! taux_ice
+            do_rotate_ice_wind=.true.
+        elseif (i.eq.4) then
+            if (.not. action) cycle	
+            stress_atmice_y(:) =  exchange(:)                    ! tauy_ice
+            do_rotate_ice_wind=.true.	     
+        elseif (i.eq.5) then
+            if (action) then 
+                prec_rain(:)    =  exchange(:)	                  ! tot_prec
+                mask=1.
+                call force_flux_consv(prec_rain, mask, i, 0,action, mesh)
+            end if
+        elseif (i.eq.6) then 
+            if (action) then
+                prec_snow(:)    =  exchange(:)                    ! snowfall
+                mask=1.
+                call force_flux_consv(prec_snow, mask,i,1,action, mesh) ! Northern hemisphere
+                call force_flux_consv(prec_snow, mask,i,2,action, mesh) ! Southern Hemisphere
+            end if
          elseif (i.eq.7) then
-             if (action) then
-	     evap_no_ifrac(:)     =  exchange(:)        	  ! tot_evap
-	     tmp_evap_no_ifrac(:) =  exchange(:) 		  ! to reset for flux 
-	     							  ! correction
-	     end if
- 	     mask=1.-a_ice
-	     evap_no_ifrac(:)     =  tmp_evap_no_ifrac(:)
-	     call force_flux_consv(evap_no_ifrac,mask,i,0,action, mesh)
-         elseif (i.eq.8) then
-             if (action) then
-	     sublimation(:)       =  exchange(:)        	  ! tot_subl
-	     tmp_sublimation(:)   =  exchange(:) 		  ! to reset for flux 
-	     							  ! correction
-	     end if
-	     mask=a_ice 
-             sublimation(:)       =  tmp_sublimation(:)
-	     call force_flux_consv(sublimation,mask,i,1,action, mesh) ! Northern hemisphere
-	     call force_flux_consv(sublimation,mask,i,2,action, mesh) ! Southern Hemisphere
-         elseif (i.eq.9) then
-             if (action) then
-	     oce_heat_flux(:)     =  exchange(:)        	  ! heat_oce
-	     tmp_oce_heat_flux(:) =  exchange(:) 		  ! to reset for flux 
-	     							  ! correction
-	     end if
-	     mask=1.-a_ice
-	     oce_heat_flux(:)     =  tmp_oce_heat_flux(:)
-	     call force_flux_consv(oce_heat_flux, mask, i, 0,action, mesh)
-         elseif (i.eq.10) then
-             if (action) then
-	     ice_heat_flux(:)     =  exchange(:)        	  ! heat_ice
-	     tmp_ice_heat_flux(:) =  exchange(:) 		  ! to reset for flux 
-	     							  ! correction
-	     end if
-	     mask=a_ice
-	     ice_heat_flux(:)     =  tmp_ice_heat_flux(:)
-	     call force_flux_consv(ice_heat_flux, mask, i, 1,action, mesh) ! Northern hemisphere
-	     call force_flux_consv(ice_heat_flux, mask, i, 2,action, mesh) ! Southern Hemisphere	     
-         elseif (i.eq.11) then
-             if (action) then
-	     shortwave(:)         =  exchange(:)		  ! heat_swr
-	     tmp_shortwave(:)     =  exchange(:) 		  ! to reset for flux 
-	     							  ! correction
-	     end if
-	     mask=1.-a_ice
-	     shortwave(:)   =  tmp_shortwave(:)
-	     call force_flux_consv(shortwave, mask, i, 0,action, mesh)
-         elseif (i.eq.12) then
-             if (action) then
-	     runoff(:)            =  exchange(:)        ! AWI-CM2: runoff, AWI-CM3: runoff + excess snow on glaciers
-    	     mask=1.
-	     call force_flux_consv(runoff, mask, i, 0,action, mesh)
-             end if
+            if (action) then
+                evap_no_ifrac(:)     =  exchange(:)        	  ! tot_evap
+                tmp_evap_no_ifrac(:) =  exchange(:) 		  ! to reset for flux 
+                                                              ! correction
+            end if
+            mask=1.-a_ice
+            evap_no_ifrac(:)     =  tmp_evap_no_ifrac(:)
+            call force_flux_consv(evap_no_ifrac,mask,i,0,action, mesh)
+        elseif (i.eq.8) then
+            if (action) then
+                sublimation(:)       =  exchange(:)        	  ! tot_subl
+                tmp_sublimation(:)   =  exchange(:) 		  ! to reset for flux 
+                                                              ! correction
+            end if
+            mask=a_ice 
+            sublimation(:)       =  tmp_sublimation(:)
+            call force_flux_consv(sublimation,mask,i,1,action, mesh) ! Northern hemisphere
+            call force_flux_consv(sublimation,mask,i,2,action, mesh) ! Southern Hemisphere
+        elseif (i.eq.9) then
+            if (action) then
+                oce_heat_flux(:)     =  exchange(:)        	  ! heat_oce
+                tmp_oce_heat_flux(:) =  exchange(:) 		  ! to reset for flux 
+                                                              ! correction
+            end if
+            mask=1.-a_ice
+            oce_heat_flux(:)     =  tmp_oce_heat_flux(:)
+            call force_flux_consv(oce_heat_flux, mask, i, 0,action, mesh)
+        elseif (i.eq.10) then
+            if (action) then
+                ice_heat_flux(:)     =  exchange(:)        	  ! heat_ice
+                tmp_ice_heat_flux(:) =  exchange(:) 		  ! to reset for flux 
+                                                              ! correction
+            end if
+            mask=a_ice
+            ice_heat_flux(:)     =  tmp_ice_heat_flux(:)
+            call force_flux_consv(ice_heat_flux, mask, i, 1,action, mesh) ! Northern hemisphere
+            call force_flux_consv(ice_heat_flux, mask, i, 2,action, mesh) ! Southern Hemisphere	     
+        elseif (i.eq.11) then
+            if (action) then
+                shortwave(:)         =  exchange(:)		        ! heat_swr
+                tmp_shortwave(:)     =  exchange(:) 		    ! to reset for flux 
+                                                                ! correction
+            end if
+            mask=1.-a_ice
+            shortwave(:)   =  tmp_shortwave(:)
+            call force_flux_consv(shortwave, mask, i, 0,action, mesh)
+        elseif (i.eq.12) then
+            if (action) then
+                runoff(:)            =  exchange(:)        ! AWI-CM2: runoff, AWI-CM3: runoff + excess snow on glaciers
+                mask=1.
+                call force_flux_consv(runoff, mask, i, 0,action, mesh)
+            end if
 #if defined (__oifs)
+
          elseif (i.eq.13) then
              if (action) then
-	     enthalpyoffuse(:)            =  exchange(:)        ! enthalpy of fusion via solid water discharge from glaciers
-    	     mask=1.
-	     call force_flux_consv(enthalpyoffuse, mask, i, 0,action, mesh)
+	              enthalpyoffuse(:)            =  exchange(:)        ! enthalpy of fusion via solid water discharge from glaciers
+    	          mask=1.
+	              call force_flux_consv(enthalpyoffuse, mask, i, 0,action, mesh)
              end if
-	 end if  
-#endif	  
+#endif	 
+         end if
+
 #ifdef VERBOSE
 	  if (mype==0) then
 		write(*,*) 'FESOM RECV: flux ', i, ', max val: ', maxval(exchange)
 	  end if
 #endif
       end do
 
-      if ((do_rotate_oce_wind .AND. do_rotate_ice_wind) .AND. rotated_grid) then
-         do n=1, myDim_nod2D+eDim_nod2D
-	    call vector_g2r(stress_atmoce_x(n), stress_atmoce_y(n), coord_nod2D(1, n), coord_nod2D(2, n), 0)
-	    call vector_g2r(stress_atmice_x(n), stress_atmice_y(n), coord_nod2D(1, n), coord_nod2D(2, n), 0)
-	 end do
-	 do_rotate_oce_wind=.false.
-         do_rotate_ice_wind=.false.
-      end if
+    if ((do_rotate_oce_wind .AND. do_rotate_ice_wind) .AND. rotated_grid) then
+        do n=1, myDim_nod2D+eDim_nod2D
+            call vector_g2r(stress_atmoce_x(n), stress_atmoce_y(n), coord_nod2D(1, n), coord_nod2D(2, n), 0)
+            call vector_g2r(stress_atmice_x(n), stress_atmice_y(n), coord_nod2D(1, n), coord_nod2D(2, n), 0)
+        end do
+        do_rotate_oce_wind=.false.
+        do_rotate_ice_wind=.false.
+    end if
 #else
   call sbc_do(mesh)
   u_wind    = atmdata(i_xwind,:)
@@ -258,20 +260,25 @@ subroutine update_atm_forcing(istep, mesh)
   Tair      = atmdata(i_tair ,:)-273.15_WP
   prec_rain = atmdata(i_prec ,:)/1000._WP
   prec_snow = atmdata(i_snow ,:)/1000._WP
-  press_air = atmdata(i_mslp ,:) ! unit should be Pa
+  if (l_mslp) then
+     press_air = atmdata(i_mslp ,:) ! unit should be Pa
+  end if
 
 
   if (use_cavity) then 
     do i=1,myDim_nod2d+eDim_nod2d
         if (ulevels_nod2d(i)>1) then
-            u_wind(i)=0.0_WP
-            v_wind(i)=0.0_WP
-            shum(i)=0.0_WP
-            longwave(i)=0.0_WP
-            Tair(i)=0.0_WP
-            prec_rain(i)=0.0_WP
-            prec_snow(i)=0.0_WP
-            press_air(i)=0.0_WP            
+            u_wind(i)   = 0.0_WP
+            v_wind(i)   = 0.0_WP
+            shum(i)     = 0.0_WP
+            longwave(i) = 0.0_WP
+            Tair(i)     = 0.0_WP
+            prec_rain(i)= 0.0_WP
+            prec_snow(i)= 0.0_WP
+            if (l_mslp) then
+               press_air(i)= 0.0_WP 
+            end if
+            runoff(i)   = 0.0_WP
         end if 
     end do
   endif