From a7bb8760919edafc7c76d044f2c043897e89b5f2 Mon Sep 17 00:00:00 2001 From: Dmitry Sidorenko Date: Mon, 4 May 2020 16:18:29 +0200 Subject: [PATCH 01/26] =?UTF-8?q?Routines=20for=20MOC=20density=20diagnost?= =?UTF-8?q?ic=20have=20been=20updated.=201.=20The=20standart=20=CF=832=20b?= =?UTF-8?q?ins=20are=20chosen=20according=20to=20Megan=20(2018)=20(72=20le?= =?UTF-8?q?vels=20for=20a=20good=20representation=20of=20deep=20and=20bott?= =?UTF-8?q?om=20waters)=20and=20augmented=20with=20density=20levels=20to?= =?UTF-8?q?=20match=20those=20presented=20in=20Xu=20et=20al.=202018.=20Alt?= =?UTF-8?q?ogether=20we=20use=2085=20density=20bins=20spanning=20the=20ran?= =?UTF-8?q?ge=20of=2030.0=E2=80=AF<=E2=80=AF=CF=832=E2=80=AF<=E2=80=AF37.2?= =?UTF-8?q?=E2=80=AF=20kg=20m=E2=88=923.=202.=20New=20fields=20which=20all?= =?UTF-8?q?ow=20to=20compute=20surfaxe=20bouyancy=20flux=20etc.have=20been?= =?UTF-8?q?=20added?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- src/gen_modules_diag.F90 | 246 ++++++++++++++++++++++++++++++--------- src/io_meandata.F90 | 18 ++- 2 files changed, 200 insertions(+), 64 deletions(-) diff --git a/src/gen_modules_diag.F90 b/src/gen_modules_diag.F90 index 7edb88a98..4c912347c 100755 --- a/src/gen_modules_diag.F90 +++ b/src/gen_modules_diag.F90 @@ -18,7 +18,7 @@ module diagnostics public :: ldiag_solver, lcurt_stress_surf, ldiag_energy, ldiag_dMOC, ldiag_DVD, ldiag_forc, ldiag_salt3D, ldiag_curl_vel3, diag_list, & compute_diagnostics, rhs_diag, curl_stress_surf, curl_vel3, wrhof, rhof, & u_x_u, u_x_v, v_x_v, v_x_w, u_x_w, dudx, dudy, dvdx, dvdy, dudz, dvdz, utau_surf, utau_bott, av_dudz_sq, av_dudz, av_dvdz, stress_bott, u_surf, v_surf, u_bott, v_bott, & - std_dens_min, std_dens_max, std_dens_N, std_dens, std_dens_UVDZ, std_dens_RHOZ, & + std_dens_min, std_dens_max, std_dens_N, std_dens, std_dens_UVDZ, std_dens_DIV, std_dens_Z, std_dens_RHOZ, std_dens_dVdT, std_dens_flux, dens_flux, & compute_diag_dvd_2ndmoment_klingbeil_etal_2014, compute_diag_dvd_2ndmoment_burchard_etal_2008, compute_diag_dvd ! Arrays used for diagnostics, some shall be accessible to the I/O ! 1. solver diagnostics: A*x=rhs? @@ -35,23 +35,22 @@ module diagnostics ! defining a set of standard density bins which will be used for computing densMOC ! integer, parameter :: std_dens_N = 100 ! real(kind=WP), save, target :: std_dens(std_dens_N) - integer, parameter :: std_dens_N =72 + integer, parameter :: std_dens_N =89 real(kind=WP), save, target :: std_dens(std_dens_N)=(/ & - 0.0000, 30.00000, 30.55556, 31.11111, 31.66667, 32.22222, & - 32.77778, 33.33333, 33.88889, 34.44444, 35.00000, 35.10622, & - 35.20319, 35.29239, 35.37498, 35.45187, 35.52380, 35.59136, & - 35.65506, 35.71531, 35.77247, 35.82685, 35.87869, 35.92823, & - 35.97566, 36.02115, 36.06487, 36.10692, 36.14746, 36.18656, & - 36.22434, 36.26089, 36.29626, 36.33056, 36.36383, 36.39613, & - 36.42753, 36.45806, 36.48778, 36.51674, 36.54495, 36.57246, & - 36.59932, 36.62555, 36.65117, 36.67621, 36.70071, 36.72467, & - 36.74813, 36.77111, 36.79363, 36.81570, 36.83733, 36.85857, & - 36.87940, 36.89985, 36.91993, 36.93965, 36.95904, 36.97808, & - 36.99682, 37.01524, 37.03336, 37.05119, 37.06874, 37.08602, & - 37.10303, 37.11979, 37.13630, 37.15257, 37.16861, 37.18441/) + 0.0000, 30.00000, 30.55556, 31.11111, 31.36000, 31.66667, 31.91000, 32.22222, 32.46000, & + 32.77778, 33.01000, 33.33333, 33.56000, 33.88889, 34.11000, 34.44444, 34.62000, 35.00000, & + 35.05000, 35.10622, 35.20319, 35.29239, 35.37498, 35.41300, 35.45187, 35.52380, 35.59136, & + 35.65506, 35.71531, 35.77247, 35.82685, 35.87869, 35.92823, 35.97566, 35.98000, 36.02115, & + 36.06487, 36.10692, 36.14746, 36.18656, 36.22434, 36.26089, 36.29626, 36.33056, 36.36383, & + 36.39613, 36.42753, 36.45806, 36.48778, 36.51674, 36.54495, 36.57246, 36.59500, 36.59932, & + 36.62555, 36.65117, 36.67621, 36.68000, 36.70071, 36.72467, 36.74813, 36.75200, 36.77111, & + 36.79363, 36.81570, 36.83733, 36.85857, 36.87500, 36.87940, 36.89985, 36.91993, 36.93965, & + 36.95904, 36.97808, 36.99682, 37.01524, 37.03336, 37.05119, 37.06874, 37.08602, 37.10303, & + 37.11979, 37.13630, 37.15257, 37.16861, 37.18441, 37.50000, 37.75000, 40.00000/) real(kind=WP), save, target :: std_dd(std_dens_N-1) - real(kind=WP), save, target :: std_dens_min=1012., std_dens_max=1039. - real(kind=WP), save, allocatable, target :: std_dens_UVDZ(:,:,:), std_dens_RHOZ(:,:) + real(kind=WP), save, target :: std_dens_min=1030., std_dens_max=1040. + real(kind=WP), save, allocatable, target :: std_dens_UVDZ(:,:,:), std_dens_RHOZ(:,:), std_dens_flux(:,:,:), std_dens_dVdT(:,:), std_dens_DIV(:,:), std_dens_Z(:,:) + real(kind=WP), save, allocatable, target :: dens_flux(:) logical :: ldiag_solver =.false. logical :: lcurt_stress_surf=.false. @@ -346,109 +345,240 @@ end subroutine diag_energy subroutine diag_densMOC(mode, mesh) implicit none integer, intent(in) :: mode + type(t_mesh), intent(in) , target :: mesh integer :: nz, snz, elem, nzmax, elnodes(3), is, ie, pos + integer :: e, edge, enodes(2), eelems(2) + real(kind=WP) :: div, deltaX, deltaY, locz integer :: jj - real(kind=WP), save :: dd - real(kind=WP) :: uvdz_el(2), rhoz_el, dz, weight, dmin, dmax, ddiff, test, test1, test2, test3 - real(kind=WP), save, allocatable :: dens(:), aux(:) - real(kind=WP), save, allocatable :: std_dens_w(:,:) - logical, save :: firstcall=.true. - type(t_mesh), intent(in) , target :: mesh + real(kind=WP), save :: dd, wdiff + real(kind=WP) :: uvdz_el(2), rhoz_el, vol_el, dz, weight, dmin, dmax, ddiff, test, test1, test2, test3 + real(kind=WP), save, allocatable :: dens(:), aux(:), el_depth(:) + real(kind=WP), save, allocatable :: std_dens_w(:,:), std_dens_VOL1(:,:), std_dens_VOL2(:,:) + logical, save :: firstcall_s=.true., firstcall_e=.true. #include "associate_mesh.h" -!===================== - - if (firstcall) then !allocate the stuff at the first call + if (firstcall_s) then !allocate the stuff at the first call allocate(std_dens_UVDZ(2,std_dens_N, myDim_elem2D)) allocate(std_dens_RHOZ( std_dens_N, myDim_elem2D)) allocate(std_dens_w ( std_dens_N, myDim_elem2D)) - allocate(aux(nl-1)) - allocate(dens(nl)) + allocate(std_dens_dVdT( std_dens_N, myDim_elem2D)) + allocate(std_dens_DIV ( std_dens_N, myDim_nod2D+eDim_nod2D)) + allocate(std_dens_VOL1( std_dens_N, myDim_elem2D)) + allocate(std_dens_VOL2( std_dens_N, myDim_elem2D)) + allocate(std_dens_flux(3,std_dens_N, myDim_elem2D)) + allocate(std_dens_Z ( std_dens_N, myDim_elem2D)) + allocate(dens_flux(elem2D)) + allocate(aux (nl-1)) + allocate(dens (nl)) + allocate(el_depth(nl)) ! -!std_dens(1)=27.5 -!do nz=2, std_dens_N -!std_dens(nz)=std_dens(nz-1)+10.5/real(std_dens_N) +!std_dens(1)=0. +!std_dens(2)=30. +!do nz=3, std_dens_N-1 +!std_dens(nz)=std_dens(nz-1)+10.5/real(std_dens_N-2) !end do +!std_dens(std_dens_N)=40. ! std_dd(:)=std_dens(2:)-std_dens(:std_dens_N-1) dens =0. std_dens_UVDZ=0. std_dens_RHOZ=0. - firstcall=.false. + std_dens_dVdT=0. + std_dens_DIV =0. + std_dens_VOL1=0. + std_dens_VOL2=0. + std_dens_Z =0. + depth =0. + el_depth =0. + firstcall_s=.false. if (mode==0) return end if std_dens_UVDZ=0. std_dens_RHOZ=0. std_dens_w =0. + std_dens_flux=0. + dens_flux =0. + std_dens_VOL2=0. + std_dens_DIV =0. + std_dens_Z =0. do elem=1, myDim_elem2D - elnodes=elem2D_nodes(:,elem) + elnodes=elem2D_nodes(:,elem) nzmax =nlevels(elem) + + dens_flux(elem)=sum(sw_alpha(1,elnodes) * heat_flux(elnodes) / vcpw + sw_beta(1,elnodes) * (relax_salt(elnodes) + water_flux(elnodes) * tr_arr(1,elnodes,2)))/3. + do nz=1, nzmax-1 - aux(nz)=sum(density_dmoc(nz, elnodes))/3.-1000. + aux(nz)=sum(density_dmoc(nz, elnodes))/3.-1000. end do + + el_depth(nzmax)=zbar_e_bot(elem) do nz=nzmax-1,2,-1 - dens(nz) = (aux(nz) * helem(nz-1,elem)+& - aux(nz-1) * helem(nz, elem))/sum(helem(nz-1:nz,elem)) + dens(nz) = (aux(nz) * helem(nz-1,elem)+& + aux(nz-1) * helem(nz, elem))/sum(helem(nz-1:nz,elem)) + el_depth(nz) = el_depth(nz+1) + helem(nz, elem) end do dens(nzmax)=dens(nzmax-1)+(dens(nzmax-1)-dens(nzmax-2))*helem(nzmax-1,elem)/helem(nzmax-2,elem) dens(1) =dens(2) +(dens(2)-dens(3)) *helem(1, elem)/helem(2,elem) + el_depth(1)=0. + is=minloc(abs(std_dens-dens(1)),1) +! std_dens_flux(is,elem)=std_dens_flux(is,elem)+dens_flux + std_dens_flux(1, is,elem)=std_dens_flux(1, is,elem)+elem_area(elem)*sum(sw_alpha(1,elnodes) * heat_flux (elnodes ))/3./vcpw + std_dens_flux(2, is,elem)=std_dens_flux(2, is,elem)+elem_area(elem)*sum(sw_beta (1,elnodes) * relax_salt(elnodes ))/3. + std_dens_flux(3, is,elem)=std_dens_flux(3, is,elem)+elem_area(elem)*sum(sw_beta (1,elnodes) * water_flux(elnodes) * tr_arr(1, elnodes, 2))/3. + do nz=nzmax-1,1,-1 - dmin=minval(dens(nz:nz+1)) - dmax=maxval(dens(nz:nz+1)) -! is=findloc(std_dens > dmin, value=.true., dim=1) - is=1 - do jj = 1, std_dens_N - if (std_dens(jj) > dmin) then - is = jj - exit - endif - end do - -! ie=findloc(std_dens < dmax, value=.true., back=.true., dim=1) - ie=std_dens_N - do jj = std_dens_N,1,-1 - if (std_dens(jj) < dmin) then - ie = jj - exit - endif - end do + dmin =minval(dens(nz:nz+1)) + dmax =maxval(dens(nz:nz+1)) + ddiff =abs(dens(nz)-dens(nz+1)) + is=findloc(std_dens > dmin, value=.true., dim=1) +! is=1 +! do jj = 1, std_dens_N +! if (std_dens(jj) > dmin) then +! is = jj +! exit +! endif +! end do + + ie=findloc(std_dens < dmax, value=.true., back=.true., dim=1) +! ie=std_dens_N +! do jj = std_dens_N,1,-1 +! if (std_dens(jj) < dmin) then +! ie = jj +! exit +! endif +! end do if (std_dens(is)>=dmax) is=ie if (std_dens(ie)<=dmin) ie=is uvdz_el=(UV(:,nz,elem)+fer_uv(:,nz,elem))*helem(nz,elem) rhoz_el=(dens(nz)-dens(nz+1))/helem(nz,elem) - ddiff=abs(dens(nz)-dens(nz+1)) + vol_el =helem(nz,elem)*elem_area(elem) + wdiff =sum(wvel(elnodes, nz)-wvel(elnodes, nz+1))/3.*elem_area(elem) if (ie-is > 0) then weight=(std_dens(is)-dmin)+std_dd(is)/2. weight=max(weight, 0.)/ddiff std_dens_UVDZ(:, is, elem)=std_dens_UVDZ(:, is, elem)+weight*uvdz_el std_dens_RHOZ( is, elem)=std_dens_RHOZ( is, elem)+weight*rhoz_el - std_dens_w( is, elem) =std_dens_w( is, elem) +weight + std_dens_VOL2( is, elem)=std_dens_VOL2( is, elem)+weight*vol_el + locz=el_depth(nz+1)+weight*helem(nz,elem) + std_dens_Z ( is, elem)=std_dens_Z ( is, elem)+locz*weight + std_dens_w( is, elem) =std_dens_w( is, elem)+weight do snz=is+1, ie-1 weight=(sum(std_dd(snz-1:snz))/2.)/ddiff std_dens_UVDZ(:, snz, elem)=std_dens_UVDZ(:, snz, elem)+weight*uvdz_el std_dens_RHOZ( snz, elem)=std_dens_RHOZ( snz, elem)+weight*rhoz_el - std_dens_w ( snz, elem)=std_dens_w ( snz, elem)+weight + std_dens_VOL2( snz, elem)=std_dens_VOL2( snz, elem)+weight*vol_el + locz=locz+weight*helem(nz,elem) + std_dens_Z ( snz, elem) =std_dens_Z ( snz, elem)+locz*weight + std_dens_w ( snz, elem) =std_dens_w ( snz, elem)+weight end do weight=(dmax-std_dens(ie))+std_dd(ie-1)/2. weight=max(weight, 0.)/ddiff std_dens_UVDZ(:, ie, elem)=std_dens_UVDZ(:, ie, elem)+weight*uvdz_el std_dens_RHOZ( ie, elem)=std_dens_RHOZ( ie, elem)+weight*rhoz_el - std_dens_w ( ie, elem)=std_dens_w ( ie, elem)+weight + std_dens_VOL2( ie, elem)=std_dens_VOL2( ie, elem)+weight*vol_el + locz=locz+weight*helem(nz,elem) + std_dens_Z ( ie, elem) =std_dens_Z ( ie, elem)+locz*weight + std_dens_w ( ie, elem) =std_dens_w ( ie, elem)+weight else std_dens_UVDZ(:, is, elem)=std_dens_UVDZ(:, is, elem)+uvdz_el std_dens_RHOZ( is, elem)=std_dens_RHOZ( is, elem)+rhoz_el + std_dens_VOL2( is, elem)=std_dens_VOL2( is, elem)+vol_el + std_dens_Z ( is, elem)=std_dens_Z ( is, elem)+el_depth(nz+1)+helem(nz,elem)/2. std_dens_w ( is, elem)=std_dens_w ( is, elem)+1._wp end if end do end do + + do edge=1, myDim_edge2D + if (myList_edge2D(edge) > edge2D_in) cycle + enodes=edges(:,edge) + eelems=edge_tri(:,edge) + nzmax =nlevels(eelems(1)) + if (eelems(2)>0) nzmax=max(nzmax, nlevels(eelems(2))) + do nz=1, nzmax-1 + aux(nz)=sum(density_dmoc(nz, enodes))/2.-1000. + end do + + do e=1,2 + elem=eelems(e) + if (elem<=0) CYCLE + deltaX=edge_cross_dxdy(1+(e-1)*2,edge) + deltaY=edge_cross_dxdy(2+(e-1)*2,edge) + nzmax =nlevels(elem) + + do nz=nzmax-1,2,-1 + dens(nz) = (aux(nz) * helem(nz-1,elem)+& + aux(nz-1) * helem(nz, elem))/sum(helem(nz-1:nz,elem)) + end do + dens(nzmax)=dens(nzmax-1)+(dens(nzmax-1)-dens(nzmax-2))*helem(nzmax-1,elem)/helem(nzmax-2,elem) + dens(1) =dens(2) +(dens(2)-dens(3)) *helem(1, elem) /helem(2,elem) + + is=minloc(abs(std_dens-dens(1)),1) + + do nz=nzmax-1,1,-1 + div=(UV(2,nz,elem)*deltaX-UV(1,nz,elem)*deltaY)*helem(nz,elem) + if (e==2) div=-div + dmin =minval(dens(nz:nz+1)) + dmax =maxval(dens(nz:nz+1)) + ddiff=abs(dens(nz)-dens(nz+1)) + is=findloc(std_dens > dmin, value=.true., dim=1) +! is=1 +! do jj = 1, std_dens_N +! if (std_dens(jj) > dmin) then +! is = jj +! exit +! endif +! end do + + ie=findloc(std_dens < dmax, value=.true., back=.true., dim=1) +! ie=std_dens_N +! do jj = std_dens_N,1,-1 +! if (std_dens(jj) < dmin) then +! ie = jj +! exit +! endif +! end do + + if (std_dens(is)>=dmax) is=ie + if (std_dens(ie)<=dmin) ie=is + if (ie-is > 0) then + weight=(std_dens(is)-dmin)+std_dd(is)/2. + weight=max(weight, 0.)/ddiff + std_dens_DIV(is, enodes(1))=std_dens_DIV(is, enodes(1))+weight*div + std_dens_DIV(is, enodes(2))=std_dens_DIV(is, enodes(2))-weight*div + do snz=is+1, ie-1 + weight=(sum(std_dd(snz-1:snz))/2.)/ddiff + std_dens_DIV(snz, enodes(1))=std_dens_DIV(snz, enodes(1))+weight*div + std_dens_DIV(snz, enodes(2))=std_dens_DIV(snz, enodes(2))-weight*div + end do + weight=(dmax-std_dens(ie))+std_dd(ie-1)/2. + weight=max(weight, 0.)/ddiff + std_dens_DIV(ie, enodes(1))=std_dens_DIV(ie, enodes(1))+weight*div + std_dens_DIV(ie, enodes(2))=std_dens_DIV(ie, enodes(2))-weight*div + else + std_dens_DIV(is, enodes(1))=std_dens_DIV(is, enodes(1))+div + std_dens_DIV(is, enodes(2))=std_dens_DIV(is, enodes(2))-div + end if + end do + end do + end do + where (std_dens_w > 0.) std_dens_RHOZ=std_dens_RHOZ/std_dens_w + std_dens_Z =std_dens_Z /std_dens_w end where + + if (.not. firstcall_e) then + std_dens_dVdT=(std_dens_VOL2-std_dens_VOL1)/dt + end if + std_dens_VOL1=std_dens_VOL2 + firstcall_e=.false. end subroutine diag_densMOC ! ============================================================== + subroutine compute_diagnostics(mode, mesh) implicit none integer, intent(in) :: mode !constructor mode (0=only allocation; any other=do diagnostic) diff --git a/src/io_meandata.F90 b/src/io_meandata.F90 index d8999c612..77fd4b963 100644 --- a/src/io_meandata.F90 +++ b/src/io_meandata.F90 @@ -319,12 +319,18 @@ subroutine ini_mean_io(mesh) end if CASE ('dMOC ') if (ldiag_dMOC) then - call def_stream((/std_dens_N, elem2D/), (/std_dens_N, myDim_elem2D/), 'U_rho_x_DZ', 'fluxes for density MOC', 'fluxes', std_dens_UVDZ(1,:,:), 1, 'm', i_real4, mesh) - call def_stream((/std_dens_N, elem2D/), (/std_dens_N, myDim_elem2D/), 'V_rho_x_DZ', 'fluxes for density MOC', 'fluxes', std_dens_UVDZ(2,:,:), 1, 'm', i_real4, mesh) - call def_stream((/std_dens_N, elem2D/), (/std_dens_N, myDim_elem2D/), 'RHO_Z', 'drho/dz', 'kg/m4' , std_dens_RHOZ(:,:), 1, 'm', i_real4, mesh) - call def_stream((/nl-1, nod2D /), (/nl-1, myDim_nod2D /), 'density_dMOC', 'density' , 'm', density_dmoc(:,:), 1, 'm', i_real4, mesh) - end if -!___________________________________________________________________________________________________________________________________ + call def_stream((/std_dens_N, elem2D/), (/std_dens_N, myDim_elem2D/), 'U_rho_x_DZ', 'fluxes for density MOC', 'fluxes', std_dens_UVDZ(1,:,:), 1, 'y', i_real4, mesh) + call def_stream((/std_dens_N, elem2D/), (/std_dens_N, myDim_elem2D/), 'V_rho_x_DZ', 'fluxes for density MOC', 'fluxes', std_dens_UVDZ(2,:,:), 1, 'y', i_real4, mesh) + call def_stream((/std_dens_N, elem2D/), (/std_dens_N, myDim_elem2D/), 'RHO_Z', 'drho/dz', 'kg/m4' , std_dens_RHOZ(:,:), 1, 'y', i_real4, mesh) + call def_stream((/std_dens_N, elem2D/), (/std_dens_N, myDim_elem2D/), 'std_heat_flux', 'drho/dz', 'kg*m/s' ,std_dens_flux(1,:,:), 1, 'y', i_real4, mesh) + call def_stream((/std_dens_N, elem2D/), (/std_dens_N, myDim_elem2D/), 'std_rest_flux', 'drho/dz', 'kg*m/s' ,std_dens_flux(2,:,:), 1, 'y', i_real4, mesh) + call def_stream((/std_dens_N, elem2D/), (/std_dens_N, myDim_elem2D/), 'std_frwt_flux', 'drho/dz', 'kg*m/s' ,std_dens_flux(3,:,:), 1, 'y', i_real4, mesh) + call def_stream((/std_dens_N, elem2D/), (/std_dens_N, myDim_elem2D/), 'std_dens_dVdT', 'dV/dT', 'm3/s' ,std_dens_dVdT(:,:), 1, 'y', i_real4, mesh) + call def_stream((/std_dens_N, nod2D /), (/std_dens_N, myDim_nod2D/), 'std_dens_DIV', 'm3/s', 'm3/s' ,std_dens_DIV(:,:), 1, 'y', i_real4, mesh) + call def_stream((/std_dens_N, elem2D/), (/std_dens_N, myDim_elem2D/), 'std_dens_Z', 'm', 'm' ,std_dens_Z(:,:), 1, 'y', i_real4, mesh) + call def_stream((/nl-1, nod2D /), (/nl-1, myDim_nod2D /), 'density_dMOC', 'density' , 'm', density_dmoc(:,:), 1, 'y', i_real4, mesh) + call def_stream(elem2D, myDim_elem2D , 'density_flux', 'density' , 'm', dens_flux(:), 1, 'y', i_real4, mesh) + end if______________________________________________________________________________________________________________________ CASE ('pgf_x ') call def_stream((/nl-1, elem2D/), (/nl-1, myDim_elem2D/), 'pgf_x', 'zonal pressure gradient force' , 'm/s^2', pgf_x(:,:), io_list(i)%freq, io_list(i)%unit, io_list(i)%precision, mesh) CASE ('pgf_y ') From ea9432009e7dc3598bdd9b00e04510e5b808d631 Mon Sep 17 00:00:00 2001 From: Jan Streffing Date: Mon, 4 May 2020 21:00:37 +0200 Subject: [PATCH 02/26] typofix --- src/io_meandata.F90 | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/io_meandata.F90 b/src/io_meandata.F90 index 77fd4b963..1a6538dd0 100644 --- a/src/io_meandata.F90 +++ b/src/io_meandata.F90 @@ -330,7 +330,8 @@ subroutine ini_mean_io(mesh) call def_stream((/std_dens_N, elem2D/), (/std_dens_N, myDim_elem2D/), 'std_dens_Z', 'm', 'm' ,std_dens_Z(:,:), 1, 'y', i_real4, mesh) call def_stream((/nl-1, nod2D /), (/nl-1, myDim_nod2D /), 'density_dMOC', 'density' , 'm', density_dmoc(:,:), 1, 'y', i_real4, mesh) call def_stream(elem2D, myDim_elem2D , 'density_flux', 'density' , 'm', dens_flux(:), 1, 'y', i_real4, mesh) - end if______________________________________________________________________________________________________________________ + end if +!___________________________________________________________________________________________________________________________________ CASE ('pgf_x ') call def_stream((/nl-1, elem2D/), (/nl-1, myDim_elem2D/), 'pgf_x', 'zonal pressure gradient force' , 'm/s^2', pgf_x(:,:), io_list(i)%freq, io_list(i)%unit, io_list(i)%precision, mesh) CASE ('pgf_y ') From 8d0fa25e14d51b81a87368e6b0ad05f7129f87d2 Mon Sep 17 00:00:00 2001 From: JanStreffing Date: Fri, 8 May 2020 09:07:56 +0200 Subject: [PATCH 03/26] wip --- src/ice_thermo_cpl.F90 | 74 +++++++++++++++++++++++++++++++++++++++--- 1 file changed, 69 insertions(+), 5 deletions(-) diff --git a/src/ice_thermo_cpl.F90 b/src/ice_thermo_cpl.F90 index 49ab6b215..7c3082299 100644 --- a/src/ice_thermo_cpl.F90 +++ b/src/ice_thermo_cpl.F90 @@ -130,6 +130,7 @@ subroutine thermodynamics(mesh) call ice_growth #if defined (__oifs) call ice_albedo(hsn,t,alb) + call ice_surftemp(hice,alb,net_heat,t) ice_alb(inod) = alb ice_temp(inod) = t+tmelt #endif @@ -423,21 +424,84 @@ end subroutine ice_growth subroutine ice_albedo (hsn,t,alb) ! INPUT: ! hsn - snow thickness, used for albedo parameterization [m] + ! t - temperature of snow/ice surface [C] + ! + ! OUTPUT: + ! alb - snow albedo use i_therm_param implicit none real(kind=WP) hsn real(kind=WP) t - real(kind=WP) alb ! Albedo of sea ice + real(kind=WP) alb ! set albedo ! ice and snow, freezing and melting conditions are distinguished. - if (hsn.gt.0.0) then ! snow cover present - alb=albsn - else ! no snow cover - alb=albi + if (t<0.0_WP) then ! freezing condition + if (hsn.gt.0.0_WP) then ! snow cover present + alb=albsn + else ! no snow cover + alb=albi + endif + else ! melting condition + if (hsn.gt.0.0_WP) then ! snow cover present + alb=albsnm + else ! no snow cover + alb=albim + endif endif end subroutine ice_albedo + + subroutine ice_surftemp (hsn,t,alb) + ! INPUT: + ! A1 - Total atmospheric heat flux + ! t - Ice surface temperature (last step) + ! + ! OUTPUT: + ! alb - snow albedo + use i_therm_param + implicit none + + !---- atmospheric heat fluxes (provided by OpenIFS) + real(kind=WP) :: a2ohf, a2ihf + !---- evaporation and sublimation (provided by OpenIFS) + real(kind=WP) :: evap, subli + !---- precipitation and runoff (provided by OpenIFS) + real(kind=WP) :: rain, snow, runo + !---- ocean variables (provided by FESOM) + real(kind=WP) :: T_oc, S_oc, ustar + real(kind=WP) hsn + real(kind=WP) t + real(kind=WP) alb + + q1 = 11637800.0_WP + q2 = -5897.8_WP + imax = 5 + + d1=rhoair*cpair*Ch_i + d2=rhoair*Ce_i + d3=d2*clhi + + ! total incoming atmospheric heat flux + A1=a2ihf + ! NEWTON-RHAPSON TO GET TEMPERATURE AT THE TOP OF THE ICE LAYER + + do iter=1,imax + + B=q1*inv_rhoair*exp(q2/(t+tmelt)) ! (saturated) specific humidity over ice + A2=-d1*ug*t-d3*ug*B & + -emiss_ice*boltzmann*((t+tmelt)**4) ! sensible and latent heat,and outward radiation + A3=-d3*ug*B*q2/((t+tmelt)**2) ! gradient coefficient for the latent heat part + C=con/hice ! gradient coefficient for downward heat conductivity + A3=A3+C+d1*ug & ! gradient coefficient for sensible heat and radiation + +4.0_WP*emiss_ice*boltzmann*((t+tmelt)**3) + C=C*(TFrez(S_oc)-t) ! downward conductivity term + t=t+(A1+A2+C)/A3 ! NEW ICE TEMPERATURE AS THE SUM OF ALL COMPONENTS + done + + t=min(0.0_WP,t) + + end subroutine thermodynamics #endif From 9c0827b87cf28ff7ed20023493d12d667dcdb4ad Mon Sep 17 00:00:00 2001 From: JanStreffing Date: Fri, 8 May 2020 17:46:41 +0200 Subject: [PATCH 04/26] compiles, not tested yet --- src/.nfs0000000004f37b9900009651 | Bin 0 -> 24576 bytes src/gen_forcing_couple.F90 | 24 +++++-- src/ice_modules.F90 | 4 +- src/ice_setup_step.F90 | 4 +- src/ice_thermo_cpl.F90 | 110 +++++++++++++++---------------- src/io_meandata.F90 | 3 +- src/io_restart.F90 | 1 + 7 files changed, 79 insertions(+), 67 deletions(-) create mode 100644 src/.nfs0000000004f37b9900009651 diff --git a/src/.nfs0000000004f37b9900009651 b/src/.nfs0000000004f37b9900009651 new file mode 100644 index 0000000000000000000000000000000000000000..79352ea962fd791b8570bd2bd75be9ca0417c078 GIT binary patch literal 24576 zcmeHPTZ|-C8Lm}821F552+?>76lP|-t1q)NyE~gU>7DJxVdutNb|sskt*JWQU2NT` zs-D?d5=qeTq9M@-O^g^KKI#MBNCXpzBt|g!fOv^9CK?oSGT-`-|5Oed5d^?U6W%c&+Bp zquq|cbse7HEk*k5^DTin9(!5Oj1wN^%q_?AvJ*xjSL1QuWs62W7e`jkv3S)B?54|O zcXy z8mI!NfxCe>KcH!E0IvWq0p9_h0lo@60Xznnz$$PX@Y`E7?VG@7fsX>edB3JT1)K!l zxmnYm1s()$1YW;M)4m84fxCb|VZ-2YfCF~}zr0b?UIV@cECD0H3-8mkF98>UP2fJ@ zQ@}@n+kxNRfHHwsffs@2fhT}j;0E9a*xWbn`%hMBG#SY~2wxIDvb$K!tN1pess z7-Q4ZEb+o>0I!DzHh<+-`VneuA+ zJ_sb;Q0`m6crfWJYN0GT5vu2FWCJ9)a1Q92d*`^se)CI@P z^lcVIb)$;Ok`JojOS31Z80xJurSi6{1X+fvydFYNT3>I{2 zJGB}{fd)_Trh*Lpgo;-Mk=jCezA`tvvdV<7oW1Sm*jB@_8kMt48MYK-(7_s0a)K2~9@IX;iX{_gNEz5B(08Lq*9l8T&2(dKXvo#mio0t-LkA5keK)P2_IcLXC>l)`{F?q4`kqbGap^nM29VHd z{-7aWJ-5_Vs?MSy5knSZ{Sh;Ji55B3);(7X(Rz zns8`FrpMvHdo}7=+Z)!lMN>o3r_ew>yRcfhzp~szTY3ZX==F7VVWmPbFB9t2G%?>n z_ZHN*A!pZiaEaPEI>pivBhL>X>ncoV7wS3EW7qsMuH%XCWD4J9&mr=A4;KeXX!pyd z0#wKr3T-JmJ0Cwhnk|g250YV3L(E^}CWf_Av8OfZgR*jXY=<%ByX=!1lvNBY=s+j( zIC^sH=NEFQI_a&m<;{VxJf z0iOd}>;DR_PXUzv9m=B)#eiZ!F`yVw3@8Q^1BwB~fMP%~pcqgLyh|B4*v%e({Z22$1e*Y7o2HXkU0sIQRjhE)xk6CV(X8k+Q z^2-sKuQA1M8;*504^yDm^n}IV3gBerCsFDgq3b)U(jXbMPUsRc(@e5KttN!PGHG4# zVPZ~M)^}E=ktZ9r@Qood=D2Ay%MIY>)#W9?&U@>4m^UW1%lw0cSA=WEF^@}KtS&UY z4VX*8D99o!snmr$&P26R4^w(rkAkhF!Qew?#m3HB(h>?*YZ&}RWyE^}zbrQ4en4+- z!YF35lXx74NxzIt$KovfNWgCFLt5ULkkoaWU`u@Ke0fvcK^hX5SkWjm z5{HuPCLBx9xEvk}m%9_}T+Yj#>8obH(2c-NE_tT!mnMy=qm}G2u)kukldYdga-10B zZs6CqjC5Tons9GR=@*;O7gj1B#s6!Qo)37bu!q;#40#`sKLQNYuxg)S8rxlz;5|MAaE1#Z|pt%9r!))1aKI53qJi{2j+l7z>na=?*I+p!vN{}GvH15>^}+w zz`ek);J5!ea0$o&hk>`@yZ<-fN#G<<2L1#8{pW!n1CIh3;4tuyLz?z|U<)__Oar?d zR~TQ^m`JTsjfrYZ?CHI`Z`7FBiVW(XSg0{^e`DfuPbb_zTd_in%O(8B%PB)PKTLvJ z!@n3U<(Q0Na`Zx-jk7aj>w|{i&~zer@8&0Wk%=TUrRb8m`#c3c$#859@LgbSQZRBw z4DLc>=!=yACRm6ROSUttyf{C*w7aafoVJio-S;a?>P0LHy^h0cS?WcT=Y6V`ifPMe ztF>QIQY*}0N<>y6UzG0xZlXRpM{;(KW;ZSU2J)3kYKv(rwht*%Cu$plse&;fi7KNl zw26{t7pBwmy}itq$W&fiw?SeNws6Jpm>o@32V)h}80!gbkr6a~4%KMLdx`N&1EV~| z3&2b!BBRo485%`5D6g*z*puWm5G$3IMq&wJSVTi6Gezi-rtM!eg@?4FD5^q?3|kbT za^ChSIF6IcaPib(-nkyKpmy{|cS&6ZA&}53o*Gg({yjIs2u*k;Uu0Oz?VDs7C=-Tm zGZw*2A$A*@aIN2?q^?*qu|dEoA{MUJd(Gd8R3lAn8SrL{W`ww?Rv6zNwds{6gt@X! z;+zUY44lqR&#d-)VYB_K42_d&jQHUpn>YRAcDz0_@hC8B`#_|i8OYg0HsLE zK`8}U2f&n+a@9($ATTferZa4{eXVKtz|R07)Vc9a^fVjK$w;jceUu`u=`+FRP86l_ VNMr4%6rdeMCV~q@k~BzL`ybX7gVX>3 literal 0 HcmV?d00001 diff --git a/src/gen_forcing_couple.F90 b/src/gen_forcing_couple.F90 index b7d4e86d9..c8dafd78b 100755 --- a/src/gen_forcing_couple.F90 +++ b/src/gen_forcing_couple.F90 @@ -58,12 +58,25 @@ subroutine update_atm_forcing(istep, mesh) do i=1,nsend exchange =0. if (i.eq.1) then +#if defined (__oifs) + ! AWI-CM3 outgoing state vectors do n=1,myDim_nod2D+eDim_nod2D -#if defined (__oifs) exchange(n)=tr_arr(1, n, 1)+tmelt ! sea surface temperature [K] + end do + elseif (i.eq.2) then + exchange(:) = a_ice(:) ! ice concentation [%] + elseif (i.eq.3) then + exchange(:) = m_snow(:) ! snow thickness + elseif (i.eq.4) then + exchange(:) = ice_alb(:) ! ice albedo + elseif (i.eq.5) then + exchange(:) = ice_surf_temp(:) ! ice surface temperature + else + print *, 'not installed yet or error in cpl_oasis3mct_send', mype #else - exchange(n)=tr_arr(1, n, 1) ! sea surface temperature [°C] -#endif + ! AWI-CM2 outgoing state vectors + do n=1,myDim_nod2D+eDim_nod2D + exchange(n)=tr_arr(1, n, 1) ! sea surface temperature [°C] end do elseif (i.eq.2) then exchange(:) = m_ice(:) ! ice thickness [m] @@ -71,12 +84,9 @@ subroutine update_atm_forcing(istep, mesh) exchange(:) = a_ice(:) ! ice concentation [%] elseif (i.eq.4) then exchange(:) = m_snow(:) ! snow thickness -#if defined (__oifs) - elseif (i.eq.5) then - exchange(:) = ice_alb(:) ! ice albedo -#endif else print *, 'not installed yet or error in cpl_oasis3mct_send', mype +#endif endif call cpl_oasis3mct_send(i, exchange, action) enddo diff --git a/src/ice_modules.F90 b/src/ice_modules.F90 index 1e5467d18..f99802cc1 100755 --- a/src/ice_modules.F90 +++ b/src/ice_modules.F90 @@ -69,7 +69,7 @@ MODULE i_ARRAYS REAL(kind=WP), ALLOCATABLE, DIMENSION(:) :: fresh_wa_flux REAL(kind=WP), ALLOCATABLE, DIMENSION(:) :: net_heat_flux #if defined (__oasis) - real(kind=WP),target, allocatable, dimension(:) :: ice_alb, ice_temp ! new fields for OIFS coupling + real(kind=WP),target, allocatable, dimension(:) :: ice_alb, ice_surf_temp ! new fields for OIFS coupling real(kind=WP),target, allocatable, dimension(:) :: oce_heat_flux, ice_heat_flux real(kind=WP),target, allocatable, dimension(:) :: tmp_oce_heat_flux, tmp_ice_heat_flux !temporary flux fields @@ -111,6 +111,8 @@ module i_therm_param REAL(kind=WP), parameter :: inv_rhosno= 1./290. ! Snow density, AOMIP REAL(kind=WP), parameter :: cpair=1005. ! Specific heat of air [J/(kg * K)] +REAL(kind=WP), parameter :: cpice=2106. ! Specific heat of ice [J/(kg * K)] +REAL(kind=WP), parameter :: cpsno=2090. ! Specific heat of snow [J/(kg * K)] REAL(kind=WP), parameter :: cc=rhowat*4190.0 ! Volumetr. heat cap. of water [J/m**3/K](cc = rhowat*cp_water) REAL(kind=WP), parameter :: cl=rhoice*3.34e5 ! Volumetr. latent heat of ice fusion [J/m**3](cl=rhoice*Lf) REAL(kind=WP), parameter :: clhw=2.501e6 ! Specific latent heat [J/kg]: water -> water vapor diff --git a/src/ice_setup_step.F90 b/src/ice_setup_step.F90 index a21efdce9..fc82468c3 100755 --- a/src/ice_setup_step.F90 +++ b/src/ice_setup_step.F90 @@ -132,9 +132,9 @@ subroutine ice_array_setup(mesh) allocate(oce_heat_flux(n_size), ice_heat_flux(n_size)) allocate(tmp_oce_heat_flux(n_size), tmp_ice_heat_flux(n_size)) #if defined (__oifs) - allocate(ice_alb(n_size), ice_temp(n_size)) + allocate(ice_alb(n_size), ice_surf_temp(n_size)) ice_alb=0._WP - ice_temp=0._WP + ice_surf_temp=0._WP #endif /* (__oifs) */ oce_heat_flux=0._WP ice_heat_flux=0._WP diff --git a/src/ice_thermo_cpl.F90 b/src/ice_thermo_cpl.F90 index 7c3082299..ebb57d8c6 100644 --- a/src/ice_thermo_cpl.F90 +++ b/src/ice_thermo_cpl.F90 @@ -4,7 +4,7 @@ subroutine thermodynamics(mesh) !=================================================================== ! ! This subroutine computes thermodynamic changes of ice and snow - ! for coupled simulations of FESOM and ECHAM. + ! for coupled simulations of FESOM2. ! It replaces the original FESOM scheme for uncoupled simulations. ! Note that atmospheric fluxes already need to be available. ! @@ -129,10 +129,14 @@ subroutine thermodynamics(mesh) call ice_growth #if defined (__oifs) + !---- For AWI-CM3 we calculate ice surface temp and albedo in fesom, + ! then send those to OpenIFS where they are used to calucate the + ! energy fluxes ---! + t = ice_surf_temp(inod) + call ice_surftemp(h,hsn,a2ihf,t) call ice_albedo(hsn,t,alb) - call ice_surftemp(hice,alb,net_heat,t) - ice_alb(inod) = alb - ice_temp(inod) = t+tmelt + ice_alb(inod) = alb + ice_surf_temp(inod) = t #endif @@ -421,7 +425,52 @@ subroutine ice_growth end subroutine ice_growth - subroutine ice_albedo (hsn,t,alb) + + + subroutine ice_surftemp(h,hsn,a2ihf,t) + ! INPUT: + ! a2ihf - Total atmo heat flux to ice + ! h - Ice thickness + ! hsn - Snow thickness + ! + ! INPUT/OUTPUT: + ! t - Ice surface temperature + + use i_therm_param + implicit none + + !---- atmospheric heat net flux into to ice (provided by OpenIFS) + real(kind=WP) a2ihf + !---- ocean variables (provided by FESOM) + real(kind=WP) h + real(kind=WP) hsn + real(kind=WP) t + !---- local variables + real(kind=WP) snicecond + real(kind=WP) zsniced + real(kind=WP) zicefl + real(kind=WP) hcapice + real(kind=WP) zcpdt + real(kind=WP) zcpdte + real(kind=WP) zcprosn + !---- local parameters + real(kind=WP), parameter :: dice = 0.05_WP ! ECHAM6's thickness for top ice "layer" + real(kind=WP), parameter :: ctfreez = 271.38_WP ! ECHAM6's temperature at which sea starts freezing/melting + + + snicecond = con/consn*rhowat/rhosno ! equivalence fraction thickness of ice/snow + zsniced=MAX(h,hmin)+snicecond*hsn ! Ice + Snow-Ice-equivalent thickness [m] + zicefl=con*ctfreez/zsniced ! Conductive heat flux through sea ice [W/m²] + hcapice=rhoice*cpice*dice ! heat capacity of upper 0.05 cm sea ice layer [J/(m²K)] + zcpdt=hcapice/dt ! Energy required to change temperature of top ice "layer" [J/(sm²K)] + zcprosn=rhowat*cpsno/dt ! Specific Energy required to change temperature of 1m snow on ice [J/(sm³K)] + zcpdte=zcpdt+zcprosn*hsn ! Combined Energy required to change temperature of snow + 0.05m of upper ice + t=(zcpdte*t+a2ihf+zicefl)/(zcpdte+con/zsniced) ! New sea ice surf temp [K] + t=min(ctfreez,t) ! Not warmer than freezing please! + end subroutine ice_surftemp + + + subroutine ice_albedo(hsn,t,alb) ! INPUT: ! hsn - snow thickness, used for albedo parameterization [m] ! t - temperature of snow/ice surface [C] @@ -452,56 +501,5 @@ subroutine ice_albedo (hsn,t,alb) endif end subroutine ice_albedo - - subroutine ice_surftemp (hsn,t,alb) - ! INPUT: - ! A1 - Total atmospheric heat flux - ! t - Ice surface temperature (last step) - ! - ! OUTPUT: - ! alb - snow albedo - use i_therm_param - implicit none - - !---- atmospheric heat fluxes (provided by OpenIFS) - real(kind=WP) :: a2ohf, a2ihf - !---- evaporation and sublimation (provided by OpenIFS) - real(kind=WP) :: evap, subli - !---- precipitation and runoff (provided by OpenIFS) - real(kind=WP) :: rain, snow, runo - !---- ocean variables (provided by FESOM) - real(kind=WP) :: T_oc, S_oc, ustar - real(kind=WP) hsn - real(kind=WP) t - real(kind=WP) alb - - q1 = 11637800.0_WP - q2 = -5897.8_WP - imax = 5 - - d1=rhoair*cpair*Ch_i - d2=rhoair*Ce_i - d3=d2*clhi - - ! total incoming atmospheric heat flux - A1=a2ihf - ! NEWTON-RHAPSON TO GET TEMPERATURE AT THE TOP OF THE ICE LAYER - - do iter=1,imax - - B=q1*inv_rhoair*exp(q2/(t+tmelt)) ! (saturated) specific humidity over ice - A2=-d1*ug*t-d3*ug*B & - -emiss_ice*boltzmann*((t+tmelt)**4) ! sensible and latent heat,and outward radiation - A3=-d3*ug*B*q2/((t+tmelt)**2) ! gradient coefficient for the latent heat part - C=con/hice ! gradient coefficient for downward heat conductivity - A3=A3+C+d1*ug & ! gradient coefficient for sensible heat and radiation - +4.0_WP*emiss_ice*boltzmann*((t+tmelt)**3) - C=C*(TFrez(S_oc)-t) ! downward conductivity term - t=t+(A1+A2+C)/A3 ! NEW ICE TEMPERATURE AS THE SUM OF ALL COMPONENTS - done - - t=min(0.0_WP,t) - - end subroutine thermodynamics #endif diff --git a/src/io_meandata.F90 b/src/io_meandata.F90 index 1a6538dd0..a462a6b5f 100644 --- a/src/io_meandata.F90 +++ b/src/io_meandata.F90 @@ -382,7 +382,8 @@ subroutine ini_mean_io(mesh) end if #if defined (__oifs) - call def_stream(nod2D, myDim_nod2D, 'alb', 'ice albedo', 'none', ice_alb(:), 1, 'm', i_real4, mesh) + call def_stream(nod2D, myDim_nod2D, 'alb', 'ice albedo', 'none', ice_alb(:), 1, 'm', i_real4, mesh) + call def_stream(nod2D, myDim_nod2D, 'alb', 'ice surface temperature', 'none', ice_surf_temp(:), 1, 'm', i_real4, mesh) #endif if (mix_scheme_nmb==5 .or. mix_scheme_nmb==56) then diff --git a/src/io_restart.F90 b/src/io_restart.F90 index dd4e0554b..1ac4d9d5b 100644 --- a/src/io_restart.F90 +++ b/src/io_restart.F90 @@ -181,6 +181,7 @@ subroutine ini_ice_io(year, mesh) call def_variable(iid, 'hsnow', (/nod2D/), 'effective snow thickness', 'm', m_snow); call def_variable(iid, 'uice', (/nod2D/), 'zonal velocity', 'm/s', u_ice); call def_variable(iid, 'vice', (/nod2D/), 'meridional velocity', 'm', v_ice); + call def_variable(iid, 'ice_surf_temp',(/nod2D/), 'ice surface temperature', 'K', ice_surf_temp); end subroutine ini_ice_io ! From fe60931ecfe5faf79d3f5693f4b94e4e61550a7d Mon Sep 17 00:00:00 2001 From: JanStreffing Date: Sat, 9 May 2020 03:13:15 +0200 Subject: [PATCH 05/26] runs for a while but looks strange and crashes --- src/.nfs0000000004f37b9900009651 | Bin 24576 -> 0 bytes src/cpl_driver.F90 | 9 ++++++-- src/gen_forcing_couple.F90 | 4 ++-- src/ice_setup_step.F90 | 4 ++-- src/ice_thermo_cpl.F90 | 34 +++++++++++++++++-------------- src/io_meandata.F90 | 2 +- src/io_restart.F90 | 1 + 7 files changed, 32 insertions(+), 22 deletions(-) delete mode 100644 src/.nfs0000000004f37b9900009651 diff --git a/src/.nfs0000000004f37b9900009651 b/src/.nfs0000000004f37b9900009651 deleted file mode 100644 index 79352ea962fd791b8570bd2bd75be9ca0417c078..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 24576 zcmeHPTZ|-C8Lm}821F552+?>76lP|-t1q)NyE~gU>7DJxVdutNb|sskt*JWQU2NT` zs-D?d5=qeTq9M@-O^g^KKI#MBNCXpzBt|g!fOv^9CK?oSGT-`-|5Oed5d^?U6W%c&+Bp zquq|cbse7HEk*k5^DTin9(!5Oj1wN^%q_?AvJ*xjSL1QuWs62W7e`jkv3S)B?54|O zcXy z8mI!NfxCe>KcH!E0IvWq0p9_h0lo@60Xznnz$$PX@Y`E7?VG@7fsX>edB3JT1)K!l zxmnYm1s()$1YW;M)4m84fxCb|VZ-2YfCF~}zr0b?UIV@cECD0H3-8mkF98>UP2fJ@ zQ@}@n+kxNRfHHwsffs@2fhT}j;0E9a*xWbn`%hMBG#SY~2wxIDvb$K!tN1pess z7-Q4ZEb+o>0I!DzHh<+-`VneuA+ zJ_sb;Q0`m6crfWJYN0GT5vu2FWCJ9)a1Q92d*`^se)CI@P z^lcVIb)$;Ok`JojOS31Z80xJurSi6{1X+fvydFYNT3>I{2 zJGB}{fd)_Trh*Lpgo;-Mk=jCezA`tvvdV<7oW1Sm*jB@_8kMt48MYK-(7_s0a)K2~9@IX;iX{_gNEz5B(08Lq*9l8T&2(dKXvo#mio0t-LkA5keK)P2_IcLXC>l)`{F?q4`kqbGap^nM29VHd z{-7aWJ-5_Vs?MSy5knSZ{Sh;Ji55B3);(7X(Rz zns8`FrpMvHdo}7=+Z)!lMN>o3r_ew>yRcfhzp~szTY3ZX==F7VVWmPbFB9t2G%?>n z_ZHN*A!pZiaEaPEI>pivBhL>X>ncoV7wS3EW7qsMuH%XCWD4J9&mr=A4;KeXX!pyd z0#wKr3T-JmJ0Cwhnk|g250YV3L(E^}CWf_Av8OfZgR*jXY=<%ByX=!1lvNBY=s+j( zIC^sH=NEFQI_a&m<;{VxJf z0iOd}>;DR_PXUzv9m=B)#eiZ!F`yVw3@8Q^1BwB~fMP%~pcqgLyh|B4*v%e({Z22$1e*Y7o2HXkU0sIQRjhE)xk6CV(X8k+Q z^2-sKuQA1M8;*504^yDm^n}IV3gBerCsFDgq3b)U(jXbMPUsRc(@e5KttN!PGHG4# zVPZ~M)^}E=ktZ9r@Qood=D2Ay%MIY>)#W9?&U@>4m^UW1%lw0cSA=WEF^@}KtS&UY z4VX*8D99o!snmr$&P26R4^w(rkAkhF!Qew?#m3HB(h>?*YZ&}RWyE^}zbrQ4en4+- z!YF35lXx74NxzIt$KovfNWgCFLt5ULkkoaWU`u@Ke0fvcK^hX5SkWjm z5{HuPCLBx9xEvk}m%9_}T+Yj#>8obH(2c-NE_tT!mnMy=qm}G2u)kukldYdga-10B zZs6CqjC5Tons9GR=@*;O7gj1B#s6!Qo)37bu!q;#40#`sKLQNYuxg)S8rxlz;5|MAaE1#Z|pt%9r!))1aKI53qJi{2j+l7z>na=?*I+p!vN{}GvH15>^}+w zz`ek);J5!ea0$o&hk>`@yZ<-fN#G<<2L1#8{pW!n1CIh3;4tuyLz?z|U<)__Oar?d zR~TQ^m`JTsjfrYZ?CHI`Z`7FBiVW(XSg0{^e`DfuPbb_zTd_in%O(8B%PB)PKTLvJ z!@n3U<(Q0Na`Zx-jk7aj>w|{i&~zer@8&0Wk%=TUrRb8m`#c3c$#859@LgbSQZRBw z4DLc>=!=yACRm6ROSUttyf{C*w7aafoVJio-S;a?>P0LHy^h0cS?WcT=Y6V`ifPMe ztF>QIQY*}0N<>y6UzG0xZlXRpM{;(KW;ZSU2J)3kYKv(rwht*%Cu$plse&;fi7KNl zw26{t7pBwmy}itq$W&fiw?SeNws6Jpm>o@32V)h}80!gbkr6a~4%KMLdx`N&1EV~| z3&2b!BBRo485%`5D6g*z*puWm5G$3IMq&wJSVTi6Gezi-rtM!eg@?4FD5^q?3|kbT za^ChSIF6IcaPib(-nkyKpmy{|cS&6ZA&}53o*Gg({yjIs2u*k;Uu0Oz?VDs7C=-Tm zGZw*2A$A*@aIN2?q^?*qu|dEoA{MUJd(Gd8R3lAn8SrL{W`ww?Rv6zNwds{6gt@X! z;+zUY44lqR&#d-)VYB_K42_d&jQHUpn>YRAcDz0_@hC8B`#_|i8OYg0HsLE zK`8}U2f&n+a@9($ATTferZa4{eXVKtz|R07)Vc9a^fVjK$w;jceUu`u=`+FRP86l_ VNMr4%6rdeMCV~q@k~BzL`ybX7gVX>3 diff --git a/src/cpl_driver.F90 b/src/cpl_driver.F90 index 3db2769d7..a2936341e 100755 --- a/src/cpl_driver.F90 +++ b/src/cpl_driver.F90 @@ -377,12 +377,17 @@ subroutine cpl_oasis3mct_define_unstr(mesh) ! ... Define symbolic names for the transient fields send by the ocean ! These must be identical to the names specified in the SMIOC file. ! +#if defined (__oifs) cpl_send( 1)='sst_feom' ! 1. sea surface temperature [K] -> + cpl_send( 2)='sie_feom' ! 2. sea ice extent [%-100] -> + cpl_send( 3)='snt_feom' ! 3. snow thickness [m] -> + cpl_send( 4)='ist_feom' ! 4. sea ice surface temperature [K] -> + cpl_send( 5)='sia_feom' ! 5. sea ice albedo [%-100] -> +#else + cpl_send( 1)='sst_feom' ! 1. sea surface temperature [°C] -> cpl_send( 2)='sit_feom' ! 2. sea ice thickness [m] -> cpl_send( 3)='sie_feom' ! 3. sea ice extent [%-100] -> cpl_send( 4)='snt_feom' ! 4. snow thickness [m] -> -#if defined (__oifs) - cpl_send( 5)='sia_feom' ! 5. sea ice albedo [%-100] -> #endif diff --git a/src/gen_forcing_couple.F90 b/src/gen_forcing_couple.F90 index c8dafd78b..5a19d444d 100755 --- a/src/gen_forcing_couple.F90 +++ b/src/gen_forcing_couple.F90 @@ -68,9 +68,9 @@ subroutine update_atm_forcing(istep, mesh) elseif (i.eq.3) then exchange(:) = m_snow(:) ! snow thickness elseif (i.eq.4) then - exchange(:) = ice_alb(:) ! ice albedo - elseif (i.eq.5) then exchange(:) = ice_surf_temp(:) ! ice surface temperature + elseif (i.eq.5) then + exchange(:) = ice_alb(:) ! ice albedo else print *, 'not installed yet or error in cpl_oasis3mct_send', mype #else diff --git a/src/ice_setup_step.F90 b/src/ice_setup_step.F90 index fc82468c3..d7afac3e8 100755 --- a/src/ice_setup_step.F90 +++ b/src/ice_setup_step.F90 @@ -133,8 +133,8 @@ subroutine ice_array_setup(mesh) allocate(tmp_oce_heat_flux(n_size), tmp_ice_heat_flux(n_size)) #if defined (__oifs) allocate(ice_alb(n_size), ice_surf_temp(n_size)) - ice_alb=0._WP - ice_surf_temp=0._WP + ice_alb=0.6_WP + ice_surf_temp=275.15_WP #endif /* (__oifs) */ oce_heat_flux=0._WP ice_heat_flux=0._WP diff --git a/src/ice_thermo_cpl.F90 b/src/ice_thermo_cpl.F90 index ebb57d8c6..4fc6d231d 100644 --- a/src/ice_thermo_cpl.F90 +++ b/src/ice_thermo_cpl.F90 @@ -134,7 +134,7 @@ subroutine thermodynamics(mesh) ! energy fluxes ---! t = ice_surf_temp(inod) call ice_surftemp(h,hsn,a2ihf,t) - call ice_albedo(hsn,t,alb) + call ice_albedo(h,hsn,t,alb) ice_alb(inod) = alb ice_surf_temp(inod) = t #endif @@ -470,7 +470,7 @@ subroutine ice_surftemp(h,hsn,a2ihf,t) end subroutine ice_surftemp - subroutine ice_albedo(hsn,t,alb) + subroutine ice_albedo(h,hsn,t,alb) ! INPUT: ! hsn - snow thickness, used for albedo parameterization [m] ! t - temperature of snow/ice surface [C] @@ -485,20 +485,24 @@ subroutine ice_albedo(hsn,t,alb) real(kind=WP) alb ! set albedo - ! ice and snow, freezing and melting conditions are distinguished. - if (t<0.0_WP) then ! freezing condition - if (hsn.gt.0.0_WP) then ! snow cover present - alb=albsn - else ! no snow cover - alb=albi + ! ice and snow, freezing and melting conditions are distinguished + if (h>0.0_WP) + if (t<0.0_WP) then ! freezing condition + if (hsn.gt.0.0_WP) then ! snow cover present + alb=albsn + else ! no snow cover + alb=albi + endif + else ! melting condition + if (hsn.gt.0.0_WP) then ! snow cover present + alb=albsnm + else ! no snow cover + alb=albim + endif endif - else ! melting condition - if (hsn.gt.0.0_WP) then ! snow cover present - alb=albsnm - else ! no snow cover - alb=albim - endif - endif + else + alb = 0.066_WP ! ocean albedo + endif end subroutine ice_albedo end subroutine thermodynamics diff --git a/src/io_meandata.F90 b/src/io_meandata.F90 index a462a6b5f..c9ff08ec3 100644 --- a/src/io_meandata.F90 +++ b/src/io_meandata.F90 @@ -383,7 +383,7 @@ subroutine ini_mean_io(mesh) #if defined (__oifs) call def_stream(nod2D, myDim_nod2D, 'alb', 'ice albedo', 'none', ice_alb(:), 1, 'm', i_real4, mesh) - call def_stream(nod2D, myDim_nod2D, 'alb', 'ice surface temperature', 'none', ice_surf_temp(:), 1, 'm', i_real4, mesh) + call def_stream(nod2D, myDim_nod2D, 'ist', 'ice surface temperature', 'K', ice_surf_temp(:), 1, 'm', i_real4, mesh) #endif if (mix_scheme_nmb==5 .or. mix_scheme_nmb==56) then diff --git a/src/io_restart.F90 b/src/io_restart.F90 index 1ac4d9d5b..4114e842e 100644 --- a/src/io_restart.F90 +++ b/src/io_restart.F90 @@ -181,6 +181,7 @@ subroutine ini_ice_io(year, mesh) call def_variable(iid, 'hsnow', (/nod2D/), 'effective snow thickness', 'm', m_snow); call def_variable(iid, 'uice', (/nod2D/), 'zonal velocity', 'm/s', u_ice); call def_variable(iid, 'vice', (/nod2D/), 'meridional velocity', 'm', v_ice); + call def_variable(iid, 'ice_albedo', (/nod2D/), 'ice albedo', '-', ice_alb); call def_variable(iid, 'ice_surf_temp',(/nod2D/), 'ice surface temperature', 'K', ice_surf_temp); end subroutine ini_ice_io From 8b71ad23e0232af4252bfc36a49a820dc3335a68 Mon Sep 17 00:00:00 2001 From: JanStreffing Date: Sat, 9 May 2020 16:43:41 +0200 Subject: [PATCH 06/26] moved ice heat flux output into section that is only used by oifs --- src/ice_thermo_cpl.F90 | 5 +++-- src/io_meandata.F90 | 3 ++- 2 files changed, 5 insertions(+), 3 deletions(-) diff --git a/src/ice_thermo_cpl.F90 b/src/ice_thermo_cpl.F90 index 4fc6d231d..431f49661 100644 --- a/src/ice_thermo_cpl.F90 +++ b/src/ice_thermo_cpl.F90 @@ -480,13 +480,14 @@ subroutine ice_albedo(h,hsn,t,alb) use i_therm_param implicit none + real(kind=WP) h real(kind=WP) hsn real(kind=WP) t real(kind=WP) alb ! set albedo ! ice and snow, freezing and melting conditions are distinguished - if (h>0.0_WP) + if (h>0.0_WP) then if (t<0.0_WP) then ! freezing condition if (hsn.gt.0.0_WP) then ! snow cover present alb=albsn @@ -501,7 +502,7 @@ subroutine ice_albedo(h,hsn,t,alb) endif endif else - alb = 0.066_WP ! ocean albedo + alb=0.066_WP ! ocean albedo endif end subroutine ice_albedo diff --git a/src/io_meandata.F90 b/src/io_meandata.F90 index c9ff08ec3..dc3e65649 100644 --- a/src/io_meandata.F90 +++ b/src/io_meandata.F90 @@ -383,7 +383,8 @@ subroutine ini_mean_io(mesh) #if defined (__oifs) call def_stream(nod2D, myDim_nod2D, 'alb', 'ice albedo', 'none', ice_alb(:), 1, 'm', i_real4, mesh) - call def_stream(nod2D, myDim_nod2D, 'ist', 'ice surface temperature', 'K', ice_surf_temp(:), 1, 'm', i_real4, mesh) + call def_stream(nod2D, myDim_nod2D, 'ist', 'ice surface temperature', 'K', ice_surf_temp(:), 1, 'd', i_real4, mesh) + call def_stream(nod2D, myDim_nod2D , 'qsi' , 'ice heat flux', 'W/m^2', ice_heat_flux(:), 1, 'm', i_real4, mesh) #endif if (mix_scheme_nmb==5 .or. mix_scheme_nmb==56) then From 8cd095538fc7445aca39e3e0730a80ec2b002d03 Mon Sep 17 00:00:00 2001 From: JanStreffing Date: Sat, 9 May 2020 18:55:40 +0200 Subject: [PATCH 07/26] added div by A and heff --- src/.nfs0000000004f806ef0000967e | Bin 0 -> 16384 bytes src/ice_thermo_cpl.F90 | 16 ++++++++++------ src/io_meandata.F90 | 2 +- 3 files changed, 11 insertions(+), 7 deletions(-) create mode 100644 src/.nfs0000000004f806ef0000967e diff --git a/src/.nfs0000000004f806ef0000967e b/src/.nfs0000000004f806ef0000967e new file mode 100644 index 0000000000000000000000000000000000000000..8f3d236c5e9f4beedc223930abf6c858939b02ae GIT binary patch literal 16384 zcmeHOOKc=Z8Lp6JAuI_Yt0WQ!66b?q($P900{}n_g7cj<0jbaF|0sa zrTskJ^;iA%|6jeUdes!U>UvEwwTKYr84uWmYw_l5z(fMLKeU>GnA7zPXjhJmY$foyOS z`#gMpllFUA|6Wu1eNBf_es4woPc3ij_-|F@uWETq$N#b-e?`m3b^H-+KxxN6wEPY& zZ&lZrOIm(X%eN}>ziaskEq|yY|68T}rz-NlYWatB{@+*RFKYR$mXB$| zbu#vXmVZ>sgNpnwT3+)1L`D8*Ef+ff(QAu`QvaW{d{N8$75N{vJl68>Rpk1Wy;IBU zM|SJi(dNf6U>GnA7zPXjh5^HXVZbn87%&VN2L4AGa6HC72HU<(8^s$LEBBXwcn4#b zfaihdfTw{+fo}j`0k(kqfK{LaG=O&jM}Ze^VC>hxuYgB^M}ThvUjn`W*uZCi*RE&m zRp7_KkANQnj{^?_J>YiWcW-Cx>%b~-9q{Ba#x{U4;K8>sb_4LW>lix@{QRwqeHFL| zI0iiP7RDX`I>0MO8G8)40E_|8y_vBu0-Hb%+zQ+b{1*G#hkyxS9QYYFzCQ&X1ilO` z0k;8{u>bx6@Cfi7AO$V}#{uH<;d0~sAFB5n?}|ty(&aM3XA{2T3cl&2JP{XinFx<( zF?Z55cBPXEzTqU&Sr0_YW1r_yKZ-Xao;XpTZ_2F4v!39CSSeB|=sk5pp^~uL_hQg_J(M1CkSE}pL@XW2G5+ND>Bac%nV1$K&DxJW{e@Th~O*~ z;IRnvi)(WqnOiAecYcv|kO)_%GLE=4$x6xW-JO-r?ApI(L54h%R4RgX>J(EC__J3L zKB}Nd_eMq30})xRJI6{*JGV-iXEN^?YW7E(R%yTn|R#kKJB!N#{y~&8D!i^&@<$e-}sxhvMu8blX zb*a5lYupk}!uBO#cVQrIadmBFuCuVJYC^TL1^o@to~^7PlD>9wLTY7Kg2!FWIEcH_ zb%G&(_I)IqQ(H!Mo(ik!J3%U%jB@Elv#)}o1kfdtBu?x!bFw_GwFC~1~sq2F?gz_lN*PF91Evm>aZf& zu20TLjk2bk`x@mEiX{IwpoWhx>9b_2KGVPiL&-<$?l4J>86w#|YSf{wbr%ZZ zgM2-ZZe=-zZZ0CKa#cQcd|GP-1xuJl+RCJ1wnx;IOz+lo*Yv11C1;eT5ovO<=5XAt zL#XHmw5lox#?j0O+{5dLPxW}_vCzdN4Y%sPP6WSPN=uoIB8kH<ArbQbNo#q&B5{9MLDe*tAh!PX4l z>`AvbcYe80ftkW#ih5^HXVZbn87%&VN1`Gp+0mFb{ zz%cOt!2n%Tyur4SQGA-OC)mczNL@wR>b3xTP7-;vQAMC<=@-v Date: Sat, 9 May 2020 19:20:10 +0200 Subject: [PATCH 08/26] wip --- src/ice_thermo_cpl.F90 | 18 ++++++++---------- 1 file changed, 8 insertions(+), 10 deletions(-) diff --git a/src/ice_thermo_cpl.F90 b/src/ice_thermo_cpl.F90 index 412efbd84..cada04019 100644 --- a/src/ice_thermo_cpl.F90 +++ b/src/ice_thermo_cpl.F90 @@ -79,6 +79,7 @@ subroutine thermodynamics(mesh) real(kind=WP) :: h0min = 0.50, h0max = 1.5 type(t_mesh), intent(in) , target :: mesh + real(kind=WP), parameter :: Aimin = 0.001, himin = 0.005 #include "associate_mesh.h" @@ -133,7 +134,7 @@ subroutine thermodynamics(mesh) ! then send those to OpenIFS where they are used to calucate the ! energy fluxes ---! t = ice_surf_temp(inod) - call ice_surftemp(A,h,hsn,a2ihf,t) + call ice_surftemp(h/(min(A,Aimin)),hsn/(min(A,Aimin)),a2ihf,t) call ice_albedo(h,hsn,t,alb) ice_alb(inod) = alb ice_surf_temp(inod) = t @@ -427,7 +428,7 @@ end subroutine ice_growth - subroutine ice_surftemp(A,h,hsn,a2ihf,t) + subroutine ice_surftemp(h,hsn,a2ihf,t) ! INPUT: ! a2ihf - Total atmo heat flux to ice ! A - Ice fraction @@ -443,10 +444,8 @@ subroutine ice_surftemp(A,h,hsn,a2ihf,t) !---- atmospheric heat net flux into to ice (provided by OpenIFS) real(kind=WP) a2ihf !---- ocean variables (provided by FESOM) - real(kind=WP) A real(kind=WP) h real(kind=WP) hsn - real(kind=WP) heff real(kind=WP) t !---- local variables real(kind=WP) snicecond @@ -459,20 +458,19 @@ subroutine ice_surftemp(A,h,hsn,a2ihf,t) !---- local parameters real(kind=WP), parameter :: dice = 0.05_WP ! ECHAM6's thickness for top ice "layer" real(kind=WP), parameter :: ctfreez = 271.38_WP ! ECHAM6's temperature at which sea starts freezing/melting - real(kind=WP), parameter :: Aimin = 0.001_WP ! ECHAM6's temperature at which sea starts freezing/melting - heff = h + hsn*con/consn - heff = heff/max(A,Aimin) - zicefl=con*ctfreez/heff ! Conductive heat flux through sea ice [W/m²] + snicecond = con/consn*rhowat/rhosno ! equivalence fraction thickness of ice/snow + zsniced=MAX(h,hmin)+snicecond*hsn ! Ice + Snow-Ice-equivalent thickness [m] + zicefl=con*ctfreez/zsniced ! Conductive heat flux through sea ice [W/m²] hcapice=rhoice*cpice*dice ! heat capacity of upper 0.05 cm sea ice layer [J/(m²K)] zcpdt=hcapice/dt ! Energy required to change temperature of top ice "layer" [J/(sm²K)] zcprosn=rhowat*cpsno/dt ! Specific Energy required to change temperature of 1m snow on ice [J/(sm³K)] - zcpdte=zcpdt+zcprosn*hsn/max(A,Aimin) ! Combined Energy required to change temperature of snow + 0.05m of upper ice + zcpdte=zcpdt+zcprosn*hsn ! Combined Energy required to change temperature of snow + 0.05m of upper ice t=(zcpdte*t+a2ihf+zicefl)/(zcpdte+con/zsniced) ! New sea ice surf temp [K] t=min(ctfreez,t) ! Not warmer than freezing please! - end subroutine ice_surftemp + end subroutine ice_surftemp subroutine ice_albedo(h,hsn,t,alb) ! INPUT: From b9678469f6db9898452926bed881b5ede3755232 Mon Sep 17 00:00:00 2001 From: JanStreffing Date: Sat, 9 May 2020 20:10:43 +0200 Subject: [PATCH 09/26] wip --- src/.nfs0000000004f806ef0000967e | Bin 16384 -> 0 bytes src/ice_thermo_cpl.F90 | 9 ++++++--- src/io_meandata.F90 | 4 ++-- 3 files changed, 8 insertions(+), 5 deletions(-) delete mode 100644 src/.nfs0000000004f806ef0000967e diff --git a/src/.nfs0000000004f806ef0000967e b/src/.nfs0000000004f806ef0000967e deleted file mode 100644 index 8f3d236c5e9f4beedc223930abf6c858939b02ae..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 16384 zcmeHOOKc=Z8Lp6JAuI_Yt0WQ!66b?q($P900{}n_g7cj<0jbaF|0sa zrTskJ^;iA%|6jeUdes!U>UvEwwTKYr84uWmYw_l5z(fMLKeU>GnA7zPXjhJmY$foyOS z`#gMpllFUA|6Wu1eNBf_es4woPc3ij_-|F@uWETq$N#b-e?`m3b^H-+KxxN6wEPY& zZ&lZrOIm(X%eN}>ziaskEq|yY|68T}rz-NlYWatB{@+*RFKYR$mXB$| zbu#vXmVZ>sgNpnwT3+)1L`D8*Ef+ff(QAu`QvaW{d{N8$75N{vJl68>Rpk1Wy;IBU zM|SJi(dNf6U>GnA7zPXjh5^HXVZbn87%&VN2L4AGa6HC72HU<(8^s$LEBBXwcn4#b zfaihdfTw{+fo}j`0k(kqfK{LaG=O&jM}Ze^VC>hxuYgB^M}ThvUjn`W*uZCi*RE&m zRp7_KkANQnj{^?_J>YiWcW-Cx>%b~-9q{Ba#x{U4;K8>sb_4LW>lix@{QRwqeHFL| zI0iiP7RDX`I>0MO8G8)40E_|8y_vBu0-Hb%+zQ+b{1*G#hkyxS9QYYFzCQ&X1ilO` z0k;8{u>bx6@Cfi7AO$V}#{uH<;d0~sAFB5n?}|ty(&aM3XA{2T3cl&2JP{XinFx<( zF?Z55cBPXEzTqU&Sr0_YW1r_yKZ-Xao;XpTZ_2F4v!39CSSeB|=sk5pp^~uL_hQg_J(M1CkSE}pL@XW2G5+ND>Bac%nV1$K&DxJW{e@Th~O*~ z;IRnvi)(WqnOiAecYcv|kO)_%GLE=4$x6xW-JO-r?ApI(L54h%R4RgX>J(EC__J3L zKB}Nd_eMq30})xRJI6{*JGV-iXEN^?YW7E(R%yTn|R#kKJB!N#{y~&8D!i^&@<$e-}sxhvMu8blX zb*a5lYupk}!uBO#cVQrIadmBFuCuVJYC^TL1^o@to~^7PlD>9wLTY7Kg2!FWIEcH_ zb%G&(_I)IqQ(H!Mo(ik!J3%U%jB@Elv#)}o1kfdtBu?x!bFw_GwFC~1~sq2F?gz_lN*PF91Evm>aZf& zu20TLjk2bk`x@mEiX{IwpoWhx>9b_2KGVPiL&-<$?l4J>86w#|YSf{wbr%ZZ zgM2-ZZe=-zZZ0CKa#cQcd|GP-1xuJl+RCJ1wnx;IOz+lo*Yv11C1;eT5ovO<=5XAt zL#XHmw5lox#?j0O+{5dLPxW}_vCzdN4Y%sPP6WSPN=uoIB8kH<ArbQbNo#q&B5{9MLDe*tAh!PX4l z>`AvbcYe80ftkW#ih5^HXVZbn87%&VN1`Gp+0mFb{ zz%cOt!2n%Tyur4SQGA-OC)mczNL@wR>b3xTP7-;vQAMC<=@-vAimin) then + call ice_surftemp(max(h,0.05)/(max(A,Aimin)),hsn/(max(A,Aimin)),a2ihf,t) + ice_surf_temp(inod) = t + else + ice_surf_temp(inod) = 275.15_WP + endif call ice_albedo(h,hsn,t,alb) ice_alb(inod) = alb - ice_surf_temp(inod) = t #endif @@ -469,7 +473,6 @@ subroutine ice_surftemp(h,hsn,a2ihf,t) zcpdte=zcpdt+zcprosn*hsn ! Combined Energy required to change temperature of snow + 0.05m of upper ice t=(zcpdte*t+a2ihf+zicefl)/(zcpdte+con/zsniced) ! New sea ice surf temp [K] t=min(ctfreez,t) ! Not warmer than freezing please! - end subroutine ice_surftemp subroutine ice_albedo(h,hsn,t,alb) diff --git a/src/io_meandata.F90 b/src/io_meandata.F90 index ff317c797..d9910238a 100644 --- a/src/io_meandata.F90 +++ b/src/io_meandata.F90 @@ -383,8 +383,8 @@ subroutine ini_mean_io(mesh) #if defined (__oifs) call def_stream(nod2D, myDim_nod2D, 'alb', 'ice albedo', 'none', ice_alb(:), 1, 'm', i_real4, mesh) - call def_stream(nod2D, myDim_nod2D, 'ist', 'ice surface temperature', 'K', ice_surf_temp(:), 1, 'd', i_real4, mesh) - call def_stream(nod2D, myDim_nod2D, 'qsi', 'ice heat flux', 'W/m^2', ice_heat_flux(:), 1, 'd', i_real4, mesh) + call def_stream(nod2D, myDim_nod2D, 'ist', 'ice surface temperature', 'K', ice_surf_temp(:), 1, 'm', i_real4, mesh) + call def_stream(nod2D, myDim_nod2D, 'qsi', 'ice heat flux', 'W/m^2', ice_heat_flux(:), 1, 'm', i_real4, mesh) #endif if (mix_scheme_nmb==5 .or. mix_scheme_nmb==56) then From e8174848380cf960558226061447893664c2036d Mon Sep 17 00:00:00 2001 From: Lorenzo Zampieri Date: Sun, 10 May 2020 17:32:34 +0200 Subject: [PATCH 10/26] Added surface temperature advection for OpenIFS coupling --- src/ice_fct.F90 | 129 +++++++++++++++++++++++++++++++++++++++-- src/ice_modules.F90 | 2 + src/ice_setup_step.F90 | 15 ++++- src/ice_thermo_oce.F90 | 16 +++++ 4 files changed, 156 insertions(+), 6 deletions(-) diff --git a/src/ice_fct.F90 b/src/ice_fct.F90 index 461619e56..af6f9f66e 100755 --- a/src/ice_fct.F90 +++ b/src/ice_fct.F90 @@ -55,6 +55,9 @@ subroutine ice_TG_rhs(mesh) rhs_m(row)=0._WP rhs_a(row)=0._WP rhs_ms(row)=0._WP +#if defined (__oifs) + ths_temp(row)=0._WP +#endif /* (__oifs) */ END DO ! Velocities at nodes do elem=1,myDim_elem2D !assembling rhs over elements @@ -89,6 +92,9 @@ subroutine ice_TG_rhs(mesh) rhs_m(row)=rhs_m(row)+sum(entries*m_ice(elnodes)) rhs_a(row)=rhs_a(row)+sum(entries*a_ice(elnodes)) rhs_ms(row)=rhs_ms(row)+sum(entries*m_snow(elnodes)) +#if defined (__oifs) + rhs_temp(row)=rhs_temp(row)+sum(entries*ice_temp(elnodes)) +#endif /* (__oifs) */ END DO end do end subroutine ice_TG_rhs @@ -113,12 +119,19 @@ subroutine ice_fct_init(mesh) ! Initialization of arrays necessary to implement FCT algorithm allocate(m_icel(n_size), a_icel(n_size), m_snowl(n_size)) ! low-order solutions allocate(dm_ice(n_size), da_ice(n_size), dm_snow(n_size)) ! increments of high - ! order solutions + ! order solutions +#if defined (__oifs) + allocate(m_templ(n_size)) + allocate(dm_temp(n_size)) +#endif /* (__oifs) */ allocate(icefluxes(myDim_elem2D,3)) allocate(icepplus(n_size), icepminus(n_size)) m_icel=0.0_WP a_icel=0.0_WP m_snowl=0.0_WP +#if defined (__oifs) + m_templ=0.0_WP +#endif /* (__oifs) */ ! Fill in the mass matrix call ice_mass_matrix_fill(mesh) @@ -141,6 +154,10 @@ subroutine ice_fct_solve(mesh) call ice_fem_fct(1, mesh) ! m_ice call ice_fem_fct(2, mesh) ! a_ice call ice_fem_fct(3, mesh) ! m_snow +#if defined (__oifs) + call ice_fem_fct(4, mesh) ! ice_temp +#endif /* (__oifs) */ + end subroutine ice_fct_solve ! !---------------------------------------------------------------------------- @@ -188,9 +205,18 @@ subroutine ice_solve_low_order(mesh) m_snowl(row)=(rhs_ms(row)+gamma*sum(mass_matrix(clo:clo2)* & m_snow(location(1:cn))))/area(1,row) + & (1.0_WP-gamma)*m_snow(row) +#if defined (__oifs) + m_templ(row)=(rhs_temp(row)+gamma*sum(mass_matrix(clo:clo2)* & + ice_temp(location(1:cn))))/area(1,row) + & + (1.0_WP-gamma)*ice_temp(row) +#endif /* (__oifs) */ end do call exchange_nod(m_icel,a_icel,m_snowl) +#if defined (__oifs) + call exchange_nod(m_templ) +#endif /* (__oifs) */ + ! Low-order solution must be known to neighbours end subroutine ice_solve_low_order ! @@ -219,9 +245,16 @@ subroutine ice_solve_high_order(mesh) dm_ice(row)=rhs_m(row)/area(1,row) da_ice(row)=rhs_a(row)/area(1,row) dm_snow(row)=rhs_ms(row)/area(1,row) +#if defined (__oifs) + dm_temp(row)=rhs_temp(row)/area(1,row) +#endif /* (__oifs) */ end do + call exchange_nod(dm_ice, da_ice, dm_snow) +#if defined (__oifs) + call exchange_nod(dm_temp) +#endif /* (__oifs) */ !iterate do n=1,num_iter_solve-1 do row=1,myDim_nod2D @@ -234,15 +267,26 @@ subroutine ice_solve_high_order(mesh) rhs_new=rhs_a(row) - sum(mass_matrix(clo:clo2)*da_ice(location(1:cn))) a_icel(row)=da_ice(row)+rhs_new/area(1,row) rhs_new=rhs_ms(row) - sum(mass_matrix(clo:clo2)*dm_snow(location(1:cn))) - m_snowl(row)=dm_snow(row)+rhs_new/area(1,row) + m_snowl(row)=dm_snow(row)+rhs_new/area(1,row) +#if defined (__oifs) + rhs_new=rhs_temp(row) - sum(mass_matrix(clo:clo2)*dm_temp(location(1:cn))) + m_templ(row)=dm_temp(row)+rhs_new/area(1,row) +#endif /* (__oifs) */ end do do row=1,myDim_nod2D dm_ice(row)=m_icel(row) da_ice(row)=a_icel(row) dm_snow(row)=m_snowl(row) +#if defined (__oifs) + dm_temp(row)=m_templ(row) +#endif /* (__oifs) */ end do call exchange_nod(dm_ice, da_ice, dm_snow) +#if defined (__oifs) + call exchange_nod(dm_temp) +#endif /* (__oifs) */ + end do end subroutine ice_solve_high_order ! @@ -317,6 +361,14 @@ subroutine ice_fem_fct(tr_array_id, mesh) dm_snow(elnodes)))*(vol/area(1,elnodes(q)))/12.0_WP end do end if +#if defined (__oifs) + if (tr_array_id==4) then + do q=1,3 + icefluxes(elem,q)=-sum(icoef(:,q)*(gamma*ice_temp(elnodes) + & + dm_temp(elnodes)))*(vol/area(1,elnodes(q)))/12.0_WP + end do + end if +#endif /* (__oifs) */ end do !========================== @@ -363,7 +415,19 @@ subroutine ice_fem_fct(tr_array_id, mesh) end do end if - +#if defined (__oifs) + if (tr_array_id==4) then + do row=1, myDim_nod2D + n=nn_num(row) + tmax(row)=maxval(m_templ(nn_pos(1:n,row))) + tmin(row)=minval(m_templ(nn_pos(1:n,row))) + ! Admissible increments + tmax(row)=tmax(row)-m_templ(row) + tmin(row)=tmin(row)-m_templ(row) + end do + end if +#endif /* (__oifs) */ + !========================= ! Sums of positive/negative fluxes to node row !========================= @@ -448,6 +512,7 @@ subroutine ice_fem_fct(tr_array_id, mesh) end do end do end if + if(tr_array_id==3) then do n=1,myDim_nod2D m_snow(n)=m_snowl(n) @@ -460,9 +525,28 @@ subroutine ice_fem_fct(tr_array_id, mesh) end do end do end if + +#if defined (__oifs) + if(tr_array_id==4) then + do n=1,myDim_nod2D + ice_temp(n)=m_templ(n) + end do + do elem=1, myDim_elem2D + elnodes=elem2D_nodes(:,elem) + do q=1,3 + n=elnodes(q) + ice_temp(n)=ice_temp(n)+icefluxes(elem,q) + end do + end do + end if +#endif /* (__oifs) */ call exchange_nod(m_ice, a_ice, m_snow) +#if defined (__oifs) + call exchange_nod(ice_temp) +#endif /* (__oifs) */ + deallocate(tmin, tmax) end subroutine ice_fem_fct ! @@ -545,7 +629,7 @@ subroutine ice_TG_rhs_div(mesh) implicit none real(kind=WP) :: diff, entries(3), um, vm, vol, dx(3), dy(3) integer :: n, q, row, elem, elnodes(3) - real(kind=WP) :: c1, c2, c3, c4, cx1, cx2, cx3, entries2(3) + real(kind=WP) :: c1, c2, c3, c4, cx1, cx2, cx3, cx4, entries2(3) type(t_mesh), intent(in) , target :: mesh #include "associate_mesh.h" @@ -559,9 +643,15 @@ subroutine ice_TG_rhs_div(mesh) rhs_m(row)=0.0_WP rhs_a(row)=0.0_WP rhs_ms(row)=0.0_WP +#if defined (__oifs) + rhs_temp(row)=0.0_WP +#endif /* (__oifs) */ rhs_mdiv(row)=0.0_WP rhs_adiv(row)=0.0_WP - rhs_msdiv(row)=0.0_WP + rhs_msdiv(row)=0.0_WP +#if defined (__oifs) + rhs_tempdiv(row)=0.0_WP +#endif /* (__oifs) */ END DO do elem=1,myDim_elem2D !assembling rhs over elements !! elem=myList_elem2D(m) @@ -592,12 +682,21 @@ subroutine ice_TG_rhs_div(mesh) cx1=vol*ice_dt*c4*(sum(m_ice(elnodes))+m_ice(elnodes(n))+sum(entries2*m_ice(elnodes)))/12.0_WP cx2=vol*ice_dt*c4*(sum(a_ice(elnodes))+a_ice(elnodes(n))+sum(entries2*a_ice(elnodes)))/12.0_WP cx3=vol*ice_dt*c4*(sum(m_snow(elnodes))+m_snow(elnodes(n))+sum(entries2*m_snow(elnodes)))/12.0_WP +#if defined (__oifs) + cx4=vol*ice_dt*c4*(sum(ice_temp(elnodes))+ice_temp(elnodes(n))+sum(entries2*ice_temp(elnodes)))/12.0_WP +#endif /* (__oifs) */ rhs_m(row)=rhs_m(row)+sum(entries*m_ice(elnodes))+cx1 rhs_a(row)=rhs_a(row)+sum(entries*a_ice(elnodes))+cx2 rhs_ms(row)=rhs_ms(row)+sum(entries*m_snow(elnodes))+cx3 +#if defined (__oifs) + rhs_temp(row)=rhs_temp(row)+sum(entries*ice_temp(elnodes))+cx4 +#endif /* (__oifs) */ rhs_mdiv(row)=rhs_mdiv(row)-cx1 rhs_adiv(row)=rhs_adiv(row)-cx2 rhs_msdiv(row)=rhs_msdiv(row)-cx3 +#if defined (__oifs) + rhs_tempdiv(row)=rhs_tempdiv(row)-cx4 +#endif /* (__oifs) */ END DO end do @@ -631,10 +730,16 @@ subroutine ice_update_for_div(mesh) dm_ice(row) =rhs_mdiv(row) /area(1,row) da_ice(row) =rhs_adiv(row) /area(1,row) dm_snow(row)=rhs_msdiv(row)/area(1,row) +#if defined (__oifs) + dm_temp(row)=rhs_tempdiv(row)/area(1,row) +#endif /* (__oifs) */ end do call exchange_nod(dm_ice) call exchange_nod(da_ice) call exchange_nod(dm_snow) +#if defined (__oifs) + call exchange_nod(dm_temp) +#endif /* (__oifs) */ !iterate do n=1,num_iter_solve-1 do row=1,myDim_nod2D @@ -649,19 +754,33 @@ subroutine ice_update_for_div(mesh) a_icel(row)=da_ice(row)+rhs_new/area(1,row) rhs_new=rhs_msdiv(row) - sum(mass_matrix(clo:clo2)*dm_snow(location(1:cn))) m_snowl(row)=dm_snow(row)+rhs_new/area(1,row) +#if defined (__oifs) + rhs_new=rhs_tempdiv(row) - sum(mass_matrix(clo:clo2)*dm_temp(location(1:cn))) + m_templ(row)=dm_temp(row)+rhs_new/area(1,row) +#endif /* (__oifs) */ end do do row=1,myDim_nod2D !! row=myList_nod2D(m) dm_ice(row)=m_icel(row) da_ice(row)=a_icel(row) dm_snow(row)=m_snowl(row) +#if defined (__oifs) + dm_temp(row)=m_templ(row) +#endif /* (__oifs) */ end do call exchange_nod(dm_ice) call exchange_nod(da_ice) call exchange_nod(dm_snow) +#if defined (__oifs) + call exchange_nod(dm_temp) +#endif /* (__oifs) */ end do m_ice=m_ice+dm_ice a_ice=a_ice+da_ice m_snow=m_snow+dm_snow +#if defined (__oifs) + ice_temp=ice_temp+dm_temp +#endif /* (__oifs) */ + end subroutine ice_update_for_div ! ============================================================= diff --git a/src/ice_modules.F90 b/src/ice_modules.F90 index 1e5467d18..2b0c1a508 100755 --- a/src/ice_modules.F90 +++ b/src/ice_modules.F90 @@ -74,6 +74,8 @@ MODULE i_ARRAYS real(kind=WP),target, allocatable, dimension(:) :: tmp_oce_heat_flux, tmp_ice_heat_flux !temporary flux fields !(for flux correction) + REAL(kind=WP), ALLOCATABLE, DIMENSION(:) :: rhs_temp, m_templ, dm_temp, rhs_tempdiv + #endif /* (__oasis) */ REAL(kind=WP), ALLOCATABLE, DIMENSION(:) :: S_oc_array, T_oc_array diff --git a/src/ice_setup_step.F90 b/src/ice_setup_step.F90 index a21efdce9..5d9a48f0b 100755 --- a/src/ice_setup_step.F90 +++ b/src/ice_setup_step.F90 @@ -133,8 +133,11 @@ subroutine ice_array_setup(mesh) allocate(tmp_oce_heat_flux(n_size), tmp_ice_heat_flux(n_size)) #if defined (__oifs) allocate(ice_alb(n_size), ice_temp(n_size)) + allocate(rhs_tempdiv(n_size), rhs_temp(n_size)) ice_alb=0._WP ice_temp=0._WP + rhs_tempdiv=0._WP + rhs_temp=0._WP #endif /* (__oifs) */ oce_heat_flux=0._WP ice_heat_flux=0._WP @@ -179,12 +182,22 @@ subroutine ice_timestep(step, mesh) ! call ice_fct_solve ! call cut_off ! new FCT routines from Sergey Danilov 08.05.2018 +#if defined (__oifs) + do i=1,myDim_nod2D+eDim_nod2D + ice_temp(i) = ice_temp(i)*a_ice(i) + end do +#endif /* (__oifs) */ if (flag_debug .and. mype==0) print *, achar(27)//'[36m'//' --> call ice_TG_rhs_div...'//achar(27)//'[0m' call ice_TG_rhs_div(mesh) if (flag_debug .and. mype==0) print *, achar(27)//'[36m'//' --> call ice_fct_solve...'//achar(27)//'[0m' call ice_fct_solve(mesh) if (flag_debug .and. mype==0) print *, achar(27)//'[36m'//' --> call ice_update_for_div...'//achar(27)//'[0m' call ice_update_for_div(mesh) +#if defined (__oifs) + do i=1,myDim_nod2D+eDim_nod2D + if (a_ice(i)>0.0_WP) ice_temp(i) = ice_temp(i)/a_ice(i) + end do +#endif /* (__oifs) */ if (flag_debug .and. mype==0) print *, achar(27)//'[36m'//' --> call cut_off...'//achar(27)//'[0m' call cut_off t2=MPI_Wtime() @@ -198,7 +211,7 @@ subroutine ice_timestep(step, mesh) write(*,*) '___ICE STEP EXECUTION TIMES____________________________' write(*,"(A, ES10.3)") ' Ice Dyn. :', t1-t0 write(*,"(A, ES10.3)") ' Ice Advect. :', t2-t1 - write(*,"(A, ES10.3)") ' Ice Thermodyn. :', t3-t2 + write(*,"(A, ES10.3)") ' Ice Thermodyn. :', t3-t2 write(*,*) ' _______________________________' write(*,"(A, ES10.3)") ' Ice TOTAL :', t3-t0 write(*,*) diff --git a/src/ice_thermo_oce.F90 b/src/ice_thermo_oce.F90 index 872440879..53d6207c9 100755 --- a/src/ice_thermo_oce.F90 +++ b/src/ice_thermo_oce.F90 @@ -10,11 +10,27 @@ subroutine cut_off() where(a_ice<0.1e-8_WP) a_ice=0.0_WP + m_ice=0.0_WP + m_snow=0.0_WP +#if defined (__oifs) + ice_temp=273.15_WP +#endif /* (__oifs) */ end where where(m_ice<0.1e-8_WP) m_ice=0.0_WP + m_snow=0.0_WP + a_ice=0.0_WP +#if defined (__oifs) + ice_temp=273.15_WP +#endif /* (__oifs) */ +end where + +#if defined (__oifs) +where(ice_temp>273.15_WP) + ice_temp=273.15_WP end where +#endif /* (__oifs) */ end subroutine cut_off #if !defined (__oasis) From 8cd5433be1db35e952c0db267841823d4e369fe1 Mon Sep 17 00:00:00 2001 From: JanStreffing Date: Fri, 15 May 2020 22:48:52 +0200 Subject: [PATCH 11/26] correct area and some small min hight fixes --- src/cpl_driver.F90 | 10 +++++----- src/ice_thermo_cpl.F90 | 4 ++-- src/io_meandata.F90 | 5 +++-- src/oce_mesh.F90 | 1 + 4 files changed, 11 insertions(+), 9 deletions(-) diff --git a/src/cpl_driver.F90 b/src/cpl_driver.F90 index a2936341e..0f807a538 100755 --- a/src/cpl_driver.F90 +++ b/src/cpl_driver.F90 @@ -222,7 +222,7 @@ subroutine cpl_oasis3mct_define_unstr(mesh) real(kind=WP), allocatable :: all_x_coords(:, :) ! longitude coordinates real(kind=WP), allocatable :: all_y_coords(:, :) ! latitude coordinates - real(kind=WP), allocatable :: all_elem_area(:,:) + real(kind=WP), allocatable :: all_area(:,:) #include "associate_mesh.h" @@ -304,11 +304,11 @@ subroutine cpl_oasis3mct_define_unstr(mesh) if (mype .eq. localroot) then ALLOCATE(all_x_coords(number_of_all_points, 1)) ALLOCATE(all_y_coords(number_of_all_points, 1)) - ALLOCATE(all_elem_area(number_of_all_points, 1)) + ALLOCATE(all_area(number_of_all_points, 1)) else ALLOCATE(all_x_coords(1, 1)) ALLOCATE(all_y_coords(1, 1)) - ALLOCATE(all_elem_area(1, 1)) + ALLOCATE(all_area(1, 1)) endif displs_from_all_pes(1) = 0 @@ -331,7 +331,7 @@ subroutine cpl_oasis3mct_define_unstr(mesh) if (mype .eq. 0) then print *, 'FESOM before 3rd GatherV' endif - CALL MPI_GATHERV(elem_area, my_number_of_points, MPI_DOUBLE_PRECISION, all_elem_area, & + CALL MPI_GATHERV(area, my_number_of_points, MPI_DOUBLE_PRECISION, all_area, & counts_from_all_pes, displs_from_all_pes, MPI_DOUBLE_PRECISION, localroot, MPI_COMM_FESOM, ierror) if (mype .eq. 0) then @@ -358,7 +358,7 @@ subroutine cpl_oasis3mct_define_unstr(mesh) DEALLOCATE(unstr_mask) print *, 'FESOM before write area' - CALL oasis_write_area(grid_name, number_of_all_points, 1, all_elem_area) + CALL oasis_write_area(grid_name, number_of_all_points, 1, all_area) end if print *, 'FESOM before terminate_grids_writing' diff --git a/src/ice_thermo_cpl.F90 b/src/ice_thermo_cpl.F90 index 52947d06d..923b363da 100644 --- a/src/ice_thermo_cpl.F90 +++ b/src/ice_thermo_cpl.F90 @@ -135,7 +135,7 @@ subroutine thermodynamics(mesh) ! energy fluxes ---! t = ice_temp(inod) if(A>Aimin) then - call ice_surftemp(max(h,0.05)/(max(A,Aimin)),hsn/(max(A,Aimin)),a2ihf,t) + call ice_surftemp(max(h/(max(A,Aimin)),0.05),hsn/(max(A,Aimin)),a2ihf,t) ice_temp(inod) = t else ice_temp(inod) = 275.15_WP @@ -465,7 +465,7 @@ subroutine ice_surftemp(h,hsn,a2ihf,t) snicecond = con/consn*rhowat/rhosno ! equivalence fraction thickness of ice/snow - zsniced=MAX(h,hmin)+snicecond*hsn ! Ice + Snow-Ice-equivalent thickness [m] + zsniced=h+snicecond*hsn ! Ice + Snow-Ice-equivalent thickness [m] zicefl=con*ctfreez/zsniced ! Conductive heat flux through sea ice [W/m²] hcapice=rhoice*cpice*dice ! heat capacity of upper 0.05 cm sea ice layer [J/(m²K)] zcpdt=hcapice/dt ! Energy required to change temperature of top ice "layer" [J/(sm²K)] diff --git a/src/io_meandata.F90 b/src/io_meandata.F90 index 664237fd2..6c8b229dd 100644 --- a/src/io_meandata.F90 +++ b/src/io_meandata.F90 @@ -383,8 +383,9 @@ subroutine ini_mean_io(mesh) #if defined (__oifs) call def_stream(nod2D, myDim_nod2D, 'alb', 'ice albedo', 'none', ice_alb(:), 1, 'm', i_real4, mesh) - call def_stream(nod2D, myDim_nod2D, 'ist', 'ice surface temperature', 'K', ice_temp(:), 1, 'd', i_real4, mesh) - call def_stream(nod2D, myDim_nod2D, 'qsi', 'ice heat flux', 'W/m^2', ice_heat_flux(:), 1, 'd', i_real4, mesh) + call def_stream(nod2D, myDim_nod2D, 'ist', 'ice surface temperature', 'K', ice_temp(:), 1, 'm', i_real4, mesh) + call def_stream(nod2D, myDim_nod2D, 'qsi', 'ice heat flux', 'W/m^2', ice_heat_flux(:), 1, 's', i_real4, mesh) + call def_stream(nod2D, myDim_nod2D, 'qso', 'oce heat flux', 'W/m^2', oce_heat_flux(:), 1, 's', i_real4, mesh) #endif if (mix_scheme_nmb==5 .or. mix_scheme_nmb==56) then diff --git a/src/oce_mesh.F90 b/src/oce_mesh.F90 index 106536ec0..182f0d576 100755 --- a/src/oce_mesh.F90 +++ b/src/oce_mesh.F90 @@ -1298,6 +1298,7 @@ SUBROUTINE mesh_areas(mesh) END DO END DO END DO +write(*,*)(mesh%elem_area) ! Only areas through which there is exchange are counted From dc85588e975734ea2dec2cb3957d49ad36dea887 Mon Sep 17 00:00:00 2001 From: JanStreffing Date: Mon, 18 May 2020 12:40:31 +0200 Subject: [PATCH 12/26] debug output --- src/io_meandata.F90 | 4 ++-- src/oce_mesh.F90 | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/io_meandata.F90 b/src/io_meandata.F90 index 6c8b229dd..5a4fc31b1 100644 --- a/src/io_meandata.F90 +++ b/src/io_meandata.F90 @@ -384,8 +384,8 @@ subroutine ini_mean_io(mesh) #if defined (__oifs) call def_stream(nod2D, myDim_nod2D, 'alb', 'ice albedo', 'none', ice_alb(:), 1, 'm', i_real4, mesh) call def_stream(nod2D, myDim_nod2D, 'ist', 'ice surface temperature', 'K', ice_temp(:), 1, 'm', i_real4, mesh) - call def_stream(nod2D, myDim_nod2D, 'qsi', 'ice heat flux', 'W/m^2', ice_heat_flux(:), 1, 's', i_real4, mesh) - call def_stream(nod2D, myDim_nod2D, 'qso', 'oce heat flux', 'W/m^2', oce_heat_flux(:), 1, 's', i_real4, mesh) + call def_stream(nod2D, myDim_nod2D, 'qsi', 'ice heat flux', 'W/m^2', ice_heat_flux(:), 1, 'm', i_real4, mesh) + call def_stream(nod2D, myDim_nod2D, 'qso', 'oce heat flux', 'W/m^2', oce_heat_flux(:), 1, 'm', i_real4, mesh) #endif if (mix_scheme_nmb==5 .or. mix_scheme_nmb==56) then diff --git a/src/oce_mesh.F90 b/src/oce_mesh.F90 index 182f0d576..0bc08bd34 100755 --- a/src/oce_mesh.F90 +++ b/src/oce_mesh.F90 @@ -1298,7 +1298,7 @@ SUBROUTINE mesh_areas(mesh) END DO END DO END DO -write(*,*)(mesh%elem_area) +write(*,*)(mesh%area) ! Only areas through which there is exchange are counted From 19b410210397176c4171637a291a7dd2bd9ac788 Mon Sep 17 00:00:00 2001 From: Dmitri Sidorenko Date: Mon, 18 May 2020 12:45:47 +0200 Subject: [PATCH 13/26] change in area computation for oasis --- src/cpl_driver.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/cpl_driver.F90 b/src/cpl_driver.F90 index 0f807a538..6107cda98 100755 --- a/src/cpl_driver.F90 +++ b/src/cpl_driver.F90 @@ -331,7 +331,7 @@ subroutine cpl_oasis3mct_define_unstr(mesh) if (mype .eq. 0) then print *, 'FESOM before 3rd GatherV' endif - CALL MPI_GATHERV(area, my_number_of_points, MPI_DOUBLE_PRECISION, all_area, & + CALL MPI_GATHERV(area(1,:), my_number_of_points, MPI_DOUBLE_PRECISION, all_area, & counts_from_all_pes, displs_from_all_pes, MPI_DOUBLE_PRECISION, localroot, MPI_COMM_FESOM, ierror) if (mype .eq. 0) then From 6cf760c384b8a8cad6c89a640c0f742c4bdcec02 Mon Sep 17 00:00:00 2001 From: JanStreffing Date: Sun, 24 May 2020 11:35:22 +0200 Subject: [PATCH 14/26] using different lead closing parameters for nh/sh --- src/ice_setup_step.F90 | 2 +- src/ice_thermo_cpl.F90 | 16 +++++++++------- src/io_meandata.F90 | 21 +++++++++++++++------ src/oce_mesh.F90 | 1 - 4 files changed, 25 insertions(+), 15 deletions(-) diff --git a/src/ice_setup_step.F90 b/src/ice_setup_step.F90 index 48c5ed59f..5fc87edee 100755 --- a/src/ice_setup_step.F90 +++ b/src/ice_setup_step.F90 @@ -135,7 +135,7 @@ subroutine ice_array_setup(mesh) allocate(ice_alb(n_size), ice_temp(n_size)) allocate(rhs_tempdiv(n_size), rhs_temp(n_size)) ice_alb=0.6_WP - ice_temp=275.15_WP + ice_temp=265.15_WP rhs_tempdiv=0._WP rhs_temp=0._WP #endif /* (__oifs) */ diff --git a/src/ice_thermo_cpl.F90 b/src/ice_thermo_cpl.F90 index 923b363da..079b0dd8a 100644 --- a/src/ice_thermo_cpl.F90 +++ b/src/ice_thermo_cpl.F90 @@ -118,15 +118,17 @@ subroutine thermodynamics(mesh) rsf = 0._WP #endif +#if defined (__oifs) !---- different lead closing parameter for NH and SH call r2g(geolon, geolat, coord_nod2d(1,inod), coord_nod2d(2,inod)) -! if (geolat.lt.0.) then -! h0min = 1.0 -! h0max = 1.5 -! else -! h0min = 0.75 -! h0max = 1.0 -! endif + if (geolat.lt.0.) then + h0min = 0.75 + h0max = 1.0 + else + h0min = 0.5 + h0max = 0.75 + endif +#endif /* (__oifs) */ call ice_growth #if defined (__oifs) diff --git a/src/io_meandata.F90 b/src/io_meandata.F90 index 5a4fc31b1..dc385568e 100644 --- a/src/io_meandata.F90 +++ b/src/io_meandata.F90 @@ -337,6 +337,20 @@ subroutine ini_mean_io(mesh) CASE ('pgf_y ') call def_stream((/nl-1, elem2D/), (/nl-1, myDim_elem2D/), 'pgf_y', 'meridional pressure gradient force', 'm/s^2', pgf_y(:,:), io_list(i)%freq, io_list(i)%unit, io_list(i)%precision, mesh) !___________________________________________________________________________________________________________________________________ + +#if defined (__oifs) +CASE ('alb ') + call def_stream(nod2D, myDim_nod2D, 'alb', 'ice albedo', 'none', ice_alb(:), io_list(i)%freq, io_list(i)%unit, io_list(i)%precision, mesh) +CASE ('ist ') + call def_stream(nod2D, myDim_nod2D, 'ist', 'ice surface temperature', 'K', ice_temp(:), io_list(i)%freq, io_list(i)%unit, io_list(i)%precision, mesh) +CASE ('qsi ') + call def_stream(nod2D, myDim_nod2D, 'qsi', 'ice heat flux', 'W/m^2', ice_heat_flux(:), io_list(i)%freq, io_list(i)%unit, io_list(i)%precision, mesh) +CASE ('qso ') + call def_stream(nod2D, myDim_nod2D, 'qso', 'oce heat flux', 'W/m^2', oce_heat_flux(:), io_list(i)%freq, io_list(i)%unit, io_list(i)%precision, mesh) +#endif +!___________________________________________________________________________________________________________________________________ + + CASE DEFAULT if (mype==0) write(*,*) 'stream ', io_list(i)%id, ' is not defined !' END SELECT @@ -381,12 +395,6 @@ subroutine ini_mean_io(mesh) if (sel_forcvar(12)==0) call def_stream(elem2D, myDim_elem2D, 'ty_sur', 'meridional wind stress to ocean','m/s2', stress_surf(2, 1:myDim_elem2D),1, 'm', i_real4, mesh) ; sel_forcvar(12)=1 end if -#if defined (__oifs) - call def_stream(nod2D, myDim_nod2D, 'alb', 'ice albedo', 'none', ice_alb(:), 1, 'm', i_real4, mesh) - call def_stream(nod2D, myDim_nod2D, 'ist', 'ice surface temperature', 'K', ice_temp(:), 1, 'm', i_real4, mesh) - call def_stream(nod2D, myDim_nod2D, 'qsi', 'ice heat flux', 'W/m^2', ice_heat_flux(:), 1, 'm', i_real4, mesh) - call def_stream(nod2D, myDim_nod2D, 'qso', 'oce heat flux', 'W/m^2', oce_heat_flux(:), 1, 'm', i_real4, mesh) -#endif if (mix_scheme_nmb==5 .or. mix_scheme_nmb==56) then ! TKE diagnostic @@ -502,6 +510,7 @@ subroutine ini_mean_io(mesh) if (sel_forcvar( 5)==0) call def_stream(nod2D , myDim_nod2D , 'prec' , 'precicipation rain' , 'm/s' , prec_rain(:) , 1, 'm', i_real4, mesh) if (sel_forcvar( 6)==0) call def_stream(nod2D , myDim_nod2D , 'snow' , 'precicipation snow' , 'm/s' , prec_snow(:) , 1, 'm', i_real4, mesh) if (sel_forcvar( 7)==0) call def_stream(nod2D , myDim_nod2D , 'evap' , 'evaporation' , 'm/s' , evaporation(:) , 1, 'm', i_real4, mesh) + if (sel_forcvar( 7)==0) call def_stream(nod2D , myDim_nod2D , 'subl' , 'sublimation' , 'm/s' , sublimation(:) , 1, 'm', i_real4, mesh) if (sel_forcvar( 8)==0) call def_stream(nod2D , myDim_nod2D , 'swr' , 'short wave radiation' , 'W/m^2', shortwave(:) , 1, 'm', i_real4, mesh) if (sel_forcvar( 9)==0) call def_stream(nod2D , myDim_nod2D , 'lwr' , 'long wave radiation' , 'W/m^2', longwave(:) , 1, 'm', i_real4, mesh) if (sel_forcvar(10)==0) call def_stream(nod2D , myDim_nod2D , 'runoff', 'river runoff' , 'none' , runoff(:) , 1, 'm', i_real4, mesh) diff --git a/src/oce_mesh.F90 b/src/oce_mesh.F90 index 0bc08bd34..106536ec0 100755 --- a/src/oce_mesh.F90 +++ b/src/oce_mesh.F90 @@ -1298,7 +1298,6 @@ SUBROUTINE mesh_areas(mesh) END DO END DO END DO -write(*,*)(mesh%area) ! Only areas through which there is exchange are counted From e43f75baccf4e16861771cb99209dde3dded3c4c Mon Sep 17 00:00:00 2001 From: JanStreffing Date: Fri, 29 May 2020 18:20:50 +0200 Subject: [PATCH 15/26] separate heat to ice and sublimation into north and south hemisphere fluxes --- src/cpl_driver.F90 | 20 +++++++++++++++++++- src/gen_forcing_couple.F90 | 14 ++++++++++++-- 2 files changed, 31 insertions(+), 3 deletions(-) diff --git a/src/cpl_driver.F90 b/src/cpl_driver.F90 index 6107cda98..135381c8f 100755 --- a/src/cpl_driver.F90 +++ b/src/cpl_driver.F90 @@ -24,10 +24,11 @@ module cpl_driver #if defined (__oifs) integer, parameter :: nsend = 5 + integer, parameter :: nrecv = 14 #else integer, parameter :: nsend = 4 -#endif integer, parameter :: nrecv = 12 +#endif integer, dimension(nsend) :: send_id integer, dimension(nrecv) :: recv_id @@ -396,6 +397,22 @@ subroutine cpl_oasis3mct_define_unstr(mesh) ! ... Define symbolic names for transient fields received by the ocean. ! These must be identical to the names specified in the SMIOC file. ! +#if defined (__oifs) + cpl_recv(1) = 'taux_oce' + cpl_recv(2) = 'tauy_oce' + cpl_recv(3) = 'taux_ico' + cpl_recv(4) = 'tauy_ico' + cpl_recv(5) = 'prec_oce' + cpl_recv(6) = 'snow_oce' + cpl_recv(7) = 'evap_oce' + cpl_recv(8) = 'subl_oce_nh' + cpl_recv(9) = 'heat_oce' + cpl_recv(10) = 'heat_ico_nh' + cpl_recv(11) = 'heat_swo' + cpl_recv(12) = 'hydr_oce' + cpl_recv(13) = 'subl_oce_sh' + cpl_recv(14) = 'heat_ico_sh' +#else cpl_recv(1) = 'taux_oce' cpl_recv(2) = 'tauy_oce' cpl_recv(3) = 'taux_ico' @@ -408,6 +425,7 @@ subroutine cpl_oasis3mct_define_unstr(mesh) cpl_recv(10) = 'heat_ico' cpl_recv(11) = 'heat_swo' cpl_recv(12) = 'hydr_oce' +#endif if (mype .eq. 0) then print *, 'FESOM after declaring the transient variables' diff --git a/src/gen_forcing_couple.F90 b/src/gen_forcing_couple.F90 index 3d47500e6..fbffdc0d8 100755 --- a/src/gen_forcing_couple.F90 +++ b/src/gen_forcing_couple.F90 @@ -68,7 +68,7 @@ subroutine update_atm_forcing(istep, mesh) elseif (i.eq.3) then exchange(:) = m_snow(:) ! snow thickness elseif (i.eq.4) then - exchange(:) = ice_temp(:) ! ice surface temperature + exchange(:) = ice_temp(:) ! ice surface temperature elseif (i.eq.5) then exchange(:) = ice_alb(:) ! ice albedo else @@ -180,10 +180,20 @@ subroutine update_atm_forcing(istep, mesh) call force_flux_consv(shortwave, mask, i, 0,action, mesh) elseif (i.eq.12) then if (action) then - runoff(:) = exchange(:) ! runoff + calving + runoff(:) = exchange(:) ! runoff + calving mask=1. call force_flux_consv(runoff, mask, i, 0,action, mesh) end if +#if defined (__oifs) + elseif (i.eq.13) then + if (action) then + sublimation(:) = sublimation(:)+exchange(:) ! Adding sh sublimation on top of nh + end if + elseif (i.eq.14) then + if (action) then + ice_heat_flux(:) = ice_heat_flux+exchange(:) ! Adding sh heat-to-ice on top of nh + end if +#endif end if #ifdef VERBOSE if (mype==0) then From fc7f9296412dfd288b29fadaac678d4362456a61 Mon Sep 17 00:00:00 2001 From: JanStreffing Date: Mon, 1 Jun 2020 22:50:44 +0200 Subject: [PATCH 16/26] adding dynamic limiter, fixing alb select, recv plus minus heat-to-ice flux --- src/cpl_driver.F90 | 9 ++++----- src/gen_forcing_couple.F90 | 8 ++------ src/ice_thermo_cpl.F90 | 2 +- src/ice_thermo_oce.F90 | 6 ++++++ 4 files changed, 13 insertions(+), 12 deletions(-) diff --git a/src/cpl_driver.F90 b/src/cpl_driver.F90 index 135381c8f..bffde7333 100755 --- a/src/cpl_driver.F90 +++ b/src/cpl_driver.F90 @@ -24,7 +24,7 @@ module cpl_driver #if defined (__oifs) integer, parameter :: nsend = 5 - integer, parameter :: nrecv = 14 + integer, parameter :: nrecv = 13 #else integer, parameter :: nsend = 4 integer, parameter :: nrecv = 12 @@ -405,13 +405,12 @@ subroutine cpl_oasis3mct_define_unstr(mesh) cpl_recv(5) = 'prec_oce' cpl_recv(6) = 'snow_oce' cpl_recv(7) = 'evap_oce' - cpl_recv(8) = 'subl_oce_nh' + cpl_recv(8) = 'subl_oce' cpl_recv(9) = 'heat_oce' - cpl_recv(10) = 'heat_ico_nh' + cpl_recv(10) = 'heat_ico_plus' cpl_recv(11) = 'heat_swo' cpl_recv(12) = 'hydr_oce' - cpl_recv(13) = 'subl_oce_sh' - cpl_recv(14) = 'heat_ico_sh' + cpl_recv(13) = 'heat_ico_minus' #else cpl_recv(1) = 'taux_oce' cpl_recv(2) = 'tauy_oce' diff --git a/src/gen_forcing_couple.F90 b/src/gen_forcing_couple.F90 index fbffdc0d8..ad384a37b 100755 --- a/src/gen_forcing_couple.F90 +++ b/src/gen_forcing_couple.F90 @@ -186,12 +186,8 @@ subroutine update_atm_forcing(istep, mesh) end if #if defined (__oifs) elseif (i.eq.13) then - if (action) then - sublimation(:) = sublimation(:)+exchange(:) ! Adding sh sublimation on top of nh - end if - elseif (i.eq.14) then - if (action) then - ice_heat_flux(:) = ice_heat_flux+exchange(:) ! Adding sh heat-to-ice on top of nh + if (action) then + ice_heat_flux(:) = ice_heat_flux(:)+exchange(:) end if #endif end if diff --git a/src/ice_thermo_cpl.F90 b/src/ice_thermo_cpl.F90 index 079b0dd8a..61eadafe8 100644 --- a/src/ice_thermo_cpl.F90 +++ b/src/ice_thermo_cpl.F90 @@ -495,7 +495,7 @@ subroutine ice_albedo(h,hsn,t,alb) ! set albedo ! ice and snow, freezing and melting conditions are distinguished if (h>0.0_WP) then - if (t<0.0_WP) then ! freezing condition + if (t<273.14_WP) then ! freezing condition if (hsn.gt.0.0_WP) then ! snow cover present alb=albsn else ! no snow cover diff --git a/src/ice_thermo_oce.F90 b/src/ice_thermo_oce.F90 index 53d6207c9..a41f4245b 100755 --- a/src/ice_thermo_oce.F90 +++ b/src/ice_thermo_oce.F90 @@ -32,6 +32,12 @@ subroutine cut_off() end where #endif /* (__oifs) */ +#if defined (__oifs) +where(ice_temp < 173.15_WP .and. a_ice >= 0.1e-8_WP) + ice_temp=271.35_WP +end where +#endif /* (__oifs) */ + end subroutine cut_off #if !defined (__oasis) !=================================================================== From fca16e4bffa4719288a6cd5eadacc1df4e35c933 Mon Sep 17 00:00:00 2001 From: JanStreffing Date: Fri, 5 Jun 2020 14:54:48 +0200 Subject: [PATCH 17/26] ready-for-deck --- cmake/FindOASIS.cmake | 2 +- src/cpl_driver.F90 | 5 ++--- src/gen_forcing_couple.F90 | 8 +------- 3 files changed, 4 insertions(+), 11 deletions(-) diff --git a/cmake/FindOASIS.cmake b/cmake/FindOASIS.cmake index acac55f21..cb669ea1a 100644 --- a/cmake/FindOASIS.cmake +++ b/cmake/FindOASIS.cmake @@ -7,5 +7,5 @@ find_library(MCT_Fortran_LIBRARIES mct HINTS ${TOPLEVEL_DIR}/../oasis/build/lib/ find_path(MPEU_Fortran_INCLUDE_DIRECTORIES m_mpout.mod HINTS ${TOPLEVEL_DIR}/../oasis/build/lib/psmile/mct/mpeu) find_library(MPEU_Fortran_LIBRARIES mpeu HINTS ${TOPLEVEL_DIR}/../oasis/build/lib/psmile/mct/mpeu) -find_path(SCRIP_Fortran_INCLUDE_DIRECTORIES remap_bicubic.mod HINTS ${TOPLEVEL_DIR}/../oasis/build/lib/psmile/scrip) +find_path(SCRIP_Fortran_INCLUDE_DIRECTORIES remap_bicubic_reduced.mod HINTS ${TOPLEVEL_DIR}/../oasis/build/lib/psmile/scrip) find_library(SCRIP_Fortran_LIBRARIES scrip HINTS ${TOPLEVEL_DIR}/../oasis/build/lib/psmile/scrip) diff --git a/src/cpl_driver.F90 b/src/cpl_driver.F90 index bffde7333..1c7eac457 100755 --- a/src/cpl_driver.F90 +++ b/src/cpl_driver.F90 @@ -24,7 +24,7 @@ module cpl_driver #if defined (__oifs) integer, parameter :: nsend = 5 - integer, parameter :: nrecv = 13 + integer, parameter :: nrecv = 12 #else integer, parameter :: nsend = 4 integer, parameter :: nrecv = 12 @@ -407,10 +407,9 @@ subroutine cpl_oasis3mct_define_unstr(mesh) cpl_recv(7) = 'evap_oce' cpl_recv(8) = 'subl_oce' cpl_recv(9) = 'heat_oce' - cpl_recv(10) = 'heat_ico_plus' + cpl_recv(10) = 'heat_ico' cpl_recv(11) = 'heat_swo' cpl_recv(12) = 'hydr_oce' - cpl_recv(13) = 'heat_ico_minus' #else cpl_recv(1) = 'taux_oce' cpl_recv(2) = 'tauy_oce' diff --git a/src/gen_forcing_couple.F90 b/src/gen_forcing_couple.F90 index ad384a37b..752641738 100755 --- a/src/gen_forcing_couple.F90 +++ b/src/gen_forcing_couple.F90 @@ -184,13 +184,7 @@ subroutine update_atm_forcing(istep, mesh) 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 - ice_heat_flux(:) = ice_heat_flux(:)+exchange(:) - end if -#endif - end if + end if #ifdef VERBOSE if (mype==0) then write(*,*) 'FESOM RECV: flux ', i, ', max val: ', maxval(exchange) From d4a668f13882b901d27aa4ed79d07ed0c7797a15 Mon Sep 17 00:00:00 2001 From: JanStreffing Date: Wed, 1 Jul 2020 01:06:55 +0200 Subject: [PATCH 18/26] faster io on juwels through chunking --- env.sh | 2 +- src/ice_thermo_cpl.F90 | 22 +++++++++++----------- src/io_meandata.F90 | 7 +++++++ src/io_restart.F90 | 17 ++++++++++++++--- 4 files changed, 33 insertions(+), 15 deletions(-) diff --git a/env.sh b/env.sh index 9530ae699..c7f4a9394 100755 --- a/env.sh +++ b/env.sh @@ -41,7 +41,7 @@ elif [[ $LOGINHOST = ubuntu ]]; then STRATEGY="ubuntu" elif [[ $LOGINHOST = bsc ]]; then STRATEGY="bsc" -elif [[ $LOGINHOST =~ ^juwels[0-9]+\.fz\-juelich\.de$ ]]; then +elif [[ $LOGINHOST =~ ^juwels[0-9]+\.ib\.juwels\.fzj\.de$ ]]; then STRATEGY="juwels" else echo "can not determine environment for host: "$LOGINHOST diff --git a/src/ice_thermo_cpl.F90 b/src/ice_thermo_cpl.F90 index 61eadafe8..eb294944e 100644 --- a/src/ice_thermo_cpl.F90 +++ b/src/ice_thermo_cpl.F90 @@ -118,17 +118,17 @@ subroutine thermodynamics(mesh) rsf = 0._WP #endif -#if defined (__oifs) - !---- different lead closing parameter for NH and SH - call r2g(geolon, geolat, coord_nod2d(1,inod), coord_nod2d(2,inod)) - if (geolat.lt.0.) then - h0min = 0.75 - h0max = 1.0 - else - h0min = 0.5 - h0max = 0.75 - endif -#endif /* (__oifs) */ +!#if defined (__oifs) +! !---- different lead closing parameter for NH and SH +! call r2g(geolon, geolat, coord_nod2d(1,inod), coord_nod2d(2,inod)) +! if (geolat.lt.0.) then +! h0min = 0.75 +! h0max = 1.0 +! else +! h0min = 0.5 +! h0max = 0.75 +! endif +!#endif /* (__oifs) */ call ice_growth #if defined (__oifs) diff --git a/src/io_meandata.F90 b/src/io_meandata.F90 index dc385568e..809627dd9 100644 --- a/src/io_meandata.F90 +++ b/src/io_meandata.F90 @@ -591,6 +591,13 @@ subroutine create_new_file(entry) entry%error_status(c) = nf_put_att_real(entry%ncid, entry%varID, 'scale_factor', NF_REAL, 1, entry%scale_factor); c=c+1 entry%error_status(c) = nf_put_att_real(entry%ncid, entry%varID, 'add_offset', NF_REAL, 1, entry%add_offset); c=c+1 endif + + if (entry%ndim==1) then + entry%error_status(c) = nf_def_var_chunking(entry%ncid, entry%varID, NF_CHUNKED, (/1/)); c=c+1 + elseif (entry%ndim==2) then + entry%error_status(c) = nf_def_var_chunking(entry%ncid, entry%varID, NF_CHUNKED, (/1, entry%glsize(1)/)); c=c+1 + endif + entry%error_status(c)=nf_put_att_text(entry%ncid, entry%varID, 'description', len_trim(entry%description), entry%description); c=c+1 entry%error_status(c)=nf_put_att_text(entry%ncid, entry%varID, 'units', len_trim(entry%units), entry%units); c=c+1 entry%error_status(c)=nf_close(entry%ncid); c=c+1 diff --git a/src/io_restart.F90 b/src/io_restart.F90 index be4d518e9..915da0a4b 100644 --- a/src/io_restart.F90 +++ b/src/io_restart.F90 @@ -307,10 +307,14 @@ subroutine create_new_file(id) do l=1, id%ndim ! list all defined dimensions if (kdim==id%dim(l)%size) dimid(k)=id%dim(l)%code end do -!________write(*,*) kdim, ' -> ', dimid(k)__________________________________ + !write(*,*) "j",j,kdim, ' -> ', dimid(k) end do - id%error_status(c) = nf_def_var(id%ncid, trim(id%var(j)%name), NF_DOUBLE, id%var(j)%ndim+1, & - (/dimid(1:n), id%rec/), id%var(j)%code); c=c+1 + id%error_status(c) = nf_def_var(id%ncid, trim(id%var(j)%name), NF_DOUBLE, id%var(j)%ndim+1, (/dimid(1:n), id%rec/), id%var(j)%code); c=c+1 + if (n==1) then + id%error_status(c)=nf_def_var_chunking(id%ncid, id%var(j)%code, NF_CHUNKED, (/1/)); c=c+1 + elseif (n==2) then + id%error_status(c)=nf_def_var_chunking(id%ncid, id%var(j)%code, NF_CHUNKED, (/1, id%dim(1)%size/)); c=c+1 + end if id%error_status(c)=nf_put_att_text(id%ncid, id%var(j)%code, 'description', len_trim(id%var(j)%longname), id%var(j)%longname); c=c+1 id%error_status(c)=nf_put_att_text(id%ncid, id%var(j)%code, 'units', len_trim(id%var(j)%units), id%var(j)%units); c=c+1 end do @@ -430,6 +434,7 @@ subroutine write_restart(id, istep, mesh) real(kind=WP), allocatable :: aux(:), laux(:) integer :: i, lev, size1, size2, shape integer :: c + real(kind=WP) :: t0, t1, t2, t3 #include "associate_mesh.h" @@ -451,6 +456,7 @@ subroutine write_restart(id, istep, mesh) if (shape==1) then size1=id%var(i)%dims(1) if (mype==0) allocate(aux(size1)) + t0=MPI_Wtime() if (size1==nod2D) call gather_nod (id%var(i)%pt1, aux) if (size1==elem2D) call gather_elem(id%var(i)%pt1, aux) if (mype==0) then @@ -468,11 +474,16 @@ subroutine write_restart(id, istep, mesh) laux=id%var(i)%pt2(lev,:) ! if (size1==nod2D .or. size2==nod2D) call gather_nod (id%var(i)%pt2(lev,:), aux) ! if (size1==elem2D .or. size2==elem2D) call gather_elem(id%var(i)%pt2(lev,:), aux) + t0=MPI_Wtime() if (size1==nod2D .or. size2==nod2D) call gather_nod (laux, aux) if (size1==elem2D .or. size2==elem2D) call gather_elem(laux, aux) + t1=MPI_Wtime() if (mype==0) then id%error_status(c)=nf_put_vara_double(id%ncid, id%var(i)%code, (/lev, 1, id%rec_count/), (/1, size2, 1/), aux, 1); c=c+1 end if + t2=MPI_Wtime() + if (mype==0 .and. size2==nod2D) write(*,*) 'nvar: ', i, 'lev: ', lev, 'gather_nod: ', t1-t0 + if (mype==0 .and. size2==nod2D) write(*,*) 'nvar: ', i, 'lev: ', lev, 'nf_put_var: ', t2-t1 end do deallocate(laux) if (mype==0) deallocate(aux) From a3d6d473e28259e5726ca88fa7ae6a2e466593eb Mon Sep 17 00:00:00 2001 From: JanStreffing Date: Fri, 3 Jul 2020 11:22:53 +0200 Subject: [PATCH 19/26] improvements and more output --- src/io_meandata.F90 | 22 ++++++++++++++++++++-- src/io_restart.F90 | 14 +++++++++----- 2 files changed, 29 insertions(+), 7 deletions(-) diff --git a/src/io_meandata.F90 b/src/io_meandata.F90 index 809627dd9..c9e6039fe 100644 --- a/src/io_meandata.F90 +++ b/src/io_meandata.F90 @@ -593,8 +593,9 @@ subroutine create_new_file(entry) endif if (entry%ndim==1) then - entry%error_status(c) = nf_def_var_chunking(entry%ncid, entry%varID, NF_CHUNKED, (/1/)); c=c+1 + entry%error_status(c) = nf_def_var_chunking(entry%ncid, entry%varID, NF_CONTIGUOUS, (/1/)); c=c+1 elseif (entry%ndim==2) then + write(*,*) "size: ", entry%glsize(1) entry%error_status(c) = nf_def_var_chunking(entry%ncid, entry%varID, NF_CHUNKED, (/1, entry%glsize(1)/)); c=c+1 endif @@ -672,13 +673,14 @@ subroutine write_mean(entry, mesh) type(t_mesh), intent(in) , target :: mesh integer :: i, size1, size2 integer :: c, lev + real(kind=WP) :: t0, t1, t2, t3 #include "associate_mesh.h" ! Serial output implemented so far if (mype==0) then c=1 - write(*,*) 'writing mean record for ', trim(entry%name), '; rec. count = ', entry%rec_count + write(*,*) 'writing mean record for ', trim(entry%name), '; rec. count = ', entry%rec_count, '; bytes = ', entry%accuracy entry%error_status(c)=nf_open(entry%filename, nf_write, entry%ncid); c=c+1 entry%error_status(c)=nf_put_vara_double(entry%ncid, entry%Tid, entry%rec_count, 1, ctime, 1); c=c+1 end if @@ -693,16 +695,22 @@ subroutine write_mean(entry, mesh) if (mype==0) then entry%error_status(c)=nf_put_vara_double(entry%ncid, entry%varID, (/1, entry%rec_count/), (/size1, 1/), aux_r8, 1); c=c+1 end if + if (mype==0) deallocate(aux_r8) !___________writing real 4 byte real _________________________________________ elseif (entry%accuracy == i_real4) then if (mype==0) allocate(aux_r4(size1)) + t0=MPI_Wtime() if (size1==nod2D) call gather_nod (entry%local_values_r4(1:entry%lcsize(1),1), aux_r4) if (size1==elem2D) call gather_elem(entry%local_values_r4(1:entry%lcsize(1),1), aux_r4) + t1=MPI_Wtime() if (mype==0) then entry%error_status(c)=nf_put_vara_real(entry%ncid, entry%varID, (/1, entry%rec_count/), (/size1, 1/), aux_r4, 1); c=c+1 end if + t2=MPI_Wtime() + if (mype==0) write(*,*) 'size: ', size1, 'gather_nod: ', t1-t0 + if (mype==0) write(*,*) 'size: ', size1, 'nf_put_var: ', t2-t1 if (mype==0) deallocate(aux_r4) !___________writing real as 2 byte integer _________________________________________ @@ -729,22 +737,32 @@ subroutine write_mean(entry, mesh) if (entry%accuracy == i_real8) then if (mype==0) allocate(aux_r8(size2)) do lev=1, size1 + t0=MPI_Wtime() if (size1==nod2D .or. size2==nod2D) call gather_nod (entry%local_values_r8(lev,1:entry%lcsize(2)), aux_r8) if (size1==elem2D .or. size2==elem2D) call gather_elem(entry%local_values_r8(lev,1:entry%lcsize(2)), aux_r8) + t1=MPI_Wtime() if (mype==0) then entry%error_status(c)=nf_put_vara_double(entry%ncid, entry%varID, (/lev, 1, entry%rec_count/), (/1, size2, 1/), aux_r8, 1); c=c+1 end if + t2=MPI_Wtime() + if (mype==0) write(*,*) 'size: ', size2, 'lev: ', lev, 'gather_nod: ',t1-t0 + if (mype==0) write(*,*) 'size: ', size2, 'lev: ', lev, 'nf_put_var: ',t2-t1 end do if (mype==0) deallocate(aux_r8) !___________writing real 4 byte real _________________________________________ elseif (entry%accuracy == i_real4) then if (mype==0) allocate(aux_r4(size2)) do lev=1, size1 + t0=MPI_Wtime() if (size1==nod2D .or. size2==nod2D) call gather_nod (entry%local_values_r4(lev,1:entry%lcsize(2)), aux_r4) if (size1==elem2D .or. size2==elem2D) call gather_elem(entry%local_values_r4(lev,1:entry%lcsize(2)), aux_r4) + t1=MPI_Wtime() if (mype==0) then entry%error_status(c)=nf_put_vara_real(entry%ncid, entry%varID, (/lev, 1, entry%rec_count/), (/1, size2, 1/), aux_r4, 1); c=c+1 end if + t2=MPI_Wtime() + if (mype==0) write(*,*) 'size: ', size2, 'lev: ', lev, 'gather_nod: ',t1-t0 + if (mype==0) write(*,*) 'size: ', size2, 'lev: ', lev, 'nf_put_var: ',t2-t1 end do if (mype==0) deallocate(aux_r4) !___________writing real as 2 byte integer _________________________________________ diff --git a/src/io_restart.F90 b/src/io_restart.F90 index 915da0a4b..811b37756 100644 --- a/src/io_restart.F90 +++ b/src/io_restart.F90 @@ -310,9 +310,9 @@ subroutine create_new_file(id) !write(*,*) "j",j,kdim, ' -> ', dimid(k) end do id%error_status(c) = nf_def_var(id%ncid, trim(id%var(j)%name), NF_DOUBLE, id%var(j)%ndim+1, (/dimid(1:n), id%rec/), id%var(j)%code); c=c+1 - if (n==1) then - id%error_status(c)=nf_def_var_chunking(id%ncid, id%var(j)%code, NF_CHUNKED, (/1/)); c=c+1 - elseif (n==2) then + !if (n==1) then + ! id%error_status(c)=nf_def_var_chunking(id%ncid, id%var(j)%code, NF_CHUNKED, (/1/)); c=c+1 + if (n==2) then id%error_status(c)=nf_def_var_chunking(id%ncid, id%var(j)%code, NF_CHUNKED, (/1, id%dim(1)%size/)); c=c+1 end if id%error_status(c)=nf_put_att_text(id%ncid, id%var(j)%code, 'description', len_trim(id%var(j)%longname), id%var(j)%longname); c=c+1 @@ -459,9 +459,13 @@ subroutine write_restart(id, istep, mesh) t0=MPI_Wtime() if (size1==nod2D) call gather_nod (id%var(i)%pt1, aux) if (size1==elem2D) call gather_elem(id%var(i)%pt1, aux) + t1=MPI_Wtime() if (mype==0) then id%error_status(c)=nf_put_vara_double(id%ncid, id%var(i)%code, (/1, id%rec_count/), (/size1, 1/), aux, 1); c=c+1 end if + t2=MPI_Wtime() + if (mype==0) write(*,*) 'nvar: ', i, 'size: ', size1, 'gather_nod: ', t1-t0 + if (mype==0) write(*,*) 'nvar: ', i, 'size: ', size1, 'nf_put_var: ', t2-t1 if (mype==0) deallocate(aux) !_______writing 3D fields________________________________________________ elseif (shape==2) then @@ -482,8 +486,8 @@ subroutine write_restart(id, istep, mesh) id%error_status(c)=nf_put_vara_double(id%ncid, id%var(i)%code, (/lev, 1, id%rec_count/), (/1, size2, 1/), aux, 1); c=c+1 end if t2=MPI_Wtime() - if (mype==0 .and. size2==nod2D) write(*,*) 'nvar: ', i, 'lev: ', lev, 'gather_nod: ', t1-t0 - if (mype==0 .and. size2==nod2D) write(*,*) 'nvar: ', i, 'lev: ', lev, 'nf_put_var: ', t2-t1 + if (mype==0) write(*,*) 'nvar: ', i, 'size: ', size2, 'lev: ', lev, 'gather_nod: ', t1-t0 + if (mype==0) write(*,*) 'nvar: ', i, 'size: ', size2, 'lev: ', lev, 'nf_put_var: ', t2-t1 end do deallocate(laux) if (mype==0) deallocate(aux) From ed0fb746149d98973b93bc22f633423521d393c9 Mon Sep 17 00:00:00 2001 From: JanStreffing Date: Sat, 26 Sep 2020 11:14:08 +0200 Subject: [PATCH 20/26] additional restart fields only for oifs coupled setup --- src/io_restart.F90 | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/io_restart.F90 b/src/io_restart.F90 index 811b37756..a6b647e6f 100644 --- a/src/io_restart.F90 +++ b/src/io_restart.F90 @@ -181,8 +181,10 @@ subroutine ini_ice_io(year, mesh) call def_variable(iid, 'hsnow', (/nod2D/), 'effective snow thickness', 'm', m_snow); call def_variable(iid, 'uice', (/nod2D/), 'zonal velocity', 'm/s', u_ice); call def_variable(iid, 'vice', (/nod2D/), 'meridional velocity', 'm', v_ice); +#if defined (__oifs) call def_variable(iid, 'ice_albedo', (/nod2D/), 'ice albedo', '-', ice_alb); call def_variable(iid, 'ice_temp',(/nod2D/), 'ice surface temperature', 'K', ice_temp); +#endif /* (__oifs) */ end subroutine ini_ice_io ! From 9bef98cad955c92a28c9a7fcf28a1848dc4cbc87 Mon Sep 17 00:00:00 2001 From: JanStreffing Date: Wed, 30 Sep 2020 11:45:23 +0200 Subject: [PATCH 21/26] allocating sublimation in standalone to avoid error in io_meandata --- src/gen_forcing_init.F90 | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/src/gen_forcing_init.F90 b/src/gen_forcing_init.F90 index 3bbb288f9..fe5db75c8 100755 --- a/src/gen_forcing_init.F90 +++ b/src/gen_forcing_init.F90 @@ -52,16 +52,15 @@ subroutine forcing_array_setup(mesh) allocate(u_wind(n2), v_wind(n2)) allocate(Tair(n2), shum(n2)) allocate(runoff(n2), evaporation(n2),ice_sublimation(n2)) + allocate(sublimation(n2)) #if defined (__oasis) - allocate(sublimation(n2), evap_no_ifrac(n2)) allocate(tmp_sublimation(n2),tmp_evap_no_ifrac(n2), tmp_shortwave(n2)) allocate(atm_net_fluxes_north(nrecv), atm_net_fluxes_south(nrecv)) allocate(oce_net_fluxes_north(nrecv), oce_net_fluxes_south(nrecv)) allocate(flux_correction_north(nrecv), flux_correction_south(nrecv)) allocate(flux_correction_total(nrecv)) - sublimation=0.0_WP - evap_no_ifrac=0.0_WP + allocate(evap_no_ifrac(n2)) tmp_sublimation = 0.0_WP tmp_evap_no_ifrac = 0.0_WP tmp_shortwave = 0.0_WP @@ -72,6 +71,7 @@ subroutine forcing_array_setup(mesh) flux_correction_north=0.0_WP flux_correction_south=0.0_WP flux_correction_total=0.0_WP + evap_no_ifrac=0.0_WP #endif @@ -86,6 +86,7 @@ subroutine forcing_array_setup(mesh) Tair=0.0_WP shum=0.0_WP runoff=0.0_WP + sublimation=0.0_WP !!PS allocate(Tair_mo(n2),shum_mo(n2)) !!PS Tair_mo=0.0_WP !!PS From ed0b4c87f6f49738a5738edf74b50e8a00ce110a Mon Sep 17 00:00:00 2001 From: JanStreffing Date: Wed, 30 Sep 2020 11:55:07 +0200 Subject: [PATCH 22/26] rm sub output --- src/gen_forcing_init.F90 | 5 ++--- src/io_meandata.F90 | 1 - 2 files changed, 2 insertions(+), 4 deletions(-) diff --git a/src/gen_forcing_init.F90 b/src/gen_forcing_init.F90 index fe5db75c8..43e9b8516 100755 --- a/src/gen_forcing_init.F90 +++ b/src/gen_forcing_init.F90 @@ -52,15 +52,14 @@ subroutine forcing_array_setup(mesh) allocate(u_wind(n2), v_wind(n2)) allocate(Tair(n2), shum(n2)) allocate(runoff(n2), evaporation(n2),ice_sublimation(n2)) - allocate(sublimation(n2)) #if defined (__oasis) allocate(tmp_sublimation(n2),tmp_evap_no_ifrac(n2), tmp_shortwave(n2)) + allocate(sublimation(n2),evap_no_ifrac(n2)) allocate(atm_net_fluxes_north(nrecv), atm_net_fluxes_south(nrecv)) allocate(oce_net_fluxes_north(nrecv), oce_net_fluxes_south(nrecv)) allocate(flux_correction_north(nrecv), flux_correction_south(nrecv)) allocate(flux_correction_total(nrecv)) - allocate(evap_no_ifrac(n2)) tmp_sublimation = 0.0_WP tmp_evap_no_ifrac = 0.0_WP tmp_shortwave = 0.0_WP @@ -72,6 +71,7 @@ subroutine forcing_array_setup(mesh) flux_correction_south=0.0_WP flux_correction_total=0.0_WP evap_no_ifrac=0.0_WP + sublimation=0.0_WP #endif @@ -86,7 +86,6 @@ subroutine forcing_array_setup(mesh) Tair=0.0_WP shum=0.0_WP runoff=0.0_WP - sublimation=0.0_WP !!PS allocate(Tair_mo(n2),shum_mo(n2)) !!PS Tair_mo=0.0_WP !!PS diff --git a/src/io_meandata.F90 b/src/io_meandata.F90 index 3c504fbe5..4934ffdc8 100644 --- a/src/io_meandata.F90 +++ b/src/io_meandata.F90 @@ -522,7 +522,6 @@ subroutine ini_mean_io(mesh) if (sel_forcvar( 5)==0) call def_stream(nod2D , myDim_nod2D , 'prec' , 'precicipation rain' , 'm/s' , prec_rain(:) , 1, 'm', i_real4, mesh) if (sel_forcvar( 6)==0) call def_stream(nod2D , myDim_nod2D , 'snow' , 'precicipation snow' , 'm/s' , prec_snow(:) , 1, 'm', i_real4, mesh) if (sel_forcvar( 7)==0) call def_stream(nod2D , myDim_nod2D , 'evap' , 'evaporation' , 'm/s' , evaporation(:) , 1, 'm', i_real4, mesh) - if (sel_forcvar( 7)==0) call def_stream(nod2D , myDim_nod2D , 'subl' , 'sublimation' , 'm/s' , sublimation(:) , 1, 'm', i_real4, mesh) if (sel_forcvar( 8)==0) call def_stream(nod2D , myDim_nod2D , 'swr' , 'short wave radiation' , 'W/m^2', shortwave(:) , 1, 'm', i_real4, mesh) if (sel_forcvar( 9)==0) call def_stream(nod2D , myDim_nod2D , 'lwr' , 'long wave radiation' , 'W/m^2', longwave(:) , 1, 'm', i_real4, mesh) if (sel_forcvar(10)==0) call def_stream(nod2D , myDim_nod2D , 'runoff', 'river runoff' , 'none' , runoff(:) , 1, 'm', i_real4, mesh) From 4d2f07354d8e3ab67726fc1fbfb9e5bb1dbd5808 Mon Sep 17 00:00:00 2001 From: Jan Streffing Date: Wed, 30 Sep 2020 16:21:23 +0200 Subject: [PATCH 23/26] moving new awicm3 cutoffs into ifdef_oifs --- src/ice_thermo_oce.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/ice_thermo_oce.F90 b/src/ice_thermo_oce.F90 index a41f4245b..b2b9c2d31 100755 --- a/src/ice_thermo_oce.F90 +++ b/src/ice_thermo_oce.F90 @@ -10,9 +10,9 @@ subroutine cut_off() where(a_ice<0.1e-8_WP) a_ice=0.0_WP +#if defined (__oifs) m_ice=0.0_WP m_snow=0.0_WP -#if defined (__oifs) ice_temp=273.15_WP #endif /* (__oifs) */ end where From eededdc314c7d4d18a743e304cf1de9643175dd3 Mon Sep 17 00:00:00 2001 From: Jan Streffing Date: Wed, 30 Sep 2020 21:54:38 +0200 Subject: [PATCH 24/26] moving new awicm3 cutoffs into ifdef_oifs p2 --- src/ice_thermo_oce.F90 | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/ice_thermo_oce.F90 b/src/ice_thermo_oce.F90 index b2b9c2d31..d78c6d839 100755 --- a/src/ice_thermo_oce.F90 +++ b/src/ice_thermo_oce.F90 @@ -18,10 +18,10 @@ subroutine cut_off() end where where(m_ice<0.1e-8_WP) - m_ice=0.0_WP - m_snow=0.0_WP - a_ice=0.0_WP + m_ice=0.0_WP #if defined (__oifs) + m_snow=0.0_WP + a_ice=0.0_WP ice_temp=273.15_WP #endif /* (__oifs) */ end where From af7fbc36ac52b13ccb0e8b0a3bff04556f44a7f4 Mon Sep 17 00:00:00 2001 From: JanStreffing Date: Thu, 1 Oct 2020 11:28:33 +0200 Subject: [PATCH 25/26] output and restart file timing info only with DEBUG settings --- src/io_meandata.F90 | 9 +++++++++ src/io_restart.F90 | 6 ++++++ 2 files changed, 15 insertions(+) diff --git a/src/io_meandata.F90 b/src/io_meandata.F90 index 4934ffdc8..e4d8a9c3f 100644 --- a/src/io_meandata.F90 +++ b/src/io_meandata.F90 @@ -731,8 +731,11 @@ subroutine write_mean(entry, mesh) entry%error_status(c)=nf_put_vara_real(entry%ncid, entry%varID, (/1, entry%rec_count/), (/size1, 1/), aux_r4, 1); c=c+1 end if t2=MPI_Wtime() +#ifdef DEBUG + ! Timeing information for collecting and writing output file if (mype==0) write(*,*) 'size: ', size1, 'gather_nod: ', t1-t0 if (mype==0) write(*,*) 'size: ', size1, 'nf_put_var: ', t2-t1 +#endif if (mype==0) deallocate(aux_r4) !___________writing real as 2 byte integer _________________________________________ @@ -767,8 +770,11 @@ subroutine write_mean(entry, mesh) entry%error_status(c)=nf_put_vara_double(entry%ncid, entry%varID, (/lev, 1, entry%rec_count/), (/1, size2, 1/), aux_r8, 1); c=c+1 end if t2=MPI_Wtime() +#ifdef DEBUG + ! Timeing information for collecting and writing output file if (mype==0) write(*,*) 'size: ', size2, 'lev: ', lev, 'gather_nod: ',t1-t0 if (mype==0) write(*,*) 'size: ', size2, 'lev: ', lev, 'nf_put_var: ',t2-t1 +#endif end do if (mype==0) deallocate(aux_r8) !___________writing real 4 byte real _________________________________________ @@ -783,8 +789,11 @@ subroutine write_mean(entry, mesh) entry%error_status(c)=nf_put_vara_real(entry%ncid, entry%varID, (/lev, 1, entry%rec_count/), (/1, size2, 1/), aux_r4, 1); c=c+1 end if t2=MPI_Wtime() +#ifdef DEBUG + ! Timeing information for collecting and writing output file if (mype==0) write(*,*) 'size: ', size2, 'lev: ', lev, 'gather_nod: ',t1-t0 if (mype==0) write(*,*) 'size: ', size2, 'lev: ', lev, 'nf_put_var: ',t2-t1 +#endif end do if (mype==0) deallocate(aux_r4) !___________writing real as 2 byte integer _________________________________________ diff --git a/src/io_restart.F90 b/src/io_restart.F90 index bc0d42d92..ff5997d95 100644 --- a/src/io_restart.F90 +++ b/src/io_restart.F90 @@ -466,8 +466,11 @@ subroutine write_restart(id, istep, mesh) id%error_status(c)=nf_put_vara_double(id%ncid, id%var(i)%code, (/1, id%rec_count/), (/size1, 1/), aux, 1); c=c+1 end if t2=MPI_Wtime() +#ifdef DEBUG + ! Timeing information for collecting and writing restart file if (mype==0) write(*,*) 'nvar: ', i, 'size: ', size1, 'gather_nod: ', t1-t0 if (mype==0) write(*,*) 'nvar: ', i, 'size: ', size1, 'nf_put_var: ', t2-t1 +#endif if (mype==0) deallocate(aux) !_______writing 3D fields________________________________________________ elseif (shape==2) then @@ -488,8 +491,11 @@ subroutine write_restart(id, istep, mesh) id%error_status(c)=nf_put_vara_double(id%ncid, id%var(i)%code, (/lev, 1, id%rec_count/), (/1, size2, 1/), aux, 1); c=c+1 end if t2=MPI_Wtime() +#ifdef DEBUG + ! Timeing information for collecting and writing output file if (mype==0) write(*,*) 'nvar: ', i, 'size: ', size2, 'lev: ', lev, 'gather_nod: ', t1-t0 if (mype==0) write(*,*) 'nvar: ', i, 'size: ', size2, 'lev: ', lev, 'nf_put_var: ', t2-t1 +#endif end do deallocate(laux) if (mype==0) deallocate(aux) From 96ae468e6a69f70e8cf6e2bb396d3a8a3783367f Mon Sep 17 00:00:00 2001 From: dsidoren Date: Mon, 5 Oct 2020 21:26:24 +0200 Subject: [PATCH 26/26] Update io_restart.F90 do not check for errors if NF_CHUNKED does not exist in the netcdf library, otherwise crashes on several machines (i.e. ollie) --- src/io_restart.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/io_restart.F90 b/src/io_restart.F90 index ff5997d95..d65cb87f0 100644 --- a/src/io_restart.F90 +++ b/src/io_restart.F90 @@ -315,7 +315,7 @@ subroutine create_new_file(id) !if (n==1) then ! id%error_status(c)=nf_def_var_chunking(id%ncid, id%var(j)%code, NF_CHUNKED, (/1/)); c=c+1 if (n==2) then - id%error_status(c)=nf_def_var_chunking(id%ncid, id%var(j)%code, NF_CHUNKED, (/1, id%dim(1)%size/)); c=c+1 + id%error_status(c)=nf_def_var_chunking(id%ncid, id%var(j)%code, NF_CHUNKED, (/1, id%dim(1)%size/)); ! c=c+1 end if id%error_status(c)=nf_put_att_text(id%ncid, id%var(j)%code, 'description', len_trim(id%var(j)%longname), id%var(j)%longname); c=c+1 id%error_status(c)=nf_put_att_text(id%ncid, id%var(j)%code, 'units', len_trim(id%var(j)%units), id%var(j)%units); c=c+1