MODULE tranpc !!============================================================================== !! *** MODULE tranpc *** !! Ocean active tracers: non penetrative convective adjustment scheme !!============================================================================== !! History : 1.0 ! 1990-09 (G. Madec) Original code !! ! 1996-01 (G. Madec) statement function for e3 !! NEMO 1.0 ! 2002-06 (G. Madec) free form F90 !! 3.0 ! 2008-06 (G. Madec) applied on ta, sa and called before tranxt in step.F90 !! 3.3 ! 2010-05 (C. Ethe, G. Madec) merge TRC-TRA !! 3.7 ! 2014-06 (L. Brodeau) new algorithm based on local Brunt-Vaisala freq. !!---------------------------------------------------------------------- !!---------------------------------------------------------------------- !! tra_npc : apply the non penetrative convection scheme !!---------------------------------------------------------------------- USE oce ! ocean dynamics and active tracers USE dom_oce ! ocean space and time domain USE phycst ! physical constants USE zdf_oce ! ocean vertical physics USE trd_oce ! ocean active tracer trends USE trdtra ! ocean active tracer trends USE eosbn2 ! equation of state (eos routine) ! USE lbclnk ! lateral boundary conditions (or mpp link) USE in_out_manager ! I/O manager USE lib_mpp ! MPP library USE wrk_nemo ! Memory Allocation USE timing ! Timing IMPLICIT NONE PRIVATE PUBLIC tra_npc ! routine called by step.F90 !! * Substitutions # include "domzgr_substitute.h90" # include "vectopt_loop_substitute.h90" !!---------------------------------------------------------------------- !! NEMO/OPA 3.6 , NEMO Consortium (2014) !! $Id$ !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) !!---------------------------------------------------------------------- CONTAINS SUBROUTINE tra_npc( kt ) !!---------------------------------------------------------------------- !! *** ROUTINE tranpc *** !! !! ** Purpose : Non-penetrative convective adjustment scheme. solve !! the static instability of the water column on after fields !! while conserving heat and salt contents. !! !! ** Method : updated algorithm able to deal with non-linear equation of state !! (i.e. static stability computed locally) !! !! ** Action : - (ta,sa) after the application od the npc scheme !! - send the associated trends for on-line diagnostics (l_trdtra=T) !! !! References : Madec, et al., 1991, JPO, 21, 9, 1349-1371. !!---------------------------------------------------------------------- INTEGER, INTENT(in) :: kt ! ocean time-step index ! INTEGER :: ji, jj, jk ! dummy loop indices INTEGER :: inpcc ! number of statically instable water column INTEGER :: jiter, ikbot, ik, ikup, ikdown, ilayer, ikm ! local integers LOGICAL :: l_bottom_reached, l_column_treated REAL(wp) :: zta, zalfa, zsum_temp, zsum_alfa, zaw, zdz, zsum_z REAL(wp) :: zsa, zbeta, zsum_sali, zsum_beta, zbw, zrw, z1_r2dt REAL(wp), POINTER, DIMENSION(:) :: zvn2 ! vertical profile of N2 at 1 given point... REAL(wp), POINTER, DIMENSION(:,:) :: zvts ! vertical profile of T and S at 1 given point... REAL(wp), POINTER, DIMENSION(:,:) :: zvab ! vertical profile of alpha and beta REAL(wp), POINTER, DIMENSION(:,:,:) :: zn2 ! N^2 REAL(wp), POINTER, DIMENSION(:,:,:,:) :: zab ! alpha and beta REAL(wp), POINTER, DIMENSION(:,:,:) :: ztrdt, ztrds ! 3D workspace ! !!LB debug: LOGICAL, PARAMETER :: l_LB_debug = .FALSE. INTEGER :: ilc1, jlc1, klc1, nncpu LOGICAL :: lp_monitor_point = .FALSE. !!LB debug. !!---------------------------------------------------------------------- ! IF( nn_timing == 1 ) CALL timing_start('tra_npc') ! IF( MOD( kt, nn_npc ) == 0 ) THEN ! CALL wrk_alloc( jpi, jpj, jpk, zn2 ) ! N2 CALL wrk_alloc( jpi, jpj, jpk, 2, zab ) ! Alpha and Beta CALL wrk_alloc( jpk, 2, zvts, zvab ) ! 1D column vector at point ji,jj CALL wrk_alloc( jpk, zvn2 ) ! 1D column vector at point ji,jj IF( l_trdtra ) THEN !* Save initial after fields CALL wrk_alloc( jpi, jpj, jpk, ztrdt, ztrds ) ztrdt(:,:,:) = tsa(:,:,:,jp_tem) ztrds(:,:,:) = tsa(:,:,:,jp_sal) ENDIF !LB debug: IF( lwp .AND. l_LB_debug ) THEN WRITE(numout,*) WRITE(numout,*) 'LOLO: entering tra_npc, kt, narea =', kt, narea ENDIF !LBdebug: Monitoring of 1 column subject to convection... IF( l_LB_debug ) THEN ! Location of 1 known convection spot to follow what's happening in the water column ilc1 = 54 ; jlc1 = 15 ; ! Labrador ORCA1 4x4 cpus: nncpu = 15 ; ! the CPU domain contains the convection spot !ilc1 = 14 ; jlc1 = 13 ; ! Labrador ORCA1 8x8 cpus: !nncpu = 54 ; ! the CPU domain contains the convection spot klc1 = mbkt(ilc1,jlc1) ! bottom of the ocean for debug point... ENDIF !LBdebug. CALL eos_rab( tsa, zab ) ! after alpha and beta CALL bn2 ( tsa, zab, zn2 ) ! after Brunt-Vaisala inpcc = 0 DO jj = 2, jpjm1 ! interior column only DO ji = fs_2, fs_jpim1 ! IF( tmask(ji,jj,2) == 1 ) THEN ! At least 2 ocean points ! ! consider one ocean column zvts(:,jp_tem) = tsa(ji,jj,:,jp_tem) ! temperature zvts(:,jp_sal) = tsa(ji,jj,:,jp_sal) ! salinity zvab(:,jp_tem) = zab(ji,jj,:,jp_tem) ! Alpha zvab(:,jp_sal) = zab(ji,jj,:,jp_sal) ! Beta zvn2(:) = zn2(ji,jj,:) ! N^2 IF( l_LB_debug ) THEN !LB debug: lp_monitor_point = .FALSE. IF( ( ji == ilc1 ).AND.( jj == jlc1 ) ) lp_monitor_point = .TRUE. ! writing only if on CPU domain where conv region is: lp_monitor_point = (narea == nncpu).AND.lp_monitor_point IF(lp_monitor_point) THEN WRITE(numout,*) '' ;WRITE(numout,*) '' ; WRITE(numout,'("Time step = ",i6.6," !!!")') kt WRITE(numout,'(" *** BEFORE anything, N^2 for point ",i3,",",i3,":" )') ji,jj DO jk = 1, klc1 WRITE(numout,*) jk, zvn2(jk) END DO WRITE(numout,*) ' ' ENDIF ENDIF !LB debug end ikbot = mbkt(ji,jj) ! ikbot: ocean bottom T-level ik = 1 ! because N2 is irrelevant at the surface level (will start at ik=2) ilayer = 0 jiter = 0 l_column_treated = .FALSE. DO WHILE ( .NOT. l_column_treated ) ! jiter = jiter + 1 IF( jiter >= 400 ) EXIT l_bottom_reached = .FALSE. DO WHILE ( .NOT. l_bottom_reached ) ik = ik + 1 !! Checking level ik for instability !! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ IF( zvn2(ik) < 0. ) THEN ! Instability found! ikm = ik ! first level whith negative N2 ilayer = ilayer + 1 ! yet another layer found.... IF(jiter == 1) inpcc = inpcc + 1 IF(l_LB_debug .AND. lp_monitor_point) & & WRITE(numout,*) 'Negative N2 at ik =', ikm, ' layer nb.', ilayer, & & ' inpcc =', inpcc !! Case we mix with upper regions where N2==0: !! All the points above ikup where N2 == 0 must also be mixed => we go !! upward to find a new ikup, where the layer doesn't have N2==0 ikup = ikm DO jk = ikm, 2, -1 ikup = ikup - 1 IF( (zvn2(jk-1) > 0.).OR.(ikup == 1) ) EXIT END DO ! adjusting ikup if the upper part of the unstable column was neutral (N2=0) IF((zvn2(ikup+1) == 0.).AND.(ikup /= 1)) ikup = ikup+1 ; IF( lp_monitor_point ) WRITE(numout,*) ' => ikup is =', ikup, ' layer nb.', ilayer zsum_temp = 0._wp zsum_sali = 0._wp zsum_alfa = 0._wp zsum_beta = 0._wp zsum_z = 0._wp DO jk = ikup, ikbot+1 ! Inside the instable (and overlying neutral) portion of the column ! IF(l_LB_debug .AND. lp_monitor_point) WRITE(numout,*) ' -> summing for jk =', jk ! zdz = fse3t(ji,jj,jk) zsum_temp = zsum_temp + zvts(jk,jp_tem)*zdz zsum_sali = zsum_sali + zvts(jk,jp_sal)*zdz zsum_alfa = zsum_alfa + zvab(jk,jp_tem)*zdz zsum_beta = zsum_beta + zvab(jk,jp_sal)*zdz zsum_z = zsum_z + zdz ! !! EXIT if we found the bottom of the unstable portion of the water column IF( (zvn2(jk+1) > 0.).OR.(jk == ikbot ).OR.((jk==ikm).AND.(zvn2(jk+1) == 0.)) ) EXIT END DO !ik = jk !LB remove? ikdown = jk ! for the current unstable layer, ikdown is the deepest point with a negative N2 IF(l_LB_debug .AND. lp_monitor_point) & & WRITE(numout,*) ' => ikdown =', ikdown, ' layer nb.', ilayer ! Mixing Temperature and salinity between ikup and ikdown: zta = zsum_temp/zsum_z zsa = zsum_sali/zsum_z zalfa = zsum_alfa/zsum_z zbeta = zsum_beta/zsum_z IF(l_LB_debug .AND. lp_monitor_point) THEN WRITE(numout,*) ' => Mean temp. in that portion =', zta WRITE(numout,*) ' => Mean sali. in that portion =', zsa WRITE(numout,*) ' => Mean Alpha in that portion =', zalfa WRITE(numout,*) ' => Mean Beta in that portion =', zbeta ENDIF !! Homogenaizing the temperature, salinity, alpha and beta in this portion of the column DO jk = ikup, ikdown zvts(jk,jp_tem) = zta zvts(jk,jp_sal) = zsa zvab(jk,jp_tem) = zalfa zvab(jk,jp_sal) = zbeta END DO ! !! Before updating N2, it is possible that another unstable !! layer exists underneath the one we just homogeneized! ik = ikdown ! ENDIF ! IF( zvn2(ik+1) < 0. ) THEN ! IF( ik == ikbot ) l_bottom_reached = .TRUE. ! END DO ! DO WHILE ( .NOT. l_bottom_reached ) IF( ik /= ikbot ) STOP 'ERROR: tranpc.F90 => PROBLEM #1' ! ******* At this stage ik == ikbot ! ******* IF( ilayer > 0 ) THEN !! least an unstable layer has been found !! Temperature, Salinity, Alpha and Beta have been homogenized in the unstable portion !! => Need to re-compute N2! will use Alpha and Beta! ! DO jk = ikup+1, ikdown+1 ! we must go 1 point deeper than ikdown! !! Doing exactly as in eosbn2.F90: !! * Except that we only are interested in the sign of N2 !!! !! => just considering the vertical gradient of density zrw = (fsdepw(ji,jj,jk ) - fsdept(ji,jj,jk)) & & / (fsdept(ji,jj,jk-1) - fsdept(ji,jj,jk)) zaw = zvab(jk,jp_tem) * (1._wp - zrw) + zvab(jk-1,jp_tem) * zrw zbw = zvab(jk,jp_sal) * (1._wp - zrw) + zvab(jk-1,jp_sal) * zrw !zvn2(jk) = grav*( zaw * ( zvts(jk-1,jp_tem) - zvts(jk,jp_tem) ) & ! & - zbw * ( zvts(jk-1,jp_sal) - zvts(jk,jp_sal) ) ) & ! & / fse3w(ji,jj,jk) * tmask(ji,jj,jk) zvn2(jk) = ( zaw * ( zvts(jk-1,jp_tem) - zvts(jk,jp_tem) ) & & - zbw * ( zvts(jk-1,jp_sal) - zvts(jk,jp_sal) ) ) END DO IF(l_LB_debug .AND. lp_monitor_point) THEN WRITE(numout, '(" *** After iteration #",i3.3,", N^2 for point ",i3,",",i3,":" )') & & jiter, ji,jj DO jk = 1, klc1 WRITE(numout,*) jk, zvn2(jk) END DO WRITE(numout,*) ' ' ENDIF ik = 1 ! starting again at the surface for the next iteration ilayer = 0 ENDIF ! IF( ik >= ikbot ) THEN IF(l_LB_debug .AND. lp_monitor_point) WRITE(numout,*) ' --- exiting jiter loop ---' l_column_treated = .TRUE. ENDIF ! END DO ! DO WHILE ( .NOT. l_column_treated ) !! Updating tsa: tsa(ji,jj,:,jp_tem) = zvts(:,jp_tem) tsa(ji,jj,:,jp_sal) = zvts(:,jp_sal) !! lolo: Should we update something else???? !! => like alpha and beta? IF(l_LB_debug .AND. lp_monitor_point) WRITE(numout,*) '' ENDIF ! IF( tmask(ji,jj,3) == 1 ) THEN END DO ! ji END DO ! jj ! IF( l_trdtra ) THEN ! send the Non penetrative mixing trends for diagnostic z1_r2dt = 1._wp / (2._wp * rdt) ztrdt(:,:,:) = ( tsa(:,:,:,jp_tem) - ztrdt(:,:,:) ) * z1_r2dt ztrds(:,:,:) = ( tsa(:,:,:,jp_sal) - ztrds(:,:,:) ) * z1_r2dt CALL trd_tra( kt, 'TRA', jp_tem, jptra_npc, ztrdt ) CALL trd_tra( kt, 'TRA', jp_sal, jptra_npc, ztrds ) CALL wrk_dealloc( jpi, jpj, jpk, ztrdt, ztrds ) ENDIF ! CALL lbc_lnk( tsa(:,:,:,jp_tem), 'T', 1. ) ; CALL lbc_lnk( tsa(:,:,:,jp_sal), 'T', 1. ) ! IF(lwp) THEN WRITE(numout,*) 'LOLO: exiting tra_npc, kt =', kt WRITE(numout,*)' => number of statically instable water column : ',inpcc WRITE(numout,*) '' ; WRITE(numout,*) '' ENDIF ! CALL wrk_dealloc(jpi, jpj, jpk, zn2 ) CALL wrk_dealloc(jpi, jpj, jpk, 2, zab ) CALL wrk_dealloc(jpk, zvn2 ) CALL wrk_dealloc(jpk, 2, zvts, zvab ) ! ENDIF ! IF( MOD( kt, nn_npc ) == 0 ) THEN ! IF( nn_timing == 1 ) CALL timing_stop('tra_npc') ! END SUBROUTINE tra_npc !!====================================================================== END MODULE tranpc