Changeset 7978 for branches/UKMO/dev_r5518_medusa_chg_trc_bio_medusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcbio_medusa.F90
- Timestamp:
- 2017-04-27T13:13:32+02:00 (7 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/UKMO/dev_r5518_medusa_chg_trc_bio_medusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcbio_medusa.F90
r7975 r7978 101 101 USE air_sea_mod, ONLY: air_sea 102 102 USE plankton_mod, ONLY: plankton 103 USE iron_chem_scav_mod, ONLY: iron_chem_scav 103 104 USE bio_medusa_diag_slice_mod, ONLY: bio_medusa_diag_slice 104 105 USE bio_medusa_fin_mod, ONLY: bio_medusa_fin … … 207 208 !! 208 209 !! iron cycle; includes parameters for Parekh et al. (2005) iron scheme 209 REAL(wp), DIMENSION(jpi,jpj) :: ffetop,ffebot,ffescav210 ! REAL(wp), DIMENSION(jpi,jpj) :: ffetop,ffebot,ffescav 210 211 REAL(wp) :: xLgF, xFeT, xFeF, xFeL !! state variables for iron-ligand system 211 212 ! REAL(wp), DIMENSION(jpi,jpj) :: xFree !! state variables for iron-ligand system … … 234 235 ! REAL(wp), DIMENSION(jpi,jpj) :: fdpn,fdpd,fdpds,fdzmi,fdzme,fdd 235 236 # if defined key_roam 236 REAL(wp), DIMENSION(jpi,jpj) :: fddc237 ! REAL(wp), DIMENSION(jpi,jpj) :: fddc 237 238 # endif 238 239 ! REAL(wp), DIMENSION(jpi,jpj) :: fdpn2,fdpd2,fdpds2,fdzmi2,fdzme2 … … 249 250 !! particle flux 250 251 ! REAL(WP), DIMENSION(jpi,jpj) :: fdep1,fcaco3 251 REAL(WP), DIMENSION(jpi,jpj) :: ftempn,ftempsi,ftempfe,ftempc,ftempca252 REAL(wp), DIMENSION(jpi,jpj) :: freminn,freminsi,freminfe,freminc,freminca252 ! REAL(WP), DIMENSION(jpi,jpj) :: ftempn,ftempsi,ftempfe,ftempc,ftempca 253 ! REAL(wp), DIMENSION(jpi,jpj) :: freminn,freminsi,freminfe,freminc,freminca 253 254 ! REAL(wp), DIMENSION(jpi,jpj) :: ffastn,ffastsi,ffastfe,ffastc,ffastca 254 REAL(wp), DIMENSION(jpi,jpj) :: fprotf255 ! REAL(wp), DIMENSION(jpi,jpj) :: fprotf 255 256 ! REAL(wp), DIMENSION(jpi,jpj) :: fsedn,fsedsi,fsedfe,fsedc,fsedca 256 257 ! REAL(wp), DIMENSION(jpi,jpj) :: fccd 257 REAL(wp), DIMENSION(jpi,jpj) :: fccd_dep258 ! REAL(wp), DIMENSION(jpi,jpj) :: fccd_dep 258 259 !! 259 260 !! AXY (06/07/11): alternative fast detritus schemes … … 937 938 ENDDO 938 939 939 DO jj = 2,jpjm1 940 DO ji = 2,jpim1 941 !! OPEN wet point IF..THEN loop 942 if (tmask(ji,jj,jk) == 1) then 943 944 !!---------------------------------------------------------------------- 945 !! Iron chemistry and fractionation 946 !! following the Parekh et al. (2004) scheme adopted by the Met. 947 !! Office, Medusa models total iron but considers "free" and 948 !! ligand-bound forms for the purposes of scavenging (only "free" 949 !! iron can be scavenged 950 !!---------------------------------------------------------------------- 951 !! 952 !! total iron concentration (mmol Fe / m3 -> umol Fe / m3) 953 xFeT = zfer(ji,jj) * 1.e3 954 !! 955 !! calculate fractionation (based on Diat-HadOCC; in turn based on Parekh et al., 2004) 956 xb_coef_tmp = xk_FeL * (xLgT - xFeT) - 1.0 957 xb2M4ac = max(((xb_coef_tmp * xb_coef_tmp) + (4.0 * xk_FeL * xLgT)), 0.0) 958 !! 959 !! "free" ligand concentration 960 xLgF = 0.5 * (xb_coef_tmp + (xb2M4ac**0.5)) / xk_FeL 961 !! 962 !! ligand-bound iron concentration 963 xFeL = xLgT - xLgF 964 !! 965 !! "free" iron concentration (and convert to mmol Fe / m3) 966 xFeF = (xFeT - xFeL) * 1.e-3 967 xFree(ji,jj)= xFeF / (zfer(ji,jj) + tiny(zfer(ji,jj))) 968 !! 969 !! scavenging of iron (multiple schemes); I'm only really happy with the 970 !! first one at the moment - the others involve assumptions (sometimes 971 !! guessed at by me) that are potentially questionable 972 !! 973 if (jiron.eq.1) then 974 !!---------------------------------------------------------------------- 975 !! Scheme 1: Dutkiewicz et al. (2005) 976 !! This scheme includes a single scavenging term based solely on a 977 !! fixed rate and the availablility of "free" iron 978 !!---------------------------------------------------------------------- 979 !! 980 ffescav(ji,jj) = xk_sc_Fe * xFeF ! = mmol/m3/d 981 !! 982 !!---------------------------------------------------------------------- 983 !! 984 !! Mick's code contains a further (optional) implicit "scavenging" of 985 !! iron that sets an upper bound on "free" iron concentration, and 986 !! essentially caps the concentration of total iron as xFeL + "free" 987 !! iron; since the former is constrained by a fixed total ligand 988 !! concentration (= 1.0 umol/m3), and the latter isn't allowed above 989 !! this upper bound, total iron is constrained to a maximum of ... 990 !! 991 !! xFeL + min(xFeF, 0.3 umol/m3) = 1.0 + 0.3 = 1.3 umol / m3 992 !! 993 !! In Mick's code, the actual value of total iron is reset to this 994 !! sum (i.e. TFe = FeL + Fe'; but Fe' <= 0.3 umol/m3); this isn't 995 !! our favoured approach to tracer updating here (not least because 996 !! of the leapfrog), so here the amount scavenged is augmented by an 997 !! additional amount that serves to drag total iron back towards that 998 !! expected from this limitation on iron concentration ... 999 !! 1000 xmaxFeF = min((xFeF * 1.e3), 0.3) ! = umol/m3 1001 !! 1002 !! Here, the difference between current total Fe and (FeL + Fe') is 1003 !! calculated and added to the scavenging flux already calculated 1004 !! above ... 1005 !! 1006 fdeltaFe = (xFeT - (xFeL + xmaxFeF)) * 1.e-3 ! = mmol/m3 1007 !! 1008 !! This assumes that the "excess" iron is dissipated with a time- 1009 !! scale of 1 day; seems reasonable to me ... (famous last words) 1010 !! 1011 ffescav(ji,jj) = ffescav(ji,jj) + fdeltaFe ! = mmol/m3/d 1012 !! 1013 # if defined key_deep_fe_fix 1014 !! AXY (17/01/13) 1015 !! stop scavenging for iron concentrations below 0.5 umol / m3 1016 !! at depths greater than 1000 m; this aims to end MEDUSA's 1017 !! continual loss of iron at depth without impacting things 1018 !! at the surface too much; the justification for this is that 1019 !! it appears to be what Mick Follows et al. do in their work 1020 !! (as evidenced by the iron initial condition they supplied 1021 !! me with); to be honest, it looks like Follow et al. do this 1022 !! at shallower depths than 1000 m, but I'll stick with this 1023 !! for now; I suspect that this seemingly arbitrary approach 1024 !! effectively "parameterises" the particle-based scavenging 1025 !! rates that other models use (i.e. at depth there are no 1026 !! sinking particles, so scavenging stops); it might be fun 1027 !! justifying this in a paper though! 1028 !! 1029 if ((fsdepw(ji,jj,jk).gt.1000.) .and. (xFeT.lt.0.5)) then 1030 ffescav(ji,jj) = 0. 1031 endif 1032 # endif 1033 !! 1034 elseif (jiron.eq.2) then 1035 !!---------------------------------------------------------------------- 1036 !! Scheme 2: Moore et al. (2004) 1037 !! This scheme includes a single scavenging term that accounts for 1038 !! both suspended and sinking particles in the water column; this 1039 !! term scavenges total iron rather than "free" iron 1040 !!---------------------------------------------------------------------- 1041 !! 1042 !! total iron concentration (mmol Fe / m3 -> umol Fe / m3) 1043 xFeT = zfer(ji,jj) * 1.e3 1044 !! 1045 !! this has a base scavenging rate (12% / y) which is modified by local 1046 !! particle concentration and sinking flux (and dust - but I'm ignoring 1047 !! that here for now) and which is accelerated when Fe concentration gets 1048 !! 0.6 nM (= 0.6 umol/m3 = 0.0006 mmol/m3), and decreased as concentrations 1049 !! below 0.4 nM (= 0.4 umol/m3 = 0.0004 mmol/m3) 1050 !! 1051 !! base scavenging rate (0.12 / y) 1052 fbase_scav = 0.12 / 365.25 1053 !! 1054 !! calculate sinking particle part of scaling factor 1055 !! this takes local fast sinking carbon (mmol C / m2 / d) and 1056 !! gets it into nmol C / cm3 / s ("rdt" below is the number of seconds in 1057 !! a model timestep) 1058 !! 1059 !! fscal_sink = ffastc(ji,jj) * 1.e2 / (86400.) 1060 !! 1061 !! ... actually, re-reading Moore et al.'s equations, it looks like he uses 1062 !! his sinking flux directly, without scaling it by time-step or anything, 1063 !! so I'll copy this here ... 1064 !! 1065 fscal_sink = ffastc(ji,jj) * 1.e2 1066 !! 1067 !! calculate particle part of scaling factor 1068 !! this totals up the carbon in suspended particles (Pn, Pd, Zmi, Zme, D), 1069 !! which comes out in mmol C / m3 (= nmol C / cm3), and then multiplies it 1070 !! by a magic factor, 0.002, to get it into nmol C / cm2 / s 1071 !! 1072 fscal_part = ((xthetapn * zphn(ji,jj)) + (xthetapd * zphd(ji,jj)) + (xthetazmi * zzmi(ji,jj)) + & 1073 (xthetazme * zzme(ji,jj)) + (xthetad * zdet(ji,jj))) * 0.002 1074 !! 1075 !! calculate scaling factor for base scavenging rate 1076 !! this uses the (now correctly scaled) sinking flux and standing 1077 !! particle concentration, divides through by some sort of reference 1078 !! value (= 0.0066 nmol C / cm2 / s) and then uses this, or not if its 1079 !! too high, to rescale the base scavenging rate 1080 !! 1081 fscal_scav = fbase_scav * min(((fscal_sink + fscal_part) / 0.0066), 4.0) 1082 !! 1083 !! the resulting scavenging rate is then scaled further according to the 1084 !! local iron concentration (i.e. diminished in low iron regions; enhanced 1085 !! in high iron regions; less alone in intermediate iron regions) 1086 !! 1087 if (xFeT.lt.0.4) then 1088 !! 1089 !! low iron region 1090 !! 1091 fscal_scav = fscal_scav * (xFeT / 0.4) 1092 !! 1093 elseif (xFeT.gt.0.6) then 1094 !! 1095 !! high iron region 1096 !! 1097 fscal_scav = fscal_scav + ((xFeT / 0.6) * (6.0 / 1.4)) 1098 !! 1099 else 1100 !! 1101 !! intermediate iron region: do nothing 1102 !! 1103 endif 1104 !! 1105 !! apply the calculated scavenging rate ... 1106 !! 1107 ffescav(ji,jj) = fscal_scav * zfer(ji,jj) 1108 !! 1109 elseif (jiron.eq.3) then 1110 !!---------------------------------------------------------------------- 1111 !! Scheme 3: Moore et al. (2008) 1112 !! This scheme includes a single scavenging term that accounts for 1113 !! sinking particles in the water column, and includes organic C, 1114 !! biogenic opal, calcium carbonate and dust in this (though the 1115 !! latter is ignored here until I work out what units the incoming 1116 !! "dust" flux is in); this term scavenges total iron rather than 1117 !! "free" iron 1118 !!---------------------------------------------------------------------- 1119 !! 1120 !! total iron concentration (mmol Fe / m3 -> umol Fe / m3) 1121 xFeT = zfer(ji,jj) * 1.e3 1122 !! 1123 !! this has a base scavenging rate which is modified by local 1124 !! particle sinking flux (including dust - but I'm ignoring that 1125 !! here for now) and which is accelerated when Fe concentration 1126 !! is > 0.6 nM (= 0.6 umol/m3 = 0.0006 mmol/m3), and decreased as 1127 !! concentrations < 0.5 nM (= 0.5 umol/m3 = 0.0005 mmol/m3) 1128 !! 1129 !! base scavenging rate (Fe_b in paper; units may be wrong there) 1130 fbase_scav = 0.00384 ! (ng)^-1 cm 1131 !! 1132 !! calculate sinking particle part of scaling factor; this converts 1133 !! mmol / m2 / d fluxes of organic carbon, silicon and calcium 1134 !! carbonate into ng / cm2 / s fluxes; it is assumed here that the 1135 !! mass conversions simply consider the mass of the main element 1136 !! (C, Si and Ca) and *not* the mass of the molecules that they are 1137 !! part of; Moore et al. (2008) is unclear on the conversion that 1138 !! should be used 1139 !! 1140 !! milli -> nano; mol -> gram; /m2 -> /cm2; /d -> /s 1141 fscal_csink = (ffastc(ji,jj) * 1.e6 * xmassc * 1.e-4 / 86400.) ! ng C / cm2 / s 1142 fscal_sisink = (ffastsi(ji,jj) * 1.e6 * xmasssi * 1.e-4 / 86400.) ! ng Si / cm2 / s 1143 fscal_casink = (ffastca(ji,jj) * 1.e6 * xmassca * 1.e-4 / 86400.) ! ng Ca / cm2 / s 1144 !! 1145 !! sum up these sinking fluxes and convert to ng / cm by dividing 1146 !! through by a sinking rate of 100 m / d = 1.157 cm / s 1147 fscal_sink = ((fscal_csink * 6.) + fscal_sisink + fscal_casink) / & 1148 (100. * 1.e3 / 86400) ! ng / cm 1149 !! 1150 !! now calculate the scavenging rate based upon the base rate and 1151 !! this particle flux scaling; according to the published units, 1152 !! the result actually has *no* units, but as it must be expressed 1153 !! per unit time for it to make any sense, I'm assuming a missing 1154 !! "per second" 1155 fscal_scav = fbase_scav * fscal_sink ! / s 1156 !! 1157 !! the resulting scavenging rate is then scaled further according to the 1158 !! local iron concentration (i.e. diminished in low iron regions; enhanced 1159 !! in high iron regions; less alone in intermediate iron regions) 1160 !! 1161 if (xFeT.lt.0.5) then 1162 !! 1163 !! low iron region (0.5 instead of the 0.4 in Moore et al., 2004) 1164 !! 1165 fscal_scav = fscal_scav * (xFeT / 0.5) 1166 !! 1167 elseif (xFeT.gt.0.6) then 1168 !! 1169 !! high iron region (functional form different in Moore et al., 2004) 1170 !! 1171 fscal_scav = fscal_scav + ((xFeT - 0.6) * 0.00904) 1172 !! 1173 else 1174 !! 1175 !! intermediate iron region: do nothing 1176 !! 1177 endif 1178 !! 1179 !! apply the calculated scavenging rate ... 1180 !! 1181 ffescav(ji,jj) = fscal_scav * zfer(ji,jj) 1182 !! 1183 elseif (jiron.eq.4) then 1184 !!---------------------------------------------------------------------- 1185 !! Scheme 4: Galbraith et al. (2010) 1186 !! This scheme includes two scavenging terms, one for organic, 1187 !! particle-based scavenging, and another for inorganic scavenging; 1188 !! both terms scavenge "free" iron only 1189 !!---------------------------------------------------------------------- 1190 !! 1191 !! Galbraith et al. (2010) present a more straightforward outline of 1192 !! the scheme in Parekh et al. (2005) ... 1193 !! 1194 !! sinking particulate carbon available for scavenging 1195 !! this assumes a sinking rate of 100 m / d (Moore & Braucher, 2008), 1196 xCscav1 = (ffastc(ji,jj) * xmassc) / 100. ! = mg C / m3 1197 !! 1198 !! scale by Honeyman et al. (1981) exponent coefficient 1199 !! multiply by 1.e-3 to express C flux in g C rather than mg C 1200 xCscav2 = (xCscav1 * 1.e-3)**0.58 1201 !! 1202 !! multiply by Galbraith et al. (2010) scavenging rate 1203 xk_org = 0.5 ! ((g C m/3)^-1) / d 1204 xORGscav = xk_org * xCscav2 * xFeF 1205 !! 1206 !! Galbraith et al. (2010) also include an inorganic bit ... 1207 !! 1208 !! this occurs at a fixed rate, again based on the availability of 1209 !! "free" iron 1210 !! 1211 !! k_inorg = 1000 d**-1 nmol Fe**-0.5 kg**-0.5 1212 !! 1213 !! to implement this here, scale xFeF by 1026 to put in units of 1214 !! umol/m3 which approximately equal nmol/kg 1215 !! 1216 xk_inorg = 1000. ! ((nmol Fe / kg)^1.5) 1217 xINORGscav = (xk_inorg * (xFeF * 1026.)**1.5) * 1.e-3 1218 !! 1219 !! sum these two terms together 1220 ffescav(ji,jj) = xORGscav + xINORGscav 1221 else 1222 !!---------------------------------------------------------------------- 1223 !! No Scheme: you coward! 1224 !! This scheme puts its head in the sand and eskews any decision about 1225 !! how iron is removed from the ocean; prepare to get deluged in iron 1226 !! you fool! 1227 !!---------------------------------------------------------------------- 1228 ffescav(ji,jj) = 0. 1229 endif 1230 1231 !!---------------------------------------------------------------------- 1232 !! Other iron cycle processes 1233 !!---------------------------------------------------------------------- 1234 !! 1235 !! aeolian iron deposition 1236 if (jk.eq.1) then 1237 !! zirondep is in mmol-Fe / m2 / day 1238 !! ffetop(ji,jj) is in mmol-dissolved-Fe / m3 / day 1239 ffetop(ji,jj) = zirondep(ji,jj) * xfe_sol / fse3t(ji,jj,jk) 1240 else 1241 ffetop(ji,jj) = 0.0 1242 endif 1243 !! 1244 !! seafloor iron addition 1245 !! AXY (10/07/12): amended to only apply sedimentary flux up to ~500 m down 1246 !! if (jk.eq.(mbathy(ji,jj)-1).AND.jk.lt.i1100) then 1247 if ((jk.eq.mbathy(ji,jj)).AND.jk.le.i0500) then 1248 !! Moore et al. (2004) cite a coastal California value of 5 umol/m2/d, but adopt a 1249 !! global value of 2 umol/m2/d for all areas < 1100 m; here we use this latter value 1250 !! but apply it everywhere 1251 !! AXY (21/07/09): actually, let's just apply it below 1100 m (levels 1-37) 1252 ffebot(ji,jj) = (xfe_sed / fse3t(ji,jj,jk)) 1253 else 1254 ffebot(ji,jj) = 0.0 1255 endif 1256 1257 !! AXY (16/12/09): remove iron addition/removal processes 1258 !! For the purposes of the quarter degree run, the iron cycle is being pegged to the 1259 !! initial condition supplied by Mick Follows via restoration with a 30 day period; 1260 !! iron addition at the seafloor is still permitted with the idea that this extra 1261 !! iron will be removed by the restoration away from the source 1262 !! ffescav(ji,jj) = 0.0 1263 !! ffetop(ji,jj) = 0.0 1264 !! ffebot(ji,jj) = 0.0 1265 1266 # if defined key_debug_medusa 1267 !! report miscellaneous calculations 1268 if (idf.eq.1.AND.idfval.eq.1) then 1269 IF (lwp) write (numout,*) '------------------------------' 1270 IF (lwp) write (numout,*) 'xfe_sol = ', xfe_sol 1271 IF (lwp) write (numout,*) 'xfe_mass = ', xfe_mass 1272 IF (lwp) write (numout,*) 'ffetop(',jk,') = ', ffetop(ji,jj) 1273 IF (lwp) write (numout,*) 'ffebot(',jk,') = ', ffebot(ji,jj) 1274 IF (lwp) write (numout,*) 'xFree(',jk,') = ', xFree(ji,jj) 1275 IF (lwp) write (numout,*) 'ffescav(',jk,') = ', ffescav(ji,jj) 1276 endif 1277 # endif 1278 1279 ! MAYBE BUT A BREAK IN HERE, MISCELLANEOUS PROCESSES - marc 20/4/17 1280 ! (iron chemistry and scavenging is 340 lines) 1281 ENDIF 1282 ENDDO 1283 ENDDO 940 !!---------------------------------------------------------------- 941 !! Iron chemistry and scavenging 942 !!---------------------------------------------------------------- 943 CALL iron_chem_scav( jk ) 944 945 ! Miscellaneous processes - marc 1284 946 1285 947 DO jj = 2,jpjm1
Note: See TracChangeset
for help on using the changeset viewer.
