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 14993 for NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE/TRA/trasbc.F90 – NEMO

Ignore:
Timestamp:
2021-06-15T00:35:18+02:00 (18 months ago)
Author:
techene
Message:

#2605 adapt and optimize mass flux forcing for RK3

File:
1 edited

Legend:

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

    r14943 r14993  
    3939   PRIVATE 
    4040 
    41    PUBLIC   tra_sbc   ! routine called by step.F90 
     41   PUBLIC   tra_sbc       ! routine called by step.F90 
     42   PUBLIC   tra_sbc_RK3   ! routine called by stprk3_.F90 
    4243 
    4344   !! * Substitutions 
     
    260261   END SUBROUTINE tra_sbc 
    261262 
     263 
     264   SUBROUTINE tra_sbc_RK3 ( kt, Kmm, pts, Krhs, kstg ) 
     265      !!---------------------------------------------------------------------- 
     266      !!                  ***  ROUTINE tra_sbc_RK3  *** 
     267      !! 
     268      !! ** Purpose :   Compute the tracer surface boundary condition trend of 
     269      !!      (flux through the interface, concentration/dilution effect) 
     270      !!      and add it to the general trend of tracer equations. 
     271      !! 
     272      !! ** Method :   The (air+ice)-sea flux has two components: 
     273      !!      (1) Fext, external forcing (i.e. flux through the (air+ice)-sea interface); 
     274      !!      (2) Fwe , tracer carried with the water that is exchanged with air+ice. 
     275      !!               The input forcing fields (emp, rnf, sfx) contain Fext+Fwe, 
     276      !!             they are simply added to the tracer trend (ts(Krhs)). 
     277      !!               In linear free surface case (ln_linssh=T), the volume of the 
     278      !!             ocean does not change with the water exchanges at the (air+ice)-sea 
     279      !!             interface. Therefore another term has to be added, to mimic the 
     280      !!             concentration/dilution effect associated with water exchanges. 
     281      !! 
     282      !! ** Action  : - Update ts(Krhs) with the surface boundary condition trend 
     283      !!              - send trends to trdtra module for further diagnostics(l_trdtra=T) 
     284      !!---------------------------------------------------------------------- 
     285      INTEGER                                  , INTENT(in   ) ::   kt, Kmm, Krhs   ! ocean time-step and time-level indices 
     286      INTEGER                                  , INTENT(in   ) ::   kstg            ! RK3 stage index 
     287      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts,jpt), INTENT(inout) ::   pts             ! active tracers and RHS of tracer Eq. 
     288      ! 
     289      INTEGER  ::   ji, jj, jk, jn               ! dummy loop indices 
     290      REAL(wp) ::   z1_rho0_e3t, zdep, ztim    ! local scalar 
     291      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::  ztrdt, ztrds 
     292      !!---------------------------------------------------------------------- 
     293      ! 
     294      IF( ln_timing )   CALL timing_start('tra_sbc_RK3') 
     295      ! 
     296      IF( ntile == 0 .OR. ntile == 1 )  THEN                       ! Do only on the first tile 
     297         ! 
     298         IF( kt == nit000 ) THEN 
     299            IF(lwp) WRITE(numout,*) 
     300            IF(lwp) WRITE(numout,*) 'tra_sbc_RK3 : TRAcer Surface Boundary Condition' 
     301            IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ ' 
     302         ENDIF 
     303         ! 
     304         IF( l_trdtra ) THEN                    !* Save ta and sa trends 
     305            ALLOCATE( ztrdt(jpi,jpj,jpk), ztrds(jpi,jpj,jpk) ) 
     306            ztrdt(:,:,:) = pts(:,:,:,jp_tem,Krhs) 
     307            ztrds(:,:,:) = pts(:,:,:,jp_sal,Krhs) 
     308         ENDIF 
     309         ! 
     310      ENDIF 
     311      ! 
     312             
     313!!gm  This should be moved into sbcmod.F90 module ? (especially now that ln_traqsr is read in namsbc namelist) 
     314      IF( .NOT.ln_traqsr  .AND. kstg == 1) THEN     ! no solar radiation penetration 
     315         DO_2D( 0, 0, 0, 0 ) 
     316            qns(ji,jj) = qns(ji,jj) + qsr(ji,jj)         ! total heat flux in qns 
     317            qsr(ji,jj) = 0._wp                           ! qsr set to zero 
     318         END_2D 
     319      ENDIF 
     320 
     321      !---------------------------------------- 
     322      !        EMP, SFX and QNS effects 
     323      !---------------------------------------- 
     324      !                             !==  update tracer trend  ==! 
     325      SELECT CASE( kstg ) 
     326      ! 
     327      CASE( 1 , 2 )                       !=  stage 1 and 2  =!   only in non linear ssh 
     328         ! 
     329         IF( .NOT.ln_linssh ) THEN           !* only heat and salt fluxes associated with mass fluxes 
     330            DO_2D( 0, 0, 0, 0 ) 
     331               z1_rho0_e3t = r1_rho0 / e3t(ji,jj,1,Kmm) 
     332               pts(ji,jj,1,jp_tem,Krhs) = pts(ji,jj,1,jp_tem,Krhs) - emp(ji,jj)*pts(ji,jj,1,jp_tem,Kmm) * z1_rho0_e3t 
     333               pts(ji,jj,1,jp_sal,Krhs) = pts(ji,jj,1,jp_sal,Krhs) - emp(ji,jj)*pts(ji,jj,1,jp_sal,Kmm) * z1_rho0_e3t 
     334            END_2D 
     335         ENDIF 
     336         ! 
     337      CASE( 3 ) 
     338         ! 
     339         IF( ln_linssh ) THEN                !* linear free surface 
     340            DO_2D( 0, 0, 0, 0 ) 
     341               z1_rho0_e3t = r1_rho0 / e3t(ji,jj,1,Kmm) 
     342               pts(ji,jj,1,jp_tem,Krhs) = pts(ji,jj,1,jp_tem,Krhs) + (  r1_rcp * qns(ji,jj)   &                                ! non solar heat flux 
     343                  &                                                   +          emp(ji,jj)*pts(ji,jj,1,jp_tem,Kmm)  ) * z1_rho0_e3t  ! add concentration/dilution effect due to constant volume cell 
     344               pts(ji,jj,1,jp_sal,Krhs) = pts(ji,jj,1,jp_sal,Krhs) + (           sfx(ji,jj)    &                               ! salt flux due to freezing/melting 
     345                  &                                                   +          emp(ji,jj)*pts(ji,jj,1,jp_sal,Kmm)  ) * z1_rho0_e3t  ! add concentration/dilution effect due to constant volume cell 
     346            END_2D 
     347            IF( ntile == 0 .OR. ntile == nijtile ) THEN             ! Do only on the last tile 
     348               IF( iom_use('emp_x_sst') )   CALL iom_put( "emp_x_sst", emp (:,:) * pts(:,:,1,jp_tem,Kmm) ) 
     349               IF( iom_use('emp_x_sss') )   CALL iom_put( "emp_x_sss", emp (:,:) * pts(:,:,1,jp_sal,Kmm) ) 
     350            ENDIF 
     351         ELSE 
     352            DO_2D( 0, 0, 0, 0 ) 
     353               z1_rho0_e3t = r1_rho0 / e3t(ji,jj,1,Kmm) 
     354               pts(ji,jj,1,jp_tem,Krhs) = pts(ji,jj,1,jp_tem,Krhs) +  r1_rcp * qns(ji,jj) * z1_rho0_e3t 
     355               pts(ji,jj,1,jp_sal,Krhs) = pts(ji,jj,1,jp_sal,Krhs) +           sfx(ji,jj) * z1_rho0_e3t 
     356            END_2D 
     357         ENDIF 
     358      END SELECT 
     359      ! 
     360      ! 
     361      !---------------------------------------- 
     362      !        River Runoff effects 
     363      !---------------------------------------- 
     364      ! 
     365      IF( ln_rnf ) THEN         ! input of heat and salt due to river runoff 
     366         DO_2D( 0, 0, 0, 0 ) 
     367            IF( rnf(ji,jj) /= 0._wp ) THEN 
     368               zdep = 1._wp / h_rnf(ji,jj) 
     369               DO jk = 1, nk_rnf(ji,jj) 
     370                                     pts(ji,jj,jk,jp_tem,Krhs) = pts(ji,jj,jk,jp_tem,Krhs)  + rnf_tsc(ji,jj,jp_tem) * zdep 
     371                  IF( ln_rnf_sal )   pts(ji,jj,jk,jp_sal,Krhs) = pts(ji,jj,jk,jp_sal,Krhs)  + rnf_tsc(ji,jj,jp_sal) * zdep 
     372               END DO 
     373            ENDIF 
     374         END_2D 
     375      ENDIF 
     376      ! 
     377      IF( ntile == 0 .OR. ntile == nijtile )  THEN                ! Do only on the last tile 
     378         IF( iom_use('rnf_x_sst') )   CALL iom_put( "rnf_x_sst", rnf*pts(:,:,1,jp_tem,Kmm) )   ! runoff term on sst 
     379         IF( iom_use('rnf_x_sss') )   CALL iom_put( "rnf_x_sss", rnf*pts(:,:,1,jp_sal,Kmm) )   ! runoff term on sss 
     380      ENDIF 
     381 
     382#if defined key_asminc 
     383      ! 
     384      !---------------------------------------- 
     385      !        Assmilation effects 
     386      !---------------------------------------- 
     387      ! 
     388      IF( ln_sshinc .AND. kstg == 3 ) THEN         ! input of heat and salt due to assimilation 
     389         ! 
     390         IF( ln_linssh ) THEN 
     391            DO_2D( 0, 0, 0, 0 ) 
     392               ztim = ssh_iau(ji,jj) / e3t(ji,jj,1,Kmm) 
     393               pts(ji,jj,1,jp_tem,Krhs) = pts(ji,jj,1,jp_tem,Krhs) + pts(ji,jj,1,jp_tem,Kmm) * ztim 
     394               pts(ji,jj,1,jp_sal,Krhs) = pts(ji,jj,1,jp_sal,Krhs) + pts(ji,jj,1,jp_sal,Kmm) * ztim 
     395            END_2D 
     396         ELSE 
     397            DO_2D( 0, 0, 0, 0 ) 
     398               ztim = ssh_iau(ji,jj) / ( ht(ji,jj) + 1. - ssmask(ji, jj) ) 
     399               pts(ji,jj,:,jp_tem,Krhs) = pts(ji,jj,:,jp_tem,Krhs) + pts(ji,jj,:,jp_tem,Kmm) * ztim 
     400               pts(ji,jj,:,jp_sal,Krhs) = pts(ji,jj,:,jp_sal,Krhs) + pts(ji,jj,:,jp_sal,Kmm) * ztim 
     401            END_2D 
     402         ENDIF 
     403         ! 
     404      ENDIF 
     405      ! 
     406#endif 
     407      ! 
     408      IF( l_trdtra )   THEN                      ! save the horizontal diffusive trends for further diagnostics 
     409         IF( ntile == 0 .OR. ntile == nijtile )  THEN 
     410            ztrdt(:,:,:) = pts(:,:,:,jp_tem,Krhs) - ztrdt(:,:,:) 
     411            ztrds(:,:,:) = pts(:,:,:,jp_sal,Krhs) - ztrds(:,:,:) 
     412            CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_tem, jptra_nsr, ztrdt ) 
     413            CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_sal, jptra_nsr, ztrds ) 
     414            DEALLOCATE( ztrdt , ztrds ) 
     415         ENDIF 
     416      ENDIF 
     417      ! 
     418      IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab3d_1=pts(:,:,:,jp_tem,Krhs), clinfo1=' sbc  - Ta: ', mask1=tmask,   & 
     419         &                                  tab3d_2=pts(:,:,:,jp_sal,Krhs), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' ) 
     420      ! 
     421      IF( ln_timing )   CALL timing_stop('tra_sbc_RK3') 
     422      ! 
     423   END SUBROUTINE tra_sbc_RK3 
     424 
    262425   !!====================================================================== 
    263426END MODULE trasbc 
Note: See TracChangeset for help on using the changeset viewer.