Transient tracers bugfixed

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,
 /

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

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

src/gen_model_setup.F90

@@ -125,18 +125,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'

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
@@ -257,17 +257,16 @@ 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
 CASE ('h_snow    ')
     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
@@ -873,13 +872,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)
@@ -915,17 +934,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        ')

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'