Mercurial > hgweb > burst-coupling
changeset 89:85cec98f272f
Several bug fixes and improvements:
+ correct sinking speeds (per timestep) wpoco etc. for off-line (sediment only);
+ SED_OFFLINE_DEBUG saves extra bottom water output;
+ start off-line (sediment only) not before nday==2 of MICOM (otherwise last
output is stored after spin-up);
+ rewind date more consistenly (and only for full MICOM/HAMOCC);
+ use INTENT(IN) for all *_ variables (and so use original vars for full);
+ sedfluxo may be updated also in off-line case; not sure if it makes a
difference.
author | Marco van Hulten <marco@hulten.org> |
---|---|
date | Fri, 26 Jul 2019 16:03:00 +0200 |
parents | 4d0b1c84ebea |
children | 918a4a892a21 |
files | beleg_bgc.F90 mo_biomod.F90 mo_sedmnt_offline.F90 mod_dia.F powach.F90 sediment_step.F90 |
diffstat | 6 files changed, 114 insertions(+), 31 deletions(-) [+] |
line wrap: on
line diff
--- a/beleg_bgc.F90 Tue Jul 02 11:02:37 2019 +0200 +++ b/beleg_bgc.F90 Fri Jul 26 16:03:00 2019 +0200 @@ -191,10 +191,20 @@ wpoc = 5.*dtb !m/d iris : 5. wcal = 30.*dtb !m/d wopal = 30.*dtb !m/d iris : 60 +# ifdef SED_OFFLINE + wpoco = 5.*dto !m/d iris : 5. + wcalo = 30.*dto !m/d + wopalo= 30.*dto !m/d iris : 60 +# endif #ifdef WLIN wmin = 1.*dtb !m/d minimum sinking speed wmax = 60.*dtb !m/d maximum sinking speed wlin = 60./2400.*dtb !m/d/m constant describing incr. with depth, r/a=1.0 +# ifdef SED_OFFLINE + wmino = 1.*dto !m/d minimum sinking speed + wmaxo = 60.*dto !m/d maximum sinking speed + wlino = 60./2400.*dto !m/d/m constant describing incr. with depth, r/a=1.0 +# endif #endif @@ -410,6 +420,9 @@ & * (claydens - 1025.) / 1.567 * 1000. & !excess density / dyn. visc. & * dustd2 * 1.e-4)*dtb wdust = dustsink +#ifdef SED_OFFLINE + wdusto= dustsink * dto/dtb +#endif IF (mnproc.eq.1) THEN WRITE(io_stdo_bgc,*) &
--- a/mo_biomod.F90 Tue Jul 02 11:02:37 2019 +0200 +++ b/mo_biomod.F90 Fri Jul 26 16:03:00 2019 +0200 @@ -49,6 +49,9 @@ REAL :: bluefix,tf2,tf1,tf0,tff REAL :: bkphy,bkzoo,bkopal,bifr13,bifr14 REAL :: wpoc,wcal,wopal,drempoc,dremdoc,dremn2o +# if defined(SED_OFFLINE) + REAL :: wpoco,wcalo,wopalo,wdusto +# endif REAL :: dphymor,dzoomor,dremopal REAL :: dremsul REAL :: psedi,csedi,ssedi @@ -63,6 +66,9 @@ REAL :: dustd1,dustd2,dustd3,dustsink,calmax #elif defined(WLIN) REAL :: wmin,wmax,wlin +# if defined(SED_OFFLINE) + REAL :: wmino,wmaxo,wlino +# endif #endif #ifdef __c_isotopes REAL :: factor_13c, factor_14c, atm_c14_cal, atm_dc14_cal
--- a/mo_sedmnt_offline.F90 Tue Jul 02 11:02:37 2019 +0200 +++ b/mo_sedmnt_offline.F90 Fri Jul 26 16:03:00 2019 +0200 @@ -75,8 +75,8 @@ use mo_carbch, only: ocetra, keqb, co3 use mo_param1_bgc, only: nocetra, idet, icalc, iopal, ifdust & & , ks, nsedtra, npowtra -use mo_biomod, only: kbo, bolay, wpoc, wmin, wmax, wdust & - & , wlin, wopal, wcal +use mo_biomod, only: kbo, bolay, wpoco, wmino, wmaxo, wdusto & + & , wlino, wopalo, wcalo use mod_xc use mo_common_bgc use mod_dia, only: ddm, depthslev, depthslev_bnds & @@ -291,11 +291,11 @@ co3_kbo_avg = co3_kbo_avg / nstep_in_month do j = 1, jdm ! Here we assume the WLIN case. do i = 1, idm - wpoc = min(wmin+wlin*bolay_avg(i,j),wmax) - prorca_avg(i,j)=ocetra_kbo_avg(i,j,idet )*wpoc - prcaca_avg(i,j)=ocetra_kbo_avg(i,j,icalc )*wcal - silpro_avg(i,j)=ocetra_kbo_avg(i,j,iopal )*wopal - produs_avg(i,j)=ocetra_kbo_avg(i,j,ifdust)*wdust + wpoco = min(wmino+wlino*bolay_avg(i,j),wmaxo) + prorca_avg(i,j)=ocetra_kbo_avg(i,j,idet )*wpoco + prcaca_avg(i,j)=ocetra_kbo_avg(i,j,icalc )*wcalo + silpro_avg(i,j)=ocetra_kbo_avg(i,j,iopal )*wopalo + produs_avg(i,j)=ocetra_kbo_avg(i,j,ifdust)*wdusto enddo enddo nstep_in_month = 0 @@ -391,6 +391,9 @@ integer :: nday_save, nmonth_save ! calendar variables outside sedmnt_offline() integer :: nyear_save, nday_of_year_save, nd_in_m_2, nday_in_year_save integer :: nacc_bgc_1, nacc_bgc_2 +#if defined(SED_OFFLINE_DEBUG) + character(len= 2) seqstring_burst, seqstring_year +#endif ! Read the bottom water climatology ! @@ -399,8 +402,11 @@ lread_clim = .false. ! but only one time. endif - if ( lsed_spinup .and. lcompleted_clim ) then - if (.not. (nmonth==1 .and. nday==1) .and. mnproc==1) then + ! Apparently we need to wait until nday = 2 before starting the off-line + ! spin-up, because otherwise data won't be saved for the last period. + ! + if ( lsed_spinup .and. lcompleted_clim .and. nday==2 ) then + if (.not. (nmonth==1 .and. nday<=2) .and. mnproc==1) then write(io_stdo_bgc,*) 'sedmnt_offline(): WARNING: Not at start of year!' endif @@ -435,6 +441,39 @@ & 'sedmnt_offline(): nyear_global = ', nyear_global do nmonth = 1, 12 +#if defined(SED_OFFLINE_DEBUG) + ! We misuse the *_avg array variables to write current bottom water tracers. + ! An array of 12 timeslices is reserved, but only the first one stored (FIXME). + ! + do iocetra = 1, nocetra + do j = 1, jdm + do i = 1, idm + ocetra_kbo_avg(i,j,iocetra) = ocetra(i,j,kbo(i,j),iocetra) + enddo + enddo + enddo + do j = 1, jdm + do i = 1, idm + bgc_t_kbo_avg(i,j) = bgc_t(i,j,kbo(i,j)) + bgc_s_kbo_avg(i,j) = bgc_s(i,j,kbo(i,j)) + bgc_rho_kbo_avg(i,j) = bgc_rho(i,j,kbo(i,j)) + co3_kbo_avg(i,j) = co3(i,j,kbo(i,j)) + wpoco = min(wmino+wlino*bolay_avg(i,j),wmaxo) + prorca_avg(i,j)=ocetra_kbo_avg(i,j,idet )*wpoco + prcaca_avg(i,j)=ocetra_kbo_avg(i,j,icalc )*wcalo + silpro_avg(i,j)=ocetra_kbo_avg(i,j,iopal )*wopalo + produs_avg(i,j)=ocetra_kbo_avg(i,j,ifdust)*wdusto + enddo + enddo + bolay_avg = bolay + keqb_avg = keqb + + write (seqstring_burst,'(I0.2)') nburst + write (seqstring_year,'(I0.2)') nyear + call aufw_bgc_onlysed(idm, jdm, kdm, nyear, nmonth, nday, nstep, & + & "bottom_seawater.burst"//seqstring_burst//".year"//seqstring_year//".DEBUG.nc") +#endif + ! set up sediment layers (much higher diffusion rate; correct timestep) dtoff = 3600*24*nd_in_m(nmonth) call bodensed(kpie,kpje,kpke,pddpo) @@ -503,8 +542,24 @@ powtra2(:,:,ks+1:2*ks,:) = powtra(:,:,:,:) burial2(:,:,1,:) = burial(:,:,:) burial2(:,:,2,:) = burial(:,:,:) + +#if defined(SED_OFFLINE_DEBUG) + ! We misused the *_avg array variables to write current bottom water tracers. + ! + ocetra_kbo_avg = 0.0 + bolay_avg = 4000.0 + keqb_avg = 0.0 + prorca_avg = 0.0 + prcaca_avg = 0.0 + silpro_avg = 0.0 + produs_avg = 0.0 + bgc_t_kbo_avg = 0.0 + bgc_s_kbo_avg = 0.0 + bgc_rho_kbo_avg = 0.0 + co3_kbo_avg = 0.0 +#endif + lcompleted_clim = .false. ! don't reuse the climatology endif - lcompleted_clim = .false. ! don't reuse the climatology end subroutine sedmnt_offline subroutine alloc_mem_sedmnt_offline(kpie, kpje)
--- a/mod_dia.F Tue Jul 02 11:02:37 2019 +0200 +++ b/mod_dia.F Fri Jul 26 16:03:00 2019 +0200 @@ -228,6 +228,7 @@ integer runid_len,nstep real diagfq,filefq logical filemon,fileann,filedec,filecen,filemil + integer mymonth c real epsil parameter (epsil=1.e-11) @@ -265,10 +266,18 @@ . runid(1:runid_len),'_',ctag,'_',ny,'.',nm,'.',nd,'.nc' endif elseif (expcnf.eq.'cesm') then - if ( filemon .or. fileann ) then - call pstdat(ny,nm,nd,1) - else - call pstdat(ny,nm,nd,nint(diagfq)) + if (.not. lspinning_up_sed ) then + if (diagfq==30) then ! always rewind the actual diagfq as the + ! fname is created after the first diagfq + if (nmonth /= 1) then + mymonth = nmonth-1 + else + mymonth = 12 + endif + call pstdat(ny,nm,nd,nd_in_m(mymonth)) + else + call pstdat(ny,nm,nd,nint(diagfq)) + endif endif if (filemon) then if ( lspinning_up_sed ) then
--- a/powach.F90 Tue Jul 02 11:02:37 2019 +0200 +++ b/powach.F90 Fri Jul 26 16:03:00 2019 +0200 @@ -47,10 +47,10 @@ ! REAL bolay - thickness of the bottom gridbox [m]. ! REAL ocetra - bottom layer ocean tracer ocetra(:,:,kbo,:). ! REAL keqb - chemical equilibrium constants. -! REAL prorca - sedimentation of carbon [kmol/m^2/s]. -! REAL prcaca - sedimentation of calcium carbonate [kmol/m^2/s]. -! REAL silpro - sedimentation of biogenic silica [kmol/m^2/s]. -! REAL produs - sedimentation of lithogenic dust [kmol/m^2/s]. +! REAL prorca - sedimentation of carbon [kmol/m^2/timestep]. +! REAL prcaca - sedimentation of calcium carbonate [kmol/m^2/timestep]. +! REAL silpro - sedimentation of biogenic silica [kmol/m^2/timestep]. +! REAL produs - sedimentation of lithogenic dust [kmol/m^2/timestep]. ! REAL co3 - dissolved carbonate in the bottom gridbox [mol/kg]. ! ! The actual argument names sometimes end with an underscore, signifying @@ -86,10 +86,10 @@ real, intent(in) :: bolay_ (kpie,kpje) real, intent(in) :: ocetra_(kpie,kpje,nocetra) real, intent(in) :: keqb_ (11,kpie,kpje) -real, intent(inout) :: prorca_(kpie,kpje) -real, intent(inout) :: prcaca_(kpie,kpje) -real, intent(inout) :: silpro_(kpie,kpje) -real, intent(inout) :: produs_(kpie,kpje) +real, intent(in) :: prorca_(kpie,kpje) +real, intent(in) :: prcaca_(kpie,kpje) +real, intent(in) :: silpro_(kpie,kpje) +real, intent(in) :: produs_(kpie,kpje) real, intent(in) :: co3_ (kpie,kpje) real, intent(in) :: omask (kpie,kpje) @@ -201,10 +201,10 @@ do i = 1, kpie if(omask(i,j) > 0.5) then - if ( .not. lspinning_up_sed ) then sedfluxo(i,j,ipowasi) = sedfluxo(i,j,ipowasi) + & & (silsat-sediso(i,0) - ocetra_(i,j,isilica))*bolay_(i,j) + if ( .not. lspinning_up_sed ) then ! NOTE: If ocetra_(:,:,:) instead of ocetra(:,:,kbo,:) were to be updated, ! the latter would need to be updated by the former later on (MvH)! ! @@ -562,16 +562,16 @@ !$OMP PARALLEL DO do j = 1, kpje do i = 1, kpie - silpro_(i,j) = 0. - prorca_(i,j) = 0. + silpro(i,j) = 0. + prorca(i,j) = 0. #ifdef __c_isotopes pror13(i,j) = 0. pror14(i,j) = 0. prca13(i,j) = 0. prca14(i,j) = 0. #endif - prcaca_(i,j) = 0. - produs_(i,j) = 0. + prcaca(i,j) = 0. + produs(i,j) = 0. enddo enddo !$OMP END PARALLEL DO
--- a/sediment_step.F90 Tue Jul 02 11:02:37 2019 +0200 +++ b/sediment_step.F90 Fri Jul 26 16:03:00 2019 +0200 @@ -71,10 +71,10 @@ real, intent(in) :: ocetra_(kpie,kpje,nocetra) real, intent(in) :: bolay_ (kpie,kpje) real, intent(in) :: keqb_ (11,kpie,kpje) -real, intent(inout) :: prorca_(kpie,kpje) -real, intent(inout) :: prcaca_(kpie,kpje) -real, intent(inout) :: silpro_(kpie,kpje) -real, intent(inout) :: produs_(kpie,kpje) +real, intent(in) :: prorca_(kpie,kpje) +real, intent(in) :: prcaca_(kpie,kpje) +real, intent(in) :: silpro_(kpie,kpje) +real, intent(in) :: produs_(kpie,kpje) real, intent(in) :: co3_ (kpie,kpje) integer, intent(in), optional :: jpowaic_(nbgcmax) integer, intent(in), optional :: jpowaal_(nbgcmax)