New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 7958 for branches/UKMO/dev_r5518_medusa_chg_trc_bio_medusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcbio_medusa.F90 – NEMO

Ignore:
Timestamp:
2017-04-24T13:30:28+02:00 (7 years ago)
Author:
marc
Message:

Pulled air-sea gas exchange and river inputs from trcbio_medusa.F90 into air_sea.F90

File:
1 edited

Legend:

Unmodified
Added
Removed

  • branches/UKMO/dev_r5518_medusa_chg_trc_bio_medusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcbio_medusa.F90

    r7938 r7958  
    9999      USE bio_medusa_init_mod,        ONLY: bio_medusa_init 
    100100      USE carb_chem_mod,              ONLY: carb_chem 
     101      USE air_sea_mod,                ONLY: air_sea 
    101102      USE bio_medusa_diag_slice_mod,  ONLY: bio_medusa_diag_slice 
    102103      USE bio_medusa_fin_mod,         ONLY: bio_medusa_fin 
    … …  
    304305      !! 
    305306      !! flags to help with calculating the position of the CCD 
    306       INTEGER, DIMENSION(jpi,jpj) ::     i2_omcal,i2_omarg 
     307! Moved into carb_chem.F90 - marc 20/4/17 
     308!      INTEGER, DIMENSION(jpi,jpj) ::     i2_omcal,i2_omarg 
    307309      !! 
    308310      !! AXY (24/11/16): add xCO2 variable for atmosphere (what we actually have) 
    … …  
    339341      !! Jpalm (11-08-2014) 
    340342      !! add DMS in MEDUSA for UKESM1 model 
    341       REAL(wp), DIMENSION(jpi,jpj) ::    dms_surf 
     343!      REAL(wp), DIMENSION(jpi,jpj) ::    dms_surf 
    342344      !! AXY (13/03/15): add in other DMS calculations 
    343       REAL(wp), DIMENSION(jpi,jpj) ::    dms_andr, dms_simo, dms_aran, dms_hall 
     345!      REAL(wp), DIMENSION(jpi,jpj) ::    dms_andr, dms_simo, dms_aran, dms_hall 
    344346 
    345347# endif 
    … …  
    578580         DO ji = 2,jpim1 
    579581            !! OPEN wet point IF..THEN loop 
    580             if (tmask(ji,jj,jk).eq.1) then                
    581                !!====================================================================== 
     582            if (tmask(ji,jj,jk) == 1) then                
     583               !!=========================================================== 
    582584               !! SETUP LOCAL GRID CELL 
    583                !!====================================================================== 
    584                !! 
    585                !!--------------------------------------------------------------------- 
     585               !!=========================================================== 
     586               !! 
     587               !!----------------------------------------------------------- 
    586588               !! Some notes on grid vertical structure 
    587                !! - fsdepw(ji,jj,jk) is the depth of the upper surface of level jk 
    588                !! - fsde3w(ji,jj,jk) is *approximately* the midpoint of level jk 
     589               !! - fsdepw(ji,jj,jk) is the depth of the upper surface of  
     590               !!   level jk 
     591               !! - fsde3w(ji,jj,jk) is *approximately* the midpoint of  
     592               !!   level jk 
    589593               !! - fse3t(ji,jj,jk)  is the thickness of level jk 
    590                !!--------------------------------------------------------------------- 
     594               !!----------------------------------------------------------- 
    591595               !! 
    592596               !! AXY (01/03/10): set up level depth (bottom of level) 
    593597               fdep1(ji,jj) = fsdepw(ji,jj,jk) + fse3t(ji,jj,jk) 
    594598               !! AXY (28/11/16): local seafloor depth 
    595                !!                 previously mbathy(ji,jj) - 1, now mbathy(ji,jj) 
     599               !!                 previously mbathy(ji,jj) - 1, now  
     600               !!                 mbathy(ji,jj) 
    596601               mbathy(ji,jj) = mbathy(ji,jj) 
    597602               !! 
    … …  
    599604               !! negative values of state variables are not allowed to 
    600605               !! contribute to the calculated fluxes 
    601                zchn(ji,jj) = max(0.,trn(ji,jj,jk,jpchn)) !! non-diatom chlorophyll 
    602                zchd(ji,jj) = max(0.,trn(ji,jj,jk,jpchd)) !! diatom chlorophyll 
    603                zphn(ji,jj) = max(0.,trn(ji,jj,jk,jpphn)) !! non-diatoms 
    604                zphd(ji,jj) = max(0.,trn(ji,jj,jk,jpphd)) !! diatoms 
    605                zpds(ji,jj) = max(0.,trn(ji,jj,jk,jppds)) !! diatom silicon 
    606                !! AXY (28/01/10): probably need to take account of chl/biomass connection 
     606               !! non-diatom chlorophyll 
     607               zchn(ji,jj) = max(0.,trn(ji,jj,jk,jpchn)) 
     608               !! diatom chlorophyll 
     609               zchd(ji,jj) = max(0.,trn(ji,jj,jk,jpchd)) 
     610               !! non-diatoms 
     611               zphn(ji,jj) = max(0.,trn(ji,jj,jk,jpphn)) 
     612               !! diatoms 
     613               zphd(ji,jj) = max(0.,trn(ji,jj,jk,jpphd)) 
     614               !! diatom silicon 
     615               zpds(ji,jj) = max(0.,trn(ji,jj,jk,jppds)) 
     616               !! AXY (28/01/10): probably need to take account of  
     617               !! chl/biomass connection 
    607618               if (zchn(ji,jj).eq.0.) zphn(ji,jj) = 0. 
    608619               if (zchd(ji,jj).eq.0.) zphd(ji,jj) = 0. 
    … …  
    612623          if (zpds(ji,jj).eq.0.) zphd(ji,jj) = 0. 
    613624          if (zphd(ji,jj).eq.0.) zpds(ji,jj) = 0. 
    614                zzmi(ji,jj) = max(0.,trn(ji,jj,jk,jpzmi)) !! microzooplankton 
    615                zzme(ji,jj) = max(0.,trn(ji,jj,jk,jpzme)) !! mesozooplankton 
    616                zdet(ji,jj) = max(0.,trn(ji,jj,jk,jpdet)) !! detrital nitrogen 
    617                zdin(ji,jj) = max(0.,trn(ji,jj,jk,jpdin)) !! dissolved inorganic nitrogen 
    618                zsil(ji,jj) = max(0.,trn(ji,jj,jk,jpsil)) !! dissolved silicic acid 
    619                zfer(ji,jj) = max(0.,trn(ji,jj,jk,jpfer)) !! dissolved "iron" 
     625               !! microzooplankton 
     626               zzmi(ji,jj) = max(0.,trn(ji,jj,jk,jpzmi)) 
     627               !! mesozooplankton 
     628               zzme(ji,jj) = max(0.,trn(ji,jj,jk,jpzme)) 
     629               !! detrital nitrogen 
     630               zdet(ji,jj) = max(0.,trn(ji,jj,jk,jpdet)) 
     631               !! dissolved inorganic nitrogen 
     632               zdin(ji,jj) = max(0.,trn(ji,jj,jk,jpdin)) 
     633               !! dissolved silicic acid 
     634               zsil(ji,jj) = max(0.,trn(ji,jj,jk,jpsil)) 
     635               !! dissolved "iron" 
     636               zfer(ji,jj) = max(0.,trn(ji,jj,jk,jpfer)) 
    620637# if defined key_roam 
    621                zdtc(ji,jj) = max(0.,trn(ji,jj,jk,jpdtc)) !! detrital carbon 
    622                zdic(ji,jj) = max(0.,trn(ji,jj,jk,jpdic)) !! dissolved inorganic carbon 
    623                zalk(ji,jj) = max(0.,trn(ji,jj,jk,jpalk)) !! alkalinity 
    624                zoxy(ji,jj) = max(0.,trn(ji,jj,jk,jpoxy)) !! oxygen 
     638               !! detrital carbon 
     639               zdtc(ji,jj) = max(0.,trn(ji,jj,jk,jpdtc)) 
     640               !! dissolved inorganic carbon 
     641               zdic(ji,jj) = max(0.,trn(ji,jj,jk,jpdic)) 
     642               !! alkalinity 
     643               zalk(ji,jj) = max(0.,trn(ji,jj,jk,jpalk)) 
     644               !! oxygen 
     645               zoxy(ji,jj) = max(0.,trn(ji,jj,jk,jpoxy)) 
    625646#  if defined key_axy_carbchem && defined key_mocsy 
    626                zpho(ji,jj) = max(0.,trn(ji,jj,jk,jpdin)) / 16.0 !! phosphate via DIN and Redfield 
     647               !! phosphate via DIN and Redfield 
     648               zpho(ji,jj) = max(0.,trn(ji,jj,jk,jpdin)) / 16.0 
    627649#  endif 
    628650               !! 
    … …  
    633655          !! AXY (28/02/14): check input fields 
    634656               if (ztmp(ji,jj) .lt. -3.0 .or. ztmp(ji,jj) .gt. 40.0 ) then 
    635                   IF(lwp) WRITE(numout,*) ' trc_bio_medusa: T WARNING 2D, ', & 
    636                   tsb(ji,jj,jk,jp_tem), tsn(ji,jj,jk,jp_tem), ' at (',    & 
    637                   ji, ',', jj, ',', jk, ') at time', kt 
    638         IF(lwp) WRITE(numout,*) ' trc_bio_medusa: T SWITCHING 2D, ', & 
    639                   tsn(ji,jj,jk,jp_tem), ' -> ', tsb(ji,jj,jk,jp_tem) 
    640                   ztmp(ji,jj) = tsb(ji,jj,jk,jp_tem) !! temperature 
     657                  IF(lwp) WRITE(numout,*) ' trc_bio_medusa: T WARNING 2D, ',  & 
     658                     tsb(ji,jj,jk,jp_tem), tsn(ji,jj,jk,jp_tem), ' at (',     & 
     659                     ji, ',', jj, ',', jk, ') at time', kt 
     660        IF(lwp) WRITE(numout,*) ' trc_bio_medusa: T SWITCHING 2D, ',& 
     661                     tsn(ji,jj,jk,jp_tem), ' -> ', tsb(ji,jj,jk,jp_tem) 
     662                  !! temperatur 
     663                  ztmp(ji,jj) = tsb(ji,jj,jk,jp_tem) 
    641664               endif 
    642665               if (zsal(ji,jj) .lt. 0.0 .or. zsal(ji,jj) .gt. 45.0 ) then 
    643666                  IF(lwp) WRITE(numout,*) ' trc_bio_medusa: S WARNING 2D, ', & 
    644                   tsb(ji,jj,jk,jp_sal), tsn(ji,jj,jk,jp_sal), ' at (',    & 
    645                   ji, ',', jj, ',', jk, ') at time', kt 
     667                     tsb(ji,jj,jk,jp_sal), tsn(ji,jj,jk,jp_sal), ' at (',    & 
     668                     ji, ',', jj, ',', jk, ') at time', kt 
    646669               endif 
    647670# else 
    648                zdtc(ji,jj) = zdet(ji,jj) * xthetad              !! implicit detrital carbon 
     671               !! implicit detrital carbon 
     672               zdtc(ji,jj) = zdet(ji,jj) * xthetad 
    649673# endif 
    650674# if defined key_debug_medusa 
    651675               if (idf.eq.1) then 
    652                !! AXY (15/01/10) 
     676                  !! AXY (15/01/10) 
    653677                  if (trn(ji,jj,jk,jpdin).lt.0.) then 
    654678                     IF (lwp) write (numout,*) '------------------------------' 
    655                      IF (lwp) write (numout,*) 'NEGATIVE DIN ERROR =', trn(ji,jj,jk,jpdin) 
    656                      IF (lwp) write (numout,*) 'NEGATIVE DIN ERROR @', ji, jj, jk, kt 
     679                     IF (lwp) write (numout,*) 'NEGATIVE DIN ERROR =',        & 
     680                        trn(ji,jj,jk,jpdin) 
     681                     IF (lwp) write (numout,*) 'NEGATIVE DIN ERROR @',        & 
     682                        ji, jj, jk, kt 
    657683                  endif 
    658684                  if (trn(ji,jj,jk,jpsil).lt.0.) then 
    659685                     IF (lwp) write (numout,*) '------------------------------' 
    660                      IF (lwp) write (numout,*) 'NEGATIVE SIL ERROR =', trn(ji,jj,jk,jpsil) 
    661                      IF (lwp) write (numout,*) 'NEGATIVE SIL ERROR @', ji, jj, jk, kt 
     686                     IF (lwp) write (numout,*) 'NEGATIVE SIL ERROR =',        & 
     687                        trn(ji,jj,jk,jpsil) 
     688                     IF (lwp) write (numout,*) 'NEGATIVE SIL ERROR @',        & 
     689                        ji, jj, jk, kt 
    662690                  endif 
    663691#  if defined key_roam 
    664692                  if (trn(ji,jj,jk,jpdic).lt.0.) then 
    665693                     IF (lwp) write (numout,*) '------------------------------' 
    666                      IF (lwp) write (numout,*) 'NEGATIVE DIC ERROR =', trn(ji,jj,jk,jpdic) 
    667                      IF (lwp) write (numout,*) 'NEGATIVE DIC ERROR @', ji, jj, jk, kt 
     694                     IF (lwp) write (numout,*) 'NEGATIVE DIC ERROR =',        & 
     695                        trn(ji,jj,jk,jpdic) 
     696                     IF (lwp) write (numout,*) 'NEGATIVE DIC ERROR @',        & 
     697                        ji, jj, jk, kt 
    668698                  endif 
    669699                  if (trn(ji,jj,jk,jpalk).lt.0.) then 
    670700                     IF (lwp) write (numout,*) '------------------------------' 
    671                      IF (lwp) write (numout,*) 'NEGATIVE ALK ERROR =', trn(ji,jj,jk,jpalk) 
    672                      IF (lwp) write (numout,*) 'NEGATIVE ALK ERROR @', ji, jj, jk, kt 
     701                     IF (lwp) write (numout,*) 'NEGATIVE ALK ERROR =',        & 
     702                        trn(ji,jj,jk,jpalk) 
     703                     IF (lwp) write (numout,*) 'NEGATIVE ALK ERROR @',        & 
     704                        ji, jj, jk, kt 
    673705                  endif 
    674706                  if (trn(ji,jj,jk,jpoxy).lt.0.) then 
    675707                     IF (lwp) write (numout,*) '------------------------------' 
    676                      IF (lwp) write (numout,*) 'NEGATIVE OXY ERROR =', trn(ji,jj,jk,jpoxy) 
    677                      IF (lwp) write (numout,*) 'NEGATIVE OXY ERROR @', ji, jj, jk, kt 
     708                     IF (lwp) write (numout,*) 'NEGATIVE OXY ERROR =',        & 
     709                        trn(ji,jj,jk,jpoxy) 
     710                     IF (lwp) write (numout,*) 'NEGATIVE OXY ERROR @',        & 
     711                        ji, jj, jk, kt 
    678712                  endif 
    679713#  endif 
    … …  
    713747               IF( lk_iomput ) THEN 
    714748                  IF ( med_diag%INVTN%dgsave )   THEN 
    715                      ftot_n(ji,jj)  = ftot_n(ji,jj) + & 
    716                              (fse3t(ji,jj,jk) * ( zphn(ji,jj) + zphd(ji,jj) + zzmi(ji,jj) + zzme(ji,jj) + zdet(ji,jj) + zdin(ji,jj) ) ) 
     749                     ftot_n(ji,jj)  = ftot_n(ji,jj) +                        & 
     750                        (fse3t(ji,jj,jk) * ( zphn(ji,jj) + zphd(ji,jj) +     & 
     751                                             zzmi(ji,jj) + zzme(ji,jj) +     & 
     752                                             zdet(ji,jj) + zdin(ji,jj) ) ) 
    717753                  ENDIF 
    718754                  IF ( med_diag%INVTSI%dgsave )  THEN 
    719                      ftot_si(ji,jj) = ftot_si(ji,jj) + &  
    720                              (fse3t(ji,jj,jk) * ( zpds(ji,jj) + zsil(ji,jj) ) ) 
     755                     ftot_si(ji,jj) = ftot_si(ji,jj) +                       &  
     756                        (fse3t(ji,jj,jk) * ( zpds(ji,jj) + zsil(ji,jj) ) ) 
    721757                  ENDIF 
    722758                  IF ( med_diag%INVTFE%dgsave )  THEN 
    723                      ftot_fe(ji,jj) = ftot_fe(ji,jj) + &  
    724                              (fse3t(ji,jj,jk) * ( xrfn * ( zphn(ji,jj) + zphd(ji,jj) + zzmi(ji,jj) + zzme(ji,jj) + zdet(ji,jj) ) + zfer(ji,jj) ) ) 
     759                     ftot_fe(ji,jj) = ftot_fe(ji,jj) +                       &  
     760                        (fse3t(ji,jj,jk) * ( xrfn *                          & 
     761                                            ( zphn(ji,jj) + zphd(ji,jj) +    & 
     762                                              zzmi(ji,jj) + zzme(ji,jj) +    & 
     763                                              zdet(ji,jj) ) + zfer(ji,jj) ) ) 
    725764                  ENDIF 
    726765# if defined key_roam 
    727766                  IF ( med_diag%INVTC%dgsave )  THEN 
    728767                     ftot_c(ji,jj)  = ftot_c(ji,jj) + &  
    729                              (fse3t(ji,jj,jk) * ( (xthetapn * zphn(ji,jj)) + (xthetapd * zphd(ji,jj)) + & 
    730                              (xthetazmi * zzmi(ji,jj)) + (xthetazme * zzme(ji,jj)) + zdtc(ji,jj) +   & 
    731                              zdic(ji,jj) ) ) 
     768                        (fse3t(ji,jj,jk) * ( (xthetapn * zphn(ji,jj)) +      & 
     769                                             (xthetapd * zphd(ji,jj)) +      & 
     770                                             (xthetazmi * zzmi(ji,jj)) +     & 
     771                                             (xthetazme * zzme(ji,jj)) +     & 
     772                                             zdtc(ji,jj) + zdic(ji,jj) ) ) 
    732773                  ENDIF 
    733774                  IF ( med_diag%INVTALK%dgsave ) THEN 
    734                      ftot_a(ji,jj)  = ftot_a(ji,jj) + (fse3t(ji,jj,jk) * ( zalk(ji,jj) ) ) 
     775                     ftot_a(ji,jj)  = ftot_a(ji,jj) + (fse3t(ji,jj,jk) *     & 
     776                                                       zalk(ji,jj)) 
    735777                  ENDIF 
    736778                  IF ( med_diag%INVTO2%dgsave )  THEN 
    737                      ftot_o2(ji,jj) = ftot_o2(ji,jj) + (fse3t(ji,jj,jk) * ( zoxy(ji,jj) ) ) 
     779                     ftot_o2(ji,jj) = ftot_o2(ji,jj) + (fse3t(ji,jj,jk) *    & 
     780                                                        zoxy(ji,jj)) 
    738781                  ENDIF 
    739782                  !! 
    740783                  !! AXY (10/11/16): CMIP6 diagnostics 
    741784                  IF ( med_diag%INTDISSIC%dgsave ) THEN 
    742                      intdissic(ji,jj) = intdissic(ji,jj) + (fse3t(ji,jj,jk) * zdic(ji,jj)) 
     785                     intdissic(ji,jj) = intdissic(ji,jj) +                   & 
     786                                        (fse3t(ji,jj,jk) * zdic(ji,jj)) 
    743787                  ENDIF 
    744788                  IF ( med_diag%INTDISSIN%dgsave ) THEN 
    745                      intdissin(ji,jj) = intdissin(ji,jj) + (fse3t(ji,jj,jk) * zdin(ji,jj)) 
     789                     intdissin(ji,jj) = intdissin(ji,jj) +                   & 
     790                                        (fse3t(ji,jj,jk) * zdin(ji,jj)) 
    746791                  ENDIF 
    747792                  IF ( med_diag%INTDISSISI%dgsave ) THEN 
    748                      intdissisi(ji,jj) = intdissisi(ji,jj) + (fse3t(ji,jj,jk) * zsil(ji,jj)) 
     793                     intdissisi(ji,jj) = intdissisi(ji,jj) +                 & 
     794                                         (fse3t(ji,jj,jk) * zsil(ji,jj)) 
    749795                  ENDIF 
    750796                  IF ( med_diag%INTTALK%dgsave ) THEN 
    751                      inttalk(ji,jj) = inttalk(ji,jj) + (fse3t(ji,jj,jk) * zalk(ji,jj)) 
     797                     inttalk(ji,jj) = inttalk(ji,jj) +                       & 
     798                                      (fse3t(ji,jj,jk) * zalk(ji,jj)) 
    752799                  ENDIF 
    753800                  IF ( med_diag%O2min%dgsave ) THEN 
    … …  
    755802                        o2min(ji,jj)  = zoxy(ji,jj) 
    756803                        IF ( med_diag%ZO2min%dgsave ) THEN 
    757                            zo2min(ji,jj) = (fsdepw(ji,jj,jk) + fdep1(ji,jj)) / 2. !! layer midpoint 
     804                           !! layer midpoint 
     805                           zo2min(ji,jj) = (fsdepw(ji,jj,jk) +               & 
     806                                            fdep1(ji,jj)) / 2.0 
    758807                        ENDIF 
    759808                     endif 
    … …  
    764813               CALL flush(numout) 
    765814 
    766                !!====================================================================== 
    767                !! LOCAL GRID CELL CALCULATIONS 
    768                !!====================================================================== 
    769                !! 
     815            ENDIF 
     816         ENDDO 
     817         ENDDO 
     818 
     819         !!---------------------------------------------------------------- 
     820         !! Calculate air-sea gas exchange and river inputs 
     821         !!---------------------------------------------------------------- 
     822         IF ( jk == 1 ) THEN 
     823            call air_sea( kt ) 
     824         END IF 
     825 
    770826# if defined key_roam 
    771                if ( jk .eq. 1 ) then 
    772                   !!---------------------------------------------------------------------- 
    773                   !! Air-sea gas exchange 
    774                   !!---------------------------------------------------------------------- 
    775                   IF (lk_oasis) THEN 
    776                      f_xco2a(ji,jj) = PCO2a_in_cpl(ji,jj)        !! use 2D atm xCO2 from atm coupling 
    777                   ENDIF 
    778                   !! 
    779                   !! AXY (23/06/15): as part of an effort to update the carbonate chemistry 
    780                   !!                 in MEDUSA, the gas transfer velocity used in the carbon 
    781                   !!                 and oxygen cycles has been harmonised and is calculated 
    782                   !!                 by the same function here; this harmonisation includes 
    783                   !!                 changes to the PML carbonate chemistry scheme so that 
    784                   !!                 it too makes use of the same gas transfer velocity; the 
    785                   !!                 preferred parameterisation of this is Wanninkhof (2014), 
    786                   !!                 option 7 
    787                   !! 
    788 #   if defined key_debug_medusa 
    789                      IF (lwp) write (numout,*) 'trc_bio_medusa: entering gas_transfer' 
    790                      CALL flush(numout) 
    791 #   endif 
    792                   CALL gas_transfer( wndm(ji,jj), 1, 7, &  ! inputs 
    793                                      f_kw660(ji,jj) )        ! outputs 
    794 #   if defined key_debug_medusa 
    795                      IF (lwp) write (numout,*) 'trc_bio_medusa: exiting gas_transfer' 
    796                      CALL flush(numout) 
    797 #   endif 
    798                   !! 
    799                   !! air pressure (atm); ultimately this will use air pressure at the base 
    800                   !! of the UKESM1 atmosphere  
    801                   !!                                      
    802                   f_pp0(ji,jj)   = 1.0 
    803                   !! 
    804                   !! IF(lwp) WRITE(numout,*) ' MEDUSA ztmp    =', ztmp(ji,jj) 
    805                   !! IF(lwp) WRITE(numout,*) ' MEDUSA zwind_i =', zwind_i(ji,jj) 
    806                   !! IF(lwp) WRITE(numout,*) ' MEDUSA zwind_j =', zwind_j(ji,jj) 
    807                   !! IF(lwp) WRITE(numout,*) ' MEDUSA wndm    =', wndm(ji,jj) 
    808                   !! IF(lwp) WRITE(numout,*) ' MEDUSA fr_i    =', fr_i(ji,jj) 
    809                   !! 
    810 #  if defined key_axy_carbchem 
    811 #   if defined key_mocsy 
    812                   !! 
    813                   !! AXY (22/06/15): use Orr & Epitalon (2015) MOCSY-2 carbonate 
    814                   !!                 chemistry package; note that depth is set to 
    815                   !!                 zero in this call 
    816                   CALL mocsy_interface( ztmp(ji,jj), zsal(ji,jj), zalk(ji,jj), zdic(ji,jj), zsil(ji,jj), zpho(ji,jj),        &  ! inputs 
    817                   f_pp0(ji,jj), 0.0, gphit(ji,jj), f_kw660(ji,jj), f_xco2a(ji,jj), 1,                   &  ! inputs 
    818                   f_ph(ji,jj), f_pco2w(ji,jj), f_fco2w(ji,jj), f_h2co3(ji,jj), f_hco3(ji,jj), f_co3(ji,jj), f_omarg(ji,jj),  &  ! outputs 
    819                   f_omcal(ji,jj), f_BetaD(ji,jj), f_rhosw(ji,jj), f_opres(ji,jj), f_insitut(ji,jj),            &  ! outputs 
    820                   f_pco2atm(ji,jj), f_fco2atm(ji,jj), f_schmidtco2(ji,jj), f_kwco2(ji,jj), f_K0(ji,jj),               &  ! outputs 
    821                   f_co2starair(ji,jj), f_co2flux(ji,jj), f_dpco2(ji,jj) )                                  ! outputs 
    822                   !! 
    823                   f_TDIC(ji,jj) = (zdic(ji,jj) / f_rhosw(ji,jj)) * 1000. ! mmol / m3 -> umol / kg 
    824                   f_TALK(ji,jj) = (zalk(ji,jj) / f_rhosw(ji,jj)) * 1000. !  meq / m3 ->  ueq / kg 
    825                   f_dcf(ji,jj)  = f_rhosw(ji,jj) 
    826 #   else                   
    827                   iters(ji,jj) = 0 
    828                   !! 
    829                   !! carbon dioxide (CO2); Jerry Blackford code (ostensibly OCMIP-2, but not) 
    830                   CALL trc_co2_medusa( ztmp(ji,jj), zsal(ji,jj), zdic(ji,jj), zalk(ji,jj), 0.0, f_kw660(ji,jj), f_xco2a(ji,jj),  &  ! inputs 
    831                   f_ph(ji,jj), f_pco2w(ji,jj), f_h2co3(ji,jj), f_hco3(ji,jj), f_co3(ji,jj), f_omcal(ji,jj),               &  ! outputs 
    832                   f_omarg(ji,jj), f_co2flux(ji,jj), f_TDIC(ji,jj), f_TALK(ji,jj), f_dcf(ji,jj), f_henry(ji,jj), iters(ji,jj) )      ! outputs 
    833                   !! 
    834                   !! AXY (09/01/14): removed iteration and NaN checks; these have 
    835                   !!                 been moved to trc_co2_medusa together with a 
    836                   !!                 fudge that amends erroneous values (this is 
    837                   !!                 intended to be a temporary fudge!); the 
    838                   !!                 output warnings are retained here so that 
    839                   !!                 failure position can be determined 
    840                   if (iters(ji,jj) .eq. 25) then 
    841                      IF(lwp) WRITE(numout,*) ' trc_bio_medusa: ITERS WARNING, ', & 
    842                      iters(ji,jj), ' AT (', ji, ', ', jj, ', ', jk, ') AT ', kt 
    843                   endif 
    844 #   endif 
    845 #  else 
    846                   !! AXY (18/04/13): switch off carbonate chemistry calculations; provide 
    847                   !!                 quasi-sensible alternatives 
    848                   f_ph(ji,jj)           = 8.1 
    849                   f_pco2w(ji,jj)        = f_xco2a(ji,jj) 
    850                   f_h2co3(ji,jj)        = 0.005 * zdic(ji,jj) 
    851                   f_hco3(ji,jj)         = 0.865 * zdic(ji,jj) 
    852                   f_co3(ji,jj)          = 0.130 * zdic(ji,jj) 
    853                   f_omcal(ji,jj) = 4. 
    854                   f_omarg(ji,jj) = 2. 
    855                   f_co2flux(ji,jj)      = 0. 
    856                   f_TDIC(ji,jj)         = zdic(ji,jj) 
    857                   f_TALK(ji,jj)         = zalk(ji,jj) 
    858                   f_dcf(ji,jj)          = 1.026 
    859                   f_henry(ji,jj)        = 1. 
    860                   !! AXY (23/06/15): add in some extra MOCSY diagnostics 
    861                   f_fco2w(ji,jj)        = f_xco2a(ji,jj) 
    862                   f_BetaD(ji,jj)        = 1. 
    863                   f_rhosw(ji,jj)        = 1.026 
    864                   f_opres(ji,jj)        = 0. 
    865                   f_insitut(ji,jj)      = ztmp(ji,jj) 
    866                   f_pco2atm(ji,jj)      = f_xco2a(ji,jj) 
    867                   f_fco2atm(ji,jj)      = f_xco2a(ji,jj) 
    868                   f_schmidtco2(ji,jj)   = 660. 
    869                   f_kwco2(ji,jj)        = 0. 
    870                   f_K0(ji,jj)           = 0. 
    871                   f_co2starair(ji,jj)   = f_xco2a(ji,jj) 
    872                   f_dpco2(ji,jj)        = 0. 
    873 #  endif 
    874                   !! 
    875                   !! mmol/m2/s -> mmol/m3/d; correct for sea-ice; divide through by layer thickness 
    876                   f_co2flux(ji,jj) = (1. - fr_i(ji,jj)) * f_co2flux(ji,jj) * 86400. / fse3t(ji,jj,jk) 
    877                   !! 
    878                   !! oxygen (O2); OCMIP-2 code 
    879                   !! AXY (23/06/15): amend input list for oxygen to account for common gas 
    880                   !!                 transfer velocity 
    881                   !! Note that f_kwo2 is an about from the subroutine below, 
    882                   !! which doesn't seem to be used - marc 10/4/17  
    883                   CALL trc_oxy_medusa( ztmp(ji,jj), zsal(ji,jj), f_kw660(ji,jj), f_pp0(ji,jj), zoxy(ji,jj),  &  ! inputs 
    884                   f_kwo2(ji,jj), f_o2flux(ji,jj), f_o2sat(ji,jj) )                                ! outputs 
    885                   !! 
    886                   !! mmol/m2/s -> mol/m3/d; correct for sea-ice; divide through by layer thickness 
    887                   f_o2flux(ji,jj)  = (1. - fr_i(ji,jj)) * f_o2flux(ji,jj) * 86400. / fse3t(ji,jj,jk) 
    888                   !! 
    889                   !! Jpalm (08-2014) 
    890                   !! DMS surface concentration calculation 
    891                   !! initialy added for UKESM1 model. 
    892                   !! using MET-OFFICE subroutine. 
    893                   !! DMS module only needs Chl concentration and MLD 
    894                   !! to get an aproximate value of DMS concentration. 
    895                   !! air-sea fluxes are calculated by atmospheric chemitry model 
    896                   !! from atm and oc-surface concentrations. 
    897                   !! 
    898                   !! AXY (13/03/15): this is amended to calculate all of the DMS 
    899                   !!                 estimates examined during UKESM1 (see comments 
    900                   !!                 in trcdms_medusa.F90) 
    901                   !! 
    902                   IF (jdms .eq. 1) THEN 
    903                      !! 
    904                      !! feed in correct inputs 
    905                      if (jdms_input .eq. 0) then 
    906                         !! use instantaneous inputs 
    907                         CALL trc_dms_medusa( zchn(ji,jj), zchd(ji,jj), hmld(ji,jj), qsr(ji,jj), zdin(ji,jj), &  ! inputs 
    908                         dms_andr(ji,jj), dms_simo(ji,jj), dms_aran(ji,jj), dms_hall(ji,jj) )                           ! outputs 
    909                      else 
    910                         !! use diel-average inputs 
    911                         CALL trc_dms_medusa( zn_dms_chn(ji,jj), zn_dms_chd(ji,jj), &  ! inputs 
    912                         zn_dms_mld(ji,jj), zn_dms_qsr(ji,jj), zn_dms_din(ji,jj),   &  ! inputs 
    913                         dms_andr(ji,jj), dms_simo(ji,jj), dms_aran(ji,jj), dms_hall(ji,jj) )                      ! outputs 
    914                      endif 
    915                      !! 
    916                      !! assign correct output to variable passed to atmosphere 
    917                      if     (jdms_model .eq. 1) then 
    918                         dms_surf(ji,jj) = dms_andr(ji,jj) 
    919                      elseif (jdms_model .eq. 2) then 
    920                         dms_surf(ji,jj) = dms_simo(ji,jj) 
    921                      elseif (jdms_model .eq. 3) then 
    922                         dms_surf(ji,jj) = dms_aran(ji,jj) 
    923                      elseif (jdms_model .eq. 4) then 
    924                         dms_surf(ji,jj) = dms_hall(ji,jj) 
    925                      endif 
    926                      !! 
    927                      !! 2D diag through iom_use 
    928                      IF( lk_iomput ) THEN 
    929                        IF( med_diag%DMS_SURF%dgsave ) THEN 
    930                          dms_surf2d(ji,jj) = dms_surf(ji,jj) 
    931                        ENDIF 
    932                        IF( med_diag%DMS_ANDR%dgsave ) THEN 
    933                          dms_andr2d(ji,jj) = dms_andr(ji,jj) 
    934                        ENDIF 
    935                        IF( med_diag%DMS_SIMO%dgsave ) THEN 
    936                          dms_simo2d(ji,jj) = dms_simo(ji,jj) 
    937                        ENDIF 
    938                        IF( med_diag%DMS_ARAN%dgsave ) THEN 
    939                          dms_aran2d(ji,jj) = dms_aran(ji,jj) 
    940                        ENDIF 
    941                        IF( med_diag%DMS_HALL%dgsave ) THEN 
    942                          dms_hall2d(ji,jj) = dms_hall(ji,jj) 
    943                        ENDIF 
    944 #   if defined key_debug_medusa 
    945                        IF (lwp) write (numout,*) 'trc_bio_medusa: finish calculating dms' 
    946                      CALL flush(numout) 
    947 #   endif  
    948                      ENDIF 
    949                      !! End iom 
    950                   ENDIF 
    951                   !! End DMS Loop 
    952                   !! 
    953                   !! store 2D outputs 
    954                   !! 
    955                   !! JPALM -- 17-11-16 -- put fgco2 out of diag request 
    956                   !!                    is needed for coupling; pass through restart 
    957                   !! IF( med_diag%FGCO2%dgsave ) THEN 
    958                      !! convert from  mol/m2/day to kg/m2/s 
    959                      fgco2(ji,jj) = f_co2flux(ji,jj) * fse3t(ji,jj,jk) * CO2flux_conv  !! mmol-C/m3/d -> kg-CO2/m2/s 
    960                   !! ENDIF 
    961                   IF ( lk_iomput ) THEN 
    962                       IF( med_diag%ATM_PCO2%dgsave ) THEN 
    963                          f_pco2a2d(ji,jj) = f_pco2atm(ji,jj) 
    964                       ENDIF 
    965                       IF( med_diag%OCN_PCO2%dgsave ) THEN 
    966                          f_pco2w2d(ji,jj) = f_pco2w(ji,jj) 
    967                       ENDIF 
    968                       IF( med_diag%CO2FLUX%dgsave ) THEN 
    969                          f_co2flux2d(ji,jj) = f_co2flux(ji,jj) * fse3t(ji,jj,jk)           !! mmol/m3/d -> mmol/m2/d 
    970                       ENDIF 
    971                       IF( med_diag%TCO2%dgsave ) THEN 
    972                          f_TDIC2d(ji,jj) = f_TDIC(ji,jj) 
    973                       ENDIF 
    974                       IF( med_diag%TALK%dgsave ) THEN 
    975                          f_TALK2d(ji,jj) = f_TALK(ji,jj) 
    976                       ENDIF 
    977                       IF( med_diag%KW660%dgsave ) THEN 
    978                          f_kw6602d(ji,jj) = f_kw660(ji,jj) 
    979                       ENDIF 
    980                       IF( med_diag%ATM_PP0%dgsave ) THEN 
    981                          f_pp02d(ji,jj) = f_pp0(ji,jj) 
    982                       ENDIF 
    983                       IF( med_diag%O2FLUX%dgsave ) THEN 
    984                          f_o2flux2d(ji,jj) = f_o2flux(ji,jj) 
    985                       ENDIF 
    986                       IF( med_diag%O2SAT%dgsave ) THEN 
    987                          f_o2sat2d(ji,jj) = f_o2sat(ji,jj) 
    988                       ENDIF 
    989                       !! AXY (24/11/16): add in extra MOCSY diagnostics 
    990                       IF( med_diag%ATM_XCO2%dgsave ) THEN 
    991                          f_xco2a_2d(ji,jj) = f_xco2a(ji,jj) 
    992                       ENDIF 
    993                       IF( med_diag%OCN_FCO2%dgsave ) THEN 
    994                          f_fco2w_2d(ji,jj) = f_fco2w(ji,jj) 
    995                       ENDIF 
    996                       IF( med_diag%ATM_FCO2%dgsave ) THEN 
    997                          f_fco2a_2d(ji,jj) = f_fco2atm(ji,jj) 
    998                       ENDIF 
    999                       IF( med_diag%OCN_RHOSW%dgsave ) THEN 
    1000                          f_ocnrhosw_2d(ji,jj) = f_rhosw(ji,jj) 
    1001                       ENDIF 
    1002                       IF( med_diag%OCN_SCHCO2%dgsave ) THEN 
    1003                          f_ocnschco2_2d(ji,jj) = f_schmidtco2(ji,jj) 
    1004                       ENDIF 
    1005                       IF( med_diag%OCN_KWCO2%dgsave ) THEN 
    1006                          f_ocnkwco2_2d(ji,jj) = f_kwco2(ji,jj) 
    1007                       ENDIF 
    1008                       IF( med_diag%OCN_K0%dgsave ) THEN 
    1009                          f_ocnk0_2d(ji,jj) = f_K0(ji,jj) 
    1010                       ENDIF 
    1011                       IF( med_diag%CO2STARAIR%dgsave ) THEN 
    1012                          f_co2starair_2d(ji,jj) = f_co2starair(ji,jj) 
    1013                       ENDIF 
    1014                       IF( med_diag%OCN_DPCO2%dgsave ) THEN 
    1015                          f_ocndpco2_2d(ji,jj) = f_dpco2(ji,jj) 
    1016                       ENDIF 
    1017                   ENDIF 
    1018                   !!  
    1019                endif 
    1020                !! End jk = 1 loop within ROAM key  
     827 
     828! Moved from above - marc 21/4/17 
     829! I think this should be moved into diagnostics at bottom - it 
     830! doesn't like it is used anyway else - marc 21/4/17 
     831         DO jj = 2,jpjm1 
     832         DO ji = 2,jpim1 
     833            !! OPEN wet point IF..THEN loop 
     834            if (tmask(ji,jj,jk) == 1) then 
    1021835 
    1022836               !! AXY (11/11/16): CMIP6 oxygen saturation 3D diagnostic 
    … …  
    1025839                  o2sat3(ji, jj, jk) = f_o2sat3 
    1026840               ENDIF 
    1027  
    1028 # endif 
    1029  
    1030                if ( jk .eq. 1 ) then 
    1031                   !!---------------------------------------------------------------------- 
    1032                   !! River inputs 
    1033                   !!---------------------------------------------------------------------- 
    1034                   !! 
    1035                   !! runoff comes in as        kg / m2 / s 
    1036                   !! used and written out as   m3 / m2 / d (= m / d) 
    1037                   !! where                     1000 kg / m2 / d = 1 m3 / m2 / d = 1 m / d 
    1038                   !! 
    1039                   !! AXY (17/07/14): the compiler doesn't like this line for some reason; 
    1040                   !!                 as MEDUSA doesn't even use runoff for riverine inputs, 
    1041                   !!                 a temporary solution is to switch off runoff entirely 
    1042                   !!                 here; again, this change is one of several that will  
    1043                   !!                 need revisiting once MEDUSA has bedded down in UKESM1; 
    1044                   !!                 particularly so if the land scheme provides information 
    1045                   !!                 concerning nutrient fluxes 
    1046                   !! 
    1047                   !! f_runoff(ji,jj) = sf_rnf(1)%fnow(ji,jj,1) / 1000. * 60. * 60. * 24. 
    1048                   f_runoff(ji,jj) = 0.0 
    1049                   !! 
    1050                   !! nutrients are added via rivers to the model in one of two ways: 
    1051                   !!   1. via river concentration; i.e. the average nutrient concentration 
    1052                   !!      of a river water is described by a spatial file, and this is 
    1053                   !!      multiplied by runoff to give a nutrient flux 
    1054                   !!   2. via direct river flux; i.e. the average nutrient flux due to 
    1055                   !!      rivers is described by a spatial file, and this is simply applied 
    1056                   !!      as a direct nutrient flux (i.e. it does not relate or respond to 
    1057                   !!      model runoff) 
    1058                   !! nutrient fields are derived from the GlobalNEWS 2 database; carbon and 
    1059                   !! alkalinity are derived from continent-scale DIC estimates (Huang et al.,  
    1060                   !! 2012) and some Arctic river alkalinity estimates (Katya?) 
    1061                   !!  
    1062                   !! as of 19/07/12, riverine nutrients can now be spread vertically across  
    1063                   !! several grid cells rather than just poured into the surface box; this 
    1064                   !! block of code is still executed, however, to set up the total amounts 
    1065                   !! of nutrient entering via rivers 
    1066                   !! 
    1067                   !! nitrogen 
    1068                   if (jriver_n .eq. 1) then 
    1069                      !! river concentration specified; use runoff to calculate input 
    1070                      f_riv_n(ji,jj) = f_runoff(ji,jj) * riv_n(ji,jj) 
    1071                   elseif (jriver_n .eq. 2) then 
    1072                      !! river flux specified; independent of runoff 
    1073                      f_riv_n(ji,jj) = riv_n(ji,jj) 
    1074                   endif 
    1075                   !! 
    1076                   !! silicon 
    1077                   if (jriver_si .eq. 1) then 
    1078                      !! river concentration specified; use runoff to calculate input 
    1079                      f_riv_si(ji,jj) = f_runoff(ji,jj) * riv_si(ji,jj) 
    1080                   elseif (jriver_si .eq. 2) then 
    1081                      !! river flux specified; independent of runoff 
    1082                      f_riv_si(ji,jj) = riv_si(ji,jj) 
    1083                   endif 
    1084                   !! 
    1085                   !! carbon 
    1086                   if (jriver_c .eq. 1) then 
    1087                      !! river concentration specified; use runoff to calculate input 
    1088                      f_riv_c(ji,jj) = f_runoff(ji,jj) * riv_c(ji,jj) 
    1089                   elseif (jriver_c .eq. 2) then 
    1090                      !! river flux specified; independent of runoff 
    1091                      f_riv_c(ji,jj) = riv_c(ji,jj) 
    1092                   endif 
    1093                   !! 
    1094                   !! alkalinity 
    1095                   if (jriver_alk .eq. 1) then 
    1096                      !! river concentration specified; use runoff to calculate input 
    1097                      f_riv_alk(ji,jj) = f_runoff(ji,jj) * riv_alk(ji,jj) 
    1098                   elseif (jriver_alk .eq. 2) then 
    1099                      !! river flux specified; independent of runoff 
    1100                      f_riv_alk(ji,jj) = riv_alk(ji,jj) 
    1101                   endif 
    1102  
    1103                endif 
     841            ENDIF 
     842         ENDDO 
     843         ENDDO 
     844# endif 
     845 
     846         DO jj = 2,jpjm1 
     847         DO ji = 2,jpim1 
     848            !! OPEN wet point IF..THEN loop 
     849            if (tmask(ji,jj,jk) == 1) then 
    1104850 
    1105851               !!---------------------------------------------------------------------- 
    … …  
    13701116               endif 
    13711117# endif 
     1118 
     1119! MAYBE BUT A BREAK IN HERE, ZOOPLANKTON GRAZING - marc 20/4/17  
     1120! (plankton growth is 281 lines) 
    13721121 
    13731122               !!---------------------------------------------------------------------- 
    … …  
    15201269                  fmegrow(ji,jj) + (xphi * (fgmepn(ji,jj) + fgmepd(ji,jj) + fgmezmi(ji,jj) + fgmed(ji,jj))) + fmeexcr(ji,jj) + ((1.0 - xbetan) * finme(ji,jj)) ) 
    15211270 
     1271! MAYBE BUT A BREAK IN HERE, MISCELLANEOUS PLANKTON LOSSES - marc 20/4/17  
     1272! (zoo plankton grazing is 152 lines) 
     1273 
    15221274               !!---------------------------------------------------------------------- 
    15231275               !! Plankton metabolic losses 
    … …  
    15691321               if (jmzme.eq.4) fdzme(ji,jj) = xmzme * zzme(ji,jj) * &        !! sigmoid 
    15701322                  ((zzme(ji,jj) * zzme(ji,jj)) / (xkzme + (zzme(ji,jj) * zzme(ji,jj)))) 
     1323 
     1324! MAYBE BUT A BREAK IN HERE, DETRITUS PROCESS- marc 20/4/17  
     1325! (misc plankton loses is 53 lines - maybe we can link with another section) 
    15711326 
    15721327               !!---------------------------------------------------------------------- 
    … …  
    16431398               endif 
    16441399 
     1400! MAYBE BUT A BREAK IN HERE, IRON CHEMISTRY AND SCAVENGING - marc 20/4/17  
     1401! (detritus processes is 74 lines) 
     1402 
    16451403               !!---------------------------------------------------------------------- 
    16461404               !! Iron chemistry and fractionation 
    … …  
    19781736# endif 
    19791737 
     1738! MAYBE BUT A BREAK IN HERE, MISCELLANEOUS PROCESSES - marc 20/4/17  
     1739! (iron chemistry and scavenging is 340 lines) 
     1740 
    19801741               !!---------------------------------------------------------------------- 
    19811742               !! Miscellaneous 
    … …  
    20531814!               (xthetazmi * fdzmi2(ji,jj)) + (xthetazme * fdzme2(ji,jj))) * fse3t(ji,jj,jk))        ! linear mortality 
    20541815# endif 
     1816 
     1817! MAYBE BUT A BREAK IN HERE, FAST-SINKINIG DETRITUS - marc 20/4/17  
     1818! (miscellaneous processes is 79 lines) 
    20551819 
    20561820               !!---------------------------------------------------------------------- 
    … …  
    26282392               CALL flush(numout) 
    26292393 
     2394! MAYBE BUT A BREAK IN HERE, BUSINESS AND UPDATING - marc 20/4/17  
     2395! (fast sinking detritus is 576 lines) 
     2396 
    26302397               !!====================================================================== 
    26312398               !! LOCAL GRID CELL TRENDS 
    … …  
    30692836               !!    endif 
    30702837               !! endif      
     2838 
     2839! MAYBE BUT A BREAK IN HERE - marc 20/4/17  
     2840! (this would make the previous section about 470 lines and the one below 
     2841! about 700 lines) 
    30712842 
    30722843# if defined key_trc_diabio 
Note: See TracChangeset for help on using the changeset viewer.