Changeset 5841 for branches/NERC/dev_r5518_NOC_MEDUSA_Stable/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcbio_medusa.F90
- Timestamp:
- 2015-10-30T12:48:06+01:00 (8 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/NERC/dev_r5518_NOC_MEDUSA_Stable/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcbio_medusa.F90
r5726 r5841 16 16 !! - ! 2013-05 (A. Yool) updated for v3.5 17 17 !! - ! 2014-08 (A. Yool, J. Palm) Add DMS module for UKESM1 model 18 !! - ! 2015-06 (A. Yool) Update to include MOCSY 19 !! - ! 2015-07 (A. Yool) Update for rolling averages 18 20 !!---------------------------------------------------------------------- 19 21 !! … … 36 38 #endif 37 39 !! 40 #if defined key_mocsy 41 !!---------------------------------------------------------------------- 42 !! Updates with the addition of MOCSY include: 43 !! - option to use PML or MOCSY carbonate chemistry (the latter is 44 !! preferred) 45 !! - central calculation of gas transfer velocity, f_kw660; previously 46 !! this was done separately for CO2 and O2 with predictable results 47 !! - distribution of f_kw660 to both PML and MOCSY CO2 air-sea flux 48 !! calculations and to those for O2 air-sea flux 49 !! - extra diagnostics included for MOCSY 50 !!---------------------------------------------------------------------- 51 #endif 52 !! 38 53 #if defined key_medusa 39 54 !!---------------------------------------------------------------------- … … 55 70 # endif 56 71 # if defined key_roam 72 USE gastransfer 73 # if defined key_mocsy 74 USE mocsy_wrapper 75 # else 57 76 USE trcco2_medusa 77 # endif 58 78 USE trcoxy_medusa 59 79 !! Jpalm (08/08/2014) … … 134 154 REAL(wp) :: ztmp, zsal 135 155 # endif 156 # if defined key_mocsy 157 REAL(wp) :: zpho 158 # endif 136 159 !! 137 160 !! integrated source and sink terms … … 266 289 REAL(wp) :: f_kw660, f_o2flux, f_o2sat 267 290 REAL(wp), DIMENSION(jpi,jpj) :: f_omcal, f_omarg 291 !! 292 !! AXY (23/06/15): additional diagnostics for MOCSY and oxygen 293 REAL(wp) :: f_fco2w, f_BetaD, f_rhosw, f_opres, f_insitut, f_pco2atm, f_fco2atm 294 REAL(wp) :: f_schmidtco2, f_kwco2, f_K0, f_co2starair, f_dpco2, f_kwo2 295 !! 268 296 INTEGER :: iters 269 297 REAL(wp) :: f_year … … 305 333 306 334 !!--------------------------------------------------------------------- 335 336 # if defined key_debug_medusa 337 IF (lwp) write (numout,*) 'trc_bio_medusa: variables defined' 338 CALL flush(numout) 339 # endif 307 340 308 341 !! AXY (20/11/14): alter this to report on first MEDUSA call … … 486 519 # endif 487 520 521 # if defined key_debug_medusa 522 IF (lwp) write (numout,*) 'trc_bio_medusa: variables initialised and checked' 523 CALL flush(numout) 524 # endif 525 488 526 # if defined key_roam 489 527 !!---------------------------------------------------------------------- … … 514 552 f_pco2a = fq4 515 553 endif 516 # if defined key_axy_pi_co2554 # if defined key_axy_pi_co2 517 555 f_pco2a = hist_pco2(1) 518 IF(lwp) WRITE(numout,*) ' MEDUSA atm pCO2 = FIXED' 519 # endif 556 # endif 520 557 !! IF(lwp) WRITE(numout,*) ' MEDUSA nyear =', nyear 521 558 !! IF(lwp) WRITE(numout,*) ' MEDUSA nsec_day =', real(nsec_day) … … 529 566 # endif 530 567 568 # if defined key_debug_medusa 569 IF (lwp) write (numout,*) 'trc_bio_medusa: ready for carbonate chemistry' 570 IF (lwp) write (numout,*) 'trc_bio_medusa: kt = ', kt 571 IF (lwp) write (numout,*) 'trc_bio_medusa: nittrc000 = ', nittrc000 572 CALL flush(numout) 573 # endif 574 531 575 # if defined key_roam 532 !! AXY (20/11/14): alter to call on first MEDUSA timestep 576 !! AXY (20/11/14): alter to call on first MEDUSA timestep and then every 577 !! month (this is hardwired as 960 timesteps but should 578 !! be calculated and done properly 533 579 !! IF( kt == nit000 .or. mod(kt,1920) == 0 ) THEN 534 IF( kt == nittrc000 .or. mod(kt, 1920) == 0 ) THEN580 IF( kt == nittrc000 .or. mod(kt,960) == 0 ) THEN 535 581 !!---------------------------------------------------------------------- 536 582 !! Calculate the carbonate chemistry for the whole ocean on the first … … 540 586 !! 541 587 IF(lwp) WRITE(numout,*) ' MEDUSA calculating all carbonate chemistry at kt =', kt 588 CALL flush(numout) 542 589 !! blank flags 543 590 i2_omcal(:,:) = 0 … … 549 596 !! OPEN wet point IF..THEN loop 550 597 if (tmask(ji,jj,jk).eq.1) then 551 !! carbonate chemistry 598 !! do carbonate chemistry 599 !! 552 600 fdep2 = fsdept(ji,jj,jk) !! set up level midpoint 601 !! 553 602 !! set up required state variables 554 603 zdic = max(0.,trn(ji,jj,jk,jpdic)) !! dissolved inorganic carbon … … 556 605 ztmp = tsn(ji,jj,jk,jp_tem) !! temperature 557 606 zsal = tsn(ji,jj,jk,jp_sal) !! salinity 607 # if defined key_mocsy 608 zsil = max(0.,trn(ji,jj,jk,jpsil)) !! silicic acid 609 zpho = max(0.,trn(ji,jj,jk,jpdin)) / 16.0 !! phosphate via DIN and Redfield 610 # endif 558 611 !! 559 612 !! AXY (28/02/14): check input fields … … 571 624 ji, ',', jj, ',', jk, ') at time', kt 572 625 endif 626 !! 627 !! blank input variables not used at this stage (they relate to air-sea flux) 628 f_kw660 = 1.0 629 f_pp0 = 1.0 630 !! 573 631 !! calculate carbonate chemistry at grid cell midpoint 574 CALL trc_co2_medusa( ztmp, zsal, zdic, zalk, fdep2, 5.0, f_pco2a, & ! inputs 575 f_ph, f_pco2w, f_h2co3, f_hco3, f_co3, f_omcal(ji,jj), & ! outputs 632 # if defined key_mocsy 633 !! AXY (22/06/15): use Orr & Epitalon (2015) MOCSY-2 carbonate 634 !! chemistry package 635 CALL mocsy_interface( ztmp, zsal, zalk, zdic, zsil, zpho, & ! inputs 636 f_pp0, fdep2, flatx, f_kw660, f_pco2a, 1, & ! inputs 637 f_ph, f_pco2w, f_fco2w, f_h2co3, f_hco3, f_co3, f_omarg(ji,jj), & ! outputs 638 f_omcal(ji,jj), f_BetaD, f_rhosw, f_opres, f_insitut, & ! outputs 639 f_pco2atm, f_fco2atm, f_schmidtco2, f_kwco2, f_K0, & ! outputs 640 f_co2starair, f_co2flux, f_dpco2 ) ! outputs 641 !! 642 f_TDIC = (zdic / f_rhosw) * 1000. ! mmol / m3 -> umol / kg 643 f_TALK = (zalk / f_rhosw) * 1000. ! meq / m3 -> ueq / kg 644 f_dcf = f_rhosw 645 # else 646 !! AXY (22/06/15): use old PML carbonate chemistry package (the 647 !! MEDUSA-2 default) 648 CALL trc_co2_medusa( ztmp, zsal, zdic, zalk, fdep2, f_kw660, & ! inputs 649 f_pco2a, f_ph, f_pco2w, f_h2co3, f_hco3, f_co3, f_omcal(ji,jj), & ! outputs 576 650 f_omarg(ji,jj), f_co2flux, f_TDIC, f_TALK, f_dcf, f_henry, iters) ! outputs 577 651 !! … … 581 655 iters, ' AT (', ji, ', ', jj, ', ', jk, ') AT ', kt 582 656 endif 657 # endif 658 !! 583 659 !! store 3D outputs 584 660 f3_pH(ji,jj,jk) = f_ph … … 588 664 f3_omcal(ji,jj,jk) = f_omcal(ji,jj) 589 665 f3_omarg(ji,jj,jk) = f_omarg(ji,jj) 666 !! 590 667 !! CCD calculation: calcite 591 668 if (i2_omcal(ji,jj) .eq. 0 .and. f_omcal(ji,jj) .lt. 1.0) then … … 607 684 i2_omcal(ji,jj) = 1 608 685 endif 686 !! 609 687 !! CCD calculation: aragonite 610 688 if (i2_omarg(ji,jj) .eq. 0 .and. f_omarg(ji,jj) .lt. 1.0) then … … 633 711 # endif 634 712 713 # if defined key_debug_medusa 714 IF (lwp) write (numout,*) 'trc_bio_medusa: ready for full domain calculations' 715 CALL flush(numout) 716 # endif 717 635 718 !!---------------------------------------------------------------------- 636 719 !! MEDUSA has unified equation through the water column … … 656 739 !! OPEN wet point IF..THEN loop 657 740 if (tmask(ji,jj,jk).eq.1) then 658 659 741 !!====================================================================== 660 742 !! SETUP LOCAL GRID CELL … … 814 896 !!---------------------------------------------------------------------- 815 897 !! 816 !! a bit of set up ...817 ! f_uwind = zwnd_i(ji,jj)818 ! f_vwind = zwnd_j(ji,jj)819 898 !! AXY (17/07/14): zwind_i and zwind_j do not exist in this 820 899 !! version of NEMO because it does not include … … 827 906 !! revisited when MEDUSA properly interacts 828 907 !! with UKESM1 physics 829 ! f_uwind = zwind_i(ji,jj) 830 ! f_vwind = zwind_j(ji,jj) 831 ! f_wind = ((f_uwind**2.0) + (f_vwind**2.0))**0.5 908 !! 832 909 f_wind = wndm(ji,jj) 833 !! AXY (17/07/14): the current oxygen code takes in separate 834 !! U and V components of the wind; to avoid 835 !! the need to change this, calculate these 836 !! components based on wndm; again, this is 837 !! not ideal, but should suffice as a 838 !! temporary measure; in the long-term all 839 !! of MEDUSA's air-sea gas exchange terms 840 !! will be revisited to ensure that they are 841 !! valid and self-consistent; and the CO2 842 !! code will be wholly replaced with a more 843 !! up-to-date parameterisation 844 f_uwind = ((f_wind**2.0) / 2.0)**0.5 845 f_vwind = f_uwind 910 !! 911 !! AXY (23/06/15): as part of an effort to update the carbonate chemistry 912 !! in MEDUSA, the gas transfer velocity used in the carbon 913 !! and oxygen cycles has been harmonised and is calculated 914 !! by the same function here; this harmonisation includes 915 !! changes to the PML carbonate chemistry scheme so that 916 !! it too makes use of the same gas transfer velocity; the 917 !! preferred parameterisation of this is Wanninkhof (2014), 918 !! option 7 919 !! 920 # if defined key_debug_medusa 921 IF (lwp) write (numout,*) 'trc_bio_medusa: entering gas_transfer' 922 CALL flush(numout) 923 # endif 924 CALL gas_transfer( f_wind, 1, 7, & ! inputs 925 f_kw660 ) ! outputs 926 # if defined key_debug_medusa 927 IF (lwp) write (numout,*) 'trc_bio_medusa: exiting gas_transfer' 928 CALL flush(numout) 929 # endif 930 !! 931 !! air pressure (atm); ultimately this will use air pressure at the base 932 !! of the UKESM1 atmosphere 933 !! 846 934 f_pp0 = 1.0 935 !! 847 936 !! IF(lwp) WRITE(numout,*) ' MEDUSA ztmp =', ztmp 848 937 !! IF(lwp) WRITE(numout,*) ' MEDUSA zwind_i =', zwind_i(ji,jj) … … 852 941 !! 853 942 # if defined key_axy_carbchem 943 # if defined key_mocsy 944 !! 945 !! AXY (22/06/15): use Orr & Epitalon (2015) MOCSY-2 carbonate 946 !! chemistry package; note that depth is set to 947 !! zero in this call 948 CALL mocsy_interface( ztmp, zsal, zalk, zdic, zsil, zpho, & ! inputs 949 f_pp0, 0.0, flatx, f_kw660, f_pco2a, 1, & ! inputs 950 f_ph, f_pco2w, f_fco2w, f_h2co3, f_hco3, f_co3, f_omarg(ji,jj), & ! outputs 951 f_omcal(ji,jj), f_BetaD, f_rhosw, f_opres, f_insitut, & ! outputs 952 f_pco2atm, f_fco2atm, f_schmidtco2, f_kwco2, f_K0, & ! outputs 953 f_co2starair, f_co2flux, f_dpco2 ) ! outputs 954 !! 955 f_TDIC = (zdic / f_rhosw) * 1000. ! mmol / m3 -> umol / kg 956 f_TALK = (zalk / f_rhosw) * 1000. ! meq / m3 -> ueq / kg 957 f_dcf = f_rhosw 958 # else 854 959 iters = 0 855 960 !! 856 961 !! carbon dioxide (CO2); Jerry Blackford code (ostensibly OCMIP-2, but not) 857 CALL trc_co2_medusa( ztmp, zsal, zdic, zalk, 0.0, f_ wind, f_pco2a, &! inputs858 f_ph, f_pco2w, f_h2co3, f_hco3, f_co3, f_omcal(ji,jj), &! outputs859 f_omarg(ji,jj), f_co2flux, f_TDIC, f_TALK, f_dcf, f_henry, iters ) 962 CALL trc_co2_medusa( ztmp, zsal, zdic, zalk, 0.0, f_kw660, f_pco2a, & ! inputs 963 f_ph, f_pco2w, f_h2co3, f_hco3, f_co3, f_omcal(ji,jj), & ! outputs 964 f_omarg(ji,jj), f_co2flux, f_TDIC, f_TALK, f_dcf, f_henry, iters ) ! outputs 860 965 !! 861 966 !! AXY (09/01/14): removed iteration and NaN checks; these have … … 869 974 iters, ' AT (', ji, ', ', jj, ', ', jk, ') AT ', kt 870 975 endif 976 # endif 871 977 # else 872 978 !! AXY (18/04/13): switch off carbonate chemistry calculations; provide … … 882 988 f_TDIC = zdic 883 989 f_TALK = zalk 884 f_dcf = 1. 990 f_dcf = 1.026 885 991 f_henry = 1. 992 !! AXY (23/06/15): add in some extra MOCSY diagnostics 993 f_fco2w = fpco2a 994 f_BetaD = 1. 995 f_rhosw = 1.026 996 f_opres = 0. 997 f_insitut = ztmp 998 f_pco2atm = fpco2a 999 f_fco2atm = fpco2a 1000 f_schmidtco2 = 660. 1001 f_kwco2 = 0. 1002 f_K0 = 0. 1003 f_co2starair = fpco2a 1004 f_dpco2 = 0. 886 1005 # endif 887 1006 !! 888 !! already in right units; correct for sea-ice; divide through by layer thickness889 f_co2flux = (1. - fr_i(ji,jj)) * f_co2flux / fthk1007 !! mmol/m2/s -> mmol/m3/d; correct for sea-ice; divide through by layer thickness 1008 f_co2flux = (1. - fr_i(ji,jj)) * f_co2flux * 86400. / fthk 890 1009 !! 891 1010 !! oxygen (O2); OCMIP-2 code 892 CALL trc_oxy_medusa( ztmp, zsal, f_uwind, f_vwind, f_pp0, zoxy / 1000., fthk, & ! inputs 893 f_kw660, f_o2flux, f_o2sat ) ! outputs 894 !! 895 !! mol/m3/s -> mmol/m3/d; correct for sea-ice 896 f_o2flux = (1. - fr_i(ji,jj)) * f_o2flux * 1000. * 60. * 60. * 24. 897 f_o2sat = f_o2sat * 1000. 1011 !! AXY (23/06/15): amend input list for oxygen to account for common gas 1012 !! transfer velocity 1013 !! CALL trc_oxy_medusa( ztmp, zsal, f_uwind, f_vwind, f_pp0, zoxy / 1000., fthk, & ! inputs 1014 !! f_kw660, f_o2flux, f_o2sat ) ! outputs 1015 CALL trc_oxy_medusa( ztmp, zsal, f_kw660, f_pp0, zoxy, & ! inputs 1016 f_kwo2, f_o2flux, f_o2sat ) ! outputs 1017 !! 1018 !! mmol/m2/s -> mol/m3/d; correct for sea-ice; divide through by layer thickness 1019 f_o2flux = (1. - fr_i(ji,jj)) * f_o2flux * 86400. / fthk 898 1020 !! 899 1021 !! Jpalm (08-2014) … … 911 1033 !! 912 1034 IF (jdms .eq. 1) THEN 913 !! 914 !! CALL trc_dms_medusa( zchn, zchd, hmld(ji,jj), & !! inputs 915 !! dms_surf ) !! outputs 916 CALL trc_dms_medusa( zchn, zchd, hmld(ji,jj), qsr(ji,jj), zdin, & !! inputs 917 dms_surf, dms_andr, dms_simo, dms_aran, dms_hall ) !! outputs 1035 !! 1036 !! feed in correct inputs 1037 if (jdms_input .eq. 0) then 1038 !! use instantaneous inputs 1039 CALL trc_dms_medusa( zchn, zchd, hmld(ji,jj), qsr(ji,jj), zdin, & ! inputs 1040 dms_andr, dms_simo, dms_aran, dms_hall ) ! outputs 1041 else 1042 !! use diel-average inputs 1043 CALL trc_dms_medusa( zn_dms_chn(ji,jj), zn_dms_chd(ji,jj), & ! inputs 1044 zn_dms_mld(ji,jj), zn_dms_qsr(ji,jj), zn_dms_din(ji,jj), & ! inputs 1045 dms_andr, dms_simo, dms_aran, dms_hall ) ! outputs 1046 endif 1047 !! 1048 !! assign correct output to variable passed to atmosphere 1049 if (jdms_model .eq. 1) then 1050 dms_surf = dms_andr 1051 elseif (jdms_model .eq. 2) then 1052 dms_surf = dms_simo 1053 elseif (jdms_model .eq. 3) then 1054 dms_surf = dms_aran 1055 elseif (jdms_model .eq. 4) then 1056 dms_surf = dms_hall 1057 endif 918 1058 ENDIF 919 1059 !! End DMS Loop … … 953 1093 !! alkalinity are derived from continent-scale DIC estimates (Huang et al., 954 1094 !! 2012) and some Arctic river alkalinity estimates (Katya?) 955 !! 1095 !! 956 1096 !! as of 19/07/12, riverine nutrients can now be spread vertically across 957 1097 !! several grid cells rather than just poured into the surface box; this … … 3317 3457 ENDDO 3318 3458 ENDDO 3459 !! AXY (07/07/15): temporary hijacking 3460 # if defined key_roam 3461 trc2d(:,:,126) = zn_dms_chn(:,:) 3462 trc2d(:,:,127) = zn_dms_chd(:,:) 3463 trc2d(:,:,128) = zn_dms_mld(:,:) 3464 trc2d(:,:,129) = zn_dms_qsr(:,:) 3465 trc2d(:,:,130) = zn_dms_din(:,:) 3466 # endif 3319 3467 3320 3468 if (ibenthic.eq.2) then … … 3322 3470 !! effectively commented out because it does not work as 3323 3471 !! anticipated; it can be deleted at a later date 3324 if (jorgben.eq.1) then 3325 za_sed_n(:,:) = ( f_sbenin_n(:,:) + f_fbenin_n(:,:) - f_benout_n(:,:) ) * rdt 3326 za_sed_fe(:,:) = ( f_sbenin_fe(:,:) + f_fbenin_fe(:,:) - f_benout_fe(:,:) ) * rdt 3327 za_sed_c(:,:) = ( f_sbenin_c(:,:) + f_fbenin_c(:,:) - f_benout_c(:,:) ) * rdt 3328 endif 3329 if (jinorgben.eq.1) then 3330 za_sed_si(:,:) = ( f_fbenin_si(:,:) - f_benout_si(:,:) ) * rdt 3331 za_sed_ca(:,:) = ( f_fbenin_ca(:,:) - f_benout_ca(:,:) ) * rdt 3332 endif 3333 !! 3334 !! Leap-frog scheme - only in explicit case, otherwise the time stepping 3335 !! is already being done in trczdf 3336 !! IF( l_trczdf_exp .AND. (ln_trcadv_cen2 .OR. ln_trcadv_tvd) ) THEN 3337 !! zfact = 2. * rdttra(jk) * FLOAT( ndttrc ) 3338 !! IF( neuler == 0 .AND. kt == nittrc000 ) zfact = rdttra(jk) * FLOAT(ndttrc) 3339 !! if (jorgben.eq.1) then 3340 !! za_sed_n(:,:) = zb_sed_n(:,:) + ( zfact * za_sed_n(:,:) ) 3341 !! za_sed_fe(:,:) = zb_sed_fe(:,:) + ( zfact * za_sed_fe(:,:) ) 3342 !! za_sed_c(:,:) = zb_sed_c(:,:) + ( zfact * za_sed_c(:,:) ) 3343 !! endif 3344 !! if (jinorgben.eq.1) then 3345 !! za_sed_si(:,:) = zb_sed_si(:,:) + ( zfact * za_sed_si(:,:) ) 3346 !! za_sed_ca(:,:) = zb_sed_ca(:,:) + ( zfact * za_sed_ca(:,:) ) 3347 !! endif 3348 !! ENDIF 3349 !! 3350 !! Time filter and swap of arrays 3351 IF( ln_trcadv_cen2 .OR. ln_trcadv_tvd ) THEN ! centred or tvd scheme 3352 IF( neuler == 0 .AND. kt == nittrc000 ) THEN 3472 if (jorgben.eq.1) then 3473 za_sed_n(:,:) = ( f_sbenin_n(:,:) + f_fbenin_n(:,:) - f_benout_n(:,:) ) * rdt 3474 za_sed_fe(:,:) = ( f_sbenin_fe(:,:) + f_fbenin_fe(:,:) - f_benout_fe(:,:) ) * rdt 3475 za_sed_c(:,:) = ( f_sbenin_c(:,:) + f_fbenin_c(:,:) - f_benout_c(:,:) ) * rdt 3476 endif 3477 if (jinorgben.eq.1) then 3478 za_sed_si(:,:) = ( f_fbenin_si(:,:) - f_benout_si(:,:) ) * rdt 3479 za_sed_ca(:,:) = ( f_fbenin_ca(:,:) - f_benout_ca(:,:) ) * rdt 3480 endif 3481 !! 3482 !! Leap-frog scheme - only in explicit case, otherwise the time stepping 3483 !! is already being done in trczdf 3484 !! IF( l_trczdf_exp .AND. (ln_trcadv_cen2 .OR. ln_trcadv_tvd) ) THEN 3485 !! zfact = 2. * rdttra(jk) * FLOAT( ndttrc ) 3486 !! IF( neuler == 0 .AND. kt == nittrc000 ) zfact = rdttra(jk) * FLOAT(ndttrc) 3487 !! if (jorgben.eq.1) then 3488 !! za_sed_n(:,:) = zb_sed_n(:,:) + ( zfact * za_sed_n(:,:) ) 3489 !! za_sed_fe(:,:) = zb_sed_fe(:,:) + ( zfact * za_sed_fe(:,:) ) 3490 !! za_sed_c(:,:) = zb_sed_c(:,:) + ( zfact * za_sed_c(:,:) ) 3491 !! endif 3492 !! if (jinorgben.eq.1) then 3493 !! za_sed_si(:,:) = zb_sed_si(:,:) + ( zfact * za_sed_si(:,:) ) 3494 !! za_sed_ca(:,:) = zb_sed_ca(:,:) + ( zfact * za_sed_ca(:,:) ) 3495 !! endif 3496 !! ENDIF 3497 !! 3498 !! Time filter and swap of arrays 3499 IF( ln_trcadv_cen2 .OR. ln_trcadv_tvd ) THEN ! centred or tvd scheme 3500 IF( neuler == 0 .AND. kt == nittrc000 ) THEN 3501 if (jorgben.eq.1) then 3502 zb_sed_n(:,:) = zn_sed_n(:,:) 3503 zn_sed_n(:,:) = za_sed_n(:,:) 3504 za_sed_n(:,:) = 0.0 3505 !! 3506 zb_sed_fe(:,:) = zn_sed_fe(:,:) 3507 zn_sed_fe(:,:) = za_sed_fe(:,:) 3508 za_sed_fe(:,:) = 0.0 3509 !! 3510 zb_sed_c(:,:) = zn_sed_c(:,:) 3511 zn_sed_c(:,:) = za_sed_c(:,:) 3512 za_sed_c(:,:) = 0.0 3513 endif 3514 if (jinorgben.eq.1) then 3515 zb_sed_si(:,:) = zn_sed_si(:,:) 3516 zn_sed_si(:,:) = za_sed_si(:,:) 3517 za_sed_si(:,:) = 0.0 3518 !! 3519 zb_sed_ca(:,:) = zn_sed_ca(:,:) 3520 zn_sed_ca(:,:) = za_sed_ca(:,:) 3521 za_sed_ca(:,:) = 0.0 3522 endif 3523 ELSE 3524 if (jorgben.eq.1) then 3525 zb_sed_n(:,:) = (atfp * ( zb_sed_n(:,:) + za_sed_n(:,:) )) + (atfp1 * zn_sed_n(:,:) ) 3526 zn_sed_n(:,:) = za_sed_n(:,:) 3527 za_sed_n(:,:) = 0.0 3528 !! 3529 zb_sed_fe(:,:) = (atfp * ( zb_sed_fe(:,:) + za_sed_fe(:,:) )) + (atfp1 * zn_sed_fe(:,:)) 3530 zn_sed_fe(:,:) = za_sed_fe(:,:) 3531 za_sed_fe(:,:) = 0.0 3532 !! 3533 zb_sed_c(:,:) = (atfp * ( zb_sed_c(:,:) + za_sed_c(:,:) )) + (atfp1 * zn_sed_c(:,:) ) 3534 zn_sed_c(:,:) = za_sed_c(:,:) 3535 za_sed_c(:,:) = 0.0 3536 endif 3537 if (jinorgben.eq.1) then 3538 zb_sed_si(:,:) = (atfp * ( zb_sed_si(:,:) + za_sed_si(:,:) )) + (atfp1 * zn_sed_si(:,:)) 3539 zn_sed_si(:,:) = za_sed_si(:,:) 3540 za_sed_si(:,:) = 0.0 3541 !! 3542 zb_sed_ca(:,:) = (atfp * ( zb_sed_ca(:,:) + za_sed_ca(:,:) )) + (atfp1 * zn_sed_ca(:,:)) 3543 zn_sed_ca(:,:) = za_sed_ca(:,:) 3544 za_sed_ca(:,:) = 0.0 3545 endif 3546 ENDIF 3547 ELSE ! case of smolar scheme or muscl 3353 3548 if (jorgben.eq.1) then 3354 zb_sed_n(:,:) = z n_sed_n(:,:)3549 zb_sed_n(:,:) = za_sed_n(:,:) 3355 3550 zn_sed_n(:,:) = za_sed_n(:,:) 3356 3551 za_sed_n(:,:) = 0.0 3357 3552 !! 3358 zb_sed_fe(:,:) = z n_sed_fe(:,:)3553 zb_sed_fe(:,:) = za_sed_fe(:,:) 3359 3554 zn_sed_fe(:,:) = za_sed_fe(:,:) 3360 3555 za_sed_fe(:,:) = 0.0 3361 3556 !! 3362 zb_sed_c(:,:) = z n_sed_c(:,:)3557 zb_sed_c(:,:) = za_sed_c(:,:) 3363 3558 zn_sed_c(:,:) = za_sed_c(:,:) 3364 3559 za_sed_c(:,:) = 0.0 3365 3560 endif 3366 3561 if (jinorgben.eq.1) then 3367 zb_sed_si(:,:) = z n_sed_si(:,:)3562 zb_sed_si(:,:) = za_sed_si(:,:) 3368 3563 zn_sed_si(:,:) = za_sed_si(:,:) 3369 3564 za_sed_si(:,:) = 0.0 3370 3565 !! 3371 zb_sed_ca(:,:) = zn_sed_ca(:,:) 3372 zn_sed_ca(:,:) = za_sed_ca(:,:) 3373 za_sed_ca(:,:) = 0.0 3374 endif 3375 ELSE 3376 if (jorgben.eq.1) then 3377 zb_sed_n(:,:) = (atfp * ( zb_sed_n(:,:) + za_sed_n(:,:) )) + (atfp1 * zn_sed_n(:,:) ) 3378 zn_sed_n(:,:) = za_sed_n(:,:) 3379 za_sed_n(:,:) = 0.0 3380 !! 3381 zb_sed_fe(:,:) = (atfp * ( zb_sed_fe(:,:) + za_sed_fe(:,:) )) + (atfp1 * zn_sed_fe(:,:)) 3382 zn_sed_fe(:,:) = za_sed_fe(:,:) 3383 za_sed_fe(:,:) = 0.0 3384 !! 3385 zb_sed_c(:,:) = (atfp * ( zb_sed_c(:,:) + za_sed_c(:,:) )) + (atfp1 * zn_sed_c(:,:) ) 3386 zn_sed_c(:,:) = za_sed_c(:,:) 3387 za_sed_c(:,:) = 0.0 3388 endif 3389 if (jinorgben.eq.1) then 3390 zb_sed_si(:,:) = (atfp * ( zb_sed_si(:,:) + za_sed_si(:,:) )) + (atfp1 * zn_sed_si(:,:)) 3391 zn_sed_si(:,:) = za_sed_si(:,:) 3392 za_sed_si(:,:) = 0.0 3393 !! 3394 zb_sed_ca(:,:) = (atfp * ( zb_sed_ca(:,:) + za_sed_ca(:,:) )) + (atfp1 * zn_sed_ca(:,:)) 3566 zb_sed_ca(:,:) = za_sed_ca(:,:) 3395 3567 zn_sed_ca(:,:) = za_sed_ca(:,:) 3396 3568 za_sed_ca(:,:) = 0.0 3397 3569 endif 3398 3570 ENDIF 3399 ELSE ! case of smolar scheme or muscl3400 if (jorgben.eq.1) then3401 zb_sed_n(:,:) = za_sed_n(:,:)3402 zn_sed_n(:,:) = za_sed_n(:,:)3403 za_sed_n(:,:) = 0.03404 !!3405 zb_sed_fe(:,:) = za_sed_fe(:,:)3406 zn_sed_fe(:,:) = za_sed_fe(:,:)3407 za_sed_fe(:,:) = 0.03408 !!3409 zb_sed_c(:,:) = za_sed_c(:,:)3410 zn_sed_c(:,:) = za_sed_c(:,:)3411 za_sed_c(:,:) = 0.03412 endif3413 if (jinorgben.eq.1) then3414 zb_sed_si(:,:) = za_sed_si(:,:)3415 zn_sed_si(:,:) = za_sed_si(:,:)3416 za_sed_si(:,:) = 0.03417 !!3418 zb_sed_ca(:,:) = za_sed_ca(:,:)3419 zn_sed_ca(:,:) = za_sed_ca(:,:)3420 za_sed_ca(:,:) = 0.03421 endif3422 ENDIF3423 3571 endif 3424 3572 3425 3573 IF( ln_diatrc ) THEN 3426 3574 !!----------------------------------------------------------------------
Note: See TracChangeset
for help on using the changeset viewer.