diff --git a/config/namelist.io b/config/namelist.io index 4a83228e9..abeef5b92 100644 --- a/config/namelist.io +++ b/config/namelist.io @@ -56,4 +56,5 @@ io_list = 'sst ',1, 'm', 4, 'bolus_w ',1, 'y', 4, 'fw ',1, 'm', 4, 'fh ',1, 'm', 4, + 'otracers ',1, 'y', 4, / diff --git a/config/namelist.transit b/config/namelist.transit index 3ca5c87cf..5510ef1ad 100644 --- a/config/namelist.transit +++ b/config/namelist.transit @@ -1,11 +1,16 @@ ! The namelist file for transient tracers &transit_param +l_r14c = .false. +l_r39ar = .false. +l_f11 = .false. +l_f12 = .false. +l_sf6 = .false. anthro_transit=.false. paleo_transit=.false. length_transit=1 ! 166 for anthro_transit=.true. ti_start_transit=1 ! 1 for D14C, 80 for CFC-12 -ifile_transit='/work/ollie/mbutzin/fesom2/input/trace_gases/Table_CO2_isoC_CFC12_SF6.txt' +ifile_transit='/work/ab0246/a270108/fesom2_recom_config/input-for-awiesm/Table_CO2_isoC_CFCs1112_SF6.txt' r14c_a = 1.0000 ! atm. 14C/C ratio, global mean r39ar_a = 1.0000 ! atm. 39Ar/Ar ratio, global mean xarg_a = 9.34e-3 ! atm. Argon concn. (mole fraction), global mean diff --git a/src/fesom_module.F90 b/src/fesom_module.F90 index dffef9186..8f0113eb2 100755 --- a/src/fesom_module.F90 +++ b/src/fesom_module.F90 @@ -53,11 +53,14 @@ module fesom_main_storage_module use recom_interface #endif +! Transient tracers +use mod_transit, only: year_ce, r14c_nh, r14c_tz, r14c_sh, r14c_ti, xCO2_ti, xf11_nh, xf11_sh, xf12_nh, xf12_sh, xsf6_nh, xsf6_sh, ti_transit, anthro_transit + implicit none type :: fesom_main_storage_type - integer :: n, from_nstep, offset, row, i, provided + integer :: n, from_nstep, offset, row, i, provided, id integer :: which_readr ! read which restart files (0=netcdf, 1=core dump,2=dtype) integer :: total_nsteps integer, pointer :: mype, npes, MPIerr, MPI_COMM_FESOM, MPI_COMM_WORLD, MPI_COMM_FESOM_IB @@ -186,6 +189,21 @@ subroutine fesom_init(fesom_total_nsteps) if (f%mype==0) write(*,*) 'FESOM mesh_setup... complete' +! Transient tracers: control output of initial input values + if(use_transit .and. anthro_transit .and. f%mype==0) then + write (*,*) + write (*,*) "*** Transient tracers: Initial atmospheric input values >>>" + write (*,*) "Year CE, xCO2, D14C_NH, D14C_TZ, D14C_SH, xCFC-11_NH, xCFC-11_SH, xCFC-12_NH, xCFC-12_SH, xSF6_NH, xSF6_SH" + write (*, fmt="(2x,i4,10(2x,f6.2))"), & + year_ce(ti_transit), xCO2_ti(ti_transit) * 1.e6, & + (r14c_nh(ti_transit) - 1.) * 1000., (r14c_tz(ti_transit) - 1.) * 1000., (r14c_sh(ti_transit) - 1.) * 1000., & + xf11_nh(ti_transit) * 1.e12, xf11_sh(ti_transit) * 1.e12, & + xf12_nh(ti_transit) * 1.e12, xf12_sh(ti_transit) * 1.e12, & + xsf6_nh(ti_transit) * 1.e12, xsf6_sh(ti_transit) * 1.e12 + write (*,*) + end if + + !===================== ! Allocate field variables ! and additional arrays needed for @@ -591,6 +609,14 @@ subroutine fesom_runloop(current_nsteps) !___model ocean step____________________________________________________ if (flag_debug .and. f%mype==0) print *, achar(27)//'[34m'//' --> call oce_timestep_ale'//achar(27)//'[0m' call oce_timestep_ale(n, f%ice, f%dynamics, f%tracers, f%partit, f%mesh) + if (use_transit) then + ! Prevent negative concentrations of SF6, CFC-11 and CFC-12 during the first years (inital values are zero) + do tr_num=1, f%tracers%num_tracers + if ((f%tracers%data(tr_num)%ID==6) .or. (f%tracers%data(tr_num)%ID==11) .or. (f%tracers%data(tr_num)%ID==12)) then + f%tracers%data(tr_num)%values(:,:) = max(f%tracers%data(tr_num)%values(:,:), 0._WP) + end if + end do + end if ! use_transit f%t3 = MPI_Wtime() !___compute energy diagnostics..._______________________________________ @@ -620,6 +646,24 @@ subroutine fesom_runloop(current_nsteps) #if defined (__recom) f%rtime_compute_recom = f%rtime_compute_recom + f%t1_recom - f%t0_recom #endif + +! Transient tracers: update of input values between restarts + if(use_transit .and. anthro_transit .and. (daynew == ndpyr) .and. (timenew==86400.)) then + ti_transit = ti_transit + 1 + if (f%mype==0) then + write (*,*) + write (*,*) "*** Transient tracers: Updated atmospheric input values >>>" + write (*,*) "Year CE, xCO2, D14C_NH, D14C_TZ, D14C_SH, xCFC-11_NH, xCFC-11_SH, xCFC-12_NH, xCFC-12_SH, xSF6_NH, xSF6_SH" + write (*, fmt="(2x,i4,10(2x,f6.2))"), & + year_ce(ti_transit), xCO2_ti(ti_transit) * 1.e6, & + (r14c_nh(ti_transit) - 1.) * 1000., (r14c_tz(ti_transit) - 1.) * 1000., (r14c_sh(ti_transit) - 1.) * 1000., & + xf11_nh(ti_transit) * 1.e12, xf11_sh(ti_transit) * 1.e12, & + xf12_nh(ti_transit) * 1.e12, xf12_sh(ti_transit) * 1.e12, & + xsf6_nh(ti_transit) * 1.e12, xsf6_sh(ti_transit) * 1.e12 + write (*,*) + end if + endif + end do !call cray_acc_set_debug_global_level(3) f%from_nstep = f%from_nstep+current_nsteps diff --git a/src/gen_model_setup.F90 b/src/gen_model_setup.F90 index e469bdf32..b9ddd34ea 100755 --- a/src/gen_model_setup.F90 +++ b/src/gen_model_setup.F90 @@ -253,18 +253,24 @@ subroutine setup_model(partit) #endif if (use_transit) then -! Transient tracer input, input file names have to be specified in -! namelist.config, nml=run_config - if(partit%mype==0) print *, "Transient tracers are ON. Tracer input file: ", ifile_transit - open (20,file=ifile_transit) - if (anthro_transit .or. paleo_transit) then - call read_transit_input + if(partit%mype==0) print *, "Transient tracers are ON." + nmlfile = 'namelist.transit' ! name of transient tracers namelist file + open (newunit=fileunit, file=nmlfile) + read (fileunit, nml=transit_param) + close (fileunit) + if (anthro_transit) then +! transient values of historical CO2, bomb radiocarbon, CFC-12, and SF6 + if(partit%mype==0) print *, "Reading transient input values from file: ", ifile_transit + open (newunit=fileunit, file=ifile_transit) + call read_transit_input(fileunit) + close (fileunit) + elseif (paleo_transit) then +! transient values of atmospheric CO2 and radiocarbon reconstructions +! under construction / not yet realized else -! Spinup / equilibrium runs with constant tracer input, -! read parameter values from namelist.oce - read (20,nml=transit_param) +! Spinup / equilibrium runs with constant tracer input specified in namelist.transit + if(partit%mype==0) print *, "Reading constant input values from file: ", nmlfile end if - close (20) end if if(partit%mype==0) write(*,*) 'Namelist files are read in' diff --git a/src/io_meandata.F90 b/src/io_meandata.F90 index 82448d28a..a70a84dab 100644 --- a/src/io_meandata.F90 +++ b/src/io_meandata.F90 @@ -152,7 +152,7 @@ subroutine ini_mean_io(ice, dynamics, tracers, partit, mesh) #endif use g_forcing_param, only: use_virt_salt, use_landice_water, use_age_tracer !---fwf-code, age-code use g_config, only : lwiso !---wiso-code - use mod_transit, only : index_transit + use mod_transit, only : index_transit_r14c, index_transit_r39ar, index_transit_f11, index_transit_f12, index_transit_sf6 implicit none integer :: i, j @@ -264,9 +264,9 @@ subroutine ini_mean_io(ice, dynamics, tracers, partit, mesh) end if CASE ('m_snow ') if (use_ice) then - call def_stream(nod2D, myDim_nod2D, 'm_snow', 'snow height per unit area', 'm', ice%data(3)%values(1:myDim_nod2D), io_list(i)%freq, io_list(i)%unit, io_list(i)%precision, partit, mesh) + call def_stream(nod2D, myDim_nod2D, 'm_snow', 'snow height per unit area', 'm', ice%data(3)%values(1:myDim_nod2D), io_list(i)%freq, io_list(i)%unit, io_list(i)%precision, partit, mesh) end if -CASE ('h_ice ') +CASE ('h_ice ') if (use_ice) then call def_stream(nod2D, myDim_nod2D, 'h_ice', 'ice thickness over ice-covered fraction', 'm', ice%h_ice(1:myDim_nod2D), io_list(i)%freq, io_list(i)%unit, io_list(i)%precision, partit, mesh) end if @@ -274,7 +274,6 @@ subroutine ini_mean_io(ice, dynamics, tracers, partit, mesh) if (use_ice) then call def_stream(nod2D, myDim_nod2D, 'h_snow', 'snow thickness over ice-covered fraction', 'm', ice%h_snow(1:myDim_nod2D), io_list(i)%freq, io_list(i)%unit, io_list(i)%precision, partit, mesh) end if - ! Debug ice variables CASE ('strength_ice') if (use_ice) then @@ -880,13 +879,33 @@ subroutine ini_mean_io(ice, dynamics, tracers, partit, mesh) endif else #endif - call def_stream((/nl-1, nod2D/), (/nl-1, myDim_nod2D/), 'tra_'//id_string, 'passive tracer ID='//id_string, 'n/a', tracers%data(j)%values(:,:), io_list(i)%freq, io_list(i)%unit, io_list(i)%precision, partit, mesh) + + if (use_transit) then +! Transient tracers + select case (tracers%data(j)%ID) + case(6) + call def_stream((/nl-1, nod2D/), (/nl-1, myDim_nod2D/), 'sf6', 'sulfur hexafluoride', 'mol / m**3', tracers%data(index_transit_sf6)%values(:,:), io_list(i)%freq, io_list(i)%unit, io_list(i)%precision, partit, mesh) + case(11) + call def_stream((/nl-1, nod2D/), (/nl-1, myDim_nod2D/), 'cfc11', 'chlorofluorocarbon CFC-11', 'mol / m**3', tracers%data(index_transit_f11)%values(:,:), io_list(i)%freq, io_list(i)%unit, io_list(i)%precision, partit, mesh) + case(12) + call def_stream((/nl-1, nod2D/), (/nl-1, myDim_nod2D/), 'cfc12', 'chlorofluorocarbon CFC-12', 'mol / m**3', tracers%data(index_transit_f12)%values(:,:), io_list(i)%freq, io_list(i)%unit, io_list(i)%precision, partit, mesh) + case(14) + call def_stream((/nl-1, nod2D/), (/nl-1, myDim_nod2D/), 'r14c', '14C / C ratio of DIC', 'none', tracers%data(index_transit_r14c)%values(:,:), io_list(i)%freq, io_list(i)%unit, io_list(i)%precision, partit, mesh) + case(39) + call def_stream((/nl-1, nod2D/), (/nl-1, myDim_nod2D/), 'r39ar', '39Ar / Ar ratio', 'none', tracers%data(index_transit_r39ar)%values(:,:), io_list(i)%freq, io_list(i)%unit, io_list(i)%precision, partit, mesh) + case default + call def_stream((/nl-1, nod2D/), (/nl-1, myDim_nod2D/), 'tra_'//id_string, 'passive tracer ID='//id_string, 'n/a', tracers%data(j)%values(:,:), io_list(i)%freq, io_list(i)%unit, io_list(i)%precision, partit, mesh) + end select + else +! Other passive but not transient tracers + call def_stream((/nl-1, nod2D/), (/nl-1, myDim_nod2D/), 'tra_'//id_string, 'passive tracer ID='//id_string, 'n/a', tracers%data(j)%values(:,:), io_list(i)%freq, io_list(i)%unit, io_list(i)%precision, partit, mesh) + end if #if defined(__recom) end if #endif - end do + end do CASE ('sigma0 ') call def_stream((/nl-1, nod2D/), (/nl-1, myDim_nod2D/), 'sigma0', 'potential density', 'kg/m3', density_m_rho0(:,:), io_list(i)%freq, io_list(i)%unit, io_list(i)%precision, partit, mesh) @@ -922,17 +941,6 @@ subroutine ini_mean_io(ice, dynamics, tracers, partit, mesh) CASE ('slope_z ') call def_stream((/nl-1, nod2D/), (/nl-1, myDim_nod2D/), 'slope_z', 'neutral slope Z', 'none', neutral_slope(3,:,:), io_list(i)%freq, io_list(i)%unit, io_list(i)%precision, partit, mesh) -! Transient tracers -CASE('SF6 ') - if (use_transit) call def_stream((/nl-1, nod2D/), (/nl-1, myDim_nod2D/), 'sf6', 'sulfur hexafluoride', 'mol / m**3', tracers%data(index_transit(1))%values(:,:), io_list(i)%freq, io_list(i)%unit, io_list(i)%precision, partit, mesh) -CASE('CFC-12 ') - if (use_transit) call def_stream((/nl-1, nod2D/), (/nl-1, myDim_nod2D/), 'cfc12', 'chlorofluorocarbon CFC-12', 'mol / m**3', tracers%data(index_transit(2))%values(:,:), io_list(i)%freq, io_list(i)%unit, io_list(i)%precision, partit, mesh) -CASE('R14C ') - if (use_transit) call def_stream((/nl-1, nod2D/), (/nl-1, myDim_nod2D/), 'r14c', '14C / C ratio of DIC', 'none', tracers%data(index_transit(3))%values(:,:), io_list(i)%freq, io_list(i)%unit, io_list(i)%precision, partit, mesh) -CASE('R39Ar ') - if (use_transit) call def_stream((/nl-1, nod2D/), (/nl-1, myDim_nod2D/), 'r39ar', '39Ar / Ar ratio', 'none', tracers%data(index_transit(4))%values(:,:), io_list(i)%freq, io_list(i)%unit, io_list(i)%precision, partit, mesh) -! Transient tracers end - CASE ('N2 ') call def_stream((/nl, nod2D/), (/nl, myDim_nod2D/), 'N2', 'brunt väisälä', '1/s2', bvfreq(:,:), io_list(i)%freq, io_list(i)%unit, io_list(i)%precision, partit, mesh) CASE ('Kv ') diff --git a/src/io_restart.F90 b/src/io_restart.F90 index a51be4bd6..c23b55843 100644 --- a/src/io_restart.F90 +++ b/src/io_restart.F90 @@ -62,7 +62,7 @@ subroutine ini_ocean_io(year, dynamics, tracers, partit, mesh) use nvfortran_subarray_workaround_module #endif integer, intent(in) :: year - integer :: j + integer :: j, id character(500) :: longname character(500) :: trname, units character(4) :: cyear @@ -118,7 +118,8 @@ subroutine ini_ocean_io(year, dynamics, tracers, partit, mesh) endif do j=1,tracers%num_tracers - SELECT CASE (j) + id=tracers%data(j)%ID !MB: Avoid hard-wired tracer assignments like SELECT CASE(j) + SELECT CASE (id) CASE(1) trname='temp' longname='potential temperature' @@ -131,6 +132,10 @@ subroutine ini_ocean_io(year, dynamics, tracers, partit, mesh) trname='sf6' longname='sulfur hexafluoride' units='mol / m**3' + CASE(11) + trname='cfc11' + longname='chlorofluorocarbon CFC-11' + units='mol / m**3' CASE(12) trname='cfc12' longname='chlorofluorocarbon CFC-12' diff --git a/src/mod_transit.F90 b/src/mod_transit.F90 index 6713032eb..7fbc8e018 100644 --- a/src/mod_transit.F90 +++ b/src/mod_transit.F90 @@ -16,6 +16,7 @@ MODULE mod_transit r39ar_a = 1.0, & ! 39Ar / 40 Ar (homogeneous) xarg_a = 9.34e-3, & ! Argon (homogeneous) xCO2_a = 284.32e-6, & ! CO2 (CMIP6 & OMIP-BGC: 284.32e-6 for 1700-1850, PMIP4: 190.00e-6 for 21 ka BP) + xf11_a = 0.0, & ! CFC-11 (latitude dependent) xf12_a = 0.0, & ! CFC-12 (latitude dependent) xsf6_a = 0.0 ! SF6 (latitude dependent) @@ -23,20 +24,27 @@ MODULE mod_transit real(kind=8), allocatable, dimension(:) :: r14c_nh, r14c_tz, r14c_sh, & ! 14CO2 / 12CO2, latitude-dependent (e.g., bomb 14C) r14c_ti, & ! 14CO2 / 12CO2, homogenous (e.g., IntCal) xCO2_ti, & ! CO2 + xf11_nh, xf11_sh, & ! CFC-11, latitude-dependent xf12_nh, xf12_sh, & ! CFC-12, latitude-dependent xsf6_nh, xsf6_sh ! SF6, latitude-dependent integer, allocatable, dimension(:) :: year_ce ! current year in anthropenic runs (control output) integer :: length_transit = 1, & ! length (years) of transient tracer input ti_start_transit = 1 ! index of the first tracer input year in ifile_transit - logical :: anthro_transit = .false., & + logical :: l_r14c = .false., & ! switch on / off specific tracerss + l_r39ar = .false. , & + l_f11 = .false., & + l_f12 = .false., & + l_sf6 = .false., & + anthro_transit = .false., & paleo_transit = .false. ! specify tracer input scenario - character(300) :: ifile_transit ='Table_CO2_isoC_CFC12_SF6.txt'! tracer input file; not neccessary for steady state simulations + character(300) :: ifile_transit ='Table_CO2_isoC_CFCs_SF6.txt'! tracer input file; not neccessary for steady state simulations ! Parameters which can be changed via namelist.oce (-> transit_param) ! Global-mean concentrations of DIC and Argon in the mixed layer (mol / m**3) real(kind=8) :: dic_0 = 2.00, & ! GLODAPv2, 0-50 m: TCO2 ~ 2050 umol / kg arg_0 = 0.01 ! Hamme et al. 2019, doi:10.1146/annurev-marine-121916-063604 + ! Radioactive decay constants (1 / s; default values assume that 1 year = 365.00 days) real(kind=8) :: decay14 = 3.8561e-12 , & ! 14C; t1/2 = 5700 a following OMIP-BGC decay39 = 8.1708e-11 ! 39Ar; t1/2 = 269 a @@ -45,16 +53,20 @@ MODULE mod_transit ! Latitude of atmospheric boundary conditions and latitudinal interpolation weight real(kind=8) :: y_abc, yy_nh ! Tracer indices of transient tracers - integer :: id_r14c, id_r39ar, id_f12, id_sf6 - integer, dimension(4) :: index_transit = (/-1, -1, -1, -1/) + integer :: id_r14c, id_r39ar, id_f11, id_f12, id_sf6, index_transit_r14c, index_transit_r39ar, index_transit_f11, index_transit_f12, index_transit_sf6 ! Time index (=year) in transient simulations integer :: ti_transit ! Namelist to modify default parameter settings namelist / transit_param / & - anthro_transit, & ! anthoropogenic transient tracers + l_r14c, & ! switch on R14C + l_r39ar, & ! switch on R39Ar + l_f11, & ! switch on CFC-11 + l_f12, & ! switch on CFC-12 + l_sf6, & ! switch on SF6 + anthro_transit, & ! anthropogenic transient tracers paleo_transit, & ! paleo transient tracers - length_transit, & ! 166 for anthro_transit=.true. + length_transit, & ! 166 if (anthro_transit=.true.) ti_start_transit, & ! 1 for D14C, 80 for CFC-12 ifile_transit, & ! forcing file r14c_a, & ! atm. 14C/C ratio, global mean @@ -135,6 +147,11 @@ function solub(which_gas, temp_c, sal) a1 = -160.7333; a2 = 215.4152; a3 = 89.8920; a4 = -1.47759; pow = 2 b1 = 0.029941; b2 = -0.027455; b3 = 0.0053407; c1 = 0. con2con = 1000. ! convert to mol / (m**3 * atm) + case ("f11") +! CFC-11 in mol / (L * atm) (Warner & Weiss 1985, doi:10.1016/0198-0149(85)90099-8, Table 5) + a1 = -229.9261; a2 = 319.6552; a3 = 119.4471; a4 = -1.39165; pow = 2 + b1 = -0.142382; b2 = 0.091459; b3 = -0.0157274; c1 = 0. + con2con = 1000. ! convert to mol / (m**3 * atm) case ("f12") ! CFC-12 in mol / (L * atm) (Warner & Weiss 1985, doi:10.1016/0198-0149(85)90099-8, Table 5) a1 = -218.0971; a2 = 298.9702; a3 = 113.8049; a4 = -1.39165; pow = 2 @@ -175,6 +192,8 @@ function sc_660(which_gas, temp_c) select case (which_gas) case ("co2") ! CO2 as = 2116.8; bs = -136.25; cs = 4.7353; ds = -0.092307; es = 0.0007555 + case ("f11") ! CFC-11 + as = 3579.2; bs = -222.63; cs = 7.5749; ds = -0.145950; es = 0.0011870 case ("f12") ! CFC-12 as = 3828.1; bs = -249.86; cs = 8.7603; ds = -0.171600; es = 0.0014080 case ("sf6") ! SF6 @@ -246,11 +265,12 @@ function speed_2(windstr_x, windstr_y) end function speed_2 - subroutine read_transit_input + subroutine read_transit_input(ifileunit) ! Read atmospheric input of isoCO2 and / or other tracers implicit none ! Internal variables + integer, intent(in) :: ifileunit integer :: jj real(kind=8), allocatable, dimension(:) :: d14c_nh, d14c_tz, d14c_sh, d14c_ti, d13c_dummy @@ -263,6 +283,8 @@ subroutine read_transit_input allocate(r14c_tz(length_transit)) allocate(r14c_sh(length_transit)) allocate(xCO2_ti(length_transit)) + allocate(xf11_nh(length_transit)) + allocate(xf11_sh(length_transit)) allocate(xf12_nh(length_transit)) allocate(xf12_sh(length_transit)) allocate(xsf6_nh(length_transit)) @@ -272,14 +294,15 @@ subroutine read_transit_input ! Skip header lines do jj = 1,15 - read (20, fmt=*) + read (ifileunit, fmt=*) end do ! Read input values do jj = 1, length_transit - read (20, fmt=*) year_ce(jj), & + read (ifileunit, fmt=*) year_ce(jj), & xCO2_ti(jj), & d14c_nh(jj), d14c_tz(jj), d14c_sh(jj), & d13c_dummy(jj), & + xf11_nh(jj), xf11_sh(jj), & xf12_nh(jj), xf12_sh(jj), & xsf6_nh(jj), xsf6_sh(jj) end do @@ -290,6 +313,8 @@ subroutine read_transit_input r14c_sh = 1. + 0.001 * d14c_sh ! Convert volume mixing ratios xCO2_ti = xCO2_ti * 1.e-6 + xf11_nh = xf11_nh * 1.e-12 + xf11_sh = xf11_sh * 1.e-12 xf12_nh = xf12_nh * 1.e-12 xf12_sh = xf12_sh * 1.e-12 xsf6_nh = xsf6_nh * 1.e-12 diff --git a/src/oce_ale_tracer.F90 b/src/oce_ale_tracer.F90 index bb40ff335..0a07f26b2 100644 --- a/src/oce_ale_tracer.F90 +++ b/src/oce_ale_tracer.F90 @@ -81,29 +81,15 @@ subroutine diff_tracers_ale(tr_num, dynamics, tracer, ice, partit, mesh) module bc_surface_interface interface - function bc_surface(n, id, sval, nzmin, partit) - use mod_mesh - USE MOD_PARTIT - USE MOD_PARSUP - integer , intent(in) :: n, id, nzmin - type(t_partit), intent(inout), target :: partit - real(kind=WP) :: bc_surface - real(kind=WP), intent(in) :: sval - end function - end interface -end module - -module transit_bc_surface_interface - interface - function transit_bc_surface(n, id, sst, sss, a_ice, sval, nzmin, partit, mesh) + function bc_surface(n, id, sval, nzmin, partit, mesh, sst, sss, a_ice) use mod_mesh USE MOD_PARTIT USE MOD_PARSUP integer , intent(in) :: n, id, nzmin type(t_partit), intent(inout), target :: partit type(t_mesh), intent(in), target :: mesh - real(kind=WP) :: transit_bc_surface - real(kind=WP), intent(in) :: sst, sss, a_ice, sval + real(kind=WP) :: bc_surface + real(kind=WP), intent(in) :: sval, sst, sss, a_ice end function end interface end module @@ -571,7 +557,6 @@ subroutine diff_ver_part_impl_ale(tr_num, dynamics, tracers, ice, partit, mesh) use o_mixing_KPP_mod !for ghats _GO_ use g_cvmix_kpp, only: kpp_nonlcltranspT, kpp_nonlcltranspS, kpp_oblmixc use bc_surface_interface - use transit_bc_surface_interface use mod_ice use iceberg_params implicit none @@ -992,11 +977,10 @@ subroutine diff_ver_part_impl_ale(tr_num, dynamics, tracers, ice, partit, mesh) !_______________________________________________________________________ ! case of activated shortwave penetration into the ocean, ad 3d contribution if (use_icebergs .and. (.not. turn_off_hf) .and. tracers%data(tr_num)%ID==1) then - do nz=nzmin, nzmax - + do nz=nzmin, nzmax-1 zinv=1.0_WP*dt !/(zbar(nz)-zbar(nz+1)) ale! !!PS tr(nz)=tr(nz)+(sw_3d(nz, n)-sw_3d(nz+1, n) * ( area(nz+1,n)/areasvol(nz,n)) ) * zinv - tr(nz)=tr(nz) + ibhf_n(nz, n) * zinv / vcpw + tr(nz)=tr(nz)+(ibhf_n(nz, n)-ibhf_n(nz+1, n) * area(nz+1,n)/areasvol(nz,n)) * zinv / vcpw end do end if @@ -1012,12 +996,8 @@ subroutine diff_ver_part_impl_ale(tr_num, dynamics, tracers, ice, partit, mesh) ! (BUT CHECK!) | | | | ! v (+) v (+) ! - tr(nzmin)= tr(nzmin)+bc_surface(n, tracers%data(tr_num)%ID, trarr(nzmin,n), nzmin, partit) + tr(nzmin)= tr(nzmin)+bc_surface(n, tracers%data(tr_num)%ID, trarr(nzmin,n), nzmin, partit, mesh, sst(nzmin,n), sss(nzmin,n), a_ice(n)) - if ((tracers%data(tr_num)%ID .ge. 6) .and.(tracers%data(tr_num)%ID .le. 40)) then - tr(nzmin)= tr(nzmin)+transit_bc_surface(n, tracers%data(tr_num)%ID, sst(nzmin,n), sss(nzmin,n), a_ice(n), trarr(nzmin,n), nzmin, partit, mesh) - end if - !_______________________________________________________________________ ! The forward sweep algorithm to solve the three-diagonal matrix ! problem @@ -1460,7 +1440,8 @@ end subroutine diff_part_bh !=============================================================================== ! this function returns a boundary conditions for a specified tracer ID and surface node ! ID = 0 and 1 are reserved for temperature and salinity -FUNCTION bc_surface(n, id, sval, nzmin, partit) +! MB: mesh, sst, sss, and aice are only needed for transient tracers +FUNCTION bc_surface(n, id, sval, nzmin, partit, mesh, sst, sss, aice) use MOD_MESH USE MOD_PARTIT USE MOD_PARSUP @@ -1472,14 +1453,32 @@ FUNCTION bc_surface(n, id, sval, nzmin, partit) use recom_glovar #endif use mod_transit + use g_clock implicit none integer, intent(in) :: n, id, nzmin - real(kind=WP), intent(in) :: sval + real(kind=WP), intent(in) :: sval, sst, sss, aice type(t_partit),intent(inout), target :: partit + type(t_mesh), intent(in), target :: mesh + REAL(kind=WP) :: bc_surface character(len=10) :: id_string - real(kind=WP), dimension(:), pointer :: a_ice + real(kind=WP), dimension(:), pointer :: a_ice !! MB: where is this needed? + + if (use_transit) then +#if defined (__oasis) + ! SLP and wind speed in coupled setups. This is a makeshift solution + ! as long as the true values are not provided by the AGCM / OASIS. + press_a = mean_slp + wind_2 = speed_2(stress_atmoce_x(n), stress_atmoce_y(n)) +#else + press_a = press_air(n) + wind_2 = u_wind(n)**2 + v_wind(n)**2 +#endif + ! Atmospheric input of bomb-14C, CFC-12, and SF6 varies with latitude. To that effect specify + y_abc = mesh%geo_coord_nod2D(2,n) / rad ! latitude of atmospheric tracer input + yy_nh = (10. - y_abc) * 0.05 ! interpolation weight for tropical CFC-12 and SF6 values + end if ! use_transit SELECT CASE (id) CASE (1) @@ -1489,6 +1488,90 @@ FUNCTION bc_surface(n, id, sval, nzmin, partit) ! by forming/melting of sea ice bc_surface= dt*(virtual_salt(n) & !--> is zeros for zlevel/zstar + relax_salt(n) - real_salt_flux(n)*is_nonlinfs) +!---Transient tracers (cases ##6,12,14,39) + CASE (6) ! SF6 + if (anthro_transit) then + ! Select atmospheric input values corresponding to the latitude + ! Annual values are interpolated to monthly values, this is omitted in the last simulation year + if (y_abc > 10.) then ! Northern hemisphere + xsf6_a = xsf6_nh(ti_transit) + if (ti_transit < length_transit) xsf6_a = xsf6_a + month * (xsf6_nh(ti_transit + 1) - xsf6_a) / 12. + else if (y_abc <- 10.) then ! Southern hemisphere + xsf6_a = xsf6_sh(ti_transit) + if (ti_transit < length_transit) xsf6_a = xsf6_a + month * (xsf6_sh(ti_transit + 1) - xsf6_a) / 12. + else ! Tropical zone, interpolate between NH and SH + xsf6_a = (1 - yy_nh) * xsf6_nh(ti_transit) + yy_nh * xsf6_sh(ti_transit) + if (ti_transit < length_transit) & + xsf6_a = xsf6_a + month * ((1 - yy_nh) * xsf6_nh(ti_transit + 1) + yy_nh * xsf6_sh(ti_transit + 1) - xsf6_a) / 12. + end if + else + ! Constant (global-mean) namelist values are taken + end if + ! Local air-sea exchange gas flux of SF6 (m / s): + bc_surface = dt * (gas_flux("sf6", sst, sss, wind_2, aice, press_a, xsf6_a, sval) - sval * water_flux(n) * is_nonlinfs) + CASE (11) ! CFC-11 + if (anthro_transit) then + ! Select atmospheric input values corresponding to the latitude + ! Annual values are interpolated to monthly values, this is omitted in the last simulation year + if (y_abc > 10.) then ! Northern hemisphere + xf11_a = xf11_nh(ti_transit) + if (ti_transit < length_transit) xf11_a = xf11_a + month * (xf11_nh(ti_transit + 1) - xf11_a) / 12. + else if (y_abc <- 10.) then ! Southern hemisphere + xf11_a = xf12_sh(ti_transit) + if (ti_transit < length_transit) xf11_a = xf11_a + month * (xf11_sh(ti_transit + 1) - xf11_a) / 12. + else ! Tropical zone, interpolate between NH and SH + xf11_a = (1 - yy_nh) * xf11_nh(ti_transit) + yy_nh * xf11_sh(ti_transit) + if (ti_transit < length_transit) & + xf11_a = xf11_a + month * ((1 - yy_nh) * xf11_nh(ti_transit + 1) + yy_nh * xf11_sh(ti_transit + 1) - xf11_a) / 12. + end if + else + ! Constant (global-mean) namelist values are taken + end if + ! Local air-sea exchange gas flux of CFC-12 (m / s): + bc_surface = dt * (gas_flux("f11", sst, sss, wind_2, aice, press_a, xf11_a, sval) - sval * water_flux(n) * is_nonlinfs) + CASE (12) ! CFC-12 + if (anthro_transit) then + ! Select atmospheric input values corresponding to the latitude + ! Annual values are interpolated to monthly values, this is omitted in the last simulation year + if (y_abc > 10.) then ! Northern hemisphere + xf12_a = xf12_nh(ti_transit) + if (ti_transit < length_transit) xf12_a = xf12_a + month * (xf12_nh(ti_transit + 1) - xf12_a) / 12. + else if (y_abc <- 10.) then ! Southern hemisphere + xf12_a = xf12_sh(ti_transit) + if (ti_transit < length_transit) xf12_a = xf12_a + month * (xf12_sh(ti_transit + 1) - xf12_a) / 12. + else ! Tropical zone, interpolate between NH and SH + xf12_a = (1 - yy_nh) * xf12_nh(ti_transit) + yy_nh * xf12_sh(ti_transit) + if (ti_transit < length_transit) & + xf12_a = xf12_a + month * ((1 - yy_nh) * xf12_nh(ti_transit + 1) + yy_nh * xf12_sh(ti_transit + 1) - xf12_a) / 12. + end if + else + ! Constant (global-mean) namelist values are taken + end if + ! Local air-sea exchange gas flux of CFC-12 (m / s): + bc_surface = dt * (gas_flux("f12", sst, sss, wind_2, aice, press_a, xf12_a, sval) - sval * water_flux(n) * is_nonlinfs) + CASE (14) ! Radiocarbon (more precisely, fractionation-corrected 14C/C): + if (anthro_transit) then + ! Select atmospheric input values corresponding to the latitude + if (y_abc > 30.) then ! Northern hemisphere + r14c_a = r14c_nh(ti_transit) + else if (y_abc <- 30.) then ! Southern hemisphere + r14c_a = r14c_sh(ti_transit) + else ! Tropical zone (values prescribed in input file) + r14c_a = r14c_tz(ti_transit) + end if + xCO2_a = xCO2_ti(ti_transit) + else if (paleo_transit) then + r14c_a = r14c_ti(ti_transit) + xCO2_a = xCO2_ti(ti_transit) + else + ! Constant (global-mean) namelist values are taken + end if + ! Local isotopic 14CO2/CO2 air-sea exchange flux (m / s) + bc_surface = dt * (iso_flux("co2", sst, sss, wind_2, aice, press_a, xco2_a, r14c_a, sval, dic_0) - sval * water_flux(n) * is_nonlinfs) + CASE (39) ! Argon-39 (fractionationation-corrected 39Ar/Ar) + ! Local isotopic 39Ar/Ar air-sea exchange flux (m / s) + bc_surface = dt * (iso_flux("arg", sst, sss, wind_2, aice, press_a, xarg_a, r39ar_a, sval, arg_0) - sval * water_flux(n) * is_nonlinfs) +!--- Done with boundary conditions for transient tracers. #if defined(__recom) CASE (1001) ! DIN bc_surface= dt*(AtmNInput(n) + RiverDIN2D(n) * is_riverinput & @@ -1523,8 +1606,6 @@ FUNCTION bc_surface(n, id, sval, nzmin, partit) #endif CASE (101) ! apply boundary conditions to tracer ID=101 bc_surface= dt*(prec_rain(n))! - real_salt_flux(n)*is_nonlinfs) -!---Transient tracers (case ##6,12,14,39) need additional input parameters -! and are considered in the separate function transit_bc_surface !---wiso-code CASE (102) ! apply boundary conditions to tracer ID=101 (H218O) bc_surface = dt*wiso_flux_oce(n,1) @@ -1562,140 +1643,3 @@ FUNCTION bc_surface(n, id, sval, nzmin, partit) RETURN END FUNCTION - - -!=============================================================================== -! This function returns a boundary conditions for a specified transient tracer ID and surface node. -! Different to function bc_surface, SST, SSS, and sea ice concentrations are always needed as -! auxiliary variable -FUNCTION transit_bc_surface(n, id, sst, sss, aice, sval, nzmin, partit, mesh) - use MOD_MESH - USE MOD_PARTIT - USE MOD_PARSUP - USE o_ARRAYS - USE g_forcing_arrays - USE g_config - use g_clock - use mod_transit - implicit none - - integer, intent(in) :: n, id, nzmin - real(kind=WP), intent(in) :: sst, sss, aice, sval - type(t_partit),intent(inout), target :: partit - type(t_mesh), intent(in), target :: mesh - REAL(kind=WP) :: transit_bc_surface - character(len=10) :: id_string - - - ! --> is_nonlinfs=1.0 for zelvel,zstar .... - ! --> is_nonlinfs=0.0 for linfs - -#if defined (__oasis) -! SLP and wind speed in coupled setups. This is a makeshift solution -! as long as the true values are not provided by the AGCM / OASIS. - press_a = mean_slp - wind_2 = speed_2(stress_atmoce_x(n), stress_atmoce_y(n)) -#else - press_a = press_air(n) - wind_2 = u_wind(n)**2 + v_wind(n)**2 -#endif - -! The atmospheric input of bomb 14C, CFC-12, and SF6 depends on latitude. To that effect specify - y_abc = mesh%geo_coord_nod2D(2,n) / rad ! latitude of atmospheric tracer input - yy_nh = (10. - y_abc) * 0.05 ! interpolation weight for tropical tracer values - - - SELECT CASE (id) - -! Boundary conditions for additional (transient) tracers (14C, 39Ar, CFC-12, and SF6) - CASE (14) ! Radiocarbon (more precisely, fractionation-corrected 14C/C): - if (anthro_transit) then -! Select atmospheric input values corresponding to the latitude - if (y_abc > 30.) then -! Northern Hemisphere - r14c_a = r14c_nh(ti_transit) - else if (y_abc <- 30.) then -! Southern Hemisphere - r14c_a = r14c_sh(ti_transit) - else -! Tropical zone - r14c_a = r14c_tz(ti_transit) - end if - xCO2_a = xCO2_ti(ti_transit) - else if (paleo_transit) then - r14c_a = r14c_ti(ti_transit) - xCO2_a = xCO2_ti(ti_transit) - else -! Constant (global-mean) namelist values are taken - end if -! Local isotopic 14CO2/CO2 air-sea exchange flux (in m / s), -! since F14C is normalized to atmospheric (water) values the isotopic flux has to be -! corrected for precipitation or evaporation fluxes with different isotopic signatures. - transit_bc_surface = dt * (iso_flux("co2", sst, sss, wind_2, aice, press_a, xco2_a, r14c_a, sval, dic_0) & - - sval * water_flux(n) * is_nonlinfs) - - CASE (39) ! Argon-39 (fractionationation-corrected 39Ar/Ar) -! Local isotopic 39Ar/Ar air-sea exchange flux (in m / s), -! since F39Ar is normalized to atmospheric (water) values the isotopic flux has to be -! corrected for precipitation or evaporation fluxes with different isotopic signatures. - transit_bc_surface = dt * (iso_flux("arg", sst, sss, wind_2, aice, press_a, xarg_a, r39ar_a, sval, arg_0) & - - sval * water_flux(n) * is_nonlinfs) - - CASE (12) ! CFC-12 - if (anthro_transit) then -! Select atmospheric input values corresponding to the latitude -! Annual values are interpolated to monthly values, this is omitted in the last simulation year - if (y_abc > 10.) then ! Northern Hemisphere -! Northern Hemisphere - xf12_a = xf12_nh(ti_transit) - if (ti_transit < length_transit) xf12_a = xf12_a + month * (xf12_nh(ti_transit + 1) - xf12_a) / 12. - else if (y_abc <- 10.) then -! Southern Hemisphere - xf12_a = xf12_sh(ti_transit) - if (ti_transit < length_transit) xf12_a = xf12_a + month * (xf12_sh(ti_transit + 1) - xf12_a) / 12. - else -! Tropical zone, interpolate between NH and SH - xf12_a = (1 - yy_nh) * xf12_nh(ti_transit) + yy_nh * xf12_sh(ti_transit) - if (ti_transit < length_transit) & - xf12_a = xf12_a + month * ((1 - yy_nh) * xf12_nh(ti_transit + 1) + yy_nh * xf12_sh(ti_transit + 1) - xf12_a) / 12. - end if - else -! Constant (global-mean) namelist values are taken - end if - -! Local air-sea exchange gas flux of CFC-12 (in m / s): - transit_bc_surface = dt * (gas_flux("f12", sst, sss, wind_2, aice, press_a, xf12_a, sval) & - - sval * water_flux(n) * is_nonlinfs) - - CASE (6) ! SF6 - if (anthro_transit) then -! Select atmospheric input values corresponding to the latitude -! Annual values are interpolated to monthly values, this is omitted in the last simulation year - if (y_abc > 10.) then ! Northern Hemisphere -! Northern Hemisphere - xsf6_a = xsf6_nh(ti_transit) - if (ti_transit < length_transit) xsf6_a = xsf6_a + month * (xsf6_nh(ti_transit + 1) - xsf6_a) / 12. - else if (y_abc <- 10.) then -! Southern Hemisphere - xsf6_a = xsf6_sh(ti_transit) - if (ti_transit < length_transit) xsf6_a = xsf6_a + month * (xsf6_sh(ti_transit + 1) - xsf6_a) / 12. - else -! Tropical zone, interpolate between NH and SH - xsf6_a = (1 - yy_nh) * xsf6_nh(ti_transit) + yy_nh * xsf6_sh(ti_transit) - if (ti_transit < length_transit) & - xsf6_a = xsf6_a + month * ((1 - yy_nh) * xsf6_nh(ti_transit + 1) + yy_nh * xsf6_sh(ti_transit + 1) - xsf6_a) / 12. - end if - else -! Constant (global-mean) namelist values are taken - end if - -! Local air-sea exchange gas flux of SF6 (in m / s): - transit_bc_surface = dt * (gas_flux("sf6", sst, sss, wind_2, aice, press_a, xsf6_a, sval) & - - sval * water_flux(n) * is_nonlinfs) - -! Done with boundary conditions for (transient) tracers. - END SELECT - RETURN - -END FUNCTION - diff --git a/src/oce_setup_step.F90 b/src/oce_setup_step.F90 index 4a1e96d9d..f465155fd 100755 --- a/src/oce_setup_step.F90 +++ b/src/oce_setup_step.F90 @@ -277,7 +277,7 @@ SUBROUTINE tracer_init(tracers, partit, mesh) USE g_ic3d use g_forcing_param, only: use_age_tracer !---age-code use g_config, only : lwiso, use_transit ! add lwiso switch and switch for transient tracers - use mod_transit, only : index_transit + use mod_transit, only : index_transit_r14c, index_transit_r39ar, index_transit_f11, index_transit_f12, index_transit_sf6, l_r14c, l_r39ar, l_f11, l_f12, l_sf6 IMPLICIT NONE type(t_tracer), intent(inout), target :: tracers type(t_partit), intent(inout), target :: partit @@ -375,24 +375,42 @@ SUBROUTINE tracer_init(tracers, partit, mesh) !---age-code-end ! Transient tracers -!! UNDER CONSTRUCTION - Actually we do not want to hardwire the number of transient tracers if (use_transit) then ! add transient tracers to the model - nml_tracer_list(num_tracers+1) = nml_tracer_list(1) - nml_tracer_list(num_tracers+2) = nml_tracer_list(1) - nml_tracer_list(num_tracers+3) = nml_tracer_list(1) - nml_tracer_list(num_tracers+4) = nml_tracer_list(1) - nml_tracer_list(num_tracers+1)%id = 6 - nml_tracer_list(num_tracers+1)%id = 12 - nml_tracer_list(num_tracers+1)%id = 14 - nml_tracer_list(num_tracers+1)%id = 39 - - index_transit(1) = num_tracers+1 - index_transit(2) = num_tracers+2 - index_transit(3) = num_tracers+3 - index_transit(4) = num_tracers+4 - - num_tracers = num_tracers + 4 + if (l_sf6) then + nml_tracer_list(num_tracers+1) = nml_tracer_list(1) + nml_tracer_list(num_tracers+1)%id = 6 + index_transit_sf6 = num_tracers+1 + num_tracers = num_tracers + 1 + endif + + if (l_f11) then + nml_tracer_list(num_tracers+1) = nml_tracer_list(1) + nml_tracer_list(num_tracers+1)%id = 11 + index_transit_f11 = num_tracers+1 + num_tracers = num_tracers + 1 + endif + + if (l_f12) then + nml_tracer_list(num_tracers+1) = nml_tracer_list(1) + nml_tracer_list(num_tracers+1)%id = 12 + index_transit_f12 = num_tracers+1 + num_tracers = num_tracers + 1 + endif + + if (l_r14c) then + nml_tracer_list(num_tracers+1) = nml_tracer_list(1) + nml_tracer_list(num_tracers+1)%id = 14 + index_transit_r14c = num_tracers+1 + num_tracers = num_tracers + 1 + endif + + if (l_r39ar) then + nml_tracer_list(num_tracers+1) = nml_tracer_list(1) + nml_tracer_list(num_tracers+1)%id = 39 + index_transit_r39ar = num_tracers+1 + num_tracers = num_tracers + 1 + endif ! tracers initialised from file idlist((n_ic3d+1):(n_ic3d+1)) = (/14/) @@ -1161,19 +1179,16 @@ SUBROUTINE oce_initial_state(tracers, partit, mesh) !---wiso-code-end ! Transient tracers - CASE (14) ! initialize tracer ID=14, fractionation-corrected 14C/C -! this initialization can be overwritten by calling do_ic3d -!! if (.not. any(idlist == 14)) then ! CHECK IF THIS LINE IS STILL NECESSARY - tracers%data(i)%values(:,:) = 0.85 - if (mype==0) then - write (i_string, "(I3)") i - write (id_string, "(I3)") id - write(*,*) 'initializing '//trim(i_string)//'th tracer with ID='//trim(id_string) - write (*,*) tracers%data(i)%values(1,1) - end if -!! end if - CASE (39) ! initialize tracer ID=39, fractionation-corrected 39Ar/Ar - tracers%data(i)%values(:,:) = 0.85 + CASE (6) ! initialize tracer ID=6, SF6 + tracers%data(i)%values(:,:) = 0.0_WP + if (mype==0) then + write (i_string, "(I3)") i + write (id_string, "(I3)") id + write(*,*) 'initializing '//trim(i_string)//'th tracer with ID='//trim(id_string) + write (*,*) tracers%data(i)%values(1,1) + end if + CASE (11) ! initialize tracer ID=11, CFC-11 + tracers%data(i)%values(:,:) = 0.0_WP if (mype==0) then write (i_string, "(I3)") i write (id_string, "(I3)") id @@ -1181,15 +1196,26 @@ SUBROUTINE oce_initial_state(tracers, partit, mesh) write (*,*) tracers%data(i)%values(1,1) end if CASE (12) ! initialize tracer ID=12, CFC-12 - tracers%data(i)%values(:,:) = 0. + tracers%data(i)%values(:,:) = 0.0_WP if (mype==0) then write (i_string, "(I3)") i write (id_string, "(I3)") id write(*,*) 'initializing '//trim(i_string)//'th tracer with ID='//trim(id_string) write (*,*) tracers%data(i)%values(1,1) end if - CASE (6) ! initialize tracer ID=6, SF6 - tracers%data(i)%values(:,:) = 0. + CASE (14) ! initialize tracer ID=14, fractionation-corrected 14C/C +! this initialization can be overwritten by calling do_ic3d + if (.not. any(idlist == 14)) then ! CHECK IF THIS LINE IS STILL NECESSARY + tracers%data(i)%values(:,:) = 0.85 + if (mype==0) then + write (i_string, "(I3)") i + write (id_string, "(I3)") id + write(*,*) 'initializing '//trim(i_string)//'th tracer with ID='//trim(id_string) + write (*,*) tracers%data(i)%values(1,1) + end if + end if + CASE (39) ! initialize tracer ID=39, fractionation-corrected 39Ar/Ar + tracers%data(i)%values(:,:) = 0.50 if (mype==0) then write (i_string, "(I3)") i write (id_string, "(I3)") id