Changeset 8436 for branches/UKMO/dev_r5518_GO6_package_asm_surf_bgc/NEMOGCM/NEMO/OPA_SRC/ASM/asmlogchlbal_medusa.F90
- Timestamp:
- 2017-08-14T15:22:09+02:00 (7 years ago)
- File:
-
- 1 copied
Legend:
- Unmodified
- Added
- Removed
-
branches/UKMO/dev_r5518_GO6_package_asm_surf_bgc/NEMOGCM/NEMO/OPA_SRC/ASM/asmlogchlbal_medusa.F90
r8428 r8436 1 MODULE asmlogchlbal_ hadocc1 MODULE asmlogchlbal_medusa 2 2 !!====================================================================== 3 !! *** MODULE asmlogchlbal_ hadocc***4 !! Calculate increments to HadOCCbased on surface logchl increments3 !! *** MODULE asmlogchlbal_medusa *** 4 !! Calculate increments to MEDUSA based on surface logchl increments 5 5 !! 6 6 !! IMPORTANT NOTE: This calls the bioanalysis routine of Hemmings et al. … … 10 10 !! 11 11 !!====================================================================== 12 !! History : 3.6 ! 2017-08 (D. Ford) Adapted from bioanal.F9013 !!---------------------------------------------------------------------- 14 #if defined key_asminc && defined key_ hadocc12 !! History : 3.6 ! 2017-08 (D. Ford) Adapted from asmlogchlbal_hadocc 13 !!---------------------------------------------------------------------- 14 #if defined key_asminc && defined key_medusa && defined key_foam_medusa 15 15 !!---------------------------------------------------------------------- 16 16 !! 'key_asminc' : assimilation increment interface 17 !! 'key_hadocc' : HadOCC model 18 !!---------------------------------------------------------------------- 19 !! asm_logchl_bal_hadocc : routine to calculate increments to HadOCC 17 !! 'key_medusa' : MEDUSA model 18 !! 'key_foam_medusa' : MEDUSA extras for FOAM OBS and ASM 19 !!---------------------------------------------------------------------- 20 !! asm_logchl_bal_medusa : routine to calculate increments to MEDUSA 20 21 !!---------------------------------------------------------------------- 21 22 USE par_kind, ONLY: wp ! kind parameters … … 24 25 USE zdfmxl ! mixed layer depth 25 26 USE iom ! i/o 26 USE trc, ONLY: trn, trb, & ! HadOCC variables 27 & HADOCC_CHL, & 28 & pgrow_avg, & 29 & ploss_avg, & 30 & phyt_avg, & 31 & mld_max 32 USE par_hadocc ! HadOCC parameters 33 USE had_bgc_stnd, ONLY: kmt ! HadOCC parameters 34 USE had_bgc_const ! HadOCC parameters 27 USE trc, ONLY: trn, trb ! MEDUSA variables 28 USE sms_medusa ! MEDUSA parameters 29 USE par_medusa ! MEDUSA parameters 35 30 USE par_trc, ONLY: jptra ! Tracer parameters 36 31 USE bioanalysis ! Nitrogen balancing … … 39 34 PRIVATE 40 35 41 PUBLIC asm_logchl_bal_ hadocc36 PUBLIC asm_logchl_bal_medusa 42 37 43 38 ! Default values for biological assimilation parameters … … 74 69 CONTAINS 75 70 76 SUBROUTINE asm_logchl_bal_hadocc( logchl_bkginc, aincper, mld_choice_bgc, & 77 & k_maxchlinc, ld_logchlbal, logchl_balinc ) 71 SUBROUTINE asm_logchl_bal_medusa( logchl_bkginc, aincper, mld_choice_bgc, & 72 & k_maxchlinc, ld_logchlbal, & 73 & pgrow_avg_bkg, ploss_avg_bkg, & 74 & phyt_avg_bkg, mld_max_bkg, & 75 & logchl_balinc ) 78 76 !!--------------------------------------------------------------------------- 79 !! *** ROUTINE asm_logchl_bal_ hadocc***80 !! 81 !! ** Purpose : calculate increments to HadOCCfrom logchl increments77 !! *** ROUTINE asm_logchl_bal_medusa *** 78 !! 79 !! ** Purpose : calculate increments to MEDUSA from logchl increments 82 80 !! 83 81 !! ** Method : convert logchl increments to chl increments 82 !! average up MEDUSA to look like HadOCC 84 83 !! call nitrogen balancing scheme 84 !! separate back out to MEDUSA 85 85 !! 86 86 !! ** Action : populate logchl_balinc … … 95 95 REAL(wp), INTENT(in ) :: k_maxchlinc ! Max chl increment 96 96 LOGICAL, INTENT(in ) :: ld_logchlbal ! Balancing y/n 97 REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: pgrow_avg_bkg ! Avg phyto growth 98 REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: ploss_avg_bkg ! Avg phyto loss 99 REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: phyt_avg_bkg ! Avg phyto 100 REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: mld_max_bkg ! Max MLD 97 101 REAL(wp), INTENT( out), DIMENSION(jpi,jpj,jpk,jptra) :: logchl_balinc ! Balancing increments 98 102 !! 99 INTEGER :: ji, jj, jk 103 INTEGER :: ji, jj, jk, jn ! Loop counters 100 104 INTEGER :: jkmax ! Loop index 101 105 INTEGER, DIMENSION(6) :: i_tracer ! Tracer indices 106 REAL(wp) :: n2be_p ! N:biomass for total phy 107 REAL(wp) :: n2be_z ! N:biomass for total zoo 108 REAL(wp) :: n2be_d ! N:biomass for detritus 109 REAL(wp) :: zfrac_chn ! Fraction of jpchn 110 REAL(wp) :: zfrac_chd ! Fraction of jpchd 111 REAL(wp) :: zfrac_phn ! Fraction of jpphn 112 REAL(wp) :: zfrac_phd ! Fraction of jpphd 113 REAL(wp) :: zfrac_zmi ! Fraction of jpzmi 114 REAL(wp) :: zfrac_zme ! Fraction of jpzme 115 REAL(wp) :: zrat_pds_phd ! Ratio of jppds:jpphd 116 REAL(wp) :: zrat_chd_phd ! Ratio of jpchd:jpphd 117 REAL(wp) :: zrat_chn_phn ! Ratio of jpchn:jpphn 118 REAL(wp) :: zrat_dtc_det ! Ratio of jpdtc:jpdet 102 119 REAL(wp), DIMENSION(jpi,jpj) :: chl_inc ! Chlorophyll increments 120 REAL(wp), DIMENSION(jpi,jpj) :: medusa_chl ! MEDUSA total chlorophyll 121 REAL(wp), DIMENSION(jpi,jpj) :: cchl_p ! C:Chl for total phy 103 122 REAL(wp), DIMENSION(jpi,jpj) :: zmld ! Mixed layer depth 104 123 REAL(wp), DIMENSION(16) :: modparm ! Model parameters … … 117 136 ! 3) Subtract background from analysis to get chl incs 118 137 ! If k_maxchlinc > 0 then cap total absolute chlorophyll increment at that value 138 medusa_chl(:,:) = trb(:,:,1,jpchn) + trb(:,:,1,jpchd) 119 139 DO jj = 1, jpj 120 140 DO ji = 1, jpi 121 IF ( HADOCC_CHL(ji,jj,1) > 0.0 ) THEN122 chl_inc(ji,jj) = 10**( LOG10( HADOCC_CHL(ji,jj,1) ) + logchl_bkginc(ji,jj) ) - HADOCC_CHL(ji,jj,1)141 IF ( medusa_chl(ji,jj) > 0.0 ) THEN 142 chl_inc(ji,jj) = 10**( LOG10( medusa_chl(ji,jj) ) + logchl_bkginc(ji,jj) ) - medusa_chl(ji,jj) 123 143 IF ( k_maxchlinc > 0.0 ) THEN 124 144 chl_inc(ji,jj) = MAX( -1.0 * k_maxchlinc, MIN( chl_inc(ji,jj), k_maxchlinc ) ) … … 158 178 IF ( ld_logchlbal ) THEN ! Nitrogen balancing 159 179 160 ! Set up model parameters to be passed into Hemmings balancing routine 161 modparm(1) = grow_sat 162 modparm(2) = psmax 163 modparm(3) = par 164 modparm(4) = alpha 165 modparm(5) = resp_rate 166 modparm(6) = pmort_rate 167 modparm(7) = phyto_min 168 modparm(8) = z_mort_1 169 modparm(9) = z_mort_2 170 modparm(10) = c2n_p 171 modparm(11) = c2n_z 172 modparm(12) = c2n_d 173 modparm(13) = graze_threshold 174 modparm(14) = holling_coef 175 modparm(15) = graze_sat 176 modparm(16) = graze_max 180 ! Set up model parameters to be passed into Hemmings balancing routine. 181 ! For now these are hardwired to the standard HadOCC parameter values 182 ! (except C:N ratios) as this is what the scheme was developed for. 183 ! Obviously, HadOCC and MEDUSA are rather different models, so this 184 ! isn't ideal, but there's not always direct analogues between the two 185 ! parameter sets, so it's the easiest way to get something running. 186 ! In the longer term, some serious MarMOT-based development is required. 187 modparm(1) = 0.1 ! grow_sat 188 modparm(2) = 2.0 ! psmax 189 modparm(3) = 0.845 ! par 190 modparm(4) = 0.02 ! alpha 191 modparm(5) = 0.05 ! resp_rate 192 modparm(6) = 0.05 ! pmort_rate 193 modparm(7) = 0.01 ! phyto_min 194 modparm(8) = 0.05 ! z_mort_1 195 modparm(9) = 1.0 ! z_mort_2 196 modparm(10) = ( xthetapn + xthetapd ) / 2.0 ! c2n_p 197 modparm(11) = ( xthetazmi + xthetazme ) / 2.0 ! c2n_z 198 modparm(12) = xthetad ! c2n_d 199 modparm(13) = 0.01 ! graze_threshold 200 modparm(14) = 2.0 ! holling_coef 201 modparm(15) = 0.5 ! graze_sat 202 modparm(16) = 2.0 ! graze_max 177 203 178 204 ! Set up assimilation parameters to be passed into balancing routine … … 207 233 208 234 ! Set background state 209 bstate(:,:,:,i_tracer(1)) = trb(:,:,:,jp_had_nut) 210 bstate(:,:,:,i_tracer(2)) = trb(:,:,:,jp_had_phy) 211 bstate(:,:,:,i_tracer(3)) = trb(:,:,:,jp_had_zoo) 212 bstate(:,:,:,i_tracer(4)) = trb(:,:,:,jp_had_pdn) 213 bstate(:,:,:,i_tracer(5)) = trb(:,:,:,jp_had_dic) 214 bstate(:,:,:,i_tracer(6)) = trb(:,:,:,jp_had_alk) 235 bstate(:,:,:,i_tracer(1)) = trb(:,:,:,jpdin) 236 bstate(:,:,:,i_tracer(2)) = trb(:,:,:,jpphn) + trb(:,:,:,jpphd) 237 bstate(:,:,:,i_tracer(3)) = trb(:,:,:,jpzmi) + trb(:,:,:,jpzme) 238 bstate(:,:,:,i_tracer(4)) = trb(:,:,:,jpdet) 239 bstate(:,:,:,i_tracer(5)) = trb(:,:,:,jpdic) 240 bstate(:,:,:,i_tracer(6)) = trb(:,:,:,jpalk) 241 242 ! Calculate carbon to chlorophyll ratio for combined phytoplankton 243 ! and nitrogen to biomass equivalent for PZD 244 ! Hardwire nitrogen mass to 14.01 for now as it doesn't seem to be set in MEDUSA 245 !cchl_p(:,:) = ( trb(:,:,1,jpchn) + trb(:,:,1,jpchd ) ) / & 246 ! & ( ( trb(:,:,1,jpphn) * xthetapn ) + ( trb(:,:,1,jpphd) * xthetapd ) ) 247 cchl_p(:,:) = 0.0 248 DO jj = 1, jpj 249 DO ji = 1, jpi 250 IF ( ( trb(ji,jj,1,jpchn) + trb(ji,jj,1,jpchd ) ) .GT. 0.0 ) THEN 251 cchl_p(ji,jj) = ( ( trb(ji,jj,1,jpphn) * xthetapn ) + ( trb(ji,jj,1,jpphd) * xthetapd ) ) / & 252 & ( trb(ji,jj,1,jpchn) + trb(ji,jj,1,jpchd ) ) 253 ENDIF 254 END DO 255 END DO 256 n2be_p = 14.01 + ( xmassc * ( ( xthetapn + xthetapd ) / 2.0 ) ) 257 n2be_z = 14.01 + ( xmassc * ( ( xthetazmi + xthetazme ) / 2.0 ) ) 258 n2be_d = 14.01 + ( xmassc * xthetad ) 259 260 WRITE(numout,*) 'DAF: nproc, min/max cchl_p, min/max chl_inc = ', nproc, MINVAL(cchl_p), MAXVAL(cchl_p), MINVAL(chl_inc), MAXVAL(chl_inc) 215 261 216 262 ! Call nitrogen balancing routine 217 CALL bio_analysis( jpi, jpj, jpk, ZDZ(:,:,:), i_tracer, modparm,&218 & n2be_p, n2be_z, n2be_d, assimparm, &219 & INT(aincper), 1, kmt(:,:), tmask(:,:,:),&220 & zmld(:,:), mld_max (:,:), chl_inc(:,:), cchl_p(:,:,1),&221 & nbal_active, phyt_avg (:,:),&222 & gl_active, pgrow_avg (:,:), ploss_avg(:,:),&223 & subsurf_active, deepneg_active, &224 & deeppos_active, nutprof_active, &225 & bstate, outincs, &226 & diag_active, diag, &263 CALL bio_analysis( jpi, jpj, jpk, gdepw_n(:,:,2:jpk), i_tracer, modparm, & 264 & n2be_p, n2be_z, n2be_d, assimparm, & 265 & INT(aincper), 1, INT(SUM(tmask,3)), tmask(:,:,:), & 266 & zmld(:,:), mld_max_bkg(:,:), chl_inc(:,:), cchl_p(:,:), & 267 & nbal_active, phyt_avg_bkg(:,:), & 268 & gl_active, pgrow_avg_bkg(:,:), ploss_avg_bkg(:,:), & 269 & subsurf_active, deepneg_active, & 270 & deeppos_active, nutprof_active, & 271 & bstate, outincs, & 272 & diag_active, diag, & 227 273 & diag_fulldepth_active, diag_fulldepth ) 228 229 ! Save balancing increments 230 logchl_balinc(:,:,:,jp_had_nut) = outincs(:,:,:,i_tracer(1)) 231 logchl_balinc(:,:,:,jp_had_phy) = outincs(:,:,:,i_tracer(2)) 232 logchl_balinc(:,:,:,jp_had_zoo) = outincs(:,:,:,i_tracer(3)) 233 logchl_balinc(:,:,:,jp_had_pdn) = outincs(:,:,:,i_tracer(4)) 234 logchl_balinc(:,:,:,jp_had_dic) = outincs(:,:,:,i_tracer(5)) 235 logchl_balinc(:,:,:,jp_had_alk) = outincs(:,:,:,i_tracer(6)) 274 275 WRITE(numout,*) 'DAF: nproc, min/max outincs(phy)1,20 = ', nproc, MINVAL(outincs(:,:,1,i_tracer(2))), MAXVAL(outincs(:,:,1,i_tracer(2))), MINVAL(outincs(:,:,20,i_tracer(2))), MAXVAL(outincs(:,:,20,i_tracer(2))) 276 277 ! Loop over each grid point partioning the increments 278 logchl_balinc(:,:,:,:) = 0.0 279 DO jk = 1, jpk 280 DO jj = 1, jpj 281 DO ji = 1, jpi 282 283 IF ( ( trb(ji,jj,jk,jpphn) > 0.0 ) .AND. ( trb(ji,jj,jk,jpphd) > 0.0 ) ) THEN 284 ! Phytoplankton nitrogen and silicate split up based on existing ratios 285 zfrac_phn = trb(ji,jj,jk,jpphn) / (trb(ji,jj,jk,jpphn) + trb(ji,jj,jk,jpphd)) 286 zfrac_phd = 1.0 - zfrac_phn 287 zrat_pds_phd = trb(ji,jj,jk,jppds) / trb(ji,jj,jk,jpphd) 288 logchl_balinc(ji,jj,jk,jpphn) = outincs(ji,jj,jk,i_tracer(2)) * zfrac_phn 289 logchl_balinc(ji,jj,jk,jpphd) = outincs(ji,jj,jk,i_tracer(2)) * zfrac_phd 290 logchl_balinc(ji,jj,jk,jppds) = logchl_balinc(ji,jj,jk,jpphd) * zrat_pds_phd 291 292 ! Chlorophyll split up based on existing ratios to phytoplankton nitrogen 293 ! Not using chl_inc directly as it's only 2D 294 ! This method should give same results at surface as splitting chl_inc would 295 zrat_chn_phn = trb(ji,jj,jk,jpchn) / trb(ji,jj,jk,jpphn) 296 zrat_chd_phd = trb(ji,jj,jk,jpchd) / trb(ji,jj,jk,jpphd) 297 logchl_balinc(ji,jj,jk,jpchn) = logchl_balinc(ji,jj,jk,jpphn) * zrat_chn_phn 298 logchl_balinc(ji,jj,jk,jpchd) = logchl_balinc(ji,jj,jk,jpphd) * zrat_chd_phd 299 ENDIF 300 301 IF ( ( trb(ji,jj,jk,jpzmi) > 0.0 ) .AND. ( trb(ji,jj,jk,jpzme) > 0.0 ) ) THEN 302 ! Zooplankton nitrogen split up based on existing ratios 303 zfrac_zmi = trb(ji,jj,jk,jpzmi) / (trb(ji,jj,jk,jpzmi) + trb(ji,jj,jk,jpzme)) 304 zfrac_zme = 1.0 - zfrac_zmi 305 logchl_balinc(ji,jj,jk,jpzmi) = outincs(ji,jj,jk,i_tracer(3)) * zfrac_zmi 306 logchl_balinc(ji,jj,jk,jpzme) = outincs(ji,jj,jk,i_tracer(3)) * zfrac_zme 307 ENDIF 308 309 ! Nitrogen nutrient straight from balancing scheme 310 logchl_balinc(ji,jj,jk,jpdin) = outincs(ji,jj,jk,i_tracer(1)) 311 312 ! Nitrogen detritus straight from balancing scheme 313 logchl_balinc(ji,jj,jk,jpdet) = outincs(ji,jj,jk,i_tracer(4)) 314 315 ! DIC straight from balancing scheme 316 logchl_balinc(ji,jj,jk,jpdic) = outincs(ji,jj,jk,i_tracer(5)) 317 318 ! Alkalinity straight from balancing scheme 319 logchl_balinc(ji,jj,jk,jpalk) = outincs(ji,jj,jk,i_tracer(6)) 320 321 ! Remove diatom silicate increment from nutrient silicate to conserve mass 322 IF ( ( trb(ji,jj,jk,jpsil) - logchl_balinc(ji,jj,jk,jppds) ) > 0.0 ) THEN 323 logchl_balinc(ji,jj,jk,jpsil) = logchl_balinc(ji,jj,jk,jppds) * (-1.0) 324 ENDIF 325 326 IF ( ( trb(ji,jj,jk,jpdet) > 0.0 ) .AND. ( trb(ji,jj,jk,jpdtc) > 0.0 ) ) THEN 327 ! Carbon detritus based on existing ratios 328 zrat_dtc_det = trb(ji,jj,jk,jpdtc) / trb(ji,jj,jk,jpdet) 329 logchl_balinc(ji,jj,jk,jpdtc) = logchl_balinc(ji,jj,jk,jpdet) * zrat_dtc_det 330 ENDIF 331 332 ! Do nothing with iron or oxygen for the time being 333 logchl_balinc(ji,jj,jk,jpfer) = 0.0 334 logchl_balinc(ji,jj,jk,jpoxy) = 0.0 335 336 END DO 337 END DO 338 END DO 236 339 237 340 ELSE ! No nitrogen balancing 238 341 239 ! Initialise phytoplankton increment to zero 240 logchl_balinc(:,:,:,jp_had_phy) = 0.0 342 ! Initialise individual chlorophyll increments to zero 343 logchl_balinc(:,:,:,jpchn) = 0.0 344 logchl_balinc(:,:,:,jpchd) = 0.0 241 345 242 ! Convert surface chlorophyll increment to phytoplankton nitrogen 243 logchl_balinc(:,:,1,jp_had_phy) = ( cchl_p(:,:,1) / (mw_carbon * c2n_p) ) * chl_inc(:,:) 346 ! Split up total surface chlorophyll increments 347 DO jj = 1, jpj 348 DO ji = 1, jpi 349 IF ( medusa_chl(ji,jj) > 0.0 ) THEN 350 zfrac_chn = trb(ji,jj,1,jpchn) / medusa_chl(ji,jj) 351 zfrac_chd = 1.0 - zfrac_chn 352 logchl_balinc(ji,jj,1,jpchn) = chl_inc(ji,jj) * zfrac_chn 353 logchl_balinc(ji,jj,1,jpchd) = chl_inc(ji,jj) * zfrac_chd 354 ENDIF 355 END DO 356 END DO 244 357 245 358 ! Propagate through mixed layer … … 257 370 ! 258 371 DO jk = 2, jkmax 259 logchl_balinc(ji,jj,jk,jp_had_phy) = logchl_balinc(ji,jj,1,jp_had_phy) 372 logchl_balinc(ji,jj,jk,jpchn) = logchl_balinc(ji,jj,1,jpchn) 373 logchl_balinc(ji,jj,jk,jpchd) = logchl_balinc(ji,jj,1,jpchd) 260 374 END DO 261 375 ! … … 264 378 265 379 ! Set other balancing increments to zero 266 logchl_balinc(:,:,:,jp_had_nut) = 0.0 267 logchl_balinc(:,:,:,jp_had_zoo) = 0.0 268 logchl_balinc(:,:,:,jp_had_pdn) = 0.0 269 logchl_balinc(:,:,:,jp_had_dic) = 0.0 270 logchl_balinc(:,:,:,jp_had_alk) = 0.0 271 380 logchl_balinc(:,:,:,jpphn) = 0.0 381 logchl_balinc(:,:,:,jpphd) = 0.0 382 logchl_balinc(:,:,:,jppds) = 0.0 383 logchl_balinc(:,:,:,jpzmi) = 0.0 384 logchl_balinc(:,:,:,jpzme) = 0.0 385 logchl_balinc(:,:,:,jpdin) = 0.0 386 logchl_balinc(:,:,:,jpsil) = 0.0 387 logchl_balinc(:,:,:,jpfer) = 0.0 388 logchl_balinc(:,:,:,jpdet) = 0.0 389 logchl_balinc(:,:,:,jpdtc) = 0.0 390 logchl_balinc(:,:,:,jpdic) = 0.0 391 logchl_balinc(:,:,:,jpalk) = 0.0 392 logchl_balinc(:,:,:,jpoxy) = 0.0 393 272 394 ENDIF 273 395 274 END SUBROUTINE asm_logchl_bal_ hadocc396 END SUBROUTINE asm_logchl_bal_medusa 275 397 276 398 #else … … 279 401 !!---------------------------------------------------------------------- 280 402 CONTAINS 281 SUBROUTINE asm_logchl_bal_ hadocc( logchl_bkginc, aincper, mld_choice_bgc, &403 SUBROUTINE asm_logchl_bal_medusa( logchl_bkginc, aincper, mld_choice_bgc, & 282 404 & k_maxchlinc, logchl_balinc ) 283 405 REAL :: logchl_bkginc(:,:) … … 286 408 REAL :: k_maxchlinc 287 409 REAL( :: logchl_balinc(:,:,:,:) 288 WRITE(*,*) 'asm_logchl_bal_ hadocc: You should not have seen this print! error?'289 END SUBROUTINE asm_logchl_bal_ hadocc410 WRITE(*,*) 'asm_logchl_bal_medusa: You should not have seen this print! error?' 411 END SUBROUTINE asm_logchl_bal_medusa 290 412 #endif 291 413 292 414 !!====================================================================== 293 END MODULE asmlogchlbal_ hadocc415 END MODULE asmlogchlbal_medusa
Note: See TracChangeset
for help on using the changeset viewer.