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 15532 – NEMO

Changeset 15532


Ignore:
Timestamp:
2021-11-24T12:47:32+01:00 (6 months ago)
Author:
techene
Message:

#2605 #2715 : version in dev (still buggy)

Location:
NEMO/branches/2021/dev_r14318_RK3_stage1/src
Files:
7 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE/nemogcm.F90

    r14547 r15532  
    489489#if defined key_top 
    490490      !                                      ! Passive tracers 
     491# if defined key_RK3 
     492                           CALL     trc_init( Nbb, Nbb, Naa ) 
     493# else 
    491494                           CALL     trc_init( Nbb, Nnn, Naa ) 
     495# endif 
    492496#endif 
    493497      IF( l_ldfslp     )   CALL ldf_slp_init    ! slope of lateral mixing 
  • NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE/stprk3_stg.F90

    r15514 r15532  
    185185!===>>>>>> Modify dyn_hpg & dyn_hpg_...  routines : rhd computed in dyn_hpg and pass in argument to dyn_hpg_... 
    186186 
    187       CALL eos    ( ts, Kmm, rhd )              ! Kmm in situ density anomaly for hpg computation 
     187!!st      IF( kstg == 3 ) THEN 
     188!         CALL eos    ( ts(:,:,:,:,Kmm), rhd, rhop, gdept_0 ) ! now in situ density for hpg computation 
     189!      ELSE 
     190         CALL eos    ( ts, Kmm, rhd )              ! Kmm in situ density anomaly for hpg computation 
     191!      ENDIF 
    188192 
    189193!!gm end 
     
    388392         vv(:,:,jk,Kaa) = vv(:,:,jk,Kaa) + zvb(:,:)*vmask(:,:,jk) 
    389393      END DO 
     394!!st pourquoi ne pas mettre A2D(0) ici ?  
    390395          
    391396 
  • NEMO/branches/2021/dev_r14318_RK3_stage1/src/OFF/dtadyn.F90

    r14310 r15532  
    176176      ENDIF 
    177177      ! 
    178       CALL eos    ( ts(:,:,:,:,Kmm), rhd, rhop, gdept_0(:,:,:) ) ! In any case, we need rhop 
     178      CALL eos    ( ts(:,:,:,:,Kmm), rhd, gdept_0(:,:,:) ) ! In any case, we need rhd 
    179179      CALL eos_rab( ts(:,:,:,:,Kmm), rab_n, Kmm )       ! now    local thermal/haline expension ratio at T-points 
    180180      CALL bn2    ( ts(:,:,:,:,Kmm), rab_n, rn2, Kmm )  ! before Brunt-Vaisala frequency need for zdfmxl 
     
    193193      ! 
    194194      ! 
    195       CALL eos( ts(:,:,:,:,Kmm), rhd, rhop, gdept_0(:,:,:) ) ! In any case, we need rhop 
     195      CALL eos( ts(:,:,:,:,Kmm), rhd, gdept_0(:,:,:) ) ! In any case, we need rhd 
    196196      ! 
    197197      IF(sn_cfctl%l_prtctl) THEN                 ! print control 
     
    666666      ! 
    667667      IF( l_ldfslp .AND. .NOT.lk_c1d ) THEN    ! Computes slopes (here avt is used as workspace) 
    668          CALL eos    ( pts, rhd, rhop, gdept_0(:,:,:) ) 
     668         CALL eos    ( pts, rhd, gdept_0(:,:,:) ) 
    669669         CALL eos_rab( pts, rab_n, Kmm )       ! now local thermal/haline expension ratio at T-points 
    670670         CALL bn2    ( pts, rab_n, rn2, Kmm  ) ! now    Brunt-Vaisala 
     
    724724      ts(:,:,:,jp_sal,Kmm) = sf_dyn(jf_sal)%fnow(:,:,:)  * tmask(:,:,:)    ! salinity 
    725725      ! 
    726       CALL eos    ( ts(:,:,:,:,Kmm), rhd, rhop, gdept_0(:,:,:) ) ! In any case, we need rhop 
     726      CALL eos    ( ts(:,:,:,:,Kmm), rhd, gdept_0(:,:,:) ) ! In any case, we need rhd 
    727727 
    728728      IF(sn_cfctl%l_prtctl) THEN                     ! print control 
  • NEMO/branches/2021/dev_r14318_RK3_stage1/src/TOP/PISCES/P4Z/p4zche.F90

    r14086 r15532  
    1212   !!                  !  2011-02  (J. Simeon, J.Orr ) update O2 solubility constants 
    1313   !!             3.6  !  2016-03  (O. Aumont) Change chemistry to MOCSY standards 
     14   !!             4.2  !  2020     (J. ORR )  rhop is replaced by "in situ  density" rhd 
    1415   !!---------------------------------------------------------------------- 
    1516   !!   p4z_che      :  Sea water chemistry computed following OCMIP protocol 
     
    195196         !                             ! LN(K0) OF SOLUBILITY OF CO2 (EQ. 12, WEISS, 1980) 
    196197         !                             !     AND FOR THE ATMOSPHERE FOR NON IDEAL GAS 
    197          zcek1 = 9345.17/ztkel - 60.2409 + 23.3585 * LOG(zt) + zsal*(0.023517 - 0.00023656*ztkel    & 
    198          &       + 0.0047036e-4*ztkel**2) 
    199          chemc(ji,jj,1) = EXP( zcek1 ) * 1E-6 * rhop(ji,jj,1) / 1000. ! mol/(L atm) 
     198         ! J. ORR: The previous code has been modified. It computed CO2 solubility in mol/(kg*atm), then converted that to mol/(L*atm). 
     199         ! But Weiss (1974) provides sets of coefficients for each of those 2 units. 
     200         ! Thus I have changed the code to use the coefficients for mol*L/atm. 
     201         ! Hence I've eliminated using the conversion (which used the variable rhop) 
     202         ! OLD - Coefficients for CO2 soulbility in mol/(kg*atm) (Weiss,1974, Table 1, column 2) 
     203         !zcek1 = 9345.17/ztkel - 60.2409 + 23.3585 * LOG(zt) + zsal*(0.023517 - 0.00023656*ztkel    & 
     204         !&       + 0.0047036e-4*ztkel**2) 
     205         ! NEW - Coefficients for CO2 soulbility in mol/(L*atm) (Weiss, 1974, Table 1, column 1) 
     206         zcek1 = 9050.69/ztkel - 58.0931 + 22.2940 * LOG(zt) + zsal*(0.027766 - 0.00025888*ztkel    & 
     207                 &       + 0.0050578e-4*ztkel**2) 
     208         ! 
     209         ! OLD:  chemc(ji,jj,1) = EXP( zcek1 ) * 1E-6 * rhop(ji,jj,1) / 1000. ! mol/(L atm) 
     210         ! The units indicated in the above line are wrong. They are actually "mol/(L*uatm)" 
     211         ! NEW: 
     212         chemc(ji,jj,1) = EXP( zcek1 ) * 1E-6 ! mol/(L * uatm) 
    200213         chemc(ji,jj,2) = -1636.75 + 12.0408*ztkel - 0.0327957*ztkel**2 + 0.0000316528*ztkel**3 
    201214         chemc(ji,jj,3) = 57.7 - 0.118*ztkel 
     
    445458      REAL(wp)  ::  zca1, zba1 
    446459      REAL(wp)  ::  zd, zsqrtd, zhmin 
    447       REAL(wp)  ::  za2, za1, za0 
     460      REAL(wp)  ::  za2, za1, za0, zrhd 
    448461      REAL(wp)  ::  p_dictot, p_bortot, p_alkcb  
    449462      !!--------------------------------------------------------------------- 
     
    452465      ! 
    453466      DO_3D( 1, 1, 1, 1, 1, jpk ) 
    454       p_alkcb  = tr(ji,jj,jk,jptal,Kbb) * 1000. / (rhop(ji,jj,jk) + rtrn) 
    455       p_dictot = tr(ji,jj,jk,jpdic,Kbb) * 1000. / (rhop(ji,jj,jk) + rtrn) 
     467      zrhd = 1._wp / ( rhd(ji,jj,jk) + 1. ) 
     468      p_alkcb  = tr(ji,jj,jk,jptal,Kbb) * zrhd 
     469      p_dictot = tr(ji,jj,jk,jpdic,Kbb) * zrhd 
    456470      p_bortot = borat(ji,jj,jk) 
    457471      IF (p_alkcb <= 0.) THEN 
     
    502516   REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(OUT) :: p_alknw_sup 
    503517   INTEGER,                          INTENT(in)  ::  Kbb      ! time level indices 
    504  
    505    p_alknw_inf(:,:,:) =  -tr(:,:,:,jppo4,Kbb) * 1000. / (rhop(:,:,:) + rtrn) - sulfat(:,:,:)  & 
    506    &              - fluorid(:,:,:) 
    507    p_alknw_sup(:,:,:) =   (2. * tr(:,:,:,jpdic,Kbb) + 2. * tr(:,:,:,jppo4,Kbb) + tr(:,:,:,jpsil,Kbb) )    & 
    508    &               * 1000. / (rhop(:,:,:) + rtrn) + borat(:,:,:)  
     518   INTEGER  ::   ji, jj, jk 
     519   REAL(wp)  ::  zrhd 
     520 
     521    DO_3D( 1, 1, 1, 1, 1, jpk ) 
     522      zrhd = 1._wp / ( rhd(ji,jj,jk) + 1. ) 
     523      p_alknw_inf(ji,jj,jk) =  -tr(ji,jj,jk,jppo4,Kbb) * zrhd - sulfat(ji,jj,jk) & 
     524      &              - fluorid(ji,jj,jk) 
     525      p_alknw_sup(ji,jj,jk) =   (2. * tr(ji,jj,jk,jpdic,Kbb) + 2. * tr(ji,jj,jk,jppo4,Kbb) + tr(ji,jj,jk,jpsil,Kbb) )    & 
     526      &               * zrhd + borat(ji,jj,jk) 
     527    END_3D 
    509528 
    510529   END SUBROUTINE anw_infsup 
     
    536555   REAL(wp)  ::  znumer_flu, zdnumer_flu, zdenom_flu, zalk_flu, zdalk_flu 
    537556   REAL(wp)  ::  zalk_wat, zdalk_wat 
    538    REAL(wp)  ::  zfact, p_alktot, zdic, zbot, zpt, zst, zft, zsit 
     557   REAL(wp)  ::  zrhd, p_alktot, zdic, zbot, zpt, zst, zft, zsit 
    539558   LOGICAL   ::  l_exitnow 
    540559   REAL(wp), PARAMETER :: pz_exp_threshold = 1.0 
     
    551570   DO_3D( 1, 1, 1, 1, 1, jpk ) 
    552571      IF (rmask(ji,jj,jk) == 1.) THEN 
    553          p_alktot = tr(ji,jj,jk,jptal,Kbb) * 1000. / (rhop(ji,jj,jk) + rtrn) 
     572         zrhd = 1._wp / ( rhd(ji,jj,jk) + 1. ) 
     573         p_alktot = tr(ji,jj,jk,jptal,Kbb) * zrhd 
    554574         aphscale = 1. + sulfat(ji,jj,jk)/aks3(ji,jj,jk) 
    555575         zh_ini = p_hini(ji,jj,jk) 
     
    580600   DO_3D( 1, 1, 1, 1, 1, jpk ) 
    581601      IF (rmask(ji,jj,jk) == 1.) THEN 
    582          zfact = rhop(ji,jj,jk) / 1000. + rtrn 
    583          p_alktot = tr(ji,jj,jk,jptal,Kbb) / zfact 
    584          zdic  = tr(ji,jj,jk,jpdic,Kbb) / zfact 
     602         zrhd = 1._wp / ( rhd(ji,jj,jk) + 1. ) 
     603         p_alktot = tr(ji,jj,jk,jptal,Kbb) * zrhd 
     604         zdic  = tr(ji,jj,jk,jpdic,Kbb) * zrhd 
    585605         zbot  = borat(ji,jj,jk) 
    586          zpt = tr(ji,jj,jk,jppo4,Kbb) / zfact * po4r 
    587          zsit = tr(ji,jj,jk,jpsil,Kbb) / zfact 
     606         zpt = tr(ji,jj,jk,jppo4,Kbb) * zrhd * po4r 
     607         zsit = tr(ji,jj,jk,jpsil,Kbb) * zrhd 
    588608         zst = sulfat (ji,jj,jk) 
    589609         zft = fluorid(ji,jj,jk) 
  • NEMO/branches/2021/dev_r14318_RK3_stage1/src/TOP/PISCES/P4Z/p4zflx.F90

    r13295 r15532  
    1010   !!            2.0  !  2007-12  (C. Ethe, G. Madec)  F90 
    1111   !!                 !  2011-02  (J. Simeon, J. Orr) Include total atm P correction  
     12   !!             4.2  !  2020     (J. ORR )  rhop is replaced by "in situ density" rhd 
    1213   !!---------------------------------------------------------------------- 
    1314   !!   p4z_flx       :   CALCULATES GAS EXCHANGE AND CHEMISTRY AT SEA SURFACE 
     
    7879      INTEGER  ::   ji, jj, jm, iind, iindm1 
    7980      REAL(wp) ::   ztc, ztc2, ztc3, ztc4, zws, zkgwan 
    80       REAL(wp) ::   zfld, zflu, zfld16, zflu16, zfact 
     81      REAL(wp) ::   zfld, zflu, zfld16, zflu16, zrhd 
    8182      REAL(wp) ::   zvapsw, zsal, zfco2, zxc2, xCO2approx, ztkel, zfugcoeff 
    8283      REAL(wp) ::   zph, zdic, zsch_o2, zsch_co2 
     
    112113      DO_2D( 1, 1, 1, 1 ) 
    113114         ! DUMMY VARIABLES FOR DIC, H+, AND BORATE 
    114          zfact = rhop(ji,jj,1) / 1000. + rtrn 
     115         zrhd = rhd(ji,jj,1) + 1._wp 
    115116         zdic  = tr(ji,jj,1,jpdic,Kbb) 
    116          zph   = MAX( hi(ji,jj,1), 1.e-10 ) / zfact 
     117         zph   = MAX( hi(ji,jj,1), 1.e-10 ) / ( zrhd + rtrn ) 
    117118         ! CALCULATE [H2CO3] 
    118119         zh2co3(ji,jj) = zdic/(1. + ak13(ji,jj,1)/zph + ak13(ji,jj,1)*ak23(ji,jj,1)/zph**2) 
  • NEMO/branches/2021/dev_r14318_RK3_stage1/src/TOP/PISCES/P4Z/p4zlys.F90

    r13295 r15532  
    1212   !!             3.4  !  2011-06  (O. Aumont, C. Ethe) Improvment of calcite dissolution 
    1313   !!             3.6  !  2015-05  (O. Aumont) PISCES quota 
     14   !!             4.2  !  2020     (J. ORR )  rhop is replaced by "in situ density" rhd 
    1415   !!---------------------------------------------------------------------- 
    1516   !!   p4z_lys        :   Compute the CaCO3 dissolution  
     
    3334 
    3435   INTEGER  ::   rmtss              ! number of seconds per month  
    35    REAL(wp) ::   calcon = 1.03E-2   ! mean calcite concentration [Ca2+] in sea water [mole/kg solution] 
    36   
     36 
     37   !! * Module variables 
     38   REAL(wp) :: calcon = 1.03E-2           !: mean calcite concentration [Ca2+]  in sea water [mole/kg solution] 
     39   ! J. ORR: Made consistent with mocsy's choice based on literature review from Munhoven 
     40!   REAL(wp) :: calcon = 1.0287E-2           !: mean calcite concentration [Ca2+] in sea water [mole/kg solution] 
     41 
    3742   !! * Substitutions 
    3843#  include "do_loop_substitute.h90" 
     
    5964      ! 
    6065      INTEGER  ::   ji, jj, jk, jn 
    61       REAL(wp) ::   zdispot, zfact, zcalcon 
     66      REAL(wp) ::   zdispot, zrhd, zcalcon 
    6267      REAL(wp) ::   zomegaca, zexcess, zexcess0 
    6368      CHARACTER (len=25) ::   charout 
     
    6772      IF( ln_timing )  CALL timing_start('p4z_lys') 
    6873      ! 
    69       zhinit  (:,:,:) = hi(:,:,:) * 1000. / ( rhop(:,:,:) + rtrn ) 
     74 
     75      zhinit  (:,:,:) = hi(:,:,:) / ( rhd(:,:,:) + 1._wp ) 
    7076      ! 
    7177      !     ------------------------------------------- 
     
    7884         zco3(ji,jj,jk) = tr(ji,jj,jk,jpdic,Kbb) * ak13(ji,jj,jk) * ak23(ji,jj,jk) / (zhi(ji,jj,jk)**2   & 
    7985            &             + ak13(ji,jj,jk) * zhi(ji,jj,jk) + ak13(ji,jj,jk) * ak23(ji,jj,jk) + rtrn ) 
    80          hi  (ji,jj,jk) = zhi(ji,jj,jk) * rhop(ji,jj,jk) / 1000. 
     86         hi  (ji,jj,jk) = zhi(ji,jj,jk) * ( rhd(ji,jj,jk) + 1._wp ) 
    8187      END_3D 
    8288 
     
    9096 
    9197         ! DEVIATION OF [CO3--] FROM SATURATION VALUE 
    92          ! Salinity dependance in zomegaca and divide by rhop/1000 to have good units 
     98         ! Salinity dependance in zomegaca and divide by rhd to have good units 
    9399         zcalcon  = calcon * ( salinprac(ji,jj,jk) / 35._wp ) 
    94          zfact    = rhop(ji,jj,jk) / 1000._wp 
    95          zomegaca = ( zcalcon * zco3(ji,jj,jk) ) / ( aksp(ji,jj,jk) * zfact + rtrn ) 
    96          zco3sat(ji,jj,jk) = aksp(ji,jj,jk) * zfact / ( zcalcon + rtrn ) 
     100         zrhd    = rhd(ji,jj,jk) + 1._wp 
     101         zomegaca = ( zcalcon * zco3(ji,jj,jk) ) / ( aksp(ji,jj,jk) * zrhd + rtrn ) 
     102         zco3sat(ji,jj,jk) = aksp(ji,jj,jk) * zrhd / ( zcalcon + rtrn ) 
    97103 
    98104         ! SET DEGREE OF UNDER-/SUPERSATURATION 
  • NEMO/branches/2021/dev_r14318_RK3_stage1/src/TOP/TRP/trcsbc.F90

    r15380 r15532  
    118118         zsfx(:,:) = emp(:,:) 
    119119      ENDIF 
     120!!st DO LOOPs used to look like DO_2D( 0, 0, 0, 1 ) I changed it because it does not work in debug mode...  
    120121 
    121122      ! 0. initialization 
     
    125126         ! 
    126127         DO jn = 1, jptra 
    127             DO_2D( 0, 0, 0, 1 ) 
     128            DO_2D( 0, 0, 0, 0 ) 
    128129               sbc_trc(ji,jj,jn) = zsfx(ji,jj) * r1_rho0 * ptr(ji,jj,1,jn,Kmm) 
    129130            END_2D 
     
    133134         ! 
    134135         DO jn = 1, jptra 
    135             DO_2D( 0, 0, 0, 1 ) 
     136            DO_2D( 0, 0, 0, 0 ) 
    136137               sbc_trc(ji,jj,jn) = ( zsfx(ji,jj) + fmmflx(ji,jj) ) * r1_rho0 * ptr(ji,jj,1,jn,Kmm) 
    137138            END_2D 
     
    141142         ! 
    142143         DO jn = 1, jptra 
    143             DO_2D( 0, 0, 0, 1 ) 
     144            DO_2D( 0, 0, 0, 0 ) 
    144145               ! tracer flux at the ice/ocean interface (tracer/m2/s) 
    145146               zftra = - trc_i(ji,jj,jn) * fmmflx(ji,jj) ! uptake of tracer in the sea ice 
     
    164165         IF( l_trdtrc )   ztrtrd(:,:,:) = ptr(:,:,:,jn,Krhs)  ! save trends 
    165166         ! 
    166          DO_2D( 0, 0, 0, 1 ) 
     167         DO_2D( 0, 0, 0, 0 ) 
    167168#if defined key_RK3 
    168169            zse3t = 1._wp / e3t(ji,jj,1,Kmm) 
     
    257258!!not sure about trc_i case... (1) 
    258259            DO jn = 1, jptra 
    259                DO_2D( 0, 0, 0, 1 )              !!st WHY 1 : exterior here ?  
     260               DO_2D( 0, 0, 0, 0 )              !!st WHY 1 : exterior here ?  
    260261                  z1_rho0_e3t = r1_rho0 / e3t(ji,jj,1,Kmm) 
    261262                  ptr(ji,jj,1,jn,Krhs) = ptr(ji,jj,1,jn,Krhs) - emp(ji,jj) * ptr(ji,jj,1,jn,Kmm) * z1_rho0_e3t 
     
    281282               ! 
    282283               DO jn = 1, jptra 
    283                   DO_2D( 0, 0, 0, 1 ) 
     284                  DO_2D( 0, 0, 0, 0 ) 
    284285                     z1_rho0_e3t = r1_rho0  / e3t(ji,jj,1,Kmm) 
    285286                     ptr(ji,jj,1,jn,Krhs) = ptr(ji,jj,1,jn,Krhs) + emp(ji,jj) * r1_rho0 * ptr(ji,jj,1,jn,Kmm) 
     
    299300               ! 
    300301               DO jn = 1, jptra 
    301                   DO_2D( 0, 0, 0, 1 ) 
     302                  DO_2D( 0, 0, 0, 0 ) 
    302303                     z1_rho0_e3t = r1_rho0  / e3t(ji,jj,1,Kmm) 
    303304                     ! tracer flux at the ice/ocean interface (tracer/m2/s) 
     
    335336               ! 
    336337               DO jn = 1, jptra 
    337                   DO_2D( 0, 0, 0, 1 ) 
     338                  DO_2D( 0, 0, 0, 0 ) 
    338339                     ! tracer flux at the ice/ocean interface (tracer/m2/s) 
    339340                     zftra = - trc_i(ji,jj,jn) * fmmflx(ji,jj) ! uptake of tracer in the sea ice 
Note: See TracChangeset for help on using the changeset viewer.