From 181fe5a41c44bdada5eb6f2162b84e7f36d78324 Mon Sep 17 00:00:00 2001 From: Jan Hegewald Date: Fri, 18 Jun 2021 15:34:53 +0200 Subject: [PATCH 01/37] fix build error on macOS --- src/CMakeLists.txt | 1 - 1 file changed, 1 deletion(-) diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 28811319d..68a7993d5 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -1,5 +1,4 @@ cmake_minimum_required(VERSION 3.4) -set(CMAKE_OSX_DEPLOYMENT_TARGET "10.9") project(fesom C Fortran) From 563c6fb2fb76d7c152895d4663da068a73ae680f Mon Sep 17 00:00:00 2001 From: Paul Gierz Date: Thu, 24 Jun 2021 16:00:22 +0200 Subject: [PATCH 02/37] Update gen_modules_rotate_grid.F90 Better helpful comments, plus a typo --- src/gen_modules_rotate_grid.F90 | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/src/gen_modules_rotate_grid.F90 b/src/gen_modules_rotate_grid.F90 index 1634ee48b..797a8a74d 100755 --- a/src/gen_modules_rotate_grid.F90 +++ b/src/gen_modules_rotate_grid.F90 @@ -1,10 +1,13 @@ ! Routines needed to support displaced poles: ! The new pole position is set with -! alphaEuler, betaEuler and gammaEuler. The degfault values +! alphaEuler, betaEuler and gammaEuler. The default values ! alphaEuler=50. [degree] Euler angles, convention: ! betaEuler=15. [degree] first around z, then around new x, ! gammaEuler=-90. [degree] then around new z. ! +! A helpful animation may be found online here: +! https://en.wikipedia.org/wiki/Euler_angles +! ! The first two define the new pole position ! as phi_p=alphaEuler-90, theta_p=90-betaEuler. ! The third, gammaEuler, is in reality irrelevant and just From 0911d1cacdd4f681ef9b313042979175ec042882 Mon Sep 17 00:00:00 2001 From: Patrick Date: Mon, 28 Jun 2021 11:47:43 +0200 Subject: [PATCH 03/37] addd comments --- src/ice_oce_coupling.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/ice_oce_coupling.F90 b/src/ice_oce_coupling.F90 index 7b9555d02..5ff9c3f57 100755 --- a/src/ice_oce_coupling.F90 +++ b/src/ice_oce_coupling.F90 @@ -309,10 +309,10 @@ subroutine oce_fluxes(mesh) ! Also balance freshwater flux that come from ocean-cavity boundary if (use_cavity) then - if (.not. use_virt_salt) then + if (.not. use_virt_salt) then !zstar, zlevel ! only for full-free surface approach otherwise total ocean volume will drift where (ulevels_nod2d > 1) flux = -water_flux - else + else ! linfs where (ulevels_nod2d > 1) flux = 0.0_WP end if end if From 6fe46834bb7402c18c436bd12ade3d75d5363415 Mon Sep 17 00:00:00 2001 From: patrickscholz Date: Mon, 28 Jun 2021 14:07:11 +0200 Subject: [PATCH 04/37] Fix problem with FESOM2.0 reproducability bias on ollie - FESOM2.0 simulations on AWI ollie HPC were not always reproducible. Runs with identical namelists, init-conditions ... could lead to sligthly different results. - Natalja solved that issue by changing the compiler flag -fast-transcendentals with -fimf-use-svml - with this option loops with exp, sin, cos ... (IEEE function) will be always computed through the vector implementation otherwise the compiler is not fixed how IEEE funtions should be treated, which can be either Vector, Scalar or something else. All three can give slightly different results --- src/CMakeLists.txt | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 68a7993d5..a04db6736 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -67,8 +67,8 @@ 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 ) From db3d62737b0149337accf4978b2c29ebe5b0dd26 Mon Sep 17 00:00:00 2001 From: Patrick Date: Tue, 29 Jun 2021 14:54:18 +0200 Subject: [PATCH 05/37] sort structure of ../src/gen_forcing_couple.F90 a bit --- src/gen_forcing_couple.F90 | 211 +++++++++++++++++++------------------ 1 file changed, 106 insertions(+), 105 deletions(-) diff --git a/src/gen_forcing_couple.F90 b/src/gen_forcing_couple.F90 index 2a39b34b7..19240fc31 100755 --- a/src/gen_forcing_couple.F90 +++ b/src/gen_forcing_couple.F90 @@ -141,97 +141,97 @@ 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) - end if - end if + 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) + end if + end if #endif #ifdef VERBOSE if (mype==0) then @@ -240,14 +240,14 @@ subroutine update_atm_forcing(istep, mesh) #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,:) @@ -264,14 +264,15 @@ subroutine update_atm_forcing(istep, mesh) 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 + press_air(i)= 0.0_WP +! runoff(i) = 0.0_WP end if end do endif From 541ddc860b1ba2beb93197695eda2df5b48e1317 Mon Sep 17 00:00:00 2001 From: Patrick Date: Tue, 29 Jun 2021 15:18:35 +0200 Subject: [PATCH 06/37] set also river runoff zeros where there is cavity since runoff is only added in ice_therm when there is no cavity otherwise balancing might be corrupted --- src/gen_forcing_couple.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/gen_forcing_couple.F90 b/src/gen_forcing_couple.F90 index 19240fc31..d62207cab 100755 --- a/src/gen_forcing_couple.F90 +++ b/src/gen_forcing_couple.F90 @@ -272,7 +272,7 @@ subroutine update_atm_forcing(istep, mesh) prec_rain(i)= 0.0_WP prec_snow(i)= 0.0_WP press_air(i)= 0.0_WP -! runoff(i) = 0.0_WP + runoff(i) = 0.0_WP end if end do endif From b70ba4f5d973c0dbe40954360efe4f3cda0f7a7f Mon Sep 17 00:00:00 2001 From: Patrick Date: Wed, 30 Jun 2021 12:08:46 +0200 Subject: [PATCH 07/37] add some comments --- src/oce_ale.F90 | 17 ++++++++++++----- 1 file changed, 12 insertions(+), 5 deletions(-) diff --git a/src/oce_ale.F90 b/src/oce_ale.F90 index 447929d75..d231cbf61 100644 --- a/src/oce_ale.F90 +++ b/src/oce_ale.F90 @@ -2719,20 +2719,27 @@ subroutine oce_timestep_ale(n, mesh) call compute_hbar_ale(mesh) !___________________________________________________________________________ - ! Current dynamic elevation alpha*hbar(n+1/2)+(1-alpha)*hbar(n-1/2) - ! equation (14) Danlov et.al "the finite volume sea ice ocean model FESOM2 - ! ...if we do it here we don't need to write hbar_old into a restart file... + ! - Current dynamic elevation alpha*hbar(n+1/2)+(1-alpha)*hbar(n-1/2) + ! equation (14) Danlov et.al "the finite volume sea ice ocean model FESOM2 + ! ...if we do it here we don't need to write hbar_old into a restart file... + ! - where(ulevels_nod2D==1) is used here because of the rigid lid + ! approximation under the cavity + ! - at points in the cavity the time derivative term in ssh matrix will be + ! omitted; and (14) will not be applied at cavity points. Additionally, + ! since there is no real elevation, but only surface pressure, there is + ! no layer motion under the cavity. In this case the ice sheet acts as a + ! rigid lid. where(ulevels_nod2D==1) eta_n=alpha*hbar+(1.0_WP-alpha)*hbar_old - ! --> eta_(n) ! call zero_dynamics !DS, zeros several dynamical variables; to be used for testing new implementations! t5=MPI_Wtime() + !___________________________________________________________________________ + ! Do horizontal and vertical scaling of GM/Redi diffusivity if (Fer_GM .or. Redi) then call init_Redi_GM(mesh) end if - !___________________________________________________________________________ ! Implementation of Gent & McWiliams parameterization after R. Ferrari et al., 2010 ! does not belong directly to ALE formalism if (Fer_GM) then From a814b895223885b42fca50f5baec4272d39597d8 Mon Sep 17 00:00:00 2001 From: Patrick Date: Wed, 30 Jun 2021 12:16:52 +0200 Subject: [PATCH 08/37] in case of cavity do freshwater balancing only for the open ocean part since under the cavity there is rigid lid approximation --- src/ice_oce_coupling.F90 | 14 ++++++++++++-- 1 file changed, 12 insertions(+), 2 deletions(-) diff --git a/src/ice_oce_coupling.F90 b/src/ice_oce_coupling.F90 index 5ff9c3f57..f0efc5843 100755 --- a/src/ice_oce_coupling.F90 +++ b/src/ice_oce_coupling.F90 @@ -316,14 +316,24 @@ subroutine oce_fluxes(mesh) where (ulevels_nod2d > 1) flux = 0.0_WP end if end if - + + ! compute total global net freshwater flux into the ocean call integrate_nod(flux, net, mesh) + + !___________________________________________________________________________ ! here the + sign must be used because we switched up the sign of the ! water_flux with water_flux = -fresh_wa_flux, but evap, prec_... and runoff still ! have there original sign ! if use_cavity=.false. --> ocean_area == ocean_areawithcav !! water_flux=water_flux+net/ocean_area - water_flux=water_flux+net/ocean_areawithcav + if (use_cavity) then + ! due to rigid lid approximation under the cavity we to not add freshwater + ! under the cavity for the freshwater balancing we do this only for the open + ! ocean + where (ulevels_nod2d == 1) water_flux=water_flux+net/ocean_area + else + water_flux=water_flux+net/ocean_area + end if !___________________________________________________________________________ if (use_sw_pene) call cal_shortwave_rad(mesh) From eb1809886d884db80519519288d57ebe18f685bd Mon Sep 17 00:00:00 2001 From: Patrick Date: Wed, 30 Jun 2021 13:56:42 +0200 Subject: [PATCH 09/37] compute integrate parameter also under the cavity --- src/write_step_info.F90 | 1008 +++++++++++++++++++-------------------- 1 file changed, 504 insertions(+), 504 deletions(-) diff --git a/src/write_step_info.F90 b/src/write_step_info.F90 index 9ca952efe..9e0ad0eab 100644 --- a/src/write_step_info.F90 +++ b/src/write_step_info.F90 @@ -1,504 +1,504 @@ -module write_step_info_interface - interface - subroutine write_step_info(istep,outfreq, mesh) - use MOD_MESH - integer :: istep,outfreq - type(t_mesh), intent(in) , target :: mesh - end subroutine - end interface -end module - -! -! -!=============================================================================== -subroutine write_step_info(istep,outfreq, mesh) - use g_config, only: dt, use_ice - use MOD_MESH - use o_PARAM - use g_PARSUP - use o_ARRAYS - use i_ARRAYS - use g_comm_auto - implicit none - - integer :: n, istep,outfreq - real(kind=WP) :: int_eta, int_hbar, int_wflux, int_hflux, int_temp, int_salt - real(kind=WP) :: min_eta, min_hbar, min_wflux, min_hflux, min_temp, min_salt, & - min_wvel,min_hnode,min_deta,min_wvel2,min_hnode2, & - min_vvel, min_vvel2, min_uvel, min_uvel2 - real(kind=WP) :: max_eta, max_hbar, max_wflux, max_hflux, max_temp, max_salt, & - max_wvel, max_hnode, max_deta, max_wvel2, max_hnode2, max_m_ice, & - max_vvel, max_vvel2, max_uvel, max_uvel2, & - max_cfl_z, max_pgfx, max_pgfy, max_kv, max_av - real(kind=WP) :: int_deta , int_dhbar - real(kind=WP) :: loc, loc_eta, loc_hbar, loc_deta, loc_dhbar, loc_wflux,loc_hflux, loc_temp, loc_salt - type(t_mesh), intent(in) , target :: mesh -#include "associate_mesh.h" - if (mod(istep,outfreq)==0) then - - !_______________________________________________________________________ - int_eta =0. - int_hbar =0. - int_deta =0. - int_dhbar =0. - int_wflux =0. - int_hflux =0. - int_temp =0. - int_salt =0. - loc_eta =0. - loc_hbar =0. - loc_deta =0. - loc_dhbar =0. - loc_wflux =0. -!!PS loc_hflux =0. -!!PS loc_temp =0. -!!PS loc_salt =0. - loc =0. - !_______________________________________________________________________ - do n=1, myDim_nod2D - if (ulevels_nod2D(n)>1) cycle - loc_eta = loc_eta + area(ulevels_nod2D(n), n)*eta_n(n) - loc_hbar = loc_hbar + area(ulevels_nod2D(n), n)*hbar(n) - loc_deta = loc_deta + area(ulevels_nod2D(n), n)*d_eta(n) - loc_dhbar = loc_dhbar + area(ulevels_nod2D(n), n)*(hbar(n)-hbar_old(n)) - loc_wflux = loc_wflux + area(ulevels_nod2D(n), n)*water_flux(n) -!!PS loc_hflux = loc_hflux + area(1, n)*heat_flux(n) -!!PS loc_temp = loc_temp + area(1, n)*sum(tr_arr(:,n,1))/(nlevels_nod2D(n)-1) -!!PS loc_salt = loc_salt + area(1, n)*sum(tr_arr(:,n,2))/(nlevels_nod2D(n)-1) - end do - - !_______________________________________________________________________ - call MPI_AllREDUCE(loc_eta , int_eta , 1, MPI_DOUBLE_PRECISION, MPI_SUM, MPI_COMM_FESOM, MPIerr) - call MPI_AllREDUCE(loc_hbar , int_hbar , 1, MPI_DOUBLE_PRECISION, MPI_SUM, MPI_COMM_FESOM, MPIerr) - call MPI_AllREDUCE(loc_deta , int_deta , 1, MPI_DOUBLE_PRECISION, MPI_SUM, MPI_COMM_FESOM, MPIerr) - call MPI_AllREDUCE(loc_dhbar, int_dhbar, 1, MPI_DOUBLE_PRECISION, MPI_SUM, MPI_COMM_FESOM, MPIerr) - call MPI_AllREDUCE(loc_wflux, int_wflux, 1, MPI_DOUBLE_PRECISION, MPI_SUM, MPI_COMM_FESOM, MPIerr) -!!PS call MPI_AllREDUCE(loc_hflux, int_hflux, 1, MPI_DOUBLE_PRECISION, MPI_SUM, MPI_COMM_FESOM, MPIerr) -!!PS call MPI_AllREDUCE(loc_temp , int_temp , 1, MPI_DOUBLE_PRECISION, MPI_SUM, MPI_COMM_FESOM, MPIerr) -!!PS call MPI_AllREDUCE(loc_salt , int_salt , 1, MPI_DOUBLE_PRECISION, MPI_SUM, MPI_COMM_FESOM, MPIerr) - - int_eta = int_eta /ocean_area - int_hbar = int_hbar /ocean_area - int_deta = int_deta /ocean_area - int_dhbar= int_dhbar/ocean_area - int_wflux= int_wflux/ocean_area - -!!PS int_eta = int_eta /ocean_areawithcav -!!PS int_hbar = int_hbar /ocean_areawithcav -!!PS int_deta = int_deta /ocean_areawithcav -!!PS int_dhbar= int_dhbar/ocean_areawithcav -!!PS int_wflux= int_wflux/ocean_areawithcav - -!!PS int_hflux= int_hflux/ocean_area -!!PS int_temp = int_temp /ocean_area -!!PS int_salt = int_salt /ocean_area - - !_______________________________________________________________________ - loc = minval(eta_n(1:myDim_nod2D)) - call MPI_AllREDUCE(loc , min_eta , 1, MPI_DOUBLE_PRECISION, MPI_MIN, MPI_COMM_FESOM, MPIerr) - loc = minval(hbar(1:myDim_nod2D)) - call MPI_AllREDUCE(loc , min_hbar , 1, MPI_DOUBLE_PRECISION, MPI_MIN, MPI_COMM_FESOM, MPIerr) - loc = minval(water_flux(1:myDim_nod2D)) - call MPI_AllREDUCE(loc , min_wflux, 1, MPI_DOUBLE_PRECISION, MPI_MIN, MPI_COMM_FESOM, MPIerr) - loc = minval(heat_flux(1:myDim_nod2D)) - call MPI_AllREDUCE(loc , min_hflux, 1, MPI_DOUBLE_PRECISION, MPI_MIN, MPI_COMM_FESOM, MPIerr) - loc = minval(tr_arr(:,1:myDim_nod2D,1),MASK=(tr_arr(:,1:myDim_nod2D,2)/=0.0)) - call MPI_AllREDUCE(loc , min_temp , 1, MPI_DOUBLE_PRECISION, MPI_MIN, MPI_COMM_FESOM, MPIerr) - loc = minval(tr_arr(:,1:myDim_nod2D,2),MASK=(tr_arr(:,1:myDim_nod2D,2)/=0.0)) - call MPI_AllREDUCE(loc , min_salt , 1, MPI_DOUBLE_PRECISION, MPI_MIN, MPI_COMM_FESOM, MPIerr) - loc = minval(Wvel(1,1:myDim_nod2D)) - call MPI_AllREDUCE(loc , min_wvel , 1, MPI_DOUBLE_PRECISION, MPI_MIN, MPI_COMM_FESOM, MPIerr) - loc = minval(Wvel(2,1:myDim_nod2D)) - call MPI_AllREDUCE(loc , min_wvel2 , 1, MPI_DOUBLE_PRECISION, MPI_MIN, MPI_COMM_FESOM, MPIerr) - loc = minval(Unode(1,1,1:myDim_nod2D)) - call MPI_AllREDUCE(loc , min_uvel , 1, MPI_DOUBLE_PRECISION, MPI_MIN, MPI_COMM_FESOM, MPIerr) - loc = minval(Unode(1,2,1:myDim_nod2D)) - call MPI_AllREDUCE(loc , min_uvel2 , 1, MPI_DOUBLE_PRECISION, MPI_MIN, MPI_COMM_FESOM, MPIerr) - loc = minval(Unode(2,1,1:myDim_nod2D)) - call MPI_AllREDUCE(loc , min_vvel , 1, MPI_DOUBLE_PRECISION, MPI_MIN, MPI_COMM_FESOM, MPIerr) - loc = minval(Unode(2,2,1:myDim_nod2D)) - call MPI_AllREDUCE(loc , min_vvel2 , 1, MPI_DOUBLE_PRECISION, MPI_MIN, MPI_COMM_FESOM, MPIerr) - loc = minval(d_eta(1:myDim_nod2D)) - call MPI_AllREDUCE(loc , min_deta , 1, MPI_DOUBLE_PRECISION, MPI_MIN, MPI_COMM_FESOM, MPIerr) - loc = minval(hnode(1,1:myDim_nod2D),MASK=(hnode(1,1:myDim_nod2D)/=0.0)) - call MPI_AllREDUCE(loc , min_hnode , 1, MPI_DOUBLE_PRECISION, MPI_MIN, MPI_COMM_FESOM, MPIerr) - loc = minval(hnode(2,1:myDim_nod2D),MASK=(hnode(2,1:myDim_nod2D)/=0.0)) - call MPI_AllREDUCE(loc , min_hnode2 , 1, MPI_DOUBLE_PRECISION, MPI_MIN, MPI_COMM_FESOM, MPIerr) - - !_______________________________________________________________________ - loc = maxval(eta_n(1:myDim_nod2D)) - call MPI_AllREDUCE(loc , max_eta , 1, MPI_DOUBLE_PRECISION, MPI_MAX, MPI_COMM_FESOM, MPIerr) - loc = maxval(hbar(1:myDim_nod2D)) - call MPI_AllREDUCE(loc , max_hbar , 1, MPI_DOUBLE_PRECISION, MPI_MAX, MPI_COMM_FESOM, MPIerr) - loc = maxval(water_flux(1:myDim_nod2D)) - call MPI_AllREDUCE(loc , max_wflux, 1, MPI_DOUBLE_PRECISION, MPI_MAX, MPI_COMM_FESOM, MPIerr) - loc = maxval(heat_flux(1:myDim_nod2D)) - call MPI_AllREDUCE(loc , max_hflux, 1, MPI_DOUBLE_PRECISION, MPI_MAX, MPI_COMM_FESOM, MPIerr) - loc = maxval(tr_arr(:,1:myDim_nod2D,1),MASK=(tr_arr(:,1:myDim_nod2D,2)/=0.0)) - call MPI_AllREDUCE(loc , max_temp , 1, MPI_DOUBLE_PRECISION, MPI_MAX, MPI_COMM_FESOM, MPIerr) - loc = maxval(tr_arr(:,1:myDim_nod2D,2),MASK=(tr_arr(:,1:myDim_nod2D,2)/=0.0)) - call MPI_AllREDUCE(loc , max_salt , 1, MPI_DOUBLE_PRECISION, MPI_MAX, MPI_COMM_FESOM, MPIerr) - loc = maxval(Wvel(1,1:myDim_nod2D)) - call MPI_AllREDUCE(loc , max_wvel , 1, MPI_DOUBLE_PRECISION, MPI_MAX, MPI_COMM_FESOM, MPIerr) - loc = maxval(Wvel(2,1:myDim_nod2D)) - call MPI_AllREDUCE(loc , max_wvel2 , 1, MPI_DOUBLE_PRECISION, MPI_MAX, MPI_COMM_FESOM, MPIerr) - loc = maxval(Unode(1,1,1:myDim_nod2D)) - call MPI_AllREDUCE(loc , max_uvel , 1, MPI_DOUBLE_PRECISION, MPI_MAX, MPI_COMM_FESOM, MPIerr) - loc = maxval(Unode(1,2,1:myDim_nod2D)) - call MPI_AllREDUCE(loc , max_uvel2 , 1, MPI_DOUBLE_PRECISION, MPI_MAX, MPI_COMM_FESOM, MPIerr) - loc = maxval(Unode(2,1,1:myDim_nod2D)) - call MPI_AllREDUCE(loc , max_vvel , 1, MPI_DOUBLE_PRECISION, MPI_MAX, MPI_COMM_FESOM, MPIerr) - loc = maxval(Unode(2,2,1:myDim_nod2D)) - call MPI_AllREDUCE(loc , max_vvel2 , 1, MPI_DOUBLE_PRECISION, MPI_MAX, MPI_COMM_FESOM, MPIerr) - loc = maxval(d_eta(1:myDim_nod2D)) - call MPI_AllREDUCE(loc , max_deta , 1, MPI_DOUBLE_PRECISION, MPI_MAX, MPI_COMM_FESOM, MPIerr) - loc = maxval(hnode(1,1:myDim_nod2D),MASK=(hnode(1,1:myDim_nod2D)/=0.0)) - call MPI_AllREDUCE(loc , max_hnode , 1, MPI_DOUBLE_PRECISION, MPI_MAX, MPI_COMM_FESOM, MPIerr) - loc = maxval(hnode(2,1:myDim_nod2D),MASK=(hnode(2,1:myDim_nod2D)/=0.0)) - call MPI_AllREDUCE(loc , max_hnode2 , 1, MPI_DOUBLE_PRECISION, MPI_MAX, MPI_COMM_FESOM, MPIerr) - loc = maxval(CFL_z(:,1:myDim_nod2D)) - call MPI_AllREDUCE(loc , max_cfl_z , 1, MPI_DOUBLE_PRECISION, MPI_MAX, MPI_COMM_FESOM, MPIerr) - loc = maxval(abs(pgf_x(:,1:myDim_nod2D))) - call MPI_AllREDUCE(loc , max_pgfx , 1, MPI_DOUBLE_PRECISION, MPI_MAX, MPI_COMM_FESOM, MPIerr) - loc = maxval(abs(pgf_y(:,1:myDim_nod2D))) - call MPI_AllREDUCE(loc , max_pgfy , 1, MPI_DOUBLE_PRECISION, MPI_MAX, MPI_COMM_FESOM, MPIerr) - if (use_ice) then - loc = maxval(m_ice(1:myDim_nod2D)) - call MPI_AllREDUCE(loc , max_m_ice , 1, MPI_DOUBLE_PRECISION, MPI_MAX, MPI_COMM_FESOM, MPIerr) - end if - loc = maxval(abs(Av(:,1:myDim_nod2D))) - call MPI_AllREDUCE(loc , max_av , 1, MPI_DOUBLE_PRECISION, MPI_MAX, MPI_COMM_FESOM, MPIerr) - loc = maxval(abs(Kv(:,1:myDim_nod2D))) - call MPI_AllREDUCE(loc , max_kv , 1, MPI_DOUBLE_PRECISION, MPI_MAX, MPI_COMM_FESOM, MPIerr) - !_______________________________________________________________________ - if (mype==0) then - write(*,*) '___CHECK GLOBAL OCEAN VARIABLES --> mstep=',mstep - write(*,*) ' ___global estimat of eta & hbar____________________' - write(*,*) ' int(eta), int(hbar) =', int_eta, int_hbar - write(*,*) ' --> error(eta-hbar) =', int_eta-int_hbar - write(*,*) ' min(eta) , max(eta) =', min_eta, max_eta - write(*,*) ' max(hbar), max(hbar) =', min_hbar, max_hbar - write(*,*) - write(*,*) ' int(deta), int(dhbar) =', int_deta, int_dhbar - write(*,*) ' --> error(deta-dhbar) =', int_deta-int_dhbar - write(*,*) ' --> error(deta-wflux) =', int_deta-int_wflux - write(*,*) ' --> error(dhbar-wflux) =', int_dhbar-int_wflux - write(*,*) - write(*,*) ' -int(wflux)*dt =', int_wflux*dt*(-1.0) - write(*,*) ' int(deta )-int(wflux)*dt =', int_deta-int_wflux*dt*(-1.0) - write(*,*) ' int(dhbar)-int(wflux)*dt =', int_dhbar-int_wflux*dt*(-1.0) - write(*,*) - write(*,*) ' ___global min/max/mean --> mstep=',mstep,'____________' - write(*,"(A, ES10.3, A, ES10.3, A, A )") ' eta= ', min_eta ,' | ',max_eta ,' | ','N.A.' - write(*,"(A, ES10.3, A, ES10.3, A, A )") ' deta= ', min_deta ,' | ',max_deta ,' | ','N.A.' - write(*,"(A, ES10.3, A, ES10.3, A, A )") ' hbar= ', min_hbar ,' | ',max_hbar ,' | ','N.A.' - write(*,"(A, ES10.3, A, ES10.3, A, ES10.3)") ' wflux= ', min_wflux,' | ',max_wflux,' | ',int_wflux - write(*,"(A, ES10.3, A, ES10.3, A, ES10.3)") ' hflux= ', min_hflux,' | ',max_hflux,' | ',int_hflux - write(*,"(A, ES10.3, A, ES10.3, A, ES10.3)") ' temp= ', min_temp ,' | ',max_temp ,' | ',int_temp - write(*,"(A, ES10.3, A, ES10.3, A, ES10.3)") ' salt= ', min_salt ,' | ',max_salt ,' | ',int_salt - write(*,"(A, ES10.3, A, ES10.3, A, A )") ' wvel(1,:)= ', min_wvel ,' | ',max_wvel ,' | ','N.A.' - write(*,"(A, ES10.3, A, ES10.3, A, A )") ' wvel(2,:)= ', min_wvel2,' | ',max_wvel2,' | ','N.A.' - write(*,"(A, ES10.3, A, ES10.3, A, A )") ' uvel(1,:)= ', min_uvel ,' | ',max_uvel ,' | ','N.A.' - write(*,"(A, ES10.3, A, ES10.3, A, A )") ' uvel(2,:)= ', min_uvel2,' | ',max_uvel2,' | ','N.A.' - write(*,"(A, ES10.3, A, ES10.3, A, A )") ' vvel(1,:)= ', min_vvel ,' | ',max_vvel ,' | ','N.A.' - write(*,"(A, ES10.3, A, ES10.3, A, A )") ' vvel(2,:)= ', min_vvel2,' | ',max_vvel2,' | ','N.A.' - write(*,"(A, ES10.3, A, ES10.3, A, A )") ' hnode(1,:)= ', min_hnode,' | ',max_hnode,' | ','N.A.' - write(*,"(A, ES10.3, A, ES10.3, A, A )") ' hnode(2,:)= ', min_hnode2,' | ',max_hnode2,' | ','N.A.' - write(*,"(A, A , A, ES10.3, A, A )") ' cfl_z= ',' N.A. ',' | ',max_cfl_z ,' | ','N.A.' - write(*,"(A, A , A, ES10.3, A, A )") ' pgf_x= ',' N.A. ',' | ',max_pgfx ,' | ','N.A.' - write(*,"(A, A , A, ES10.3, A, A )") ' pgf_y= ',' N.A. ',' | ',max_pgfy ,' | ','N.A.' - write(*,"(A, A , A, ES10.3, A, A )") ' Av= ',' N.A. ',' | ',max_av ,' | ','N.A.' - write(*,"(A, A , A, ES10.3, A, A )") ' Kv= ',' N.A. ',' | ',max_kv ,' | ','N.A.' - if (use_ice) write(*,"(A, A , A, ES10.3, A, A )") ' m_ice= ',' N.A. ',' | ',max_m_ice ,' | ','N.A.' - write(*,*) - endif - endif ! --> if (mod(istep,logfile_outfreq)==0) then -end subroutine write_step_info -! -! -!=============================================================================== -subroutine check_blowup(istep, mesh) - use g_config, only: logfile_outfreq, which_ALE - use MOD_MESH - use o_PARAM - use g_PARSUP - use o_ARRAYS - use i_ARRAYS - use g_comm_auto - use io_BLOWUP - use g_forcing_arrays - use diagnostics - use write_step_info_interface - implicit none - - integer :: n, nz, istep, found_blowup_loc=0, found_blowup=0 - integer :: el, elidx - type(t_mesh), intent(in), target :: mesh -#include "associate_mesh.h" - !___________________________________________________________________________ -! ! if (mod(istep,logfile_outfreq)==0) then -! ! if (mype==0) then -! ! write(*,*) '___CHECK FOR BLOW UP___________ --> mstep=',istep -! ! write(*,*) -! ! endif - do n=1, myDim_nod2d - - !___________________________________________________________________ - ! check ssh - if ( ((eta_n(n) /= eta_n(n)) .or. & - eta_n(n)<-50.0 .or. eta_n(n)>50.0 .or. & - (d_eta(n) /= d_eta(n)) ) ) then -!!PS eta_n(n)<-10.0 .or. eta_n(n)>10.0)) then - found_blowup_loc=1 - write(*,*) '___CHECK FOR BLOW UP___________ --> mstep=',istep - write(*,*) ' --STOP--> found eta_n become NaN or <-10.0, >10.0' - write(*,*) 'mype = ',mype - write(*,*) 'mstep = ',istep - write(*,*) 'node = ',n - write(*,*) 'uln, nln = ',ulevels_nod2D(n), nlevels_nod2D(n) - write(*,*) 'glon,glat = ',geo_coord_nod2D(:,n)/rad - write(*,*) - write(*,*) 'eta_n(n) = ',eta_n(n) - write(*,*) 'd_eta(n) = ',d_eta(n) - write(*,*) - write(*,*) 'zbar_3d_n = ',zbar_3d_n(:,n) - write(*,*) 'Z_3d_n = ',Z_3d_n(:,n) - write(*,*) - write(*,*) 'ssh_rhs = ',ssh_rhs(n),', ssh_rhs_old = ',ssh_rhs_old(n) - write(*,*) - write(*,*) 'hbar = ',hbar(n),', hbar_old = ',hbar_old(n) - write(*,*) - write(*,*) 'wflux = ',water_flux(n) - write(*,*) - write(*,*) 'u_wind = ',u_wind(n),', v_wind = ',v_wind(n) - write(*,*) - do nz=1,nod_in_elem2D_num(n) - write(*,*) 'stress_surf(1:2,',nz,') = ',stress_surf(:,nod_in_elem2D(nz,n)) - end do - write(*,*) - write(*,*) 'm_ice = ',m_ice(n),', m_ice_old = ',m_ice_old(n) - write(*,*) 'a_ice = ',a_ice(n),', a_ice_old = ',a_ice_old(n) -!!PS write(*,*) 'thdgr = ',thdgr(n),', thdgr_old = ',thdgr_old(n) -!!PS write(*,*) 'thdgrsn = ',thdgrsn(n) - write(*,*) -!!PS if (lcurt_stress_surf) then -!!PS write(*,*) 'curl_stress_surf = ',curl_stress_surf(n) -!!PS write(*,*) -!!PS endif -!!PS do el=1,nod_in_elem2d_num(n) -!!PS elidx = nod_in_elem2D(el,n) -!!PS write(*,*) ' elem#=',el,', elemidx=',elidx -!!PS write(*,*) ' pgf_x =',pgf_x(:,elidx) -!!PS write(*,*) ' pgf_y =',pgf_y(:,elidx) -!!PS ! write(*,*) ' U =',UV(1,:,elidx) -!!PS ! write(*,*) ' V =',UV(2,:,elidx) -!!PS write(*,*) -!!PS enddo -!!PS write(*,*) 'Wvel(1, n) = ',Wvel(,n) - write(*,*) 'Wvel(:, n) = ',Wvel(ulevels_nod2D(n):nlevels_nod2D(n),n) - write(*,*) - write(*,*) 'CFL_z(:,n) = ',CFL_z(ulevels_nod2D(n):nlevels_nod2D(n),n) - write(*,*) -!!PS write(*,*) 'hnode(1, n) = ',hnode(1,n) - write(*,*) 'hnode(:, n) = ',hnode(ulevels_nod2D(n):nlevels_nod2D(n),n) - write(*,*) - - endif - - !___________________________________________________________________ - ! check surface vertical velocity --> in case of zlevel and zstar - ! vertical coordinate its indicator if Volume is conserved for - ! Wvel(1,n)~maschine preccision -!!PS if ( .not. trim(which_ALE)=='linfs' .and. ( Wvel(1, n) /= Wvel(1, n) .or. abs(Wvel(1,n))>1e-12 )) then - if ( .not. trim(which_ALE)=='linfs' .and. ( Wvel(1, n) /= Wvel(1, n) )) then - found_blowup_loc=1 - write(*,*) '___CHECK FOR BLOW UP___________ --> mstep=',istep - write(*,*) ' --STOP--> found surface layer vertical velocity becomes NaN or >1e-12' - write(*,*) 'mype = ',mype - write(*,*) 'mstep = ',istep - write(*,*) 'node = ',n - write(*,*) 'uln, nln = ',ulevels_nod2D(n), nlevels_nod2D(n) - write(*,*) 'glon,glat = ',geo_coord_nod2D(:,n)/rad - write(*,*) - write(*,*) 'Wvel(1, n) = ',Wvel(1,n) - write(*,*) 'Wvel(:, n) = ',Wvel(:,n) - write(*,*) - write(*,*) 'hnode(1, n) = ',hnode(1,n) - write(*,*) 'hnode(:, n) = ',hnode(:,n) - write(*,*) 'hflux = ',heat_flux(n) - write(*,*) 'wflux = ',water_flux(n) - write(*,*) - write(*,*) 'eta_n = ',eta_n(n) - write(*,*) 'd_eta(n) = ',d_eta(n) - write(*,*) 'hbar = ',hbar(n) - write(*,*) 'hbar_old = ',hbar_old(n) - write(*,*) 'ssh_rhs = ',ssh_rhs(n) - write(*,*) 'ssh_rhs_old = ',ssh_rhs_old(n) - write(*,*) - write(*,*) 'CFL_z(:,n) = ',CFL_z(:,n) - write(*,*) - - end if ! --> if ( .not. trim(which_ALE)=='linfs' .and. ... - - !___________________________________________________________________ - ! check surface layer thinknesss - if ( .not. trim(which_ALE)=='linfs' .and. ( hnode(1, n) /= hnode(1, n) .or. hnode(1,n)< 0 )) then - found_blowup_loc=1 - write(*,*) '___CHECK FOR BLOW UP___________ --> mstep=',istep - write(*,*) ' --STOP--> found surface layer thickness becomes NaN or <0' - write(*,*) 'mype = ',mype - write(*,*) 'mstep = ',istep - write(*,*) 'node = ',n - write(*,*) - write(*,*) 'hnode(1, n) = ',hnode(1,n) - write(*,*) 'hnode(:, n) = ',hnode(:,n) - write(*,*) - write(*,*) 'glon,glat = ',geo_coord_nod2D(:,n)/rad - write(*,*) - end if ! --> if ( .not. trim(which_ALE)=='linfs' .and. ... - - - do nz=1,nlevels_nod2D(n)-1 - !_______________________________________________________________ - ! check temp - if ( (tr_arr(nz, n,1) /= tr_arr(nz, n,1)) .or. & - tr_arr(nz, n,1) < -5.0 .or. tr_arr(nz, n,1)>60) then - found_blowup_loc=1 - write(*,*) '___CHECK FOR BLOW UP___________ --> mstep=',istep - write(*,*) ' --STOP--> found temperture becomes NaN or <-5.0, >60' - write(*,*) 'mype = ',mype - write(*,*) 'mstep = ',istep - write(*,*) 'node = ',n - write(*,*) 'lon,lat = ',geo_coord_nod2D(:,n)/rad - write(*,*) 'nz = ',nz - write(*,*) 'nzmin, nzmax= ',ulevels_nod2D(n),nlevels_nod2D(n) - write(*,*) 'x=', geo_coord_nod2D(1,n)/rad, ' ; ', 'y=', geo_coord_nod2D(2,n)/rad - write(*,*) 'z=', Z_n(nz) - write(*,*) 'temp(nz, n) = ',tr_arr(nz, n,1) - write(*,*) 'temp(: , n) = ',tr_arr(:, n,1) - write(*,*) 'temp_old(nz,n)= ',tr_arr_old(nz, n,1) - write(*,*) 'temp_old(: ,n)= ',tr_arr_old(:, n,1) - write(*,*) - write(*,*) 'hflux = ',heat_flux(n) - write(*,*) 'wflux = ',water_flux(n) - write(*,*) - write(*,*) 'eta_n = ',eta_n(n) - write(*,*) 'd_eta(n) = ',d_eta(n) - write(*,*) 'hbar = ',hbar(n) - write(*,*) 'hbar_old = ',hbar_old(n) - write(*,*) 'ssh_rhs = ',ssh_rhs(n) - write(*,*) 'ssh_rhs_old = ',ssh_rhs_old(n) - write(*,*) - write(*,*) 'm_ice = ',m_ice(n) - write(*,*) 'm_ice_old = ',m_ice_old(n) - write(*,*) 'm_snow = ',m_snow(n) - write(*,*) 'm_snow_old = ',m_snow_old(n) - write(*,*) - write(*,*) 'hnode = ',hnode(:,n) - write(*,*) 'hnode_new = ',hnode_new(:,n) - write(*,*) - write(*,*) 'Kv = ',Kv(:,n) - write(*,*) - write(*,*) 'W = ',Wvel(:,n) - write(*,*) - write(*,*) 'CFL_z(:,n) = ',CFL_z(:,n) - write(*,*) -! do el=1,nod_in_elem2d_num(n) -! elidx = nod_in_elem2D(el,n) -! write(*,*) ' elem#=',el,', elemidx=',elidx -! write(*,*) ' helem =',helem(:,elidx) -! write(*,*) ' U =',UV(1,:,elidx) -! write(*,*) ' V =',UV(2,:,elidx) -! enddo - write(*,*) - - endif ! --> if ( (tr_arr(nz, n,1) /= tr_arr(nz, n,1)) .or. & ... - - !_______________________________________________________________ - ! check salt - if ( (tr_arr(nz, n,2) /= tr_arr(nz, n,2)) .or. & - tr_arr(nz, n,2) < 0 .or. tr_arr(nz, n,2)>50 ) then - found_blowup_loc=1 - write(*,*) '___CHECK FOR BLOW UP___________ --> mstep=',istep - write(*,*) ' --STOP--> found salinity becomes NaN or <0, >50' - write(*,*) 'mype = ',mype - write(*,*) 'mstep = ',istep - write(*,*) 'node = ',n - write(*,*) 'nz = ',nz - write(*,*) 'nzmin, nzmax= ',ulevels_nod2D(n),nlevels_nod2D(n) - write(*,*) 'x=', geo_coord_nod2D(1,n)/rad, ' ; ', 'y=', geo_coord_nod2D(2,n)/rad - write(*,*) 'z=', Z_n(nz) - write(*,*) 'salt(nz, n) = ',tr_arr(nz, n,2) - write(*,*) 'salt(: , n) = ',tr_arr(:, n,2) - write(*,*) - write(*,*) 'temp(nz, n) = ',tr_arr(nz, n,1) - write(*,*) 'temp(: , n) = ',tr_arr(:, n,1) - write(*,*) - write(*,*) 'hflux = ',heat_flux(n) - write(*,*) - write(*,*) 'wflux = ',water_flux(n) - write(*,*) 'eta_n = ',eta_n(n) - write(*,*) 'd_eta(n) = ',d_eta(n) - write(*,*) 'hbar = ',hbar(n) - write(*,*) 'hbar_old = ',hbar_old(n) - write(*,*) 'ssh_rhs = ',ssh_rhs(n) - write(*,*) 'ssh_rhs_old = ',ssh_rhs_old(n) - write(*,*) - write(*,*) 'hnode = ',hnode(:,n) - write(*,*) 'hnode_new = ',hnode_new(:,n) - write(*,*) - write(*,*) 'zbar_3d_n = ',zbar_3d_n(:,n) - write(*,*) 'Z_3d_n = ',Z_3d_n(:,n) - write(*,*) - write(*,*) 'Kv = ',Kv(:,n) - write(*,*) - do el=1,nod_in_elem2d_num(n) - elidx = nod_in_elem2D(el,n) - write(*,*) ' elem#=',el,', elemidx=',elidx - write(*,*) ' Av =',Av(:,elidx) -! write(*,*) ' helem =',helem(:,elidx) -! write(*,*) ' U =',UV(1,:,elidx) -! write(*,*) ' V =',UV(2,:,elidx) - enddo - write(*,*) 'Wvel = ',Wvel(:,n) - write(*,*) - write(*,*) 'CFL_z(:,n) = ',CFL_z(:,n) - write(*,*) - write(*,*) 'glon,glat = ',geo_coord_nod2D(:,n)/rad - write(*,*) - endif ! --> if ( (tr_arr(nz, n,2) /= tr_arr(nz, n,2)) .or. & ... - end do ! --> do nz=1,nlevels_nod2D(n)-1 - end do ! --> do n=1, myDim_nod2d -! ! end if - - !_______________________________________________________________________ - ! check globally if one of the cpus hat a blowup situation. if its the - ! case CPU mype==0 needs to write out the stuff. Write out occurs in - ! moment only over CPU mype==0 - call MPI_AllREDUCE(found_blowup_loc , found_blowup , 1, MPI_INTEGER, MPI_MAX, MPI_COMM_FESOM, MPIerr) - if (found_blowup==1) then - call write_step_info(istep,1,mesh) - if (mype==0) then - call sleep(1) - write(*,*) - write(*,*) ' MODEL BLOW UP !!!' - write(*,*) ' ____' - write(*,*) ' __,-~~/~ `---.' - write(*,*) ' _/_,---( , )' - write(*,*) ' __ / < / ) \___' - write(*,*) '- -- ----===;;;`====------------------===;;;===---- -- -' - write(*,*) ' \/ ~"~"~"~"~"~\~"~)~"/' - write(*,*) ' (_ ( \ ( > \)' - write(*,*) ' \_( _ < >_>`' - write(*,*) ' ~ `-i` ::>|--"' - write(*,*) ' I;|.|.|' - write(*,*) ' <|i::|i|`' - write(*,*) ' (` ^`"`- ")' - write(*,*) ' _____.,-#%&$@%#&#~,._____' - write(*,*) - end if - call blowup(istep, mesh) - if (mype==0) write(*,*) ' --> finished writing blow up file' - call par_ex - endif -end subroutine +module write_step_info_interface + interface + subroutine write_step_info(istep,outfreq, mesh) + use MOD_MESH + integer :: istep,outfreq + type(t_mesh), intent(in) , target :: mesh + end subroutine + end interface +end module + +! +! +!=============================================================================== +subroutine write_step_info(istep,outfreq, mesh) + use g_config, only: dt, use_ice + use MOD_MESH + use o_PARAM + use g_PARSUP + use o_ARRAYS + use i_ARRAYS + use g_comm_auto + implicit none + + integer :: n, istep,outfreq + real(kind=WP) :: int_eta, int_hbar, int_wflux, int_hflux, int_temp, int_salt + real(kind=WP) :: min_eta, min_hbar, min_wflux, min_hflux, min_temp, min_salt, & + min_wvel,min_hnode,min_deta,min_wvel2,min_hnode2, & + min_vvel, min_vvel2, min_uvel, min_uvel2 + real(kind=WP) :: max_eta, max_hbar, max_wflux, max_hflux, max_temp, max_salt, & + max_wvel, max_hnode, max_deta, max_wvel2, max_hnode2, max_m_ice, & + max_vvel, max_vvel2, max_uvel, max_uvel2, & + max_cfl_z, max_pgfx, max_pgfy, max_kv, max_av + real(kind=WP) :: int_deta , int_dhbar + real(kind=WP) :: loc, loc_eta, loc_hbar, loc_deta, loc_dhbar, loc_wflux,loc_hflux, loc_temp, loc_salt + type(t_mesh), intent(in) , target :: mesh +#include "associate_mesh.h" + if (mod(istep,outfreq)==0) then + + !_______________________________________________________________________ + int_eta =0. + int_hbar =0. + int_deta =0. + int_dhbar =0. + int_wflux =0. + int_hflux =0. + int_temp =0. + int_salt =0. + loc_eta =0. + loc_hbar =0. + loc_deta =0. + loc_dhbar =0. + loc_wflux =0. +!!PS loc_hflux =0. +!!PS loc_temp =0. +!!PS loc_salt =0. + loc =0. + !_______________________________________________________________________ + do n=1, myDim_nod2D +!!PS if (ulevels_nod2D(n)>1) cycle + loc_eta = loc_eta + areasvol(ulevels_nod2D(n), n)*eta_n(n) + loc_hbar = loc_hbar + areasvol(ulevels_nod2D(n), n)*hbar(n) + loc_deta = loc_deta + areasvol(ulevels_nod2D(n), n)*d_eta(n) + loc_dhbar = loc_dhbar + areasvol(ulevels_nod2D(n), n)*(hbar(n)-hbar_old(n)) + loc_wflux = loc_wflux + areasvol(ulevels_nod2D(n), n)*water_flux(n) +!!PS loc_hflux = loc_hflux + area(1, n)*heat_flux(n) +!!PS loc_temp = loc_temp + area(1, n)*sum(tr_arr(:,n,1))/(nlevels_nod2D(n)-1) +!!PS loc_salt = loc_salt + area(1, n)*sum(tr_arr(:,n,2))/(nlevels_nod2D(n)-1) + end do + + !_______________________________________________________________________ + call MPI_AllREDUCE(loc_eta , int_eta , 1, MPI_DOUBLE_PRECISION, MPI_SUM, MPI_COMM_FESOM, MPIerr) + call MPI_AllREDUCE(loc_hbar , int_hbar , 1, MPI_DOUBLE_PRECISION, MPI_SUM, MPI_COMM_FESOM, MPIerr) + call MPI_AllREDUCE(loc_deta , int_deta , 1, MPI_DOUBLE_PRECISION, MPI_SUM, MPI_COMM_FESOM, MPIerr) + call MPI_AllREDUCE(loc_dhbar, int_dhbar, 1, MPI_DOUBLE_PRECISION, MPI_SUM, MPI_COMM_FESOM, MPIerr) + call MPI_AllREDUCE(loc_wflux, int_wflux, 1, MPI_DOUBLE_PRECISION, MPI_SUM, MPI_COMM_FESOM, MPIerr) +!!PS call MPI_AllREDUCE(loc_hflux, int_hflux, 1, MPI_DOUBLE_PRECISION, MPI_SUM, MPI_COMM_FESOM, MPIerr) +!!PS call MPI_AllREDUCE(loc_temp , int_temp , 1, MPI_DOUBLE_PRECISION, MPI_SUM, MPI_COMM_FESOM, MPIerr) +!!PS call MPI_AllREDUCE(loc_salt , int_salt , 1, MPI_DOUBLE_PRECISION, MPI_SUM, MPI_COMM_FESOM, MPIerr) +! +!!PS int_eta = int_eta /ocean_area +!!PS int_hbar = int_hbar /ocean_area +!!PS int_deta = int_deta /ocean_area +!!PS int_dhbar= int_dhbar/ocean_area +!!PS int_wflux= int_wflux/ocean_area + + int_eta = int_eta /ocean_areawithcav + int_hbar = int_hbar /ocean_areawithcav + int_deta = int_deta /ocean_areawithcav + int_dhbar= int_dhbar/ocean_areawithcav + int_wflux= int_wflux/ocean_areawithcav + +!!PS int_hflux= int_hflux/ocean_area +!!PS int_temp = int_temp /ocean_area +!!PS int_salt = int_salt /ocean_area + + !_______________________________________________________________________ + loc = minval(eta_n(1:myDim_nod2D)) + call MPI_AllREDUCE(loc , min_eta , 1, MPI_DOUBLE_PRECISION, MPI_MIN, MPI_COMM_FESOM, MPIerr) + loc = minval(hbar(1:myDim_nod2D)) + call MPI_AllREDUCE(loc , min_hbar , 1, MPI_DOUBLE_PRECISION, MPI_MIN, MPI_COMM_FESOM, MPIerr) + loc = minval(water_flux(1:myDim_nod2D)) + call MPI_AllREDUCE(loc , min_wflux, 1, MPI_DOUBLE_PRECISION, MPI_MIN, MPI_COMM_FESOM, MPIerr) + loc = minval(heat_flux(1:myDim_nod2D)) + call MPI_AllREDUCE(loc , min_hflux, 1, MPI_DOUBLE_PRECISION, MPI_MIN, MPI_COMM_FESOM, MPIerr) + loc = minval(tr_arr(:,1:myDim_nod2D,1),MASK=(tr_arr(:,1:myDim_nod2D,2)/=0.0)) + call MPI_AllREDUCE(loc , min_temp , 1, MPI_DOUBLE_PRECISION, MPI_MIN, MPI_COMM_FESOM, MPIerr) + loc = minval(tr_arr(:,1:myDim_nod2D,2),MASK=(tr_arr(:,1:myDim_nod2D,2)/=0.0)) + call MPI_AllREDUCE(loc , min_salt , 1, MPI_DOUBLE_PRECISION, MPI_MIN, MPI_COMM_FESOM, MPIerr) + loc = minval(Wvel(1,1:myDim_nod2D)) + call MPI_AllREDUCE(loc , min_wvel , 1, MPI_DOUBLE_PRECISION, MPI_MIN, MPI_COMM_FESOM, MPIerr) + loc = minval(Wvel(2,1:myDim_nod2D)) + call MPI_AllREDUCE(loc , min_wvel2 , 1, MPI_DOUBLE_PRECISION, MPI_MIN, MPI_COMM_FESOM, MPIerr) + loc = minval(Unode(1,1,1:myDim_nod2D)) + call MPI_AllREDUCE(loc , min_uvel , 1, MPI_DOUBLE_PRECISION, MPI_MIN, MPI_COMM_FESOM, MPIerr) + loc = minval(Unode(1,2,1:myDim_nod2D)) + call MPI_AllREDUCE(loc , min_uvel2 , 1, MPI_DOUBLE_PRECISION, MPI_MIN, MPI_COMM_FESOM, MPIerr) + loc = minval(Unode(2,1,1:myDim_nod2D)) + call MPI_AllREDUCE(loc , min_vvel , 1, MPI_DOUBLE_PRECISION, MPI_MIN, MPI_COMM_FESOM, MPIerr) + loc = minval(Unode(2,2,1:myDim_nod2D)) + call MPI_AllREDUCE(loc , min_vvel2 , 1, MPI_DOUBLE_PRECISION, MPI_MIN, MPI_COMM_FESOM, MPIerr) + loc = minval(d_eta(1:myDim_nod2D)) + call MPI_AllREDUCE(loc , min_deta , 1, MPI_DOUBLE_PRECISION, MPI_MIN, MPI_COMM_FESOM, MPIerr) + loc = minval(hnode(1,1:myDim_nod2D),MASK=(hnode(1,1:myDim_nod2D)/=0.0)) + call MPI_AllREDUCE(loc , min_hnode , 1, MPI_DOUBLE_PRECISION, MPI_MIN, MPI_COMM_FESOM, MPIerr) + loc = minval(hnode(2,1:myDim_nod2D),MASK=(hnode(2,1:myDim_nod2D)/=0.0)) + call MPI_AllREDUCE(loc , min_hnode2 , 1, MPI_DOUBLE_PRECISION, MPI_MIN, MPI_COMM_FESOM, MPIerr) + + !_______________________________________________________________________ + loc = maxval(eta_n(1:myDim_nod2D)) + call MPI_AllREDUCE(loc , max_eta , 1, MPI_DOUBLE_PRECISION, MPI_MAX, MPI_COMM_FESOM, MPIerr) + loc = maxval(hbar(1:myDim_nod2D)) + call MPI_AllREDUCE(loc , max_hbar , 1, MPI_DOUBLE_PRECISION, MPI_MAX, MPI_COMM_FESOM, MPIerr) + loc = maxval(water_flux(1:myDim_nod2D)) + call MPI_AllREDUCE(loc , max_wflux, 1, MPI_DOUBLE_PRECISION, MPI_MAX, MPI_COMM_FESOM, MPIerr) + loc = maxval(heat_flux(1:myDim_nod2D)) + call MPI_AllREDUCE(loc , max_hflux, 1, MPI_DOUBLE_PRECISION, MPI_MAX, MPI_COMM_FESOM, MPIerr) + loc = maxval(tr_arr(:,1:myDim_nod2D,1),MASK=(tr_arr(:,1:myDim_nod2D,2)/=0.0)) + call MPI_AllREDUCE(loc , max_temp , 1, MPI_DOUBLE_PRECISION, MPI_MAX, MPI_COMM_FESOM, MPIerr) + loc = maxval(tr_arr(:,1:myDim_nod2D,2),MASK=(tr_arr(:,1:myDim_nod2D,2)/=0.0)) + call MPI_AllREDUCE(loc , max_salt , 1, MPI_DOUBLE_PRECISION, MPI_MAX, MPI_COMM_FESOM, MPIerr) + loc = maxval(Wvel(1,1:myDim_nod2D)) + call MPI_AllREDUCE(loc , max_wvel , 1, MPI_DOUBLE_PRECISION, MPI_MAX, MPI_COMM_FESOM, MPIerr) + loc = maxval(Wvel(2,1:myDim_nod2D)) + call MPI_AllREDUCE(loc , max_wvel2 , 1, MPI_DOUBLE_PRECISION, MPI_MAX, MPI_COMM_FESOM, MPIerr) + loc = maxval(Unode(1,1,1:myDim_nod2D)) + call MPI_AllREDUCE(loc , max_uvel , 1, MPI_DOUBLE_PRECISION, MPI_MAX, MPI_COMM_FESOM, MPIerr) + loc = maxval(Unode(1,2,1:myDim_nod2D)) + call MPI_AllREDUCE(loc , max_uvel2 , 1, MPI_DOUBLE_PRECISION, MPI_MAX, MPI_COMM_FESOM, MPIerr) + loc = maxval(Unode(2,1,1:myDim_nod2D)) + call MPI_AllREDUCE(loc , max_vvel , 1, MPI_DOUBLE_PRECISION, MPI_MAX, MPI_COMM_FESOM, MPIerr) + loc = maxval(Unode(2,2,1:myDim_nod2D)) + call MPI_AllREDUCE(loc , max_vvel2 , 1, MPI_DOUBLE_PRECISION, MPI_MAX, MPI_COMM_FESOM, MPIerr) + loc = maxval(d_eta(1:myDim_nod2D)) + call MPI_AllREDUCE(loc , max_deta , 1, MPI_DOUBLE_PRECISION, MPI_MAX, MPI_COMM_FESOM, MPIerr) + loc = maxval(hnode(1,1:myDim_nod2D),MASK=(hnode(1,1:myDim_nod2D)/=0.0)) + call MPI_AllREDUCE(loc , max_hnode , 1, MPI_DOUBLE_PRECISION, MPI_MAX, MPI_COMM_FESOM, MPIerr) + loc = maxval(hnode(2,1:myDim_nod2D),MASK=(hnode(2,1:myDim_nod2D)/=0.0)) + call MPI_AllREDUCE(loc , max_hnode2 , 1, MPI_DOUBLE_PRECISION, MPI_MAX, MPI_COMM_FESOM, MPIerr) + loc = maxval(CFL_z(:,1:myDim_nod2D)) + call MPI_AllREDUCE(loc , max_cfl_z , 1, MPI_DOUBLE_PRECISION, MPI_MAX, MPI_COMM_FESOM, MPIerr) + loc = maxval(abs(pgf_x(:,1:myDim_nod2D))) + call MPI_AllREDUCE(loc , max_pgfx , 1, MPI_DOUBLE_PRECISION, MPI_MAX, MPI_COMM_FESOM, MPIerr) + loc = maxval(abs(pgf_y(:,1:myDim_nod2D))) + call MPI_AllREDUCE(loc , max_pgfy , 1, MPI_DOUBLE_PRECISION, MPI_MAX, MPI_COMM_FESOM, MPIerr) + if (use_ice) then + loc = maxval(m_ice(1:myDim_nod2D)) + call MPI_AllREDUCE(loc , max_m_ice , 1, MPI_DOUBLE_PRECISION, MPI_MAX, MPI_COMM_FESOM, MPIerr) + end if + loc = maxval(abs(Av(:,1:myDim_nod2D))) + call MPI_AllREDUCE(loc , max_av , 1, MPI_DOUBLE_PRECISION, MPI_MAX, MPI_COMM_FESOM, MPIerr) + loc = maxval(abs(Kv(:,1:myDim_nod2D))) + call MPI_AllREDUCE(loc , max_kv , 1, MPI_DOUBLE_PRECISION, MPI_MAX, MPI_COMM_FESOM, MPIerr) + !_______________________________________________________________________ + if (mype==0) then + write(*,*) '___CHECK GLOBAL OCEAN VARIABLES --> mstep=',mstep + write(*,*) ' ___global estimat of eta & hbar____________________' + write(*,*) ' int(eta), int(hbar) =', int_eta, int_hbar + write(*,*) ' --> error(eta-hbar) =', int_eta-int_hbar + write(*,*) ' min(eta) , max(eta) =', min_eta, max_eta + write(*,*) ' max(hbar), max(hbar) =', min_hbar, max_hbar + write(*,*) + write(*,*) ' int(deta), int(dhbar) =', int_deta, int_dhbar + write(*,*) ' --> error(deta-dhbar) =', int_deta-int_dhbar + write(*,*) ' --> error(deta-wflux) =', int_deta-int_wflux + write(*,*) ' --> error(dhbar-wflux) =', int_dhbar-int_wflux + write(*,*) + write(*,*) ' -int(wflux)*dt =', int_wflux*dt*(-1.0) + write(*,*) ' int(deta )-int(wflux)*dt =', int_deta-int_wflux*dt*(-1.0) + write(*,*) ' int(dhbar)-int(wflux)*dt =', int_dhbar-int_wflux*dt*(-1.0) + write(*,*) + write(*,*) ' ___global min/max/mean --> mstep=',mstep,'____________' + write(*,"(A, ES10.3, A, ES10.3, A, A )") ' eta= ', min_eta ,' | ',max_eta ,' | ','N.A.' + write(*,"(A, ES10.3, A, ES10.3, A, A )") ' deta= ', min_deta ,' | ',max_deta ,' | ','N.A.' + write(*,"(A, ES10.3, A, ES10.3, A, A )") ' hbar= ', min_hbar ,' | ',max_hbar ,' | ','N.A.' + write(*,"(A, ES10.3, A, ES10.3, A, ES10.3)") ' wflux= ', min_wflux,' | ',max_wflux,' | ',int_wflux + write(*,"(A, ES10.3, A, ES10.3, A, ES10.3)") ' hflux= ', min_hflux,' | ',max_hflux,' | ',int_hflux + write(*,"(A, ES10.3, A, ES10.3, A, ES10.3)") ' temp= ', min_temp ,' | ',max_temp ,' | ',int_temp + write(*,"(A, ES10.3, A, ES10.3, A, ES10.3)") ' salt= ', min_salt ,' | ',max_salt ,' | ',int_salt + write(*,"(A, ES10.3, A, ES10.3, A, A )") ' wvel(1,:)= ', min_wvel ,' | ',max_wvel ,' | ','N.A.' + write(*,"(A, ES10.3, A, ES10.3, A, A )") ' wvel(2,:)= ', min_wvel2,' | ',max_wvel2,' | ','N.A.' + write(*,"(A, ES10.3, A, ES10.3, A, A )") ' uvel(1,:)= ', min_uvel ,' | ',max_uvel ,' | ','N.A.' + write(*,"(A, ES10.3, A, ES10.3, A, A )") ' uvel(2,:)= ', min_uvel2,' | ',max_uvel2,' | ','N.A.' + write(*,"(A, ES10.3, A, ES10.3, A, A )") ' vvel(1,:)= ', min_vvel ,' | ',max_vvel ,' | ','N.A.' + write(*,"(A, ES10.3, A, ES10.3, A, A )") ' vvel(2,:)= ', min_vvel2,' | ',max_vvel2,' | ','N.A.' + write(*,"(A, ES10.3, A, ES10.3, A, A )") ' hnode(1,:)= ', min_hnode,' | ',max_hnode,' | ','N.A.' + write(*,"(A, ES10.3, A, ES10.3, A, A )") ' hnode(2,:)= ', min_hnode2,' | ',max_hnode2,' | ','N.A.' + write(*,"(A, A , A, ES10.3, A, A )") ' cfl_z= ',' N.A. ',' | ',max_cfl_z ,' | ','N.A.' + write(*,"(A, A , A, ES10.3, A, A )") ' pgf_x= ',' N.A. ',' | ',max_pgfx ,' | ','N.A.' + write(*,"(A, A , A, ES10.3, A, A )") ' pgf_y= ',' N.A. ',' | ',max_pgfy ,' | ','N.A.' + write(*,"(A, A , A, ES10.3, A, A )") ' Av= ',' N.A. ',' | ',max_av ,' | ','N.A.' + write(*,"(A, A , A, ES10.3, A, A )") ' Kv= ',' N.A. ',' | ',max_kv ,' | ','N.A.' + if (use_ice) write(*,"(A, A , A, ES10.3, A, A )") ' m_ice= ',' N.A. ',' | ',max_m_ice ,' | ','N.A.' + write(*,*) + endif + endif ! --> if (mod(istep,logfile_outfreq)==0) then +end subroutine write_step_info +! +! +!=============================================================================== +subroutine check_blowup(istep, mesh) + use g_config, only: logfile_outfreq, which_ALE + use MOD_MESH + use o_PARAM + use g_PARSUP + use o_ARRAYS + use i_ARRAYS + use g_comm_auto + use io_BLOWUP + use g_forcing_arrays + use diagnostics + use write_step_info_interface + implicit none + + integer :: n, nz, istep, found_blowup_loc=0, found_blowup=0 + integer :: el, elidx + type(t_mesh), intent(in), target :: mesh +#include "associate_mesh.h" + !___________________________________________________________________________ +! ! if (mod(istep,logfile_outfreq)==0) then +! ! if (mype==0) then +! ! write(*,*) '___CHECK FOR BLOW UP___________ --> mstep=',istep +! ! write(*,*) +! ! endif + do n=1, myDim_nod2d + + !___________________________________________________________________ + ! check ssh + if ( ((eta_n(n) /= eta_n(n)) .or. & + eta_n(n)<-50.0 .or. eta_n(n)>50.0 .or. & + (d_eta(n) /= d_eta(n)) ) ) then +!!PS eta_n(n)<-10.0 .or. eta_n(n)>10.0)) then + found_blowup_loc=1 + write(*,*) '___CHECK FOR BLOW UP___________ --> mstep=',istep + write(*,*) ' --STOP--> found eta_n become NaN or <-10.0, >10.0' + write(*,*) 'mype = ',mype + write(*,*) 'mstep = ',istep + write(*,*) 'node = ',n + write(*,*) 'uln, nln = ',ulevels_nod2D(n), nlevels_nod2D(n) + write(*,*) 'glon,glat = ',geo_coord_nod2D(:,n)/rad + write(*,*) + write(*,*) 'eta_n(n) = ',eta_n(n) + write(*,*) 'd_eta(n) = ',d_eta(n) + write(*,*) + write(*,*) 'zbar_3d_n = ',zbar_3d_n(:,n) + write(*,*) 'Z_3d_n = ',Z_3d_n(:,n) + write(*,*) + write(*,*) 'ssh_rhs = ',ssh_rhs(n),', ssh_rhs_old = ',ssh_rhs_old(n) + write(*,*) + write(*,*) 'hbar = ',hbar(n),', hbar_old = ',hbar_old(n) + write(*,*) + write(*,*) 'wflux = ',water_flux(n) + write(*,*) + write(*,*) 'u_wind = ',u_wind(n),', v_wind = ',v_wind(n) + write(*,*) + do nz=1,nod_in_elem2D_num(n) + write(*,*) 'stress_surf(1:2,',nz,') = ',stress_surf(:,nod_in_elem2D(nz,n)) + end do + write(*,*) + write(*,*) 'm_ice = ',m_ice(n),', m_ice_old = ',m_ice_old(n) + write(*,*) 'a_ice = ',a_ice(n),', a_ice_old = ',a_ice_old(n) +!!PS write(*,*) 'thdgr = ',thdgr(n),', thdgr_old = ',thdgr_old(n) +!!PS write(*,*) 'thdgrsn = ',thdgrsn(n) + write(*,*) +!!PS if (lcurt_stress_surf) then +!!PS write(*,*) 'curl_stress_surf = ',curl_stress_surf(n) +!!PS write(*,*) +!!PS endif +!!PS do el=1,nod_in_elem2d_num(n) +!!PS elidx = nod_in_elem2D(el,n) +!!PS write(*,*) ' elem#=',el,', elemidx=',elidx +!!PS write(*,*) ' pgf_x =',pgf_x(:,elidx) +!!PS write(*,*) ' pgf_y =',pgf_y(:,elidx) +!!PS ! write(*,*) ' U =',UV(1,:,elidx) +!!PS ! write(*,*) ' V =',UV(2,:,elidx) +!!PS write(*,*) +!!PS enddo +!!PS write(*,*) 'Wvel(1, n) = ',Wvel(,n) + write(*,*) 'Wvel(:, n) = ',Wvel(ulevels_nod2D(n):nlevels_nod2D(n),n) + write(*,*) + write(*,*) 'CFL_z(:,n) = ',CFL_z(ulevels_nod2D(n):nlevels_nod2D(n),n) + write(*,*) +!!PS write(*,*) 'hnode(1, n) = ',hnode(1,n) + write(*,*) 'hnode(:, n) = ',hnode(ulevels_nod2D(n):nlevels_nod2D(n),n) + write(*,*) + + endif + + !___________________________________________________________________ + ! check surface vertical velocity --> in case of zlevel and zstar + ! vertical coordinate its indicator if Volume is conserved for + ! Wvel(1,n)~maschine preccision +!!PS if ( .not. trim(which_ALE)=='linfs' .and. ( Wvel(1, n) /= Wvel(1, n) .or. abs(Wvel(1,n))>1e-12 )) then + if ( .not. trim(which_ALE)=='linfs' .and. ( Wvel(1, n) /= Wvel(1, n) )) then + found_blowup_loc=1 + write(*,*) '___CHECK FOR BLOW UP___________ --> mstep=',istep + write(*,*) ' --STOP--> found surface layer vertical velocity becomes NaN or >1e-12' + write(*,*) 'mype = ',mype + write(*,*) 'mstep = ',istep + write(*,*) 'node = ',n + write(*,*) 'uln, nln = ',ulevels_nod2D(n), nlevels_nod2D(n) + write(*,*) 'glon,glat = ',geo_coord_nod2D(:,n)/rad + write(*,*) + write(*,*) 'Wvel(1, n) = ',Wvel(1,n) + write(*,*) 'Wvel(:, n) = ',Wvel(:,n) + write(*,*) + write(*,*) 'hnode(1, n) = ',hnode(1,n) + write(*,*) 'hnode(:, n) = ',hnode(:,n) + write(*,*) 'hflux = ',heat_flux(n) + write(*,*) 'wflux = ',water_flux(n) + write(*,*) + write(*,*) 'eta_n = ',eta_n(n) + write(*,*) 'd_eta(n) = ',d_eta(n) + write(*,*) 'hbar = ',hbar(n) + write(*,*) 'hbar_old = ',hbar_old(n) + write(*,*) 'ssh_rhs = ',ssh_rhs(n) + write(*,*) 'ssh_rhs_old = ',ssh_rhs_old(n) + write(*,*) + write(*,*) 'CFL_z(:,n) = ',CFL_z(:,n) + write(*,*) + + end if ! --> if ( .not. trim(which_ALE)=='linfs' .and. ... + + !___________________________________________________________________ + ! check surface layer thinknesss + if ( .not. trim(which_ALE)=='linfs' .and. ( hnode(1, n) /= hnode(1, n) .or. hnode(1,n)< 0 )) then + found_blowup_loc=1 + write(*,*) '___CHECK FOR BLOW UP___________ --> mstep=',istep + write(*,*) ' --STOP--> found surface layer thickness becomes NaN or <0' + write(*,*) 'mype = ',mype + write(*,*) 'mstep = ',istep + write(*,*) 'node = ',n + write(*,*) + write(*,*) 'hnode(1, n) = ',hnode(1,n) + write(*,*) 'hnode(:, n) = ',hnode(:,n) + write(*,*) + write(*,*) 'glon,glat = ',geo_coord_nod2D(:,n)/rad + write(*,*) + end if ! --> if ( .not. trim(which_ALE)=='linfs' .and. ... + + + do nz=1,nlevels_nod2D(n)-1 + !_______________________________________________________________ + ! check temp + if ( (tr_arr(nz, n,1) /= tr_arr(nz, n,1)) .or. & + tr_arr(nz, n,1) < -5.0 .or. tr_arr(nz, n,1)>60) then + found_blowup_loc=1 + write(*,*) '___CHECK FOR BLOW UP___________ --> mstep=',istep + write(*,*) ' --STOP--> found temperture becomes NaN or <-5.0, >60' + write(*,*) 'mype = ',mype + write(*,*) 'mstep = ',istep + write(*,*) 'node = ',n + write(*,*) 'lon,lat = ',geo_coord_nod2D(:,n)/rad + write(*,*) 'nz = ',nz + write(*,*) 'nzmin, nzmax= ',ulevels_nod2D(n),nlevels_nod2D(n) + write(*,*) 'x=', geo_coord_nod2D(1,n)/rad, ' ; ', 'y=', geo_coord_nod2D(2,n)/rad + write(*,*) 'z=', Z_n(nz) + write(*,*) 'temp(nz, n) = ',tr_arr(nz, n,1) + write(*,*) 'temp(: , n) = ',tr_arr(:, n,1) + write(*,*) 'temp_old(nz,n)= ',tr_arr_old(nz, n,1) + write(*,*) 'temp_old(: ,n)= ',tr_arr_old(:, n,1) + write(*,*) + write(*,*) 'hflux = ',heat_flux(n) + write(*,*) 'wflux = ',water_flux(n) + write(*,*) + write(*,*) 'eta_n = ',eta_n(n) + write(*,*) 'd_eta(n) = ',d_eta(n) + write(*,*) 'hbar = ',hbar(n) + write(*,*) 'hbar_old = ',hbar_old(n) + write(*,*) 'ssh_rhs = ',ssh_rhs(n) + write(*,*) 'ssh_rhs_old = ',ssh_rhs_old(n) + write(*,*) + write(*,*) 'm_ice = ',m_ice(n) + write(*,*) 'm_ice_old = ',m_ice_old(n) + write(*,*) 'm_snow = ',m_snow(n) + write(*,*) 'm_snow_old = ',m_snow_old(n) + write(*,*) + write(*,*) 'hnode = ',hnode(:,n) + write(*,*) 'hnode_new = ',hnode_new(:,n) + write(*,*) + write(*,*) 'Kv = ',Kv(:,n) + write(*,*) + write(*,*) 'W = ',Wvel(:,n) + write(*,*) + write(*,*) 'CFL_z(:,n) = ',CFL_z(:,n) + write(*,*) +! do el=1,nod_in_elem2d_num(n) +! elidx = nod_in_elem2D(el,n) +! write(*,*) ' elem#=',el,', elemidx=',elidx +! write(*,*) ' helem =',helem(:,elidx) +! write(*,*) ' U =',UV(1,:,elidx) +! write(*,*) ' V =',UV(2,:,elidx) +! enddo + write(*,*) + + endif ! --> if ( (tr_arr(nz, n,1) /= tr_arr(nz, n,1)) .or. & ... + + !_______________________________________________________________ + ! check salt + if ( (tr_arr(nz, n,2) /= tr_arr(nz, n,2)) .or. & + tr_arr(nz, n,2) < 0 .or. tr_arr(nz, n,2)>50 ) then + found_blowup_loc=1 + write(*,*) '___CHECK FOR BLOW UP___________ --> mstep=',istep + write(*,*) ' --STOP--> found salinity becomes NaN or <0, >50' + write(*,*) 'mype = ',mype + write(*,*) 'mstep = ',istep + write(*,*) 'node = ',n + write(*,*) 'nz = ',nz + write(*,*) 'nzmin, nzmax= ',ulevels_nod2D(n),nlevels_nod2D(n) + write(*,*) 'x=', geo_coord_nod2D(1,n)/rad, ' ; ', 'y=', geo_coord_nod2D(2,n)/rad + write(*,*) 'z=', Z_n(nz) + write(*,*) 'salt(nz, n) = ',tr_arr(nz, n,2) + write(*,*) 'salt(: , n) = ',tr_arr(:, n,2) + write(*,*) + write(*,*) 'temp(nz, n) = ',tr_arr(nz, n,1) + write(*,*) 'temp(: , n) = ',tr_arr(:, n,1) + write(*,*) + write(*,*) 'hflux = ',heat_flux(n) + write(*,*) + write(*,*) 'wflux = ',water_flux(n) + write(*,*) 'eta_n = ',eta_n(n) + write(*,*) 'd_eta(n) = ',d_eta(n) + write(*,*) 'hbar = ',hbar(n) + write(*,*) 'hbar_old = ',hbar_old(n) + write(*,*) 'ssh_rhs = ',ssh_rhs(n) + write(*,*) 'ssh_rhs_old = ',ssh_rhs_old(n) + write(*,*) + write(*,*) 'hnode = ',hnode(:,n) + write(*,*) 'hnode_new = ',hnode_new(:,n) + write(*,*) + write(*,*) 'zbar_3d_n = ',zbar_3d_n(:,n) + write(*,*) 'Z_3d_n = ',Z_3d_n(:,n) + write(*,*) + write(*,*) 'Kv = ',Kv(:,n) + write(*,*) + do el=1,nod_in_elem2d_num(n) + elidx = nod_in_elem2D(el,n) + write(*,*) ' elem#=',el,', elemidx=',elidx + write(*,*) ' Av =',Av(:,elidx) +! write(*,*) ' helem =',helem(:,elidx) +! write(*,*) ' U =',UV(1,:,elidx) +! write(*,*) ' V =',UV(2,:,elidx) + enddo + write(*,*) 'Wvel = ',Wvel(:,n) + write(*,*) + write(*,*) 'CFL_z(:,n) = ',CFL_z(:,n) + write(*,*) + write(*,*) 'glon,glat = ',geo_coord_nod2D(:,n)/rad + write(*,*) + endif ! --> if ( (tr_arr(nz, n,2) /= tr_arr(nz, n,2)) .or. & ... + end do ! --> do nz=1,nlevels_nod2D(n)-1 + end do ! --> do n=1, myDim_nod2d +! ! end if + + !_______________________________________________________________________ + ! check globally if one of the cpus hat a blowup situation. if its the + ! case CPU mype==0 needs to write out the stuff. Write out occurs in + ! moment only over CPU mype==0 + call MPI_AllREDUCE(found_blowup_loc , found_blowup , 1, MPI_INTEGER, MPI_MAX, MPI_COMM_FESOM, MPIerr) + if (found_blowup==1) then + call write_step_info(istep,1,mesh) + if (mype==0) then + call sleep(1) + write(*,*) + write(*,*) ' MODEL BLOW UP !!!' + write(*,*) ' ____' + write(*,*) ' __,-~~/~ `---.' + write(*,*) ' _/_,---( , )' + write(*,*) ' __ / < / ) \___' + write(*,*) '- -- ----===;;;`====------------------===;;;===---- -- -' + write(*,*) ' \/ ~"~"~"~"~"~\~"~)~"/' + write(*,*) ' (_ ( \ ( > \)' + write(*,*) ' \_( _ < >_>`' + write(*,*) ' ~ `-i` ::>|--"' + write(*,*) ' I;|.|.|' + write(*,*) ' <|i::|i|`' + write(*,*) ' (` ^`"`- ")' + write(*,*) ' _____.,-#%&$@%#&#~,._____' + write(*,*) + end if + call blowup(istep, mesh) + if (mype==0) write(*,*) ' --> finished writing blow up file' + call par_ex + endif +end subroutine From cbacb09cbec77c728662fed4bace7d6f38b06eb2 Mon Sep 17 00:00:00 2001 From: Patrick Date: Wed, 30 Jun 2021 16:48:26 +0200 Subject: [PATCH 10/37] update my mesh rotation in my python plotting routine --- view_pscholz/sub_fesom_mesh.py | 20 +++++++++++++------- 1 file changed, 13 insertions(+), 7 deletions(-) diff --git a/view_pscholz/sub_fesom_mesh.py b/view_pscholz/sub_fesom_mesh.py index 2aefb37ff..f3fcc336a 100644 --- a/view_pscholz/sub_fesom_mesh.py +++ b/view_pscholz/sub_fesom_mesh.py @@ -130,7 +130,8 @@ def __init__(self,inputarray): #_______________________________________________________________________ # remove+augment periodic boundary - self.fesom_remove_pbnd() + if (inputarray['mesh_remove_cyc' ] == True): + self.fesom_remove_pbnd() # calculate fesom land mask interactivly #self.fesom_calc_landmask() @@ -285,9 +286,9 @@ def fesom_grid_rot_r2g(self,str_mode='r2g'): #_______________________________________________________________________ # make inverse of rotation matrix - if (str_mode == 'r2g') or (str_mode == 'focus'): - from numpy.linalg import inv - rotate_matrix=inv(rotate_matrix); +# if (str_mode == 'r2g') or (str_mode == 'focus'): +# from numpy.linalg import inv +# rotate_matrix=inv(rotate_matrix); #____3D_________________________________________________________________ # calculate Cartesian coordinates @@ -303,9 +304,14 @@ def fesom_grid_rot_r2g(self,str_mode='r2g'): #_______________________________________________________________________ # rotate to geographical cartesian coordinates: - xg=rotate_matrix[0,0]*xr + rotate_matrix[0,1]*yr + rotate_matrix[0,2]*zr; - yg=rotate_matrix[1,0]*xr + rotate_matrix[1,1]*yr + rotate_matrix[1,2]*zr; - zg=rotate_matrix[2,0]*xr + rotate_matrix[2,1]*yr + rotate_matrix[2,2]*zr; + if (str_mode == 'r2g') or (str_mode == 'focus'): + xg=rotate_matrix[0,0]*xr + rotate_matrix[1,0]*yr + rotate_matrix[2,0]*zr; + yg=rotate_matrix[0,1]*xr + rotate_matrix[1,1]*yr + rotate_matrix[2,1]*zr; + zg=rotate_matrix[0,2]*xr + rotate_matrix[1,2]*yr + rotate_matrix[2,2]*zr; + else: + xg=rotate_matrix[0,0]*xr + rotate_matrix[0,1]*yr + rotate_matrix[0,2]*zr; + yg=rotate_matrix[1,0]*xr + rotate_matrix[1,1]*yr + rotate_matrix[1,2]*zr; + zg=rotate_matrix[2,0]*xr + rotate_matrix[2,1]*yr + rotate_matrix[2,2]*zr; ##______________________________________________________________________ #self.nodes_2d_yg = np.degrees(np.arcsin(zg)); From 0d774f479a09d5d03a0a676dc357ae22276125ae Mon Sep 17 00:00:00 2001 From: dsidoren Date: Mon, 19 Jul 2021 11:32:42 +0200 Subject: [PATCH 11/37] Update gen_forcing_couple.F90 a compilation mistake for AWICM 2.0 was pointed out by Gregor --- src/gen_forcing_couple.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/gen_forcing_couple.F90 b/src/gen_forcing_couple.F90 index 2a39b34b7..d17ee5138 100755 --- a/src/gen_forcing_couple.F90 +++ b/src/gen_forcing_couple.F90 @@ -231,8 +231,8 @@ subroutine update_atm_forcing(istep, mesh) 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) From 6b0511166c03b384561f265866f79e5709720fa8 Mon Sep 17 00:00:00 2001 From: Patrick Date: Tue, 18 May 2021 16:05:36 +0200 Subject: [PATCH 12/37] include non local fluxes from KPP --- src/oce_ale_mixing_kpp.F90 | 17 ++- src/oce_ale_tracer.F90 | 278 +++++++++++++++++++++++++------------ src/oce_modules.F90 | 7 +- 3 files changed, 209 insertions(+), 93 deletions(-) diff --git a/src/oce_ale_mixing_kpp.F90 b/src/oce_ale_mixing_kpp.F90 index d08060e98..14b90ccfe 100755 --- a/src/oce_ale_mixing_kpp.F90 +++ b/src/oce_ale_mixing_kpp.F90 @@ -70,7 +70,9 @@ MODULE o_mixing_KPP_mod real(KIND=WP), dimension(0:nni+1,0:nnj+1) :: wst ! lookup table for ws, the turbulent velocity scale scalars logical :: smooth_blmc=.true. logical :: smooth_hbl =.false. - logical :: smooth_Ri =.false. + logical :: smooth_Ri_hor =.false. + logical :: smooth_Ri_ver =.false. + logical :: limit_hbl_ekmmob =.false. !.true. contains @@ -336,8 +338,10 @@ subroutine oce_mixing_KPP(viscAE, diffK, mesh) DO node=1, myDim_nod2D !+eDim_nod2D nzmin = ulevels_nod2D(node) - ustar(node) = sqrt( sqrt( stress_atmoce_x(node)**2 + stress_atmoce_y(node)**2 )*density_0_r ) ! @ the surface (eqn. 2) - +!!PS ustar(node) = sqrt( sqrt( stress_atmoce_x(node)**2 + stress_atmoce_y(node)**2 )*density_0_r ) ! @ the surface (eqn. 2) + + ustar(node) = sqrt( sqrt( stress_node_surf(1,node)**2 + stress_node_surf(2,node)**2 )*density_0_r ) ! @ the surface (eqn. 2) + ! Surface buoyancy forcing (eqns. A2c & A2d & A3b & A3d) !!PS Bo(node) = -g * ( sw_alpha(1,node) * heat_flux(node) / vcpw & !heat_flux & water_flux: positive up !!PS + sw_beta (1,node) * water_flux(node) * tr_arr(1,node,2)) @@ -397,7 +401,8 @@ subroutine oce_mixing_KPP(viscAE, diffK, mesh) ! only at the end should save some time call exchange_nod(diffK(:,:,1)) call exchange_nod(diffK(:,:,2)) - + call exchange_nod(ghats) + ! OVER ELEMENTS call exchange_nod(viscA) !Warning: don't forget to communicate before averaging on elements!!! minmix=3.0e-3_WP @@ -768,7 +773,7 @@ subroutine ri_iwmix(viscA, diffK, mesh) ! smooth Richardson number in the vertical using a 1-2-1 filter !!PS IF(smooth_richardson_number .and. nlevels_nod2d(node)>2) then - IF(smooth_richardson_number .and. nzmax>2) then + IF(smooth_Ri_ver .and. nzmax>2) then DO mr=1,num_smoothings ri_prev = 0.25_WP * diffK(1, node, 1) !!PS DO nz=2,nlevels_nod2d(node)-1 @@ -781,7 +786,7 @@ subroutine ri_iwmix(viscA, diffK, mesh) END IF END DO - if (smooth_Ri) then + if (smooth_Ri_hor) then call smooth_nod(diffK(:,:,1), 3, mesh) end if diff --git a/src/oce_ale_tracer.F90 b/src/oce_ale_tracer.F90 index 06b2ad8b4..ee689cb26 100644 --- a/src/oce_ale_tracer.F90 +++ b/src/oce_ale_tracer.F90 @@ -400,7 +400,8 @@ subroutine diff_ver_part_impl_ale(tr_num, mesh) use g_PARSUP use g_CONFIG use g_forcing_arrays - use o_mixing_KPP_mod !for ghats _GO_ + use o_mixing_KPP_mod !for ghats _GO_ + use g_cvmix_kpp, only: kpp_nonlcltranspT, kpp_nonlcltranspS, kpp_oblmixc use bc_surface_interface implicit none @@ -425,36 +426,80 @@ subroutine diff_ver_part_impl_ale(tr_num, mesh) Ty1 =0.0_WP ! solve equation diffusion equation implicite part: - ! --> h^(n+0.5)* (T^(n+0.5)-Tstar) = dt*( K_33*d/dz*(T^(n+0.5)-Tstar) + K_33*d/dz*Tstar ) - ! --> dTnew = T^(n+0.5)-Tstar - ! --> h^(n+0.5)* (dTnew) = dt*(K_33*d/dz*dTnew) + K_33*dt*d/dz*Tstar - ! --> h^(n+0.5)* (dTnew) = dt*(K_33*d/dz*dTnew) + RHS - ! --> solve for dT_new - ! - ! ----------- zbar_1, V_1 (Skalar Volume), A_1 (Area of edge), no Cavity A1==V1, yes Cavity A1 !=V1 - ! Z_1 o T_1 - ! ----------- zbar_2, V_2 - ! Z_2 o T_2 - ! ----------- zbar_3, V_3 - ! Z_3 o T_3 - ! ----------- zbar_4 - ! : - ! --> Difference Quotient at Volume _2: ddTnew_2/dt + d/dz*K_33 d/dz*dTnew_2 = 0 --> homogene solution - ! V2*dTnew_2 *h^(n+0.5) = -dt * [ (dTnew_1-dTnew_2)/(Z_1-Z_2)*A_2 + (dTnew_2-dTnew_3)/(Z_2-Z_3)*A_3 ] + RHS - ! dTnew_2 *h^(n+0.5) = -dt * [ (dTnew_1-dTnew_2)/(Z_1-Z_2)*A_2/V_2 + (dTnew_2-dTnew_3)/(Z_2-Z_3)*A_3/V_2 ] + RHS - ! | | - ! v v - ! diffusive flux towards diffusive flux towards - ! T_2 trough boundary V2 T_2 trough boundary V3 - ! - ! --> solve coefficents for homogene part - ! dTnew_2 *h^(n+0.5) = -dt * [ a*dTnew_1 + b*dTnew_2 + c*dTnew_3 ] + ! --> h^(n+0.5)* (T^(n+0.5)-Tstar) = dt*( K_33*d/dz*(T^(n+0.5)-Tstar) + K_33*d/dz*Tstar ) + ! --> Tnew = T^(n+0.5)-Tstar + ! --> h^(n+0.5)* (Tnew) = dt*(K_33*d/dz*Tnew) + K_33*dt*d/dz*Tstar + ! --> h^(n+0.5)* (Tnew) = dt*(K_33*d/dz*Tnew) + RHS + ! --> solve for T_new + ! --> V_1 (Skalar Volume), A_1 (Area of edge), . + ! no Cavity A1==V1, yes Cavity A1 !=V1 /I\ nvec_up (+1) + ! I + ! ----------- zbar_1, A_1 *----I----* + ! Z_1 o T_1, V1 |\ I ./| + ! ----------- zbar_2, A_2 | \ ./ | Gaus Theorem: + ! Z_2 o T_2, V2 | \ / | --> Flux form + ! ----------- zbar_3, A_3 | | | --> normal vec outwards facing + ! Z_3 o T_3, V3 *---|-----* + ! ----------- zbar_4 \ | I ./ + ! : \ | I/ + ! \|/I + ! * I + ! \I/ + ! * nvec_dwn (-1) + ! --> 1st. solve homogenouse part: + ! f(Tnew) = h^(n+0.5)* (Tnew) - dt*(K_33*dTnew/dz) = 0 ! - ! --> a = -dt*K_33/(Z_1-Z_2)*A_2/V_2 + ! --> 2nd. Difference Quotient at Tnew_i in flux form (Gaus Theorem, dont forget normal vectors!!!): + ! V_i*Tnew_i *h_i = -dt * [ K_33 * (Tnew_i-1 - Tnew_i)/(Z_i-1 - Z_i) * A_i * nvec_up + ! +K_33 * (Tnew_i - Tnew_i+1)/(Z_i - Z_i+1) * A_i+1 * nvec_dwn ] + ! Tnew_i *h_i = -dt * [ K_33 * (Tnew_i-1 - Tnew_i)/(Z_i-1 - Z_i) * A_i /V_i * nvec_up + ! +K_33 * (Tnew_i - Tnew_i+1)/(Z_i - Z_i+1) * A_i+1/V_i * nvec_dwn ] + ! + ! --> 3rd. solve for coefficents a, b, c: + ! f(Tnew) = [ a*dTnew_i-1 + b*dTnew_i + c*dTnew_i+1 ] + ! + ! df(Tnew)/dTnew_i-1 = a = -dt*K_33/(Z_i-1 - Z_i) * A_i/V_i * (nvec_up =1) + ! + ! df(Tnew)/dTnew_i+1 = c = dt * K_33 * 1/(Z_i - Z_i+1) * A_i+1/V_i * (nvec_dwn=-1) + ! = -dt * K_33 * 1/(Z_i - Z_i+1) * A_i+1/V_i + ! + ! df(Tnew)/dTnew_i = b = h_i + dt*K_33/(Z_i-1 - Z_i) * A_i/V_i * (nvec_up=+1) + ! - dt*K_33/(Z_i - Z_i+1) * A_i+1/V_i * (nvec_dwn=-1) + ! = h_i + dt*K_33/(Z_i-1 - Z_i) * A_i/V_i + ! + dt*K_33/(Z_i - Z_i+1) * A_i+1/V_i + ! = h_i -(a+c) + ! + ! --> 4th. solve inhomogenous part: + ! [ a*dTnew_i-1 + b*dTnew_i + c*dTnew_i+1 ] = RHS/V_i + ! + ! RHS = K_33*dt*d/dz*Tstar + ! + ! --> write as Difference Quotient in flux form + ! RHS/V_i = K_33 * dt * (Tstar_i-1 - Tstar_i)/(Z_i-1 - Z_i) * A_i/V_i * (nvec_up=1) + ! + K_33 * dt * (Tstar_i - Tstar_i+1)/(Z_i - Z_i+1) * A_i+1/V_i * (nvec_dwn=-1) + ! + ! = K_33*dt/(Z_i-1 - Z_i) * A_i/V_i * Tstar_i-1 + ! - K_33*dt/(Z_i-1 - Z_i) * A_i/V_i * Tstar_i + ! - K_33*dt/(Z_i - Z_i+1) * A_i+1/V_i * Tstar_i + ! + K_33*dt/(Z_i - Z_i+1) * A_i+1/V_i * Tstar_i+1 + ! + ! = -a*Tstar_i-1 + (a+c)*Tstar_i - c * Tstar_i+1 + ! |-> b = h_i - (a+c), a+c = h_i-b + ! + ! = -a*Tstar_i-1 - (b-h_i)*Tstar_i - c * Tstar_i+1 + ! + ! --> 5th. solve for Tnew_i --> forward sweep algorithm --> see lower + ! | b_1 c_1 ... | |dTnew_1| + ! | a_2 b_2 c_2 ... | |dTnew_2| + ! | a_3 b_3 c_3 ... | * |dTnew_3| = RHS/V_i + ! | a_4 b_4 c_4 ...| |dTnew_4| + ! | : | | : | + ! + ! --> a = -dt*K_33 / (Z_i-1 - Z_i) * A_i/V_i ! - ! --> c = -dt*K_33/(Z_2-Z_3)*A_3/V_2 + ! --> c = -dt*K_33 / (Z_i - Z_i+1) * A_i+1/V_i ! - ! --> b = h^(n+0.5) -[ dt*K_33/(Z_1-Z_2)*A_2/V_2 + dt*K_33/(Z_2-Z_3)*A_3/V_2 ] = -(a+c) + h^(n+0.5) + ! --> b = h^(n+0.5) -[ dt*K_33/(Z_i-1 - Z_i)*A_i/V_i + dt*K_33/(Z_i - Z_i+1) * A_i+1/V_i ] = -(a+c) + h^(n+0.5) !___________________________________________________________________________ ! loop over local nodes @@ -480,20 +525,16 @@ subroutine diff_ver_part_impl_ale(tr_num, mesh) ! Be carefull here vertical operation have to be done on NEW vertical mesh !!! zbar_n=0.0_WP Z_n=0.0_WP -! zbar_n(nzmax)=zbar(nzmax) zbar_n(nzmax)=zbar_n_bot(n) Z_n(nzmax-1)=zbar_n(nzmax) + hnode_new(nzmax-1,n)/2.0_WP - !!PS do nz=nzmax-1,2,-1 do nz=nzmax-1,nzmin+1,-1 zbar_n(nz) = zbar_n(nz+1) + hnode_new(nz,n) Z_n(nz-1) = zbar_n(nz) + hnode_new(nz-1,n)/2.0_WP end do - !!PS zbar_n(1) = zbar_n(2) + hnode_new(1,n) zbar_n(nzmin) = zbar_n(nzmin+1) + hnode_new(nzmin,n) !_______________________________________________________________________ ! Regular part of coefficients: --> surface layer - !!PS nz=1 nz=nzmin ! 1/dz(nz) @@ -507,25 +548,20 @@ subroutine diff_ver_part_impl_ale(tr_num, mesh) ! layer dependent coefficients for for solving dT(1)/dt+d/dz*K_33*d/dz*T(1) = ... a(nz)=0.0_WP - !!PS c(nz)=-(Kv(2,n)+Ty1)*zinv2*zinv*area(nz+1,n)/area(nz,n) - c(nz)=-(Kv(nz+1,n)+Ty1)*zinv2*zinv*area(nz+1,n)/areasvol(nz,n) + c(nz)=-(Kv(nz+1,n)+Ty1)*zinv2*zinv * (area(nz+1,n)/areasvol(nz,n)) b(nz)=-c(nz)+hnode_new(nz,n) ! ale ! update from the vertical advection --> comes from splitting of vert ! velocity into explicite and implicite contribution if (do_wimpl) then - !!PS v_adv =zinv*area(nz+1,n)/areasvol(nz,n) - !!PS b(nz) =b(nz)+Wvel_i(nz, n)*zinv-min(0._WP, Wvel_i(nz+1, n))*v_adv - !!PS c(nz) =c(nz)-max(0._WP, Wvel_i(nz+1, n))*v_adv - !___________________________________________________________________ ! use brackets when computing ( area(nz ,n)/areasvol(nz,n) ) for ! numerical reasons, to gurante that area/areasvol in case of no ! cavity is ==1.0_WP - v_adv =zinv* ( area(nz ,n)/areasvol(nz,n) ) + v_adv =zinv * ( area(nz ,n)/areasvol(nz,n) ) b(nz) =b(nz)+Wvel_i(nz, n)*v_adv - v_adv =zinv*area(nz+1,n)/areasvol(nz,n) + v_adv =zinv * ( area(nz+1,n)/areasvol(nz,n) ) b(nz) =b(nz)-min(0._WP, Wvel_i(nz+1, n))*v_adv c(nz) =c(nz)-max(0._WP, Wvel_i(nz+1, n))*v_adv end if @@ -534,7 +570,6 @@ subroutine diff_ver_part_impl_ale(tr_num, mesh) !_______________________________________________________________________ ! Regular part of coefficients: --> 2nd...nl-2 layer - !!PS do nz=2, nzmax-2 do nz=nzmin+1, nzmax-2 ! 1/dz(nz) @@ -552,8 +587,8 @@ subroutine diff_ver_part_impl_ale(tr_num, mesh) ! use brackets when computing ( area(nz ,n)/areasvol(nz,n) ) for ! numerical reasons, to gurante that area/areasvol in case of no ! cavity is ==1.0_WP - a(nz)=-(Kv(nz,n) +Ty )*zinv1*zinv* ( area(nz ,n)/areasvol(nz,n) ) - c(nz)=-(Kv(nz+1,n)+Ty1)*zinv2*zinv*area(nz+1,n)/areasvol(nz,n) + a(nz)=-(Kv(nz,n) +Ty )*zinv1*zinv * ( area(nz ,n)/areasvol(nz,n) ) + c(nz)=-(Kv(nz+1,n)+Ty1)*zinv2*zinv * ( area(nz+1,n)/areasvol(nz,n) ) b(nz)=-a(nz)-c(nz)+hnode_new(nz,n) ! backup zinv2 for next depth level @@ -565,15 +600,15 @@ subroutine diff_ver_part_impl_ale(tr_num, mesh) ! use brackets when computing ( area(nz ,n)/areasvol(nz,n) ) for ! numerical reasons, to gurante that area/areasvol in case of no ! cavity is ==1.0_WP - v_adv=zinv* ( area(nz ,n)/areasvol(nz,n) ) + v_adv=zinv * ( area(nz ,n)/areasvol(nz,n) ) a(nz)=a(nz)+min(0._WP, Wvel_i(nz, n))*v_adv b(nz)=b(nz)+max(0._WP, Wvel_i(nz, n))*v_adv !!PS v_adv=v_adv*areasvol(nz+1,n)/areasvol(nz,n) - v_adv=zinv*area(nz+1,n)/areasvol(nz,n) + v_adv=zinv * ( area(nz+1,n)/areasvol(nz,n) ) b(nz)=b(nz)-min(0._WP, Wvel_i(nz+1, n))*v_adv c(nz)=c(nz)-max(0._WP, Wvel_i(nz+1, n))*v_adv end if - end do ! --> do nz=2, nzmax-2 + end do ! --> do nz=nzmin+1, nzmax-2 !_______________________________________________________________________ ! Regular part of coefficients: --> nl-1 layer @@ -582,8 +617,8 @@ subroutine diff_ver_part_impl_ale(tr_num, mesh) zinv=1.0_WP*dt ! no ... /(zbar(nzmax-1)-zbar(nzmax)) because of ale ! calculate isoneutral diffusivity : Kd*s^2 --> K_33 = Kv + Kd*s^2 - Ty= (Z_n(nz-1)-zbar_n(nz)) *zinv1 *slope_tapered(3,nz-1,n)**2*Ki(nz-1,n) + & - (zbar_n(nz)-Z_n(nz)) *zinv1 *slope_tapered(3,nz,n)**2 *Ki(nz,n) + Ty= (Z_n(nz-1) -zbar_n(nz)) * zinv1 * slope_tapered(3,nz-1,n)**2 * Ki(nz-1,n) + & + (zbar_n(nz)-Z_n(nz) ) * zinv1 * slope_tapered(3,nz ,n)**2 * Ki(nz,n) Ty =Ty *isredi ! layer dependent coefficients for for solving dT(nz)/dt+d/dz*K_33*d/dz*T(nz) = ... @@ -616,56 +651,134 @@ subroutine diff_ver_part_impl_ale(tr_num, mesh) ! ! -+--> tr(1) =(a(1)+c(1))*tr_arr(1,n,tr_num)-c(1)*tr_arr(2,n,tr_num) ! |--> a(1)=0 - !!PS nz=1 nz=nzmin dz=hnode_new(nz,n) ! It would be (zbar(nz)-zbar(nz+1)) if not ALE tr(nz)=-(b(nz)-dz)*tr_arr(nz,n,tr_num)-c(nz)*tr_arr(nz+1,n,tr_num) !tr(nz)=c(nz)*(tr_arr(nz,n,tr_num) - tr_arr(nz+1,n,tr_num)) - - ! ******************************************************************* - ! nonlocal transport to the rhs (only T and S currently) _GO_ - ! ******************************************************************* - ! rsss will be used later to compute: - ! 1. the virtual salinity flux - ! 2. the contribution from the nonlocal term in KPP for salinity - if (tr_num==2) then - rsss=ref_sss - if (ref_sss_local) rsss=tr_arr(1,n,2) - end if - - !!PS do nz=2,nzmax-2 do nz=nzmin+1,nzmax-2 dz=hnode_new(nz,n) - tr(nz)=-a(nz)*tr_arr(nz-1,n,tr_num)-(b(nz)-dz)*tr_arr(nz,n,tr_num)-c(nz)*tr_arr(nz+1,n,tr_num) - !tr(nz)=-a(nz)*tr_arr(nz-1,n,tr_num) & - ! -c(nz)*tr_arr(nz+1,n,tr_num) & - ! +(a(nz)+c(nz))*tr_arr(nz,n,tr_num) + tr(nz)= -a(nz) * tr_arr(nz-1,n,tr_num) & + -(b(nz)-dz)* tr_arr(nz,n,tr_num) & + -c(nz) * tr_arr(nz+1,n,tr_num) + !tr(nz)=-a(nz) * tr_arr(nz-1,n,tr_num) & + ! +(a(nz)+c(nz))* tr_arr(nz,n,tr_num) & + ! -c(nz) * tr_arr(nz+1,n,tr_num) - ! ******************************************************************* - ! nonlocal transport to the rhs (only T and S currently) _GO_ - ! ******************************************************************* -!leads to non conservation in 8th digit. needs to be checked! -! if (mix_scheme_nmb==1 .or. mix_scheme_nmb==17) then -! if (tr_num==1) then ! T -! tr(nz)=tr(nz)+(MIN(ghats(nz,n)*Kv(nz,n), 1.0_WP)-MIN(ghats(nz+1,n)*Kv(nz+1,n), 1.0_WP)*area(nz+1,n)/area(nz,n))*heat_flux(n)/vcpw -! elseif (tr_num==2) then ! S -! tr(nz)=tr(nz)-(MIN(ghats(nz,n)*Kv(nz,n), 1.0_WP)-MIN(ghats(nz+1,n)*Kv(nz+1,n), 1.0_WP)*area(nz+1,n)/area(nz,n))*rsss*water_flux(n) -! end if -! end if end do + nz=nzmax-1 dz=hnode_new(nz,n) tr(nz)=-a(nz)*tr_arr(nz-1,n,tr_num)-(b(nz)-dz)*tr_arr(nz,n,tr_num) !tr(nz)=-a(nz)*tr_arr(nz-1,n,tr_num)+a(nz)*tr_arr(nz,n,tr_num) + !_______________________________________________________________________ + ! Add KPP nonlocal fluxes to the rhs (only T and S currently) + ! use here blmc or kpp_oblmixc instead of Kv, since Kv already contains + ! at this point the mixing enhancments from momix, instable + ! mixing or windmixing which are to much for nonlocal + ! transports and lead to instability of the model + if (use_kpp_nonlclflx) then + if (tr_num==2) then + rsss=ref_sss + if (ref_sss_local) rsss=tr_arr(1,n,2) + end if + + !___________________________________________________________________ + ! use fesom1.4 KPP + if (mix_scheme_nmb==1 .or. mix_scheme_nmb==17) then + if (tr_num==1) then ! temperature + ! --> no fluxes to the top out of the surface, no fluxes + ! downwards out of the bottom + !___surface_________________________________________________ + nz = nzmin + tr(nz)=tr(nz) & + +(-MIN(ghats(nz+1,n)*blmc(nz+1,n,2), 1.0_WP)*(area(nz+1,n)/areasvol(nz,n)) & + ) * heat_flux(n) / vcpw * dt + !___bulk____________________________________________________ + do nz=nzmin+1, nzmax-2 + tr(nz)=tr(nz) & + +( MIN(ghats(nz ,n)*blmc(nz ,n,2), 1.0_WP)*(area(nz ,n)/areasvol(nz,n)) & + -MIN(ghats(nz+1,n)*blmc(nz+1,n,2), 1.0_WP)*(area(nz+1,n)/areasvol(nz,n)) & + ) * heat_flux(n) / vcpw * dt + end do + !___bottom__________________________________________________ + nz = nzmax-1 + tr(nz)=tr(nz) & + +( MIN(ghats(nz ,n)*blmc(nz ,n,2), 1.0_WP)*(area(nz ,n)/areasvol(nz,n)) & + ) * heat_flux(n) / vcpw * dt + + elseif (tr_num==2) then ! salinity + ! --> no fluxes to the top out of the surface, no fluxes + ! downwards out of the bottom + !___surface_________________________________________________ + nz = nzmin + tr(nz)=tr(nz) & + -(-MIN(ghats(nz+1,n)*blmc(nz+1,n,3), 1.0_WP)*(area(nz+1,n)/areasvol(nz,n)) & + ) * rsss * water_flux(n) * dt + !___bulk____________________________________________________ + do nz=nzmin+1, nzmax-2 + tr(nz)=tr(nz) & + -( MIN(ghats(nz ,n)*blmc(nz ,n,3), 1.0_WP)*(area(nz ,n)/areasvol(nz,n)) & + -MIN(ghats(nz+1,n)*blmc(nz+1,n,3), 1.0_WP)*(area(nz+1,n)/areasvol(nz,n)) & + ) * rsss * water_flux(n) * dt + end do + !___bottom__________________________________________________ + nz = nzmax-1 + tr(nz)=tr(nz) & + -( MIN(ghats(nz ,n)*blmc(nz ,n,3), 1.0_WP)*(area(nz ,n)/areasvol(nz,n)) & + ) * rsss * water_flux(n) * dt + end if + !___________________________________________________________________ + ! use cvmix KPP + elseif (mix_scheme_nmb==3 .or. mix_scheme_nmb==37) then + if (tr_num==1) then ! temperature + !___surface_________________________________________________ + nz = nzmin + tr(nz)=tr(nz) & + +(-MIN(kpp_nonlcltranspT(nz+1,n)*kpp_oblmixc(nz+1,n,2), 1.0_WP)*(area(nz+1,n)/areasvol(nz,n)) & + ) * heat_flux(n) / vcpw * dt + !___bulk____________________________________________________ + do nz=nzmin+1, nzmax-2 + tr(nz)=tr(nz) & + +( MIN(kpp_nonlcltranspT(nz ,n)*kpp_oblmixc(nz ,n,2), 1.0_WP)*(area(nz ,n)/areasvol(nz,n)) & + -MIN(kpp_nonlcltranspT(nz+1,n)*kpp_oblmixc(nz+1,n,2), 1.0_WP)*(area(nz+1,n)/areasvol(nz,n)) & + ) * heat_flux(n) / vcpw * dt + end do + !___bottom__________________________________________________ + nz = nzmax-1 + tr(nz)=tr(nz) & + +( MIN(kpp_nonlcltranspT(nz ,n)*kpp_oblmixc(nz ,n,2), 1.0_WP)*(area(nz ,n)/areasvol(nz,n)) & + ) * heat_flux(n) / vcpw * dt + + elseif (tr_num==2) then ! salinity + !___surface_________________________________________________ + nz = nzmin + tr(nz)=tr(nz) & + -(-MIN(kpp_nonlcltranspS(nz+1,n)*kpp_oblmixc(nz+1,n,3), 1.0_WP)*(area(nz+1,n)/areasvol(nz,n)) & + ) * rsss * water_flux(n) * dt + !___bulk____________________________________________________ + do nz=nzmin+1, nzmax-2 + tr(nz)=tr(nz) & + -( MIN(kpp_nonlcltranspS(nz ,n)*kpp_oblmixc(nz ,n,3), 1.0_WP)*(area(nz ,n)/areasvol(nz,n)) & + -MIN(kpp_nonlcltranspS(nz+1,n)*kpp_oblmixc(nz+1,n,3), 1.0_WP)*(area(nz+1,n)/areasvol(nz,n)) & + ) * rsss * water_flux(n) * dt + end do + !___bottom__________________________________________________ + nz = nzmax-1 + tr(nz)=tr(nz) & + -( MIN(kpp_nonlcltranspS(nz ,n)*kpp_oblmixc(nz ,n,3), 1.0_WP)*(area(nz ,n)/areasvol(nz,n)) & + ) * rsss * water_flux(n) * dt + end if + end if + end if ! --> if (use_kpp_nonlclflx) then + !_______________________________________________________________________ ! case of activated shortwave penetration into the ocean, ad 3d contribution if (use_sw_pene .and. tr_num==1) then - !!PS do nz=1, nzmax-1 do nz=nzmin, nzmax-1 zinv=1.0_WP*dt !/(zbar(nz)-zbar(nz+1)) ale! - tr(nz)=tr(nz)+(sw_3d(nz, n)-sw_3d(nz+1, n)*area(nz+1,n)/areasvol(nz,n))*zinv + tr(nz)=tr(nz)+(sw_3d(nz, n)-sw_3d(nz+1, n) * ( area(nz+1,n)/areasvol(nz,n)) ) * zinv end do end if @@ -681,7 +794,6 @@ subroutine diff_ver_part_impl_ale(tr_num, mesh) ! (BUT CHECK!) | | | | ! v (+) v (+) ! - !!PS tr(1)= tr(1)+bc_surface(n, tracer_id(tr_num)) tr(nzmin)= tr(nzmin)+bc_surface(n, tracer_id(tr_num),mesh) !_______________________________________________________________________ @@ -705,13 +817,10 @@ subroutine diff_ver_part_impl_ale(tr_num, mesh) ! --> dTnew_i = rhs'_i-c'_i*dTnew_i+1 ; i = n-1,n-2,...,1 ! ! initialize c-prime and s,t-prime - !!PS cp(1) = c(1)/b(1) - !!PS tp(1) = tr(1)/b(1) cp(nzmin) = c(nzmin)/b(nzmin) tp(nzmin) = tr(nzmin)/b(nzmin) ! solve for vectors c-prime and t, s-prime - !!PS do nz = 2,nzmax-1 do nz = nzmin+1,nzmax-1 m = b(nz)-cp(nz-1)*a(nz) cp(nz) = c(nz)/m @@ -722,7 +831,6 @@ subroutine diff_ver_part_impl_ale(tr_num, mesh) tr(nzmax-1) = tp(nzmax-1) ! solve for x from the vectors c-prime and d-prime - !!PS do nz = nzmax-2, 1, -1 do nz = nzmax-2, nzmin, -1 tr(nz) = tp(nz)-cp(nz)*tr(nz+1) end do @@ -730,12 +838,10 @@ subroutine diff_ver_part_impl_ale(tr_num, mesh) !_______________________________________________________________________ ! update tracer ! tr ... dTnew = T^(n+0.5) - T* - !!PS do nz=1,nzmax-1 do nz=nzmin,nzmax-1 ! tr_arr - before ... T* tr_arr(nz,n,tr_num)=tr_arr(nz,n,tr_num)+tr(nz) ! tr_arr - after ... T^(n+0.5) = dTnew + T* = T^(n+0.5) - T* + T* - end do end do ! --> do n=1,myDim_nod2D diff --git a/src/oce_modules.F90 b/src/oce_modules.F90 index 75d40868c..d2e93315b 100755 --- a/src/oce_modules.F90 +++ b/src/oce_modules.F90 @@ -145,6 +145,10 @@ MODULE o_PARAM real(kind=WP) :: density_ref_T = 2.0_WP real(kind=WP) :: density_ref_S = 34.0_WP +!_______________________________________________________________________________ +! use k-profile nonlocal fluxes +logical :: use_kpp_nonlclflx = .false. + !_______________________________________________________________________________ ! *** active tracer cutoff logical :: limit_salinity=.true. !set an allowed range for salinity @@ -182,7 +186,8 @@ MODULE o_PARAM use_momix, momix_lat, momix_kv, & use_instabmix, instabmix_kv, & use_windmix, windmix_kv, windmix_nl, & - smooth_bh_tra, gamma0_tra, gamma1_tra, gamma2_tra + smooth_bh_tra, gamma0_tra, gamma1_tra, gamma2_tra, & + use_kpp_nonlclflx END MODULE o_PARAM !========================================================== From ee76e29f2f0fc160c71b9d51b3bcaf327357352c Mon Sep 17 00:00:00 2001 From: Patrick Date: Tue, 18 May 2021 16:07:44 +0200 Subject: [PATCH 13/37] revert total back total surface stress so github testcase does not fail --- src/oce_ale_mixing_kpp.F90 | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/src/oce_ale_mixing_kpp.F90 b/src/oce_ale_mixing_kpp.F90 index 14b90ccfe..fd0d3335c 100755 --- a/src/oce_ale_mixing_kpp.F90 +++ b/src/oce_ale_mixing_kpp.F90 @@ -338,9 +338,8 @@ subroutine oce_mixing_KPP(viscAE, diffK, mesh) DO node=1, myDim_nod2D !+eDim_nod2D nzmin = ulevels_nod2D(node) -!!PS ustar(node) = sqrt( sqrt( stress_atmoce_x(node)**2 + stress_atmoce_y(node)**2 )*density_0_r ) ! @ the surface (eqn. 2) - - ustar(node) = sqrt( sqrt( stress_node_surf(1,node)**2 + stress_node_surf(2,node)**2 )*density_0_r ) ! @ the surface (eqn. 2) + ustar(node) = sqrt( sqrt( stress_atmoce_x(node)**2 + stress_atmoce_y(node)**2 )*density_0_r ) ! @ the surface (eqn. 2) +!!PS ustar(node) = sqrt( sqrt( stress_node_surf(1,node)**2 + stress_node_surf(2,node)**2 )*density_0_r ) ! @ the surface (eqn. 2) ! Surface buoyancy forcing (eqns. A2c & A2d & A3b & A3d) !!PS Bo(node) = -g * ( sw_alpha(1,node) * heat_flux(node) / vcpw & !heat_flux & water_flux: positive up From 3ce19489ca404338c4fa7c120c078d0b353b7a9d Mon Sep 17 00:00:00 2001 From: Patrick Date: Wed, 19 May 2021 13:53:11 +0200 Subject: [PATCH 14/37] revert back some of the ( area(nz+1,n)/areasvol(nz,n) ) brackets for numerical reasons so that github testcase does not fail but with brackets would be most likely better. Add them add later point --- src/oce_ale_tracer.F90 | 15 ++++++++++----- 1 file changed, 10 insertions(+), 5 deletions(-) diff --git a/src/oce_ale_tracer.F90 b/src/oce_ale_tracer.F90 index ee689cb26..fd4378cbd 100644 --- a/src/oce_ale_tracer.F90 +++ b/src/oce_ale_tracer.F90 @@ -548,7 +548,8 @@ subroutine diff_ver_part_impl_ale(tr_num, mesh) ! layer dependent coefficients for for solving dT(1)/dt+d/dz*K_33*d/dz*T(1) = ... a(nz)=0.0_WP - c(nz)=-(Kv(nz+1,n)+Ty1)*zinv2*zinv * (area(nz+1,n)/areasvol(nz,n)) + !!PS c(nz)=-(Kv(nz+1,n)+Ty1)*zinv2*zinv * (area(nz+1,n)/areasvol(nz,n)) + c(nz)=-(Kv(nz+1,n)+Ty1)*zinv2*zinv * area(nz+1,n)/areasvol(nz,n) b(nz)=-c(nz)+hnode_new(nz,n) ! ale ! update from the vertical advection --> comes from splitting of vert @@ -561,7 +562,8 @@ subroutine diff_ver_part_impl_ale(tr_num, mesh) v_adv =zinv * ( area(nz ,n)/areasvol(nz,n) ) b(nz) =b(nz)+Wvel_i(nz, n)*v_adv - v_adv =zinv * ( area(nz+1,n)/areasvol(nz,n) ) + !!PS v_adv =zinv * ( area(nz+1,n)/areasvol(nz,n) ) + v_adv =zinv * area(nz+1,n)/areasvol(nz,n) b(nz) =b(nz)-min(0._WP, Wvel_i(nz+1, n))*v_adv c(nz) =c(nz)-max(0._WP, Wvel_i(nz+1, n))*v_adv end if @@ -588,7 +590,8 @@ subroutine diff_ver_part_impl_ale(tr_num, mesh) ! numerical reasons, to gurante that area/areasvol in case of no ! cavity is ==1.0_WP a(nz)=-(Kv(nz,n) +Ty )*zinv1*zinv * ( area(nz ,n)/areasvol(nz,n) ) - c(nz)=-(Kv(nz+1,n)+Ty1)*zinv2*zinv * ( area(nz+1,n)/areasvol(nz,n) ) + !!PS c(nz)=-(Kv(nz+1,n)+Ty1)*zinv2*zinv * ( area(nz+1,n)/areasvol(nz,n) ) + c(nz)=-(Kv(nz+1,n)+Ty1)*zinv2*zinv * area(nz+1,n)/areasvol(nz,n) b(nz)=-a(nz)-c(nz)+hnode_new(nz,n) ! backup zinv2 for next depth level @@ -604,7 +607,8 @@ subroutine diff_ver_part_impl_ale(tr_num, mesh) a(nz)=a(nz)+min(0._WP, Wvel_i(nz, n))*v_adv b(nz)=b(nz)+max(0._WP, Wvel_i(nz, n))*v_adv !!PS v_adv=v_adv*areasvol(nz+1,n)/areasvol(nz,n) - v_adv=zinv * ( area(nz+1,n)/areasvol(nz,n) ) + !!PS v_adv=zinv * ( area(nz+1,n)/areasvol(nz,n) ) + v_adv=zinv * area(nz+1,n)/areasvol(nz,n) b(nz)=b(nz)-min(0._WP, Wvel_i(nz+1, n))*v_adv c(nz)=c(nz)-max(0._WP, Wvel_i(nz+1, n))*v_adv end if @@ -778,7 +782,8 @@ subroutine diff_ver_part_impl_ale(tr_num, mesh) if (use_sw_pene .and. tr_num==1) then do nz=nzmin, nzmax-1 zinv=1.0_WP*dt !/(zbar(nz)-zbar(nz+1)) ale! - tr(nz)=tr(nz)+(sw_3d(nz, n)-sw_3d(nz+1, n) * ( area(nz+1,n)/areasvol(nz,n)) ) * zinv + !!PS tr(nz)=tr(nz)+(sw_3d(nz, n)-sw_3d(nz+1, n) * ( area(nz+1,n)/areasvol(nz,n)) ) * zinv + tr(nz)=tr(nz)+(sw_3d(nz, n)-sw_3d(nz+1, n) * area(nz+1,n)/areasvol(nz,n)) * zinv end do end if From 461fb1b4903220f5b836adb2fb9abff4b674e9ea Mon Sep 17 00:00:00 2001 From: Jan Hegewald Date: Mon, 28 Jun 2021 15:49:18 +0200 Subject: [PATCH 15/37] - do not create uninitialized copies of variables for threads via OpenMP - explicitly disable OpenMP compiling for the Cray ftn compiler --- src/CMakeLists.txt | 2 +- src/gen_modules_partitioning.F90 | 6 ------ 2 files changed, 1 insertion(+), 7 deletions(-) diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index a04db6736..5a0417889 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -75,7 +75,7 @@ elseif(${CMAKE_Fortran_COMPILER_ID} STREQUAL GNU ) 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}) diff --git a/src/gen_modules_partitioning.F90 b/src/gen_modules_partitioning.F90 index a914a9d51..770229964 100644 --- a/src/gen_modules_partitioning.F90 +++ b/src/gen_modules_partitioning.F90 @@ -72,12 +72,6 @@ module g_PARSUP integer, allocatable :: remPtr_elem2D(:), remList_elem2D(:) logical :: elem_full_flag -!$OMP threadprivate(com_nod2D,com_elem2D,com_elem2D_full) -!$OMP threadprivate(mype) -!$OMP threadprivate(myDim_nod2D, eDim_nod2D, myList_nod2D) -!$OMP threadprivate(myDim_elem2D, eDim_elem2D, eXDim_elem2D, myList_elem2D) -!$OMP threadprivate(myDim_edge2D, eDim_edge2D, myList_edge2D) - contains subroutine par_init ! initializes MPI From d9943489853bf51f511c944daa60aa650cb38c43 Mon Sep 17 00:00:00 2001 From: Jan Hegewald Date: Thu, 8 Jul 2021 15:08:21 +0200 Subject: [PATCH 16/37] give default values to uninitialized variables --- src/oce_ale_tracer.F90 | 3 +++ src/oce_ice_init_state.F90 | 4 ++++ src/oce_vel_rhs_vinv.F90 | 1 + 3 files changed, 8 insertions(+) diff --git a/src/oce_ale_tracer.F90 b/src/oce_ale_tracer.F90 index fd4378cbd..0df3b0eb6 100644 --- a/src/oce_ale_tracer.F90 +++ b/src/oce_ale_tracer.F90 @@ -342,6 +342,9 @@ subroutine diff_ver_part_expl_ale(tr_num, mesh) real(kind=WP) :: zinv1,Ty #include "associate_mesh.h" + + Ty = 0.0_WP + !___________________________________________________________________________ do n=1, myDim_nod2D nl1=nlevels_nod2D(n)-1 diff --git a/src/oce_ice_init_state.F90 b/src/oce_ice_init_state.F90 index ecb31fd31..a637515ad 100755 --- a/src/oce_ice_init_state.F90 +++ b/src/oce_ice_init_state.F90 @@ -372,6 +372,8 @@ subroutine init_fields_na_test(mesh) #include "associate_mesh.h" + c_status = .false. + ! =================== ! Fill the model fields with dummy values ! =================== @@ -478,6 +480,8 @@ subroutine init_fields_global_test(mesh) #include "associate_mesh.h" + c_status = .false. + ! =================== ! Fill the model fields with dummy values ! =================== diff --git a/src/oce_vel_rhs_vinv.F90 b/src/oce_vel_rhs_vinv.F90 index 849b5aea9..46881e065 100755 --- a/src/oce_vel_rhs_vinv.F90 +++ b/src/oce_vel_rhs_vinv.F90 @@ -120,6 +120,7 @@ subroutine compute_vel_rhs_vinv(mesh) !vector invariant real(kind=WP) :: density0_inv = 1./density_0 #include "associate_mesh.h" + w = 0.0_WP uvert=0.0_WP From 522d46c4931c17d124fc513878b16426a9675d7c Mon Sep 17 00:00:00 2001 From: suvarchal Date: Tue, 13 Jul 2021 03:48:11 +0200 Subject: [PATCH 17/37] fixes #152 by adding lon, lat to fesom.mesh.diag.nc --- src/io_mesh_info.F90 | 13 ++++++++++++- 1 file changed, 12 insertions(+), 1 deletion(-) diff --git a/src/io_mesh_info.F90 b/src/io_mesh_info.F90 index 551c063e7..7378332cc 100644 --- a/src/io_mesh_info.F90 +++ b/src/io_mesh_info.F90 @@ -4,6 +4,7 @@ module io_mesh_info use g_config use g_comm_auto use o_ARRAYS +use o_PARAM implicit none #include "netcdf.inc" @@ -50,6 +51,7 @@ subroutine write_mesh_info(mesh) integer :: zbar_e_bot_id,zbar_n_bot_id integer :: elem_id integer :: nod_id + integer :: lon_id, lat_id character(100) :: longname character(2000) :: filename real(kind=WP), allocatable :: rbuffer(:), lrbuffer(:) @@ -88,7 +90,10 @@ subroutine write_mesh_info(mesh) call my_def_var(ncid, 'nod_part', NF_INT, 1, (/nod_n_id/), nod_part_id, 'nodal partitioning at the cold start' ) call my_def_var(ncid, 'elem_part', NF_INT, 1, (/elem_n_id/), elem_part_id, 'element partitioning at the cold start') call my_def_var(ncid, 'zbar_e_bottom', NF_DOUBLE, 1, (/elem_n_id/), zbar_e_bot_id, 'element bottom depth') - call my_def_var(ncid, 'zbar_n_bottom', NF_DOUBLE, 1, (/nod_n_id/) , zbar_n_bot_id, 'nodal bottom depth') + call my_def_var(ncid, 'zbar_n_bottom', NF_DOUBLE, 1, (/nod_n_id/), zbar_n_bot_id, 'nodal bottom depth') + call my_def_var(ncid, 'lon', NF_DOUBLE, 1, (/nod_n_id/), lon_id, 'longitude') + call my_def_var(ncid, 'lat', NF_DOUBLE, 1, (/nod_n_id/), lat_id, 'latitude') + ! 2D call my_def_var(ncid, 'nod_area', NF_DOUBLE, 2, (/nod_n_id, nl_id/), nod_area_id, 'nodal areas' ) call my_def_var(ncid, 'elements', NF_INT, 2, (/elem_n_id, id_3/), elem_id, 'elements' ) @@ -185,7 +190,13 @@ subroutine write_mesh_info(mesh) allocate(rbuffer(nod2D)) do i=1, 2 call gather_nod(geo_coord_nod2D(i, 1:myDim_nod2D), rbuffer) + rbuffer = rbuffer/rad call my_put_vara(ncid, nod_id, (/1, i/), (/nod2D, 1/), rbuffer) + if (i == 1) then + call my_put_vara(ncid, lon_id, 1, nod2D, rbuffer) + else + call my_put_vara(ncid, lat_id, 1, nod2D, rbuffer) + endif end do deallocate(rbuffer) From c0b96bfdae9a9f33c8c17bbbcf951edcde570b08 Mon Sep 17 00:00:00 2001 From: Patrick Date: Thu, 22 Jul 2021 16:33:16 +0200 Subject: [PATCH 18/37] when no snow file is given, rewrite prec_rain and prec_snow array, so that freshwater balancing adds up --- config/namelist.forcing.ncep2 | 56 +++++++++++++++++++++++++++++++++++ src/ice_thermo_oce.F90 | 10 ++++++- 2 files changed, 65 insertions(+), 1 deletion(-) create mode 100644 config/namelist.forcing.ncep2 diff --git a/config/namelist.forcing.ncep2 b/config/namelist.forcing.ncep2 new file mode 100644 index 000000000..925dd1828 --- /dev/null +++ b/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' +/ diff --git a/src/ice_thermo_oce.F90 b/src/ice_thermo_oce.F90 index 9e12224dc..d4b6896b3 100755 --- a/src/ice_thermo_oce.F90 +++ b/src/ice_thermo_oce.F90 @@ -149,6 +149,7 @@ subroutine thermodynamics(mesh) snow=prec_rain(i) endif evap_in=evaporation(i) !evap_in: positive up +!!PS evap_in=0.0_WP else rain = prec_rain(i) snow = prec_snow(i) @@ -195,7 +196,7 @@ subroutine thermodynamics(mesh) net_heat_flux(i) = ehf !positive down evaporation(i) = evap !negative up ice_sublimation(i)= subli - + thdgr(i) = ithdgr thdgrsn(i) = ithdgrsn flice(i) = iflice @@ -206,6 +207,13 @@ subroutine thermodynamics(mesh) ! real salt flux due to salinity that is contained in the sea ice 4-5 psu real_salt_flux(i)= rsf !PS + ! if snow file is not given snow computed from prec_rain --> but prec_snow + ! array needs to be filled --> so that the freshwater balancing adds up + if (.not. l_snow) then + prec_rain(i) = rain + prec_snow(i) = snow + end if + end do deallocate(ustar_aux) end subroutine thermodynamics From 8205239a273ba536edf4148a7da9b28cabbc9d35 Mon Sep 17 00:00:00 2001 From: Patrick Date: Thu, 22 Jul 2021 16:35:08 +0200 Subject: [PATCH 19/37] add another flag to overwrite (switch off) the fleapyears calendar checking --- src/gen_modules_config.F90 | 3 ++- src/gen_surface_forcing.F90 | 2 +- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/src/gen_modules_config.F90 b/src/gen_modules_config.F90 index fa4051421..f265ea898 100755 --- a/src/gen_modules_config.F90 +++ b/src/gen_modules_config.F90 @@ -88,7 +88,8 @@ module g_config !_____________________________________________________________________________ ! *** fleap_year *** logical :: include_fleapyear=.false. - namelist /calendar/ include_fleapyear + logical :: use_flpyrcheck =.true. + namelist /calendar/ include_fleapyear, use_flpyrcheck !_____________________________________________________________________________ ! *** machine *** diff --git a/src/gen_surface_forcing.F90 b/src/gen_surface_forcing.F90 index eb9b5eb5d..bfd057638 100644 --- a/src/gen_surface_forcing.F90 +++ b/src/gen_surface_forcing.F90 @@ -354,7 +354,7 @@ SUBROUTINE nc_readTimeGrid(flf) call check_nferr(iost,flf%file_name) ! digg for calendar attribute in time axis variable - if (mype==0) then + if (mype==0 .and. use_flpyrcheck) then iost = nf_inq_attlen(ncid, id_time,'calendar',aux_len) iost = nf_get_att(ncid, id_time,'calendar',aux_calendar) aux_calendar = aux_calendar(1:aux_len) From 0a9d58269053d759641998ebf290e4bb4ff97b64 Mon Sep 17 00:00:00 2001 From: Patrick Date: Thu, 22 Jul 2021 16:36:45 +0200 Subject: [PATCH 20/37] comment NaN checking in gen_support.F90, smooth_nod3d routine --- src/gen_support.F90 | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/src/gen_support.F90 b/src/gen_support.F90 index c420f3de3..114a4db92 100644 --- a/src/gen_support.F90 +++ b/src/gen_support.F90 @@ -126,15 +126,15 @@ subroutine smooth_nod3D(arr, N_smooth, mesh) nln = min(nlev,nlevels_nod2d(n)) DO nz=uln,nln arr(nz, n) = work_array(nz, n) *vol(nz,n) - if (arr(nz,n)/=arr(nz,n)) then - write(*,*) ' --> found NaN in smoothing' - write(*,*) ' mype = ', mype - write(*,*) ' n = ', n - write(*,*) ' nz,uln,nln = ', nz,uln,nln - write(*,*) ' arr(nz,n) = ', arr(nz,n) - write(*,*) ' work_array(nz,n)= ', work_array(nz,n) - write(*,*) ' vol(nz,n) = ', vol(nz,n) - endif +!!PS if (arr(nz,n)/=arr(nz,n)) then +!!PS write(*,*) ' --> found NaN in smoothing' +!!PS write(*,*) ' mype = ', mype +!!PS write(*,*) ' n = ', n +!!PS write(*,*) ' nz,uln,nln = ', nz,uln,nln +!!PS write(*,*) ' arr(nz,n) = ', arr(nz,n) +!!PS write(*,*) ' work_array(nz,n)= ', work_array(nz,n) +!!PS write(*,*) ' vol(nz,n) = ', vol(nz,n) +!!PS endif END DO end DO From 98a69cd63f4183ce54cc2cfcb7e501dee1aa5c90 Mon Sep 17 00:00:00 2001 From: Patrick Date: Thu, 22 Jul 2021 16:37:57 +0200 Subject: [PATCH 21/37] correct namelist.forcing.ncep2 to use NCEP 2 Reanalysis forcing --- config/namelist.forcing.ncep2 | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/config/namelist.forcing.ncep2 b/config/namelist.forcing.ncep2 index 925dd1828..7dcc0ef59 100644 --- a/config/namelist.forcing.ncep2 +++ b/config/namelist.forcing.ncep2 @@ -25,22 +25,22 @@ 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_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_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 From 357894d098b571d9f2b372cdb8867bec0074cb2e Mon Sep 17 00:00:00 2001 From: dsidoren Date: Wed, 28 Jul 2021 16:14:13 +0200 Subject: [PATCH 22/37] Update ice_oce_coupling.F90 --- src/ice_oce_coupling.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/ice_oce_coupling.F90 b/src/ice_oce_coupling.F90 index 7b588503c..2e92d3d42 100755 --- a/src/ice_oce_coupling.F90 +++ b/src/ice_oce_coupling.F90 @@ -329,9 +329,9 @@ subroutine oce_fluxes(mesh) ! due to rigid lid approximation under the cavity we to not add freshwater ! under the cavity for the freshwater balancing we do this only for the open ! ocean - where (ulevels_nod2d == 1) water_flux=water_flux+net/ocean_area + where (ulevels_nod2d == 1) water_flux=water_flux-net/ocean_area else - water_flux=water_flux+net/ocean_area + water_flux=water_flux-net/ocean_area end if !___________________________________________________________________________ From 03194e6b31d8ac6401c2116a4b480ed3f69c59b5 Mon Sep 17 00:00:00 2001 From: dsidoren Date: Wed, 28 Jul 2021 16:24:25 +0200 Subject: [PATCH 23/37] Update io_meandata.F90 --- src/io_meandata.F90 | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/src/io_meandata.F90 b/src/io_meandata.F90 index 0cbc617e9..9fd53e400 100644 --- a/src/io_meandata.F90 +++ b/src/io_meandata.F90 @@ -625,9 +625,10 @@ subroutine create_new_file(entry, mesh) call assert_nf( nf_def_var(entry%ncid, trim(entry%name), entry%data_strategy%netcdf_type(), entry%ndim+1, & (/entry%dimid(1:entry%ndim), entry%recID/), entry%varID), __LINE__) - if (entry%ndim==2) then - call assert_nf( nf_def_var_chunking(entry%ncid, entry%varID, NF_CHUNKED, (/1, entry%glsize(1)/)), __LINE__); - end if +!CHUNKING stuff (netcdf libraries not always compited with it) +!if (entry%ndim==2) then +! call assert_nf( nf_def_var_chunking(entry%ncid, entry%varID, NF_CHUNKED, (/1, entry%glsize(1)/)), __LINE__); +! end if call assert_nf( nf_put_att_text(entry%ncid, entry%varID, 'description', len_trim(entry%description), entry%description), __LINE__) call assert_nf( nf_put_att_text(entry%ncid, entry%varID, 'long_name', len_trim(entry%description), entry%description), __LINE__) call assert_nf( nf_put_att_text(entry%ncid, entry%varID, 'units', len_trim(entry%units), entry%units), __LINE__) From cd083cb4e96fcd2c319eeaddb3470b2e2f071d18 Mon Sep 17 00:00:00 2001 From: dsidoren Date: Wed, 28 Jul 2021 16:40:39 +0200 Subject: [PATCH 24/37] Update ice_oce_coupling.F90 --- src/ice_oce_coupling.F90 | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/src/ice_oce_coupling.F90 b/src/ice_oce_coupling.F90 index 2e92d3d42..760d604af 100755 --- a/src/ice_oce_coupling.F90 +++ b/src/ice_oce_coupling.F90 @@ -289,10 +289,10 @@ subroutine oce_fluxes(mesh) ! enforce the total freshwater/salt flux be zero ! 1. water flux ! if (.not. use_virt_salt) can be used! ! we conserve only the fluxes from the database plus evaporation. - !flux = evaporation-ice_sublimation & ! the ice2atmos subplimation does not contribute to the freshwater flux into the ocean - ! +prec_rain & - ! +prec_snow*(1.0_WP-a_ice_old) & - ! +runoff + flux = evaporation-ice_sublimation & ! the ice2atmos subplimation does not contribute to the freshwater flux into the ocean + +prec_rain & + +prec_snow*(1.0_WP-a_ice_old) & + +runoff ! --> In case of zlevel and zstar and levitating sea ice, sea ice is just sitting ! on top of the ocean without displacement of water, there the thermodynamic ! growth rates of sea ice have to be taken into account to preserve the fresh water @@ -303,7 +303,7 @@ subroutine oce_fluxes(mesh) ! salinity flux !!PS if ( .not. use_floatice .and. .not. use_virt_salt) then if (.not. use_virt_salt) then - flux = water_flux+thdgr*rhoice*inv_rhowat+thdgrsn*rhosno*inv_rhowat + flux = flux-thdgr*rhoice*inv_rhowat-thdgrsn*rhosno*inv_rhowat end if ! Also balance freshwater flux that come from ocean-cavity boundary @@ -329,9 +329,9 @@ subroutine oce_fluxes(mesh) ! due to rigid lid approximation under the cavity we to not add freshwater ! under the cavity for the freshwater balancing we do this only for the open ! ocean - where (ulevels_nod2d == 1) water_flux=water_flux-net/ocean_area + where (ulevels_nod2d == 1) water_flux=water_flux+net/ocean_area else - water_flux=water_flux-net/ocean_area + water_flux=water_flux+net/ocean_area end if !___________________________________________________________________________ From 233d72b62585c490e27c1212437d0ba283754169 Mon Sep 17 00:00:00 2001 From: Dmitry Sidorenko Date: Wed, 28 Jul 2021 14:12:03 +0200 Subject: [PATCH 25/37] ignore the error message from reading the restart if the variable is not in the file (i.e. the model should not stop but use the values with which the variable was initialised initially) --- src/io_restart.F90 | 24 ++++++++++++++++++------ 1 file changed, 18 insertions(+), 6 deletions(-) diff --git a/src/io_restart.F90 b/src/io_restart.F90 index 5f6d0d674..02ae07791 100644 --- a/src/io_restart.F90 +++ b/src/io_restart.F90 @@ -29,6 +29,7 @@ MODULE io_RESTART integer :: ndim integer :: dims(2) !<=2; assume there are no variables with dimension more than 2xNLxT real(kind=WP), pointer :: pt1(:), pt2(:,:) + logical :: is_in_restart !if the variable is in the restart file end type nc_vars ! !-------------------------------------------------------------------------------------------- @@ -560,9 +561,9 @@ subroutine read_restart(id, mesh, arg) integer, optional, intent(in) :: arg real(kind=WP), allocatable :: aux(:), laux(:) integer :: i, lev, size1, size2, size_gen, size_lev, shape - integer :: rec2read, c, order + integer :: rec2read, c, order, ierror real(kind=WP) :: rtime !timestamp of the record - logical :: file_exist=.False. + logical :: file_exist=.False., var_exist type(t_mesh), intent(in) , target :: mesh #include "associate_mesh.h" @@ -605,10 +606,17 @@ subroutine read_restart(id, mesh, arg) end if call was_error(id); c=1 - + do i=1, id%nvar shape=id%var(i)%ndim - if (mype==0) write(*,*) 'reading restart for ', trim(id%var(i)%name) + var_exist=id%var(i)%is_in_restart + call MPI_BCast(var_exist, 1, MPI_INTEGER, 0, MPI_COMM_FESOM, ierror) + if (var_exist) then + if (mype==0) write(*,*) 'reading restart for ', trim(id%var(i)%name) + else + if (mype==0) write(*,*) '...skip reading for ', trim(id%var(i)%name) + cycle + end if !_______writing 2D fields________________________________________________ if (shape==1) then size1=id%var(i)%dims(1) @@ -681,7 +689,7 @@ subroutine assoc_ids(id) type(nc_file), intent(inout) :: id character(500) :: longname - integer :: c, j, k + integer :: c, j, k, status real(kind=WP) :: rtime !timestamp of the record ! Serial output implemented so far if (mype/=0) return @@ -730,7 +738,11 @@ subroutine assoc_ids(id) id%rec_count=max(id%rec_count, 1) !___Associate physical variables____________________________________________ do j=1, id%nvar - id%error_status(c) = nf_inq_varid(id%ncid, id%var(j)%name, id%var(j)%code); c=c+1 + status = nf_inq_varid(id%ncid, id%var(j)%name, id%var(j)%code); c=c+1 + id%var(j)%is_in_restart=(status .eq. nf_noerr) !does the requested variable exist in the file + if (.not. id%var(j)%is_in_restart) then !if not give the error message + write(*,*) 'WARNING: entry ', trim(id%var(j)%name), ' does not exist in the restart file!' + end if end do id%error_status(c)=nf_close(id%ncid); c=c+1 id%error_count=c-1 From 7a80700545dd9d90783b3e0e575e0a5aa8e796b1 Mon Sep 17 00:00:00 2001 From: patrickscholz Date: Mon, 4 Oct 2021 12:12:19 +0200 Subject: [PATCH 26/37] Update oce_ale.F90 (#170) Fix bug when making restart with cavity from partial cell off --> partial cell on --- src/oce_ale.F90 | 10 ---------- 1 file changed, 10 deletions(-) diff --git a/src/oce_ale.F90 b/src/oce_ale.F90 index d231cbf61..1bc33a1e6 100644 --- a/src/oce_ale.F90 +++ b/src/oce_ale.F90 @@ -1025,11 +1025,6 @@ subroutine restart_thickness_ale(mesh) nzmin = ulevels_nod2D(n) nzmax = nlevels_nod2D(n)-1 - !___________________________________________________________________ - ! if there is a cavity layer thickness is not updated, its - ! kept fixed - if (nzmin > 1) cycle - !___________________________________________________________________ ! be sure that bottom layerthickness uses partial cell layer thickness ! in case its activated, especially when you make a restart from a non @@ -1050,11 +1045,6 @@ subroutine restart_thickness_ale(mesh) nzmin = ulevels(elem) nzmax = nlevels(elem)-1 - !___________________________________________________________________ - ! if there is a cavity layer thickness is not updated, its - ! kept fixed - if (nzmin > 1) cycle - !___________________________________________________________________ elnodes=elem2D_nodes(:, elem) do nz=nzmin,nzmax-1 From 212952ee950cfda037682bcb3f757ab12810eedb Mon Sep 17 00:00:00 2001 From: Claudia Wekerle Date: Thu, 30 Sep 2021 13:16:01 +0200 Subject: [PATCH 27/37] fix bug with reading runoff and SSS --- src/gen_surface_forcing.F90 | 13 ++++++++----- 1 file changed, 8 insertions(+), 5 deletions(-) diff --git a/src/gen_surface_forcing.F90 b/src/gen_surface_forcing.F90 index bfd057638..a58c698ca 100644 --- a/src/gen_surface_forcing.F90 +++ b/src/gen_surface_forcing.F90 @@ -1099,13 +1099,15 @@ SUBROUTINE sbc_do(mesh) !========================================================================== ! prepare a flag which checks whether to update monthly data (SSS, river runoff) - update_monthly_flag=((day_in_month==num_day_in_month(fleapyear,month) .and. timenew==86400._WP)) + update_monthly_flag=( (day_in_month==num_day_in_month(fleapyear,month) .and. + timenew==86400._WP) .or. mstep==1 ) ! read in SSS for applying SSS restoring if (surf_relax_S > 0._WP) then if (sss_data_source=='CORE1' .or. sss_data_source=='CORE2') then if (update_monthly_flag) then - i=month+1 + i=month + if (mstep > 1) i=i+1 if (i > 12) i=1 if (mype==0) write(*,*) 'Updating SSS restoring data for month ', i call read_other_NetCDF(nm_sss_data_file, 'SALT', i, Ssurf, .true., mesh) @@ -1119,7 +1121,8 @@ SUBROUTINE sbc_do(mesh) if(update_monthly_flag) then if(runoff_climatology) then !climatology monthly mean - i=month+1 + i=month + if (mstep > 1) i=i+1 if (i > 12) i=1 if (mype==0) write(*,*) 'Updating monthly climatology runoff for month ', i filename=trim(nm_runoff_file) @@ -1130,8 +1133,8 @@ SUBROUTINE sbc_do(mesh) else !monthly data - - i=month+1 + i=month + if (mstep > 1) i=i+1 if (i > 12) i=1 if (mype==0) write(*,*) 'Updating monthly runoff for month ', i filename=trim(nm_runoff_file)//cyearnew//'.nc' From 554622f3e4efd04005bff31b4ba3565fa22cf953 Mon Sep 17 00:00:00 2001 From: dsidoren Date: Thu, 30 Sep 2021 22:06:12 +0200 Subject: [PATCH 28/37] Update gen_surface_forcing.F90 @koldunovn you were right regarding the line splitting --- src/gen_surface_forcing.F90 | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/gen_surface_forcing.F90 b/src/gen_surface_forcing.F90 index a58c698ca..ad22831e6 100644 --- a/src/gen_surface_forcing.F90 +++ b/src/gen_surface_forcing.F90 @@ -1099,8 +1099,7 @@ SUBROUTINE sbc_do(mesh) !========================================================================== ! prepare a flag which checks whether to update monthly data (SSS, river runoff) - update_monthly_flag=( (day_in_month==num_day_in_month(fleapyear,month) .and. - timenew==86400._WP) .or. mstep==1 ) + update_monthly_flag=( (day_in_month==num_day_in_month(fleapyear,month) .AND. timenew==86400._WP) .OR. mstep==1 ) ! read in SSS for applying SSS restoring if (surf_relax_S > 0._WP) then From f3dbf4830fec8bbf8d6ff2242d7563b3d14ac529 Mon Sep 17 00:00:00 2001 From: Nikolay Koldunov Date: Tue, 5 Oct 2021 16:10:21 +0200 Subject: [PATCH 29/37] Update setup.yml Update tests --- setups/test_pi/setup.yml | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/setups/test_pi/setup.yml b/setups/test_pi/setup.yml index e7f20c760..b38c480f0 100644 --- a/setups/test_pi/setup.yml +++ b/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.2691276443855294 + salt: 23.944024712806094 + temp: 1.701768707848739 + sst: 8.531522995932146 + u: -0.001407225233294229 + v: 0.00014182969591235959 From f94c17bb217bdd79c5a4f7c261f605f125fc316e Mon Sep 17 00:00:00 2001 From: Nikolay Koldunov Date: Tue, 5 Oct 2021 18:08:56 +0200 Subject: [PATCH 30/37] Update fesom2.1.yml Actions are not triggered for whatever reason, try again. --- .github/workflows/fesom2.1.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/fesom2.1.yml b/.github/workflows/fesom2.1.yml index 4facc60cc..6f50d8c61 100644 --- a/.github/workflows/fesom2.1.yml +++ b/.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] From c9ad54dcf0747629668efe6f632622ecd48b607d Mon Sep 17 00:00:00 2001 From: Nikolay Koldunov Date: Tue, 5 Oct 2021 18:18:01 +0200 Subject: [PATCH 31/37] Update of the icepack tests --- setups/test_pi_icepack/setup.yml | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/setups/test_pi_icepack/setup.yml b/setups/test_pi_icepack/setup.yml index 9180b3d59..b7a18cf82 100644 --- a/setups/test_pi_icepack/setup.yml +++ b/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 From 98b11078e28f380c2602a8f57ca74d46ca5d0030 Mon Sep 17 00:00:00 2001 From: Nikolay Koldunov Date: Tue, 5 Oct 2021 18:19:19 +0200 Subject: [PATCH 32/37] Trying to trigger actions again --- .github/workflows/fesom2_icepack.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/fesom2_icepack.yml b/.github/workflows/fesom2_icepack.yml index aeb50481f..d5681da06 100644 --- a/.github/workflows/fesom2_icepack.yml +++ b/.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] From 48d01d6bcd890d57ce2dd4425c13cc07ce32687a Mon Sep 17 00:00:00 2001 From: Julian Rodrigo Berlin Date: Tue, 19 Oct 2021 10:45:40 +0200 Subject: [PATCH 33/37] set more recent compilers and mpi configuration to bsc environment --- env/bsc/shell | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/env/bsc/shell b/env/bsc/shell index a5ccfb77a..0dd072599 100644 --- a/env/bsc/shell +++ b/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 From d3df069bd10a31be9cda497c58ff29e96e1a31ef Mon Sep 17 00:00:00 2001 From: Nikolay Koldunov Date: Wed, 17 Nov 2021 10:29:07 +0100 Subject: [PATCH 34/37] fix KPP bug discovered by Dima. Affects results, so tests are changed as well. --- setups/test_pi/setup.yml | 12 ++++++------ src/oce_ale_mixing_kpp.F90 | 3 ++- 2 files changed, 8 insertions(+), 7 deletions(-) diff --git a/setups/test_pi/setup.yml b/setups/test_pi/setup.yml index b38c480f0..c28aaf340 100644 --- a/setups/test_pi/setup.yml +++ b/setups/test_pi/setup.yml @@ -61,12 +61,12 @@ namelist.io: prec: 8 fcheck: - a_ice: 0.2691276443855294 - salt: 23.944024712806094 - temp: 1.701768707848739 - sst: 8.531522995932146 - u: -0.001407225233294229 - v: 0.00014182969591235959 + a_ice: 0.2691276598479261 + salt: 23.944024679303666 + temp: 1.701768750033021 + sst: 8.531528640978305 + u: -0.0014072137861434184 + v: 0.00014184602459601167 diff --git a/src/oce_ale_mixing_kpp.F90 b/src/oce_ale_mixing_kpp.F90 index fd0d3335c..b7b33ca85 100755 --- a/src/oce_ale_mixing_kpp.F90 +++ b/src/oce_ale_mixing_kpp.F90 @@ -624,7 +624,8 @@ SUBROUTINE bldepth(mesh) !----------------------------------------------------------------------- ! find stability and buoyancy forcing for final hbl values !----------------------------------------------------------------------- - IF (use_sw_pene) THEN + IF (use_sw_pene) THEN + coeff_sw = g * sw_alpha(nzmin,node) ! @ the surface @ Z (m/s2/K) ! Linear interpolation of sw_3d to depth of hbl bfsfc(node) = Bo(node) + & coeff_sw * & From 744f59fb302a00d378ecfe667376c5452058f463 Mon Sep 17 00:00:00 2001 From: Dmitry Sidorenko Date: Mon, 22 Nov 2021 17:11:29 +0100 Subject: [PATCH 35/37] bug fix in update_atm_forcing. press_air array was accessed even if not allocated. --- src/gen_forcing_couple.F90 | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/src/gen_forcing_couple.F90 b/src/gen_forcing_couple.F90 index 72a448b26..4d36ef385 100755 --- a/src/gen_forcing_couple.F90 +++ b/src/gen_forcing_couple.F90 @@ -260,7 +260,9 @@ 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 @@ -273,7 +275,9 @@ subroutine update_atm_forcing(istep, mesh) Tair(i) = 0.0_WP prec_rain(i)= 0.0_WP prec_snow(i)= 0.0_WP - press_air(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 From 060f5b1a0e8c595bb91ff55bfc5896b50606aabe Mon Sep 17 00:00:00 2001 From: dsidoren Date: Sun, 23 Jan 2022 18:59:57 +0100 Subject: [PATCH 36/37] Update oce_ale_tracer.F90 gamma0, gamma1 & gamma2 were used instead gamma0_tra, gamma1_tra & gamma2_tra in diff_part_bh --- src/oce_ale_tracer.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/oce_ale_tracer.F90 b/src/oce_ale_tracer.F90 index 0df3b0eb6..5cfed0a51 100644 --- a/src/oce_ale_tracer.F90 +++ b/src/oce_ale_tracer.F90 @@ -1112,7 +1112,7 @@ SUBROUTINE diff_part_bh(ttf, mesh) v1=UV(2, nz,el(1))-UV(2, nz,el(2)) vi=u1*u1+v1*v1 tt=ttf(nz,en(1))-ttf(nz,en(2)) - vi=sqrt(max(gamma0, max(gamma1*sqrt(vi), gamma2*vi))*len) + vi=sqrt(max(gamma0_tra, max(gamma1_tra*sqrt(vi), gamma2_tra*vi))*len) !vi=sqrt(max(sqrt(u1*u1+v1*v1),0.04)*le) ! 10m^2/s for 10 km (0.04 h/50) !vi=sqrt(10.*le) tt=tt*vi @@ -1136,7 +1136,7 @@ SUBROUTINE diff_part_bh(ttf, mesh) v1=UV(2, nz,el(1))-UV(2, nz,el(2)) vi=u1*u1+v1*v1 tt=temporary_ttf(nz,en(1))-temporary_ttf(nz,en(2)) - vi=sqrt(max(gamma0, max(gamma1*sqrt(vi), gamma2*vi))*len) + vi=sqrt(max(gamma0_tra, max(gamma1_tra*sqrt(vi), gamma2_tra*vi))*len) !vi=sqrt(max(sqrt(u1*u1+v1*v1),0.04)*le) ! 10m^2/s for 10 km (0.04 h/50) !vi=sqrt(10.*le) tt=-tt*vi*dt From 00df069b3baf9aa92786862e9c6b6af655075b37 Mon Sep 17 00:00:00 2001 From: dsidoren Date: Tue, 7 Jun 2022 15:14:55 +0200 Subject: [PATCH 37/37] Update io_meandata.F90 --- src/io_meandata.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/io_meandata.F90 b/src/io_meandata.F90 index 9fd53e400..e04fc6509 100644 --- a/src/io_meandata.F90 +++ b/src/io_meandata.F90 @@ -266,10 +266,10 @@ subroutine ini_mean_io(mesh) end if CASE ('tx_sur ') sel_forcvar(11) = 1 - call def_stream(elem2D, myDim_elem2D, 'tx_sur', 'zonal wind str. to ocean', 'm/s2', stress_surf(1, :), io_list(i)%freq, io_list(i)%unit, io_list(i)%precision, mesh) + call def_stream(elem2D, myDim_elem2D, 'tx_sur', 'total zonal str. to ocean', 'm/s2', stress_surf(1, :), io_list(i)%freq, io_list(i)%unit, io_list(i)%precision, mesh) CASE ('ty_sur ') sel_forcvar(12) = 1 - call def_stream(elem2D, myDim_elem2D, 'ty_sur', 'meridional wind str. to ocean', 'm/s2', stress_surf(2, :), io_list(i)%freq, io_list(i)%unit, io_list(i)%precision, mesh) + call def_stream(elem2D, myDim_elem2D, 'ty_sur', 'total meridional str. to ocean', 'm/s2', stress_surf(2, :), io_list(i)%freq, io_list(i)%unit, io_list(i)%precision, mesh) CASE ('curl_surf ') if (lcurt_stress_surf) then call def_stream(nod2D, myDim_nod2D, 'curl_surf', 'vorticity of the surface stress','none', curl_stress_surf(:), io_list(i)%freq, io_list(i)%unit, io_list(i)%precision, mesh)