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 15380 for NEMO/branches/2021/dev_r14318_RK3_stage1/src/TOP/TRP/trcsbc.F90 – NEMO

Ignore:
Timestamp:
2021-10-15T11:08:37+02:00 (12 months ago)
Author:
techene
Message:

#2715 refactoring and bug correction : sbc_trc variable no longer useful for RK3

File:
1 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2021/dev_r14318_RK3_stage1/src/TOP/TRP/trcsbc.F90

    r15373 r15380  
    232232      REAL(wp), DIMENSION(jpi,jpj,jpk,jptra,jpt), INTENT(inout) ::   ptr       ! passive tracers and RHS of tracer equation 
    233233      ! 
    234       INTEGER  ::   ji, jj, jn                      ! dummy loop indices 
    235       REAL(wp) ::   zse3t, zrtrn, zfact     ! local scalars 
    236       REAL(wp) ::   zftra, zdtra, ztfx, ztra   !   -      - 
     234      INTEGER  ::   ji, jj, jn           ! dummy loop indices 
     235      REAL(wp) ::   z1_rho0_e3t          ! local scalars 
     236      REAL(wp) ::   zftra, zdtra, ztfx   !   -      - 
    237237      CHARACTER (len=22) :: charout 
    238       REAL(wp), DIMENSION(jpi,jpj)   ::   zmfx 
    239238      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ztrtrd 
    240239      !!--------------------------------------------------------------------- 
    241240      ! 
    242241      IF( ln_timing )   CALL timing_start('trc_sbc_RK3') 
    243       ! 
    244       ! Allocate temporary workspace 
    245       IF( l_trdtrc )  ALLOCATE( ztrtrd(jpi,jpj,jpk) ) 
    246242      ! 
    247243      IF( kt == nittrc000 ) THEN 
     
    250246         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ ' 
    251247      ENDIF 
    252 !!st: not sure of the way to deal about this 
     248      ! 
     249!!st note that trc_sbc can be removed only re-use in atf (not relevant for RK3) 
    253250      SELECT CASE( kstg ) 
    254       ! 
     251         ! 
    255252      CASE( 1 , 2 )                       !=  stage 1 and 2  =!   only in non linear ssh 
    256253         ! 
    257254         IF( .NOT.ln_linssh ) THEN           !* only passive tracer fluxes associated with mass fluxes 
     255            !                                        ! no passive tracer concentration modification due to ssh variation 
     256!!st emp includes fmm see iceupdate.F90 
     257!!not sure about trc_i case... (1) 
     258            DO jn = 1, jptra 
     259               DO_2D( 0, 0, 0, 1 )              !!st WHY 1 : exterior here ?  
     260                  z1_rho0_e3t = r1_rho0 / e3t(ji,jj,1,Kmm) 
     261                  ptr(ji,jj,1,jn,Krhs) = ptr(ji,jj,1,jn,Krhs) - emp(ji,jj) * ptr(ji,jj,1,jn,Kmm) * z1_rho0_e3t 
     262               END_2D 
     263            END DO 
     264            ! 
     265         ENDIF 
     266         ! 
     267      CASE( 3 ) 
     268         ! 
     269         ! Allocate temporary workspace 
     270         IF( l_trdtrc )  ALLOCATE( ztrtrd(jpi,jpj,jpk) ) 
     271         ! 
     272         DO jn = 1, jptra 
     273            IF( l_trdtrc )   ztrtrd(:,:,:) = ptr(:,:,:,jn,Krhs)  ! save trends 
     274         END DO 
     275         ! 
     276         IF( ln_linssh ) THEN                !* linear free surface (add concentration/dilution effect artificially since no volume variation) 
    258277            ! 
    259278            SELECT CASE ( nn_ice_tr ) 
    260279               ! 
    261             CASE ( -1 , 0 , 1 )                      ! no passive tracer concentration modification due to ssh variation 
    262                ! 
    263 !!st emp includes fmm see iceupdate.F90 
     280            CASE ( -1 ) ! No tracers in sea ice (null concentration in sea ice) 
     281               ! 
    264282               DO jn = 1, jptra 
    265283                  DO_2D( 0, 0, 0, 1 ) 
    266                      sbc_trc(ji,jj,jn) =  - emp(ji,jj) * ptr(ji,jj,1,jn,Kmm) * r1_rho0 
     284                     z1_rho0_e3t = r1_rho0  / e3t(ji,jj,1,Kmm) 
     285                     ptr(ji,jj,1,jn,Krhs) = ptr(ji,jj,1,jn,Krhs) + emp(ji,jj) * r1_rho0 * ptr(ji,jj,1,jn,Kmm) 
    267286                  END_2D 
    268287               END DO 
    269288               ! 
    270             END SELECT 
    271             ! 
    272          ENDIF 
    273          ! 
    274       CASE( 3 ) 
    275 !!st copy of existing code 
    276          ! 
    277          IF( .NOT.ln_linssh ) THEN    !!st concentration/dilution effect due to volume variation  
    278             SELECT CASE ( nn_ice_tr ) 
    279             ! 
    280             CASE ( 0 )  ! Same concentration in sea ice and in the ocean 
     289            CASE ( 0 )  ! Same concentration in sea ice and in the ocean fmm contribution to concentration/dilution effect has to be removed 
    281290               ! 
    282291               DO jn = 1, jptra 
    283292                  DO_2D( 0, 0, 0, 1 ) 
    284                      sbc_trc(ji,jj,jn) =  - fmmflx(ji,jj) * r1_rho0 * ptr(ji,jj,1,jn,Kmm) 
     293                     z1_rho0_e3t = r1_rho0  / e3t(ji,jj,1,Kmm) 
     294                     ptr(ji,jj,1,jn,Krhs) = ptr(ji,jj,1,jn,Krhs) + ( emp(ji,jj) - fmmflx(ji,jj) ) * r1_rho0 * ptr(ji,jj,1,jn,Kmm) 
    285295                  END_2D 
    286296               END DO 
    287297               ! 
    288             CASE ( 1 )  ! Specific treatment of sea ice fluxes with an imposed concentration in sea ice  
     298            CASE ( 1 )  ! Specific treatment of sea ice fluxes with an imposed concentration in sea ice !!st TODO : check Christian new implementation 
    289299               ! 
    290300               DO jn = 1, jptra 
    291301                  DO_2D( 0, 0, 0, 1 ) 
    292                   ! tracer flux at the ice/ocean interface (tracer/m2/s) 
    293                   zftra = - trc_i(ji,jj,jn) * fmmflx(ji,jj) ! uptake of tracer in the sea ice 
    294                   !                                         ! only used in the levitating sea ice case 
    295                   ! tracer flux only       : add concentration dilution term in net tracer flux, no F-M in volume flux 
    296                   ! tracer and mass fluxes : no concentration dilution term in net tracer flux, F-M term in volume flux 
    297                   ztfx  = zftra                        ! net tracer flux 
    298                   ! 
    299                   zdtra = r1_rho0 * ( ztfx -  fmmflx(ji,jj) * ptr(ji,jj,1,jn,Kmm) )  
    300                   IF ( zdtra < 0. ) THEN 
    301                      zdtra  = MAX(zdtra, -ptr(ji,jj,1,jn,Kmm) * e3t(ji,jj,1,Kmm) / rDt_trc )   ! avoid negative concentrations to arise 
    302                   ENDIF 
    303                   sbc_trc(ji,jj,jn) =  zdtra  
     302                     z1_rho0_e3t = r1_rho0  / e3t(ji,jj,1,Kmm) 
     303                     ! tracer flux at the ice/ocean interface (tracer/m2/s) 
     304                     zftra = - trc_i(ji,jj,jn) * fmmflx(ji,jj) ! uptake of tracer in the sea ice 
     305                     !                                         ! only used in the levitating sea ice case 
     306                     ! tracer flux only       : add concentration dilution term in net tracer flux, no F-M in volume flux 
     307                     ! tracer and mass fluxes : no concentration dilution term in net tracer flux, F-M term in volume flux 
     308                     ztfx  = zftra                        ! net tracer flux 
     309                     ! 
     310                     zdtra = r1_rho0 * ( ztfx +  ( emp(ji,jj) - fmmflx(ji,jj) ) * ptr(ji,jj,1,jn,Kmm) )  
     311                     IF ( zdtra < 0. ) THEN 
     312                        zdtra  = MAX(zdtra, -ptr(ji,jj,1,jn,Kmm) * e3t(ji,jj,1,Kmm) / rDt_trc )   ! avoid negative concentrations to arise 
     313                     ENDIF 
     314                     ptr(ji,jj,1,jn,Krhs) = ptr(ji,jj,1,jn,Krhs) + zdtra 
    304315                  END_2D 
    305316               END DO 
    306317               ! 
    307318            END SELECT 
    308          ELSE                         !linear ssh !!st need to mimic concentration/dilution effect since no volume variation  
     319            ! 
     320         ELSE                                !* non linear free surface (concentration/dilution effect due to volume variation) 
     321            ! 
    309322            SELECT CASE ( nn_ice_tr ) 
    310             ! 
    311             CASE ( -1 ) ! No tracers in sea ice (null concentration in sea ice) 
     323            ! CASE ( -1 ) natural concentration/dilution effect due to volume variation : nothing to do 
     324            ! 
     325            CASE ( 0 )  ! Same concentration in sea ice and in the ocean : correct concentration/dilution effect due to "freezing - melting" 
    312326               ! 
    313327               DO jn = 1, jptra 
    314328                  DO_2D( 0, 0, 0, 1 ) 
    315                   sbc_trc(ji,jj,jn) = emp(ji,jj) * r1_rho0 * ptr(ji,jj,1,jn,Kmm) 
     329                     z1_rho0_e3t = r1_rho0  / e3t(ji,jj,1,Kmm) 
     330                     ptr(ji,jj,1,jn,Krhs) = ptr(ji,jj,1,jn,Krhs) - fmmflx(ji,jj) * r1_rho0 * ptr(ji,jj,1,jn,Kmm) 
    316331                  END_2D 
    317332               END DO 
    318333               ! 
    319             CASE ( 0 )  ! Same concentration in sea ice and in the ocean 
     334            CASE ( 1 )  ! Specific treatment of sea ice fluxes with an imposed concentration in sea ice  
    320335               ! 
    321336               DO jn = 1, jptra 
    322337                  DO_2D( 0, 0, 0, 1 ) 
    323                      sbc_trc(ji,jj,jn) = ( emp(ji,jj) - fmmflx(ji,jj) ) * r1_rho0 * ptr(ji,jj,1,jn,Kmm) 
     338                     ! tracer flux at the ice/ocean interface (tracer/m2/s) 
     339                     zftra = - trc_i(ji,jj,jn) * fmmflx(ji,jj) ! uptake of tracer in the sea ice 
     340                     !                                         ! only used in the levitating sea ice case 
     341                     ! tracer flux only       : add concentration dilution term in net tracer flux, no F-M in volume flux 
     342                     ! tracer and mass fluxes : no concentration dilution term in net tracer flux, F-M term in volume flux 
     343                     ztfx  = zftra                        ! net tracer flux 
     344                     ! 
     345                     zdtra = r1_rho0 * ( ztfx -  fmmflx(ji,jj) * ptr(ji,jj,1,jn,Kmm) )  
     346                     IF ( zdtra < 0. ) THEN 
     347                        zdtra  = MAX(zdtra, -ptr(ji,jj,1,jn,Kmm) * e3t(ji,jj,1,Kmm) / rDt_trc )   ! avoid negative concentrations to arise 
     348                     ENDIF 
     349                     ptr(ji,jj,1,jn,Krhs) = ptr(ji,jj,1,jn,Krhs) + zdtra 
    324350                  END_2D 
    325351               END DO 
    326352               ! 
    327             CASE ( 1 )  ! Specific treatment of sea ice fluxes with an imposed concentration in sea ice  
    328                ! 
    329                DO jn = 1, jptra 
    330                   DO_2D( 0, 0, 0, 1 ) 
    331                   ! tracer flux at the ice/ocean interface (tracer/m2/s) 
    332                   zftra = - trc_i(ji,jj,jn) * fmmflx(ji,jj) ! uptake of tracer in the sea ice 
    333                   !                                         ! only used in the levitating sea ice case 
    334                   ! tracer flux only       : add concentration dilution term in net tracer flux, no F-M in volume flux 
    335                   ! tracer and mass fluxes : no concentration dilution term in net tracer flux, F-M term in volume flux 
    336                   ztfx  = zftra                        ! net tracer flux 
    337                   ! 
    338                   zdtra = r1_rho0 * ( ztfx +  ( emp(ji,jj) - fmmflx(ji,jj) ) * ptr(ji,jj,1,jn,Kmm) )  
    339                   IF ( zdtra < 0. ) THEN 
    340                      zdtra  = MAX(zdtra, -ptr(ji,jj,1,jn,Kmm) * e3t(ji,jj,1,Kmm) / rDt_trc )   ! avoid negative concentrations to arise 
    341                   ENDIF 
    342                   sbc_trc(ji,jj,jn) =  zdtra  
    343                   END_2D 
    344                END DO 
    345                ! 
    346353            END SELECT 
    347354            ! 
     
    349356         ! 
    350357         ! 
     358!!st useless trc_sbc only in the interior even in MLF case         CALL lbc_lnk( 'trcsbc', sbc_trc(:,:,:), 'T', 1.0_wp ) 
     359         !                                       Concentration dilution effect on tracers due to evaporation & precipitation  
     360         DO jn = 1, jptra 
     361            ! 
     362            IF(lwp) WRITE(numout,*) 
     363            IF(lwp) WRITE(numout,*) 'trc_sbc_RK3 : Runge Kutta 3rd order at stage ', kstg, jn 
     364            IF(lwp) WRITE(numout,*) 
     365            ! 
     366            IF( l_trdtrc ) THEN 
     367               ztrtrd(:,:,:) = ptr(:,:,:,jn,Krhs) - ztrtrd(:,:,:) 
     368               CALL trd_tra( kt, Kmm, Krhs, 'TRC', jn, jptra_nsr, ztrtrd ) 
     369            END IF 
     370            ! 
     371         END DO 
     372         ! 
     373         IF( l_trdtrc )  DEALLOCATE( ztrtrd ) 
    351374         ! 
    352375      END SELECT 
    353       ! 
    354       CALL lbc_lnk( 'trcsbc', sbc_trc(:,:,:), 'T', 1.0_wp ) 
    355       !                                       Concentration dilution effect on tracers due to evaporation & precipitation  
    356       DO jn = 1, jptra 
    357          ! 
    358          IF( l_trdtrc )   ztrtrd(:,:,:) = ptr(:,:,:,jn,Krhs)  ! save trends 
    359          ! 
    360          IF(lwp) WRITE(numout,*) 
    361          IF(lwp) WRITE(numout,*) 'trc_sbc_RK3 : Runge Kutta 3rd order at stage ', kstg, jn 
    362          IF(lwp) WRITE(numout,*) 'emp', MAXVAL(emp(:,:)) 
    363          IF(lwp) WRITE(numout,*) 'sbc_trc', MAXVAL(sbc_trc(:,:,jn)) 
    364          IF(lwp) WRITE(numout,*) 
    365          DO_2D( 0, 0, 0, 1 ) 
    366             zse3t = 1._wp / e3t(ji,jj,1,Kmm) 
    367             ptr(ji,jj,1,jn,Krhs) = ptr(ji,jj,1,jn,Krhs) +  sbc_trc(ji,jj,jn) * zse3t 
    368          END_2D 
    369          ! 
    370          IF( l_trdtrc ) THEN 
    371             ztrtrd(:,:,:) = ptr(:,:,:,jn,Krhs) - ztrtrd(:,:,:) 
    372             CALL trd_tra( kt, Kmm, Krhs, 'TRC', jn, jptra_nsr, ztrtrd ) 
    373          END IF 
    374          ! 
    375       END DO 
    376376      ! 
    377377      IF( sn_cfctl%l_prttrc )   THEN 
     
    379379                                           CALL prt_ctl( tab4d_1=ptr(:,:,:,:,Krhs), mask1=tmask, clinfo=ctrcnm, clinfo3='trd' ) 
    380380      ENDIF 
    381       ! 
    382       IF( l_trdtrc )  DEALLOCATE( ztrtrd ) 
    383381      ! 
    384382      IF( ln_timing )   CALL timing_stop('trc_sbc_RK3') 
Note: See TracChangeset for help on using the changeset viewer.