Changeset 10874 for NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA
- Timestamp:
- 2019-04-15T15:57:37+02:00 (5 years ago)
- Location:
- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA
- Files:
-
- 16 edited
-
traadv.F90 (modified) (6 diffs)
-
traadv_cen.F90 (modified) (11 diffs)
-
traadv_fct.F90 (modified) (21 diffs)
-
traadv_mus.F90 (modified) (11 diffs)
-
traadv_qck.F90 (modified) (20 diffs)
-
traadv_ubs.F90 (modified) (19 diffs)
-
trabbc.F90 (modified) (3 diffs)
-
trabbl.F90 (modified) (17 diffs)
-
tradmp.F90 (modified) (7 diffs)
-
traldf.F90 (modified) (3 diffs)
-
traldf_iso.F90 (modified) (11 diffs)
-
traldf_lap_blp.F90 (modified) (11 diffs)
-
traldf_triad.F90 (modified) (18 diffs)
-
traqsr.F90 (modified) (8 diffs)
-
trasbc.F90 (modified) (14 diffs)
-
trazdf.F90 (modified) (13 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traadv.F90
r10806 r10874 75 75 CONTAINS 76 76 77 SUBROUTINE tra_adv( kt , ktlev1, ktlev2, ktlev3, kt2lev, pts_rhs)77 SUBROUTINE tra_adv( kt ) 78 78 !!---------------------------------------------------------------------- 79 79 !! *** ROUTINE tra_adv *** … … 81 81 !! ** Purpose : compute the ocean tracer advection trend. 82 82 !! 83 !! ** Method : - Update (pu_rhs,pv_rhs) with the advection term following nadv 84 !!---------------------------------------------------------------------- 85 INTEGER, INTENT(in) :: kt ! ocean time-step index 86 INTEGER, INTENT(in) :: ktlev1, ktlev2, ktlev3 ! time level indices for 3-time-level source terms 87 INTEGER, INTENT(in) :: kt2lev ! time level index for 2-time-level source terms 88 REAL(wp), INTENT( inout), DIMENSION(jpi,jpj,jpk,jpts) :: pts_rhs ! temperature and salinity trends 83 !! ** Method : - Update (ua,va) with the advection term following nadv 84 !!---------------------------------------------------------------------- 85 INTEGER, INTENT(in) :: kt ! ocean time-step index 89 86 ! 90 87 INTEGER :: jk ! dummy loop index … … 106 103 IF( ln_wave .AND. ln_sdw ) THEN 107 104 DO jk = 1, jpkm1 ! eulerian transport + Stokes Drift 108 zun(:,:,jk) = e2u (:,:) * e3u (:,:,jk,ktlev2) * ( uu(:,:,jk,ktlev2) + usd(:,:,jk) )109 zvn(:,:,jk) = e1v (:,:) * e3v (:,:,jk,ktlev2) * ( vv(:,:,jk,ktlev2) + vsd(:,:,jk) )110 zwn(:,:,jk) = e1e2t(:,:) * ( w w(:,:,jk) + wsd(:,:,jk) )105 zun(:,:,jk) = e2u (:,:) * e3u_n(:,:,jk) * ( un(:,:,jk) + usd(:,:,jk) ) 106 zvn(:,:,jk) = e1v (:,:) * e3v_n(:,:,jk) * ( vn(:,:,jk) + vsd(:,:,jk) ) 107 zwn(:,:,jk) = e1e2t(:,:) * ( wn(:,:,jk) + wsd(:,:,jk) ) 111 108 END DO 112 109 ELSE 113 110 DO jk = 1, jpkm1 114 zun(:,:,jk) = e2u (:,:) * e3u (:,:,jk,ktlev2) * uu(:,:,jk,ktlev2) ! eulerian transport only115 zvn(:,:,jk) = e1v (:,:) * e3v (:,:,jk,ktlev2) * vv(:,:,jk,ktlev2)116 zwn(:,:,jk) = e1e2t(:,:) * w w(:,:,jk)111 zun(:,:,jk) = e2u (:,:) * e3u_n(:,:,jk) * un(:,:,jk) ! eulerian transport only 112 zvn(:,:,jk) = e1v (:,:) * e3v_n(:,:,jk) * vn(:,:,jk) 113 zwn(:,:,jk) = e1e2t(:,:) * wn(:,:,jk) 117 114 END DO 118 115 ENDIF … … 142 139 IF( l_trdtra ) THEN !* Save ta and sa trends 143 140 ALLOCATE( ztrdt(jpi,jpj,jpk), ztrds(jpi,jpj,jpk) ) 144 ztrdt(:,:,:) = pts_rhs(:,:,:,jp_tem)145 ztrds(:,:,:) = pts_rhs(:,:,:,jp_sal)141 ztrdt(:,:,:) = tsa(:,:,:,jp_tem) 142 ztrds(:,:,:) = tsa(:,:,:,jp_sal) 146 143 ENDIF 147 144 ! … … 149 146 ! 150 147 CASE ( np_CEN ) ! Centered scheme : 2nd / 4th order 151 CALL tra_adv_cen ( kt, nit000, ktlev2, 'TRA', zun, zvn, zwn , ts(:,:,:,:,ktlev2), pts_rhs, jpts, nn_cen_h, nn_cen_v )148 CALL tra_adv_cen ( kt, nit000, 'TRA', zun, zvn, zwn , tsn, tsa, jpts, nn_cen_h, nn_cen_v ) 152 149 CASE ( np_FCT ) ! FCT scheme : 2nd / 4th order 153 CALL tra_adv_fct ( kt, nit000, ktlev1, ktlev2, ktlev3, 'TRA', r2dt, zun, zvn, zwn, & 154 & ts(:,:,:,:,ktlev1), ts(:,:,:,:,ktlev2), pts_rhs, jpts, nn_fct_h, nn_fct_v ) 150 CALL tra_adv_fct ( kt, nit000, 'TRA', r2dt, zun, zvn, zwn, tsb, tsn, tsa, jpts, nn_fct_h, nn_fct_v ) 155 151 CASE ( np_MUS ) ! MUSCL 156 CALL tra_adv_mus ( kt, nit000, ktlev2, kt2lev, 'TRA', r2dt, zun, zvn, zwn, ts(:,:,:,:,ktlev1), pts_rhs, jpts , ln_mus_ups )152 CALL tra_adv_mus ( kt, nit000, 'TRA', r2dt, zun, zvn, zwn, tsb, tsa, jpts , ln_mus_ups ) 157 153 CASE ( np_UBS ) ! UBS 158 CALL tra_adv_ubs ( kt, nit000, ktlev2, 'TRA', r2dt, zun, zvn, zwn, ts(:,:,:,:,ktlev1), ts(:,:,:,:,ktlev2), pts_rhs, jpts , nn_ubs_v )154 CALL tra_adv_ubs ( kt, nit000, 'TRA', r2dt, zun, zvn, zwn, tsb, tsn, tsa, jpts , nn_ubs_v ) 159 155 CASE ( np_QCK ) ! QUICKEST 160 CALL tra_adv_qck ( kt, nit000, ktlev2, 'TRA', r2dt, zun, zvn, zwn, ts(:,:,:,:,ktlev1), ts(:,:,:,:,ktlev2), pts_rhs, jpts )156 CALL tra_adv_qck ( kt, nit000, 'TRA', r2dt, zun, zvn, zwn, tsb, tsn, tsa, jpts ) 161 157 ! 162 158 END SELECT … … 164 160 IF( l_trdtra ) THEN ! save the advective trends for further diagnostics 165 161 DO jk = 1, jpkm1 166 ztrdt(:,:,jk) = pts_rhs(:,:,jk,jp_tem) - ztrdt(:,:,jk)167 ztrds(:,:,jk) = pts_rhs(:,:,jk,jp_sal) - ztrds(:,:,jk)162 ztrdt(:,:,jk) = tsa(:,:,jk,jp_tem) - ztrdt(:,:,jk) 163 ztrds(:,:,jk) = tsa(:,:,jk,jp_sal) - ztrds(:,:,jk) 168 164 END DO 169 165 CALL trd_tra( kt, 'TRA', jp_tem, jptra_totad, ztrdt ) -
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traadv_cen.F90
r10806 r10874 44 44 CONTAINS 45 45 46 SUBROUTINE tra_adv_cen( kt, kit000, ktlev, cdtype, pu, pv, pw, &47 & pt, pt_rhs, kjpt, kn_cen_h, kn_cen_v )46 SUBROUTINE tra_adv_cen( kt, kit000, cdtype, pun, pvn, pwn, & 47 & ptn, pta, kjpt, kn_cen_h, kn_cen_v ) 48 48 !!---------------------------------------------------------------------- 49 49 !! *** ROUTINE tra_adv_cen *** … … 59 59 !! = 4 ==>> 4th order COMPACT scheme - - 60 60 !! 61 !! ** Action : - update pt _rhswith the now advective tracer trends61 !! ** Action : - update pta with the now advective tracer trends 62 62 !! - send trends to trdtra module for further diagnostcs (l_trdtra=T) 63 63 !! - htr_adv, str_adv : poleward advective heat and salt transport (ln_diaptr=T) … … 65 65 INTEGER , INTENT(in ) :: kt ! ocean time-step index 66 66 INTEGER , INTENT(in ) :: kit000 ! first time step index 67 INTEGER , INTENT(in ) :: ktlev ! time level index for source terms68 67 CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator) 69 68 INTEGER , INTENT(in ) :: kjpt ! number of tracers 70 69 INTEGER , INTENT(in ) :: kn_cen_h ! =2/4 (2nd or 4th order scheme) 71 70 INTEGER , INTENT(in ) :: kn_cen_v ! =2/4 (2nd or 4th order scheme) 72 REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pu , pv, pw! 3 ocean velocity components73 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: pt ! now tracer fields74 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pt _rhs! tracer trend71 REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pun, pvn, pwn ! 3 ocean velocity components 72 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: ptn ! now tracer fields 73 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pta ! tracer trend 75 74 ! 76 75 INTEGER :: ji, jj, jk, jn ! dummy loop indices … … 107 106 DO jj = 1, jpjm1 108 107 DO ji = 1, fs_jpim1 ! vector opt. 109 zwx(ji,jj,jk) = 0.5_wp * pu (ji,jj,jk) * ( pt(ji,jj,jk,jn) + pt(ji+1,jj ,jk,jn) )110 zwy(ji,jj,jk) = 0.5_wp * pv (ji,jj,jk) * ( pt(ji,jj,jk,jn) + pt(ji ,jj+1,jk,jn) )108 zwx(ji,jj,jk) = 0.5_wp * pun(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji+1,jj ,jk,jn) ) 109 zwy(ji,jj,jk) = 0.5_wp * pvn(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji ,jj+1,jk,jn) ) 111 110 END DO 112 111 END DO … … 119 118 DO jj = 2, jpjm1 120 119 DO ji = fs_2, fs_jpim1 ! vector opt. 121 ztu(ji,jj,jk) = ( pt (ji+1,jj ,jk,jn) - pt(ji,jj,jk,jn) ) * umask(ji,jj,jk)122 ztv(ji,jj,jk) = ( pt (ji ,jj+1,jk,jn) - pt(ji,jj,jk,jn) ) * vmask(ji,jj,jk)120 ztu(ji,jj,jk) = ( ptn(ji+1,jj ,jk,jn) - ptn(ji,jj,jk,jn) ) * umask(ji,jj,jk) 121 ztv(ji,jj,jk) = ( ptn(ji ,jj+1,jk,jn) - ptn(ji,jj,jk,jn) ) * vmask(ji,jj,jk) 123 122 END DO 124 123 END DO … … 129 128 DO jj = 2, jpjm1 130 129 DO ji = 1, fs_jpim1 ! vector opt. 131 zC2t_u = pt (ji,jj,jk,jn) + pt(ji+1,jj ,jk,jn) ! C2 interpolation of T at u- & v-points (x2)132 zC2t_v = pt (ji,jj,jk,jn) + pt(ji ,jj+1,jk,jn)130 zC2t_u = ptn(ji,jj,jk,jn) + ptn(ji+1,jj ,jk,jn) ! C2 interpolation of T at u- & v-points (x2) 131 zC2t_v = ptn(ji,jj,jk,jn) + ptn(ji ,jj+1,jk,jn) 133 132 ! ! C4 interpolation of T at u- & v-points (x2) 134 133 zC4t_u = zC2t_u + r1_6 * ( ztu(ji-1,jj,jk) - ztu(ji+1,jj,jk) ) 135 134 zC4t_v = zC2t_v + r1_6 * ( ztv(ji,jj-1,jk) - ztv(ji,jj+1,jk) ) 136 135 ! ! C4 fluxes 137 zwx(ji,jj,jk) = 0.5_wp * pu (ji,jj,jk) * zC4t_u138 zwy(ji,jj,jk) = 0.5_wp * pv (ji,jj,jk) * zC4t_v136 zwx(ji,jj,jk) = 0.5_wp * pun(ji,jj,jk) * zC4t_u 137 zwy(ji,jj,jk) = 0.5_wp * pvn(ji,jj,jk) * zC4t_v 139 138 END DO 140 139 END DO … … 151 150 DO jj = 2, jpjm1 152 151 DO ji = fs_2, fs_jpim1 ! vector opt. 153 zwz(ji,jj,jk) = 0.5 * pw (ji,jj,jk) * ( pt(ji,jj,jk,jn) + pt(ji,jj,jk-1,jn) ) * wmask(ji,jj,jk)152 zwz(ji,jj,jk) = 0.5 * pwn(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji,jj,jk-1,jn) ) * wmask(ji,jj,jk) 154 153 END DO 155 154 END DO … … 157 156 ! 158 157 CASE( 4 ) !* 4th order compact 159 CALL interp_4th_cpt( pt (:,:,:,jn) , ztw ) ! ztw = interpolated value of T at w-point158 CALL interp_4th_cpt( ptn(:,:,:,jn) , ztw ) ! ztw = interpolated value of T at w-point 160 159 DO jk = 2, jpkm1 161 160 DO jj = 2, jpjm1 162 161 DO ji = fs_2, fs_jpim1 163 zwz(ji,jj,jk) = pw (ji,jj,jk) * ztw(ji,jj,jk) * wmask(ji,jj,jk)162 zwz(ji,jj,jk) = pwn(ji,jj,jk) * ztw(ji,jj,jk) * wmask(ji,jj,jk) 164 163 END DO 165 164 END DO … … 172 171 DO jj = 1, jpj 173 172 DO ji = 1, jpi 174 zwz(ji,jj, mikt(ji,jj) ) = pw (ji,jj,mikt(ji,jj)) * pt(ji,jj,mikt(ji,jj),jn)173 zwz(ji,jj, mikt(ji,jj) ) = pwn(ji,jj,mikt(ji,jj)) * ptn(ji,jj,mikt(ji,jj),jn) 175 174 END DO 176 175 END DO 177 176 ELSE ! no ice-shelf cavities (only ocean surface) 178 zwz(:,:,1) = pw (:,:,1) * pt(:,:,1,jn)177 zwz(:,:,1) = pwn(:,:,1) * ptn(:,:,1,jn) 179 178 ENDIF 180 179 ENDIF … … 183 182 DO jj = 2, jpjm1 184 183 DO ji = fs_2, fs_jpim1 ! vector opt. 185 pt _rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn) &184 pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) & 186 185 & - ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) & 187 186 & + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) & 188 & + zwz(ji,jj,jk) - zwz(ji ,jj ,jk+1) ) * r1_e1e2t(ji,jj) / e3t (ji,jj,jk,ktlev)187 & + zwz(ji,jj,jk) - zwz(ji ,jj ,jk+1) ) * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) 189 188 END DO 190 189 END DO … … 192 191 ! ! trend diagnostics 193 192 IF( l_trd ) THEN 194 CALL trd_tra( kt, cdtype, jn, jptra_xad, zwx, pu , pt(:,:,:,jn) )195 CALL trd_tra( kt, cdtype, jn, jptra_yad, zwy, pv , pt(:,:,:,jn) )196 CALL trd_tra( kt, cdtype, jn, jptra_zad, zwz, pw , pt(:,:,:,jn) )193 CALL trd_tra( kt, cdtype, jn, jptra_xad, zwx, pun, ptn(:,:,:,jn) ) 194 CALL trd_tra( kt, cdtype, jn, jptra_yad, zwy, pvn, ptn(:,:,:,jn) ) 195 CALL trd_tra( kt, cdtype, jn, jptra_zad, zwz, pwn, ptn(:,:,:,jn) ) 197 196 END IF 198 197 ! ! "Poleward" heat and salt transports -
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traadv_fct.F90
r10806 r10874 52 52 CONTAINS 53 53 54 SUBROUTINE tra_adv_fct( kt, kit000, ktlev1, ktlev2, ktlev3, cdtype, p2dt, pu, pv, pw, &55 & pt_lev1, pt_lev2, pt_rhs, kjpt, kn_fct_h, kn_fct_v )54 SUBROUTINE tra_adv_fct( kt, kit000, cdtype, p2dt, pun, pvn, pwn, & 55 & ptb, ptn, pta, kjpt, kn_fct_h, kn_fct_v ) 56 56 !!---------------------------------------------------------------------- 57 57 !! *** ROUTINE tra_adv_fct *** … … 65 65 !! - corrected flux (monotonic correction) 66 66 !! 67 !! ** Action : - update pt _rhswith the now advective tracer trends67 !! ** Action : - update pta with the now advective tracer trends 68 68 !! - send trends to trdtra module for further diagnostics (l_trdtra=T) 69 69 !! - htr_adv, str_adv : poleward advective heat and salt transport (ln_diaptr=T) 70 70 !!---------------------------------------------------------------------- 71 71 INTEGER , INTENT(in ) :: kt ! ocean time-step index 72 INTEGER , INTENT(in ) :: ktlev1, ktlev2, ktlev3 ! time level indices for source terms73 72 INTEGER , INTENT(in ) :: kit000 ! first time step index 74 73 CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator) … … 77 76 INTEGER , INTENT(in ) :: kn_fct_v ! order of the FCT scheme (=2 or 4) 78 77 REAL(wp) , INTENT(in ) :: p2dt ! tracer time-step 79 REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pu , pv, pw! 3 ocean velocity components80 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: pt _lev1, pt_lev2! before and now tracer fields81 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pt _rhs! tracer trend78 REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pun, pvn, pwn ! 3 ocean velocity components 79 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: ptb, ptn ! before and now tracer fields 80 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pta ! tracer trend 82 81 ! 83 82 INTEGER :: ji, jj, jk, jn ! dummy loop indices … … 126 125 DO ji = 1, fs_jpim1 ! vector opt. 127 126 ! upstream scheme 128 zfp_ui = pu (ji,jj,jk) + ABS( pu(ji,jj,jk) )129 zfm_ui = pu (ji,jj,jk) - ABS( pu(ji,jj,jk) )130 zfp_vj = pv (ji,jj,jk) + ABS( pv(ji,jj,jk) )131 zfm_vj = pv (ji,jj,jk) - ABS( pv(ji,jj,jk) )132 zwx(ji,jj,jk) = 0.5 * ( zfp_ui * pt _lev1(ji,jj,jk,jn) + zfm_ui * pt_lev1(ji+1,jj ,jk,jn) )133 zwy(ji,jj,jk) = 0.5 * ( zfp_vj * pt _lev1(ji,jj,jk,jn) + zfm_vj * pt_lev1(ji ,jj+1,jk,jn) )127 zfp_ui = pun(ji,jj,jk) + ABS( pun(ji,jj,jk) ) 128 zfm_ui = pun(ji,jj,jk) - ABS( pun(ji,jj,jk) ) 129 zfp_vj = pvn(ji,jj,jk) + ABS( pvn(ji,jj,jk) ) 130 zfm_vj = pvn(ji,jj,jk) - ABS( pvn(ji,jj,jk) ) 131 zwx(ji,jj,jk) = 0.5 * ( zfp_ui * ptb(ji,jj,jk,jn) + zfm_ui * ptb(ji+1,jj ,jk,jn) ) 132 zwy(ji,jj,jk) = 0.5 * ( zfp_vj * ptb(ji,jj,jk,jn) + zfm_vj * ptb(ji ,jj+1,jk,jn) ) 134 133 END DO 135 134 END DO … … 139 138 DO jj = 1, jpj 140 139 DO ji = 1, jpi 141 zfp_wk = pw (ji,jj,jk) + ABS( pw(ji,jj,jk) )142 zfm_wk = pw (ji,jj,jk) - ABS( pw(ji,jj,jk) )143 zwz(ji,jj,jk) = 0.5 * ( zfp_wk * pt _lev1(ji,jj,jk,jn) + zfm_wk * pt_lev1(ji,jj,jk-1,jn) ) * wmask(ji,jj,jk)140 zfp_wk = pwn(ji,jj,jk) + ABS( pwn(ji,jj,jk) ) 141 zfm_wk = pwn(ji,jj,jk) - ABS( pwn(ji,jj,jk) ) 142 zwz(ji,jj,jk) = 0.5 * ( zfp_wk * ptb(ji,jj,jk,jn) + zfm_wk * ptb(ji,jj,jk-1,jn) ) * wmask(ji,jj,jk) 144 143 END DO 145 144 END DO … … 149 148 DO jj = 1, jpj 150 149 DO ji = 1, jpi 151 zwz(ji,jj, mikt(ji,jj) ) = pw (ji,jj,mikt(ji,jj)) * pt_lev1(ji,jj,mikt(ji,jj),jn) ! linear free surface150 zwz(ji,jj, mikt(ji,jj) ) = pwn(ji,jj,mikt(ji,jj)) * ptb(ji,jj,mikt(ji,jj),jn) ! linear free surface 152 151 END DO 153 152 END DO 154 153 ELSE ! no cavities: only at the ocean surface 155 zwz(:,:,1) = pw (:,:,1) * pt_lev1(:,:,1,jn)154 zwz(:,:,1) = pwn(:,:,1) * ptb(:,:,1,jn) 156 155 ENDIF 157 156 ENDIF … … 165 164 & + zwz(ji,jj,jk) - zwz(ji ,jj ,jk+1) ) * r1_e1e2t(ji,jj) 166 165 ! ! update and guess with monotonic sheme 167 pt _rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn) + ztra / e3t(ji,jj,jk,ktlev2) * tmask(ji,jj,jk)168 zwi(ji,jj,jk) = ( e3t (ji,jj,jk,ktlev1) * pt_lev1(ji,jj,jk,jn) + p2dt * ztra ) / e3t(ji,jj,jk,ktlev3) * tmask(ji,jj,jk)166 pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra / e3t_n(ji,jj,jk) * tmask(ji,jj,jk) 167 zwi(ji,jj,jk) = ( e3t_b(ji,jj,jk) * ptb(ji,jj,jk,jn) + p2dt * ztra ) / e3t_a(ji,jj,jk) * tmask(ji,jj,jk) 169 168 END DO 170 169 END DO … … 185 184 DO jj = 1, jpjm1 186 185 DO ji = 1, fs_jpim1 ! vector opt. 187 zwx(ji,jj,jk) = 0.5_wp * pu (ji,jj,jk) * ( pt_lev2(ji,jj,jk,jn) + pt_lev2(ji+1,jj,jk,jn) ) - zwx(ji,jj,jk)188 zwy(ji,jj,jk) = 0.5_wp * pv (ji,jj,jk) * ( pt_lev2(ji,jj,jk,jn) + pt_lev2(ji,jj+1,jk,jn) ) - zwy(ji,jj,jk)186 zwx(ji,jj,jk) = 0.5_wp * pun(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji+1,jj,jk,jn) ) - zwx(ji,jj,jk) 187 zwy(ji,jj,jk) = 0.5_wp * pvn(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji,jj+1,jk,jn) ) - zwy(ji,jj,jk) 189 188 END DO 190 189 END DO … … 197 196 DO jj = 1, jpjm1 ! 1st derivative (gradient) 198 197 DO ji = 1, fs_jpim1 ! vector opt. 199 ztu(ji,jj,jk) = ( pt _lev2(ji+1,jj ,jk,jn) - pt_lev2(ji,jj,jk,jn) ) * umask(ji,jj,jk)200 ztv(ji,jj,jk) = ( pt _lev2(ji ,jj+1,jk,jn) - pt_lev2(ji,jj,jk,jn) ) * vmask(ji,jj,jk)198 ztu(ji,jj,jk) = ( ptn(ji+1,jj ,jk,jn) - ptn(ji,jj,jk,jn) ) * umask(ji,jj,jk) 199 ztv(ji,jj,jk) = ( ptn(ji ,jj+1,jk,jn) - ptn(ji,jj,jk,jn) ) * vmask(ji,jj,jk) 201 200 END DO 202 201 END DO … … 213 212 DO jj = 1, jpjm1 214 213 DO ji = 1, fs_jpim1 ! vector opt. 215 zC2t_u = pt _lev2(ji,jj,jk,jn) + pt_lev2(ji+1,jj ,jk,jn) ! 2 x C2 interpolation of T at u- & v-points216 zC2t_v = pt _lev2(ji,jj,jk,jn) + pt_lev2(ji ,jj+1,jk,jn)214 zC2t_u = ptn(ji,jj,jk,jn) + ptn(ji+1,jj ,jk,jn) ! 2 x C2 interpolation of T at u- & v-points 215 zC2t_v = ptn(ji,jj,jk,jn) + ptn(ji ,jj+1,jk,jn) 217 216 ! ! C4 minus upstream advective fluxes 218 zwx(ji,jj,jk) = 0.5_wp * pu (ji,jj,jk) * ( zC2t_u + zltu(ji,jj,jk) - zltu(ji+1,jj,jk) ) - zwx(ji,jj,jk)219 zwy(ji,jj,jk) = 0.5_wp * pv (ji,jj,jk) * ( zC2t_v + zltv(ji,jj,jk) - zltv(ji,jj+1,jk) ) - zwy(ji,jj,jk)217 zwx(ji,jj,jk) = 0.5_wp * pun(ji,jj,jk) * ( zC2t_u + zltu(ji,jj,jk) - zltu(ji+1,jj,jk) ) - zwx(ji,jj,jk) 218 zwy(ji,jj,jk) = 0.5_wp * pvn(ji,jj,jk) * ( zC2t_v + zltv(ji,jj,jk) - zltv(ji,jj+1,jk) ) - zwy(ji,jj,jk) 220 219 END DO 221 220 END DO … … 228 227 DO jj = 1, jpjm1 229 228 DO ji = 1, fs_jpim1 ! vector opt. 230 ztu(ji,jj,jk) = ( pt _lev2(ji+1,jj ,jk,jn) - pt_lev2(ji,jj,jk,jn) ) * umask(ji,jj,jk)231 ztv(ji,jj,jk) = ( pt _lev2(ji ,jj+1,jk,jn) - pt_lev2(ji,jj,jk,jn) ) * vmask(ji,jj,jk)229 ztu(ji,jj,jk) = ( ptn(ji+1,jj ,jk,jn) - ptn(ji,jj,jk,jn) ) * umask(ji,jj,jk) 230 ztv(ji,jj,jk) = ( ptn(ji ,jj+1,jk,jn) - ptn(ji,jj,jk,jn) ) * vmask(ji,jj,jk) 232 231 END DO 233 232 END DO … … 238 237 DO jj = 2, jpjm1 239 238 DO ji = 2, fs_jpim1 ! vector opt. 240 zC2t_u = pt _lev2(ji,jj,jk,jn) + pt_lev2(ji+1,jj ,jk,jn) ! 2 x C2 interpolation of T at u- & v-points (x2)241 zC2t_v = pt _lev2(ji,jj,jk,jn) + pt_lev2(ji ,jj+1,jk,jn)239 zC2t_u = ptn(ji,jj,jk,jn) + ptn(ji+1,jj ,jk,jn) ! 2 x C2 interpolation of T at u- & v-points (x2) 240 zC2t_v = ptn(ji,jj,jk,jn) + ptn(ji ,jj+1,jk,jn) 242 241 ! ! C4 interpolation of T at u- & v-points (x2) 243 242 zC4t_u = zC2t_u + r1_6 * ( ztu(ji-1,jj ,jk) - ztu(ji+1,jj ,jk) ) 244 243 zC4t_v = zC2t_v + r1_6 * ( ztv(ji ,jj-1,jk) - ztv(ji ,jj+1,jk) ) 245 244 ! ! C4 minus upstream advective fluxes 246 zwx(ji,jj,jk) = 0.5_wp * pu (ji,jj,jk) * zC4t_u - zwx(ji,jj,jk)247 zwy(ji,jj,jk) = 0.5_wp * pv (ji,jj,jk) * zC4t_v - zwy(ji,jj,jk)245 zwx(ji,jj,jk) = 0.5_wp * pun(ji,jj,jk) * zC4t_u - zwx(ji,jj,jk) 246 zwy(ji,jj,jk) = 0.5_wp * pvn(ji,jj,jk) * zC4t_v - zwy(ji,jj,jk) 248 247 END DO 249 248 END DO … … 258 257 DO jj = 2, jpjm1 259 258 DO ji = fs_2, fs_jpim1 260 zwz(ji,jj,jk) = ( pw (ji,jj,jk) * 0.5_wp * ( pt_lev2(ji,jj,jk,jn) + pt_lev2(ji,jj,jk-1,jn) ) &259 zwz(ji,jj,jk) = ( pwn(ji,jj,jk) * 0.5_wp * ( ptn(ji,jj,jk,jn) + ptn(ji,jj,jk-1,jn) ) & 261 260 & - zwz(ji,jj,jk) ) * wmask(ji,jj,jk) 262 261 END DO … … 265 264 ! 266 265 CASE( 4 ) !- 4th order COMPACT 267 CALL interp_4th_cpt( pt _lev2(:,:,:,jn) , ztw ) ! zwt = COMPACT interpolation of T at w-point266 CALL interp_4th_cpt( ptn(:,:,:,jn) , ztw ) ! zwt = COMPACT interpolation of T at w-point 268 267 DO jk = 2, jpkm1 269 268 DO jj = 2, jpjm1 270 269 DO ji = fs_2, fs_jpim1 271 zwz(ji,jj,jk) = ( pw (ji,jj,jk) * ztw(ji,jj,jk) - zwz(ji,jj,jk) ) * wmask(ji,jj,jk)270 zwz(ji,jj,jk) = ( pwn(ji,jj,jk) * ztw(ji,jj,jk) - zwz(ji,jj,jk) ) * wmask(ji,jj,jk) 272 271 END DO 273 272 END DO … … 283 282 ! !== monotonicity algorithm ==! 284 283 ! 285 CALL nonosc( pt _lev1(:,:,:,jn), zwx, zwy, zwz, zwi, p2dt, e3t(:,:,:,ktlev2))284 CALL nonosc( ptb(:,:,:,jn), zwx, zwy, zwz, zwi, p2dt ) 286 285 ! 287 286 ! !== final trend with corrected fluxes ==! … … 290 289 DO jj = 2, jpjm1 291 290 DO ji = fs_2, fs_jpim1 ! vector opt. 292 pt _rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn) - ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) &291 pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) - ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) & 293 292 & + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) & 294 293 & + zwz(ji,jj,jk) - zwz(ji ,jj ,jk+1) ) & 295 & * r1_e1e2t(ji,jj) / e3t (ji,jj,jk,ktlev2)294 & * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) 296 295 END DO 297 296 END DO … … 304 303 ! 305 304 IF( l_trd ) THEN ! trend diagnostics 306 CALL trd_tra( kt, cdtype, jn, jptra_xad, ztrdx, pu , pt_lev2(:,:,:,jn) )307 CALL trd_tra( kt, cdtype, jn, jptra_yad, ztrdy, pv , pt_lev2(:,:,:,jn) )308 CALL trd_tra( kt, cdtype, jn, jptra_zad, ztrdz, pw , pt_lev2(:,:,:,jn) )305 CALL trd_tra( kt, cdtype, jn, jptra_xad, ztrdx, pun, ptn(:,:,:,jn) ) 306 CALL trd_tra( kt, cdtype, jn, jptra_yad, ztrdy, pvn, ptn(:,:,:,jn) ) 307 CALL trd_tra( kt, cdtype, jn, jptra_zad, ztrdz, pwn, ptn(:,:,:,jn) ) 309 308 ENDIF 310 309 ! ! heat/salt transport … … 329 328 330 329 331 SUBROUTINE nonosc( pbef, paa, pbb, pcc, paft, p2dt , pe3t)330 SUBROUTINE nonosc( pbef, paa, pbb, pcc, paft, p2dt ) 332 331 !!--------------------------------------------------------------------- 333 332 !! *** ROUTINE nonosc *** … … 342 341 !! in-space based differencing for fluid 343 342 !!---------------------------------------------------------------------- 344 REAL(wp) , INTENT(in ) :: p2dt 345 REAL(wp), DIMENSION (jpi,jpj,jpk), INTENT(in ) :: pbef, paft , pe3t ! before & after field, now e3tfield346 REAL(wp), DIMENSION (jpi,jpj,jpk), INTENT(inout) :: paa, pbb, pcc 343 REAL(wp) , INTENT(in ) :: p2dt ! tracer time-step 344 REAL(wp), DIMENSION (jpi,jpj,jpk), INTENT(in ) :: pbef, paft ! before & after field 345 REAL(wp), DIMENSION (jpi,jpj,jpk), INTENT(inout) :: paa, pbb, pcc ! monotonic fluxes in the 3 directions 347 346 ! 348 347 INTEGER :: ji, jj, jk ! dummy loop indices … … 393 392 394 393 ! up & down beta terms 395 zbt = e1e2t(ji,jj) * pe3t(ji,jj,jk) / p2dt394 zbt = e1e2t(ji,jj) * e3t_n(ji,jj,jk) / p2dt 396 395 zbetup(ji,jj,jk) = ( zup - paft(ji,jj,jk) ) / ( zpos + zrtrn ) * zbt 397 396 zbetdo(ji,jj,jk) = ( paft(ji,jj,jk) - zdo ) / ( zneg + zrtrn ) * zbt … … 635 634 !! The tri-diagonals matrix is given as input 3D arrays: pD, pU, pL 636 635 !! (i.e. the Diagonal, the Upper diagonal, and the Lower diagonal). 637 !! The solution is pt _rhs.636 !! The solution is pta. 638 637 !! The 3d array zwt is used as a work space array. 639 638 !!---------------------------------------------------------------------- -
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traadv_mus.F90
r10806 r10874 54 54 CONTAINS 55 55 56 SUBROUTINE tra_adv_mus( kt, kit000, ktlev, kt2lev, cdtype, p2dt, pu, pv, pw, &57 & pt, pt_rhs, kjpt, ld_msc_ups )56 SUBROUTINE tra_adv_mus( kt, kit000, cdtype, p2dt, pun, pvn, pwn, & 57 & ptb, pta, kjpt, ld_msc_ups ) 58 58 !!---------------------------------------------------------------------- 59 59 !! *** ROUTINE tra_adv_mus *** … … 66 66 !! ld_msc_ups=T : 67 67 !! 68 !! ** Action : - update pt _rhswith the now advective tracer trends68 !! ** Action : - update pta with the now advective tracer trends 69 69 !! - send trends to trdtra module for further diagnostcs (l_trdtra=T) 70 70 !! - htr_adv, str_adv : poleward advective heat and salt transport (ln_diaptr=T) … … 75 75 INTEGER , INTENT(in ) :: kt ! ocean time-step index 76 76 INTEGER , INTENT(in ) :: kit000 ! first time step index 77 INTEGER , INTENT(in ) :: ktlev ! time level index for 3-time-level source terms78 INTEGER , INTENT(in ) :: kt2lev ! time level index for 2-time-level source terms79 77 CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator) 80 78 INTEGER , INTENT(in ) :: kjpt ! number of tracers 81 79 LOGICAL , INTENT(in ) :: ld_msc_ups ! use upstream scheme within muscl 82 80 REAL(wp) , INTENT(in ) :: p2dt ! tracer time-step 83 REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pu , pv, pw! 3 ocean velocity components84 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: pt ! before tracer field85 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pt _rhs! tracer trend81 REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pun, pvn, pwn ! 3 ocean velocity components 82 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: ptb ! before tracer field 83 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pta ! tracer trend 86 84 ! 87 85 INTEGER :: ji, jj, jk, jn ! dummy loop indices … … 136 134 DO jj = 1, jpjm1 137 135 DO ji = 1, fs_jpim1 ! vector opt. 138 zwx(ji,jj,jk) = umask(ji,jj,jk) * ( pt (ji+1,jj,jk,jn) - pt(ji,jj,jk,jn) )139 zwy(ji,jj,jk) = vmask(ji,jj,jk) * ( pt (ji,jj+1,jk,jn) - pt(ji,jj,jk,jn) )136 zwx(ji,jj,jk) = umask(ji,jj,jk) * ( ptb(ji+1,jj,jk,jn) - ptb(ji,jj,jk,jn) ) 137 zwy(ji,jj,jk) = vmask(ji,jj,jk) * ( ptb(ji,jj+1,jk,jn) - ptb(ji,jj,jk,jn) ) 140 138 END DO 141 139 END DO … … 174 172 DO ji = fs_2, fs_jpim1 ! vector opt. 175 173 ! MUSCL fluxes 176 z0u = SIGN( 0.5, pu (ji,jj,jk) )174 z0u = SIGN( 0.5, pun(ji,jj,jk) ) 177 175 zalpha = 0.5 - z0u 178 zu = z0u - 0.5 * pu (ji,jj,jk) * p2dt * r1_e1e2u(ji,jj) / e3u(ji,jj,jk,ktlev)179 zzwx = pt (ji+1,jj,jk,jn) + xind(ji,jj,jk) * zu * zslpx(ji+1,jj,jk)180 zzwy = pt (ji ,jj,jk,jn) + xind(ji,jj,jk) * zu * zslpx(ji ,jj,jk)181 zwx(ji,jj,jk) = pu (ji,jj,jk) * ( zalpha * zzwx + (1.-zalpha) * zzwy )176 zu = z0u - 0.5 * pun(ji,jj,jk) * p2dt * r1_e1e2u(ji,jj) / e3u_n(ji,jj,jk) 177 zzwx = ptb(ji+1,jj,jk,jn) + xind(ji,jj,jk) * zu * zslpx(ji+1,jj,jk) 178 zzwy = ptb(ji ,jj,jk,jn) + xind(ji,jj,jk) * zu * zslpx(ji ,jj,jk) 179 zwx(ji,jj,jk) = pun(ji,jj,jk) * ( zalpha * zzwx + (1.-zalpha) * zzwy ) 182 180 ! 183 z0v = SIGN( 0.5, pv (ji,jj,jk) )181 z0v = SIGN( 0.5, pvn(ji,jj,jk) ) 184 182 zalpha = 0.5 - z0v 185 zv = z0v - 0.5 * pv (ji,jj,jk) * p2dt * r1_e1e2v(ji,jj) / e3v(ji,jj,jk,ktlev)186 zzwx = pt (ji,jj+1,jk,jn) + xind(ji,jj,jk) * zv * zslpy(ji,jj+1,jk)187 zzwy = pt (ji,jj ,jk,jn) + xind(ji,jj,jk) * zv * zslpy(ji,jj ,jk)188 zwy(ji,jj,jk) = pv (ji,jj,jk) * ( zalpha * zzwx + (1.-zalpha) * zzwy )183 zv = z0v - 0.5 * pvn(ji,jj,jk) * p2dt * r1_e1e2v(ji,jj) / e3v_n(ji,jj,jk) 184 zzwx = ptb(ji,jj+1,jk,jn) + xind(ji,jj,jk) * zv * zslpy(ji,jj+1,jk) 185 zzwy = ptb(ji,jj ,jk,jn) + xind(ji,jj,jk) * zv * zslpy(ji,jj ,jk) 186 zwy(ji,jj,jk) = pvn(ji,jj,jk) * ( zalpha * zzwx + (1.-zalpha) * zzwy ) 189 187 END DO 190 188 END DO … … 195 193 DO jj = 2, jpjm1 196 194 DO ji = fs_2, fs_jpim1 ! vector opt. 197 pt _rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn) - ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) &195 pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) - ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) & 198 196 & + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) ) & 199 & * r1_e1e2t(ji,jj) / e3t (ji,jj,jk,ktlev)197 & * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) 200 198 END DO 201 199 END DO … … 203 201 ! ! trend diagnostics 204 202 IF( l_trd ) THEN 205 CALL trd_tra( kt, cdtype, jn, jptra_xad, zwx, pu , pt(:,:,:,jn) )206 CALL trd_tra( kt, cdtype, jn, jptra_yad, zwy, pv , pt(:,:,:,jn) )203 CALL trd_tra( kt, cdtype, jn, jptra_xad, zwx, pun, ptb(:,:,:,jn) ) 204 CALL trd_tra( kt, cdtype, jn, jptra_yad, zwy, pvn, ptb(:,:,:,jn) ) 207 205 END IF 208 206 ! ! "Poleward" heat and salt transports … … 217 215 zwx(:,:,jpk) = 0._wp 218 216 DO jk = 2, jpkm1 ! interior values 219 zwx(:,:,jk) = tmask(:,:,jk) * ( pt (:,:,jk-1,jn) - pt(:,:,jk,jn) )217 zwx(:,:,jk) = tmask(:,:,jk) * ( ptb(:,:,jk-1,jn) - ptb(:,:,jk,jn) ) 220 218 END DO 221 219 ! !-- Slopes of tracer … … 241 239 DO jj = 2, jpjm1 242 240 DO ji = fs_2, fs_jpim1 ! vector opt. 243 z0w = SIGN( 0.5, pw (ji,jj,jk+1) )241 z0w = SIGN( 0.5, pwn(ji,jj,jk+1) ) 244 242 zalpha = 0.5 + z0w 245 zw = z0w - 0.5 * pw (ji,jj,jk+1) * p2dt * r1_e1e2t(ji,jj) / e3w(ji,jj,jk+1,kt2lev)246 zzwx = pt (ji,jj,jk+1,jn) + xind(ji,jj,jk) * zw * zslpx(ji,jj,jk+1)247 zzwy = pt (ji,jj,jk ,jn) + xind(ji,jj,jk) * zw * zslpx(ji,jj,jk )248 zwx(ji,jj,jk+1) = pw (ji,jj,jk+1) * ( zalpha * zzwx + (1.-zalpha) * zzwy ) * wmask(ji,jj,jk)243 zw = z0w - 0.5 * pwn(ji,jj,jk+1) * p2dt * r1_e1e2t(ji,jj) / e3w_n(ji,jj,jk+1) 244 zzwx = ptb(ji,jj,jk+1,jn) + xind(ji,jj,jk) * zw * zslpx(ji,jj,jk+1) 245 zzwy = ptb(ji,jj,jk ,jn) + xind(ji,jj,jk) * zw * zslpx(ji,jj,jk ) 246 zwx(ji,jj,jk+1) = pwn(ji,jj,jk+1) * ( zalpha * zzwx + (1.-zalpha) * zzwy ) * wmask(ji,jj,jk) 249 247 END DO 250 248 END DO … … 254 252 DO jj = 1, jpj 255 253 DO ji = 1, jpi 256 zwx(ji,jj, mikt(ji,jj) ) = pw (ji,jj,mikt(ji,jj)) * pt(ji,jj,mikt(ji,jj),jn)254 zwx(ji,jj, mikt(ji,jj) ) = pwn(ji,jj,mikt(ji,jj)) * ptb(ji,jj,mikt(ji,jj),jn) 257 255 END DO 258 256 END DO 259 257 ELSE ! no cavities: only at the ocean surface 260 zwx(:,:,1) = pw (:,:,1) * pt(:,:,1,jn)258 zwx(:,:,1) = pwn(:,:,1) * ptb(:,:,1,jn) 261 259 ENDIF 262 260 ENDIF … … 265 263 DO jj = 2, jpjm1 266 264 DO ji = fs_2, fs_jpim1 ! vector opt. 267 pt _rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn) - ( zwx(ji,jj,jk) - zwx(ji,jj,jk+1) ) * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,ktlev)265 pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) - ( zwx(ji,jj,jk) - zwx(ji,jj,jk+1) ) * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) 268 266 END DO 269 267 END DO 270 268 END DO 271 269 ! ! send trends for diagnostic 272 IF( l_trd ) CALL trd_tra( kt, cdtype, jn, jptra_zad, zwx, pw , pt(:,:,:,jn) )270 IF( l_trd ) CALL trd_tra( kt, cdtype, jn, jptra_zad, zwx, pwn, ptb(:,:,:,jn) ) 273 271 ! 274 272 END DO ! end of tracer loop -
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traadv_qck.F90
r10806 r10874 47 47 CONTAINS 48 48 49 SUBROUTINE tra_adv_qck ( kt, kit000, ktlev, cdtype, p2dt, pu, pv, pw, &50 & pt_lev1, pt_lev2, pt_rhs, kjpt )49 SUBROUTINE tra_adv_qck ( kt, kit000, cdtype, p2dt, pun, pvn, pwn, & 50 & ptb, ptn, pta, kjpt ) 51 51 !!---------------------------------------------------------------------- 52 52 !! *** ROUTINE tra_adv_qck *** … … 72 72 !! dt = 2*rdtra and the scalar values are tb and sb 73 73 !! 74 !! On the vertical, the simple centered scheme used pt _lev274 !! On the vertical, the simple centered scheme used ptn 75 75 !! 76 76 !! The fluxes are bounded by the ULTIMATE limiter to … … 78 78 !! prevent the appearance of spurious numerical oscillations 79 79 !! 80 !! ** Action : - update pt _rhswith the now advective tracer trends80 !! ** Action : - update pta with the now advective tracer trends 81 81 !! - send trends to trdtra module for further diagnostcs (l_trdtra=T) 82 82 !! - htr_adv, str_adv : poleward advective heat and salt transport (ln_diaptr=T) … … 86 86 INTEGER , INTENT(in ) :: kt ! ocean time-step index 87 87 INTEGER , INTENT(in ) :: kit000 ! first time step index 88 INTEGER , INTENT(in ) :: ktlev ! time level index for source terms89 88 CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator) 90 89 INTEGER , INTENT(in ) :: kjpt ! number of tracers 91 90 REAL(wp) , INTENT(in ) :: p2dt ! tracer time-step 92 REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pu , pv, pw! 3 ocean velocity components93 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: pt _lev1, pt_lev2! before and now tracer fields94 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pt _rhs! tracer trend91 REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pun, pvn, pwn ! 3 ocean velocity components 92 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: ptb, ptn ! before and now tracer fields 93 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pta ! tracer trend 95 94 !!---------------------------------------------------------------------- 96 95 ! … … 109 108 ! 110 109 ! ! horizontal fluxes are computed with the QUICKEST + ULTIMATE scheme 111 CALL tra_adv_qck_i( kt, ktlev, cdtype, p2dt, pu, pt_lev1, pt_lev2, pt_rhs, kjpt )112 CALL tra_adv_qck_j( kt, ktlev, cdtype, p2dt, pv, pt_lev1, pt_lev2, pt_rhs, kjpt )110 CALL tra_adv_qck_i( kt, cdtype, p2dt, pun, ptb, ptn, pta, kjpt ) 111 CALL tra_adv_qck_j( kt, cdtype, p2dt, pvn, ptb, ptn, pta, kjpt ) 113 112 114 113 ! ! vertical fluxes are computed with the 2nd order centered scheme 115 CALL tra_adv_cen2_k( kt, ktlev, cdtype, pw, pt_lev2, pt_rhs, kjpt )114 CALL tra_adv_cen2_k( kt, cdtype, pwn, ptn, pta, kjpt ) 116 115 ! 117 116 END SUBROUTINE tra_adv_qck 118 117 119 118 120 SUBROUTINE tra_adv_qck_i( kt, ktlev, cdtype, p2dt, pu, &121 & pt _lev1, pt_lev2, pt_rhs, kjpt )119 SUBROUTINE tra_adv_qck_i( kt, cdtype, p2dt, pun, & 120 & ptb, ptn, pta, kjpt ) 122 121 !!---------------------------------------------------------------------- 123 122 !! 124 123 !!---------------------------------------------------------------------- 125 124 INTEGER , INTENT(in ) :: kt ! ocean time-step index 126 INTEGER , INTENT(in ) :: ktlev ! time level index for source terms127 125 CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator) 128 126 INTEGER , INTENT(in ) :: kjpt ! number of tracers 129 127 REAL(wp) , INTENT(in ) :: p2dt ! tracer time-step 130 REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pu ! i-velocity components131 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: pt _lev1, pt_lev2! before and now tracer fields132 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pt _rhs! tracer trend128 REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pun ! i-velocity components 129 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: ptb, ptn ! before and now tracer fields 130 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pta ! tracer trend 133 131 !! 134 132 INTEGER :: ji, jj, jk, jn ! dummy loop indices … … 147 145 DO jj = 2, jpjm1 148 146 DO ji = fs_2, fs_jpim1 ! vector opt. 149 zfc(ji,jj,jk) = pt _lev1(ji-1,jj,jk,jn) ! Upstream in the x-direction for the tracer150 zfd(ji,jj,jk) = pt _lev1(ji+1,jj,jk,jn) ! Downstream in the x-direction for the tracer147 zfc(ji,jj,jk) = ptb(ji-1,jj,jk,jn) ! Upstream in the x-direction for the tracer 148 zfd(ji,jj,jk) = ptb(ji+1,jj,jk,jn) ! Downstream in the x-direction for the tracer 151 149 END DO 152 150 END DO … … 160 158 DO jj = 2, jpjm1 161 159 DO ji = fs_2, fs_jpim1 ! vector opt. 162 zdir = 0.5 + SIGN( 0.5, pu (ji,jj,jk) ) ! if pu> 0 : zdir = 1 otherwise zdir = 0160 zdir = 0.5 + SIGN( 0.5, pun(ji,jj,jk) ) ! if pun > 0 : zdir = 1 otherwise zdir = 0 163 161 zfu(ji,jj,jk) = zdir * zfc(ji,jj,jk ) + ( 1. - zdir ) * zfd(ji+1,jj,jk) ! FU in the x-direction for T 164 162 END DO … … 169 167 DO jj = 2, jpjm1 170 168 DO ji = fs_2, fs_jpim1 ! vector opt. 171 zdir = 0.5 + SIGN( 0.5, pu (ji,jj,jk) ) ! if pu> 0 : zdir = 1 otherwise zdir = 0172 zdx = ( zdir * e1t(ji,jj) + ( 1. - zdir ) * e1t(ji+1,jj) ) * e2u(ji,jj) * e3u (ji,jj,jk,ktlev)173 zwx(ji,jj,jk) = ABS( pu (ji,jj,jk) ) * p2dt / zdx ! (0<zc_cfl<1 : Courant number on x-direction)174 zfc(ji,jj,jk) = zdir * pt _lev1(ji ,jj,jk,jn) + ( 1. - zdir ) * pt_lev1(ji+1,jj,jk,jn) ! FC in the x-direction for T175 zfd(ji,jj,jk) = zdir * pt _lev1(ji+1,jj,jk,jn) + ( 1. - zdir ) * pt_lev1(ji ,jj,jk,jn) ! FD in the x-direction for T169 zdir = 0.5 + SIGN( 0.5, pun(ji,jj,jk) ) ! if pun > 0 : zdir = 1 otherwise zdir = 0 170 zdx = ( zdir * e1t(ji,jj) + ( 1. - zdir ) * e1t(ji+1,jj) ) * e2u(ji,jj) * e3u_n(ji,jj,jk) 171 zwx(ji,jj,jk) = ABS( pun(ji,jj,jk) ) * p2dt / zdx ! (0<zc_cfl<1 : Courant number on x-direction) 172 zfc(ji,jj,jk) = zdir * ptb(ji ,jj,jk,jn) + ( 1. - zdir ) * ptb(ji+1,jj,jk,jn) ! FC in the x-direction for T 173 zfd(ji,jj,jk) = zdir * ptb(ji+1,jj,jk,jn) + ( 1. - zdir ) * ptb(ji ,jj,jk,jn) ! FD in the x-direction for T 176 174 END DO 177 175 END DO … … 199 197 DO jj = 2, jpjm1 200 198 DO ji = fs_2, fs_jpim1 ! vector opt. 201 zdir = 0.5 + SIGN( 0.5, pu (ji,jj,jk) ) ! if pu> 0 : zdir = 1 otherwise zdir = 0199 zdir = 0.5 + SIGN( 0.5, pun(ji,jj,jk) ) ! if pun > 0 : zdir = 1 otherwise zdir = 0 202 200 !--- If the second ustream point is a land point 203 201 !--- the flux is computed by the 1st order UPWIND scheme 204 202 zmsk = zdir * zfu(ji,jj,jk) + ( 1. - zdir ) * zfu(ji+1,jj,jk) 205 203 zwx(ji,jj,jk) = zmsk * zwx(ji,jj,jk) + ( 1. - zmsk ) * zfc(ji,jj,jk) 206 zwx(ji,jj,jk) = zwx(ji,jj,jk) * pu (ji,jj,jk)204 zwx(ji,jj,jk) = zwx(ji,jj,jk) * pun(ji,jj,jk) 207 205 END DO 208 206 END DO … … 215 213 DO jj = 2, jpjm1 216 214 DO ji = fs_2, fs_jpim1 ! vector opt. 217 zbtr = r1_e1e2t(ji,jj) / e3t (ji,jj,jk,ktlev)215 zbtr = r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) 218 216 ! horizontal advective trends 219 217 ztra = - zbtr * ( zwx(ji,jj,jk) - zwx(ji-1,jj,jk) ) 220 218 !--- add it to the general tracer trends 221 pt _rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn) + ztra219 pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra 222 220 END DO 223 221 END DO 224 222 END DO 225 223 ! ! trend diagnostics 226 IF( l_trd ) CALL trd_tra( kt, cdtype, jn, jptra_xad, zwx, pu , pt_lev2(:,:,:,jn) )224 IF( l_trd ) CALL trd_tra( kt, cdtype, jn, jptra_xad, zwx, pun, ptn(:,:,:,jn) ) 227 225 ! 228 226 END DO … … 231 229 232 230 233 SUBROUTINE tra_adv_qck_j( kt, ktlev, cdtype, p2dt, pv, &234 & pt _lev1, pt_lev2, pt_rhs, kjpt )231 SUBROUTINE tra_adv_qck_j( kt, cdtype, p2dt, pvn, & 232 & ptb, ptn, pta, kjpt ) 235 233 !!---------------------------------------------------------------------- 236 234 !! 237 235 !!---------------------------------------------------------------------- 238 236 INTEGER , INTENT(in ) :: kt ! ocean time-step index 239 INTEGER , INTENT(in ) :: ktlev ! time level index for source terms240 237 CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator) 241 238 INTEGER , INTENT(in ) :: kjpt ! number of tracers 242 239 REAL(wp) , INTENT(in ) :: p2dt ! tracer time-step 243 REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pv ! j-velocity components244 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: pt _lev1, pt_lev2! before and now tracer fields245 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pt _rhs! tracer trend240 REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pvn ! j-velocity components 241 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: ptb, ptn ! before and now tracer fields 242 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pta ! tracer trend 246 243 !! 247 244 INTEGER :: ji, jj, jk, jn ! dummy loop indices … … 262 259 DO ji = fs_2, fs_jpim1 ! vector opt. 263 260 ! Upstream in the x-direction for the tracer 264 zfc(ji,jj,jk) = pt _lev1(ji,jj-1,jk,jn)261 zfc(ji,jj,jk) = ptb(ji,jj-1,jk,jn) 265 262 ! Downstream in the x-direction for the tracer 266 zfd(ji,jj,jk) = pt _lev1(ji,jj+1,jk,jn)263 zfd(ji,jj,jk) = ptb(ji,jj+1,jk,jn) 267 264 END DO 268 265 END DO … … 278 275 DO jj = 2, jpjm1 279 276 DO ji = fs_2, fs_jpim1 ! vector opt. 280 zdir = 0.5 + SIGN( 0.5, pv (ji,jj,jk) ) ! if pu> 0 : zdir = 1 otherwise zdir = 0277 zdir = 0.5 + SIGN( 0.5, pvn(ji,jj,jk) ) ! if pun > 0 : zdir = 1 otherwise zdir = 0 281 278 zfu(ji,jj,jk) = zdir * zfc(ji,jj,jk ) + ( 1. - zdir ) * zfd(ji,jj+1,jk) ! FU in the x-direction for T 282 279 END DO … … 287 284 DO jj = 2, jpjm1 288 285 DO ji = fs_2, fs_jpim1 ! vector opt. 289 zdir = 0.5 + SIGN( 0.5, pv (ji,jj,jk) ) ! if pu> 0 : zdir = 1 otherwise zdir = 0290 zdx = ( zdir * e2t(ji,jj) + ( 1. - zdir ) * e2t(ji,jj+1) ) * e1v(ji,jj) * e3v (ji,jj,jk,ktlev)291 zwy(ji,jj,jk) = ABS( pv (ji,jj,jk) ) * p2dt / zdx ! (0<zc_cfl<1 : Courant number on x-direction)292 zfc(ji,jj,jk) = zdir * pt _lev1(ji,jj ,jk,jn) + ( 1. - zdir ) * pt_lev1(ji,jj+1,jk,jn) ! FC in the x-direction for T293 zfd(ji,jj,jk) = zdir * pt _lev1(ji,jj+1,jk,jn) + ( 1. - zdir ) * pt_lev1(ji,jj ,jk,jn) ! FD in the x-direction for T286 zdir = 0.5 + SIGN( 0.5, pvn(ji,jj,jk) ) ! if pun > 0 : zdir = 1 otherwise zdir = 0 287 zdx = ( zdir * e2t(ji,jj) + ( 1. - zdir ) * e2t(ji,jj+1) ) * e1v(ji,jj) * e3v_n(ji,jj,jk) 288 zwy(ji,jj,jk) = ABS( pvn(ji,jj,jk) ) * p2dt / zdx ! (0<zc_cfl<1 : Courant number on x-direction) 289 zfc(ji,jj,jk) = zdir * ptb(ji,jj ,jk,jn) + ( 1. - zdir ) * ptb(ji,jj+1,jk,jn) ! FC in the x-direction for T 290 zfd(ji,jj,jk) = zdir * ptb(ji,jj+1,jk,jn) + ( 1. - zdir ) * ptb(ji,jj ,jk,jn) ! FD in the x-direction for T 294 291 END DO 295 292 END DO … … 317 314 DO jj = 2, jpjm1 318 315 DO ji = fs_2, fs_jpim1 ! vector opt. 319 zdir = 0.5 + SIGN( 0.5, pv (ji,jj,jk) ) ! if pu> 0 : zdir = 1 otherwise zdir = 0316 zdir = 0.5 + SIGN( 0.5, pvn(ji,jj,jk) ) ! if pun > 0 : zdir = 1 otherwise zdir = 0 320 317 !--- If the second ustream point is a land point 321 318 !--- the flux is computed by the 1st order UPWIND scheme 322 319 zmsk = zdir * zfu(ji,jj,jk) + ( 1. - zdir ) * zfu(ji,jj+1,jk) 323 320 zwy(ji,jj,jk) = zmsk * zwy(ji,jj,jk) + ( 1. - zmsk ) * zfc(ji,jj,jk) 324 zwy(ji,jj,jk) = zwy(ji,jj,jk) * pv (ji,jj,jk)321 zwy(ji,jj,jk) = zwy(ji,jj,jk) * pvn(ji,jj,jk) 325 322 END DO 326 323 END DO … … 333 330 DO jj = 2, jpjm1 334 331 DO ji = fs_2, fs_jpim1 ! vector opt. 335 zbtr = r1_e1e2t(ji,jj) / e3t (ji,jj,jk,ktlev)332 zbtr = r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) 336 333 ! horizontal advective trends 337 334 ztra = - zbtr * ( zwy(ji,jj,jk) - zwy(ji,jj-1,jk) ) 338 335 !--- add it to the general tracer trends 339 pt _rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn) + ztra336 pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra 340 337 END DO 341 338 END DO 342 339 END DO 343 340 ! ! trend diagnostics 344 IF( l_trd ) CALL trd_tra( kt, cdtype, jn, jptra_yad, zwy, pv , pt_lev2(:,:,:,jn) )341 IF( l_trd ) CALL trd_tra( kt, cdtype, jn, jptra_yad, zwy, pvn, ptn(:,:,:,jn) ) 345 342 ! ! "Poleward" heat and salt transports (contribution of upstream fluxes) 346 343 IF( l_ptr ) CALL dia_ptr_hst( jn, 'adv', zwy(:,:,:) ) … … 351 348 352 349 353 SUBROUTINE tra_adv_cen2_k( kt, ktlev, cdtype, pw, &354 & pt _lev2, pt_rhs, kjpt )350 SUBROUTINE tra_adv_cen2_k( kt, cdtype, pwn, & 351 & ptn, pta, kjpt ) 355 352 !!---------------------------------------------------------------------- 356 353 !! 357 354 !!---------------------------------------------------------------------- 358 355 INTEGER , INTENT(in ) :: kt ! ocean time-step index 359 INTEGER , INTENT(in ) :: ktlev ! time level index for source terms360 356 CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator) 361 357 INTEGER , INTENT(in ) :: kjpt ! number of tracers 362 REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pw ! vertical velocity363 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: pt _lev2! before and now tracer fields364 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pt _rhs! tracer trend358 REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pwn ! vertical velocity 359 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: ptn ! before and now tracer fields 360 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pta ! tracer trend 365 361 ! 366 362 INTEGER :: ji, jj, jk, jn ! dummy loop indices … … 378 374 DO jj = 2, jpjm1 379 375 DO ji = fs_2, fs_jpim1 ! vector opt. 380 zwz(ji,jj,jk) = 0.5 * pw (ji,jj,jk) * ( pt_lev2(ji,jj,jk-1,jn) + pt_lev2(ji,jj,jk,jn) ) * wmask(ji,jj,jk)376 zwz(ji,jj,jk) = 0.5 * pwn(ji,jj,jk) * ( ptn(ji,jj,jk-1,jn) + ptn(ji,jj,jk,jn) ) * wmask(ji,jj,jk) 381 377 END DO 382 378 END DO … … 386 382 DO jj = 1, jpj 387 383 DO ji = 1, jpi 388 zwz(ji,jj, mikt(ji,jj) ) = pw (ji,jj,mikt(ji,jj)) * pt_lev2(ji,jj,mikt(ji,jj),jn) ! linear free surface384 zwz(ji,jj, mikt(ji,jj) ) = pwn(ji,jj,mikt(ji,jj)) * ptn(ji,jj,mikt(ji,jj),jn) ! linear free surface 389 385 END DO 390 386 END DO 391 387 ELSE ! no ocean cavities (only ocean surface) 392 zwz(:,:,1) = pw (:,:,1) * pt_lev2(:,:,1,jn)388 zwz(:,:,1) = pwn(:,:,1) * ptn(:,:,1,jn) 393 389 ENDIF 394 390 ENDIF … … 397 393 DO jj = 2, jpjm1 398 394 DO ji = fs_2, fs_jpim1 ! vector opt. 399 pt _rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn) - ( zwz(ji,jj,jk) - zwz(ji,jj,jk+1) ) &400 & * r1_e1e2t(ji,jj) / e3t (ji,jj,jk,ktlev)395 pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) - ( zwz(ji,jj,jk) - zwz(ji,jj,jk+1) ) & 396 & * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) 401 397 END DO 402 398 END DO 403 399 END DO 404 400 ! ! Send trends for diagnostic 405 IF( l_trd ) CALL trd_tra( kt, cdtype, jn, jptra_zad, zwz, pw , pt_lev2(:,:,:,jn) )401 IF( l_trd ) CALL trd_tra( kt, cdtype, jn, jptra_zad, zwz, pwn, ptn(:,:,:,jn) ) 406 402 ! 407 403 END DO -
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traadv_ubs.F90
r10806 r10874 46 46 CONTAINS 47 47 48 SUBROUTINE tra_adv_ubs( kt, kit000, ktlev, cdtype, p2dt, pu, pv, pw, &49 & pt _lev1, pt_lev2, pt_rhs, kjpt, kn_ubs_v )48 SUBROUTINE tra_adv_ubs( kt, kit000, cdtype, p2dt, pun, pvn, pwn, & 49 & ptb, ptn, pta, kjpt, kn_ubs_v ) 50 50 !!---------------------------------------------------------------------- 51 51 !! *** ROUTINE tra_adv_ubs *** … … 58 58 !! It is only used in the horizontal direction. 59 59 !! For example the i-component of the advective fluxes are given by : 60 !! ! e2u e3u u u ( mi(Tn) - zltu(i ) ,ktlev) if uu(i,ktlev) >= 060 !! ! e2u e3u un ( mi(Tn) - zltu(i ) ) if un(i) >= 0 61 61 !! ztu = ! or 62 !! ! e2u e3u u u ( mi(Tn) - zltu(i+1) ,ktlev) if uu(i,ktlev) < 062 !! ! e2u e3u un ( mi(Tn) - zltu(i+1) ) if un(i) < 0 63 63 !! where zltu is the second derivative of the before temperature field: 64 64 !! zltu = 1/e3t di[ e2u e3u / e1u di[Tb] ] … … 77 77 !! scheme (kn_ubs_v=4). 78 78 !! 79 !! ** Action : - update pt _rhswith the now advective tracer trends79 !! ** Action : - update pta with the now advective tracer trends 80 80 !! - send trends to trdtra module for further diagnostcs (l_trdtra=T) 81 81 !! - htr_adv, str_adv : poleward advective heat and salt transport (ln_diaptr=T) … … 86 86 INTEGER , INTENT(in ) :: kt ! ocean time-step index 87 87 INTEGER , INTENT(in ) :: kit000 ! first time step index 88 INTEGER , INTENT(in ) :: ktlev ! time level index for source terms89 88 CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator) 90 89 INTEGER , INTENT(in ) :: kjpt ! number of tracers 91 90 INTEGER , INTENT(in ) :: kn_ubs_v ! number of tracers 92 91 REAL(wp) , INTENT(in ) :: p2dt ! tracer time-step 93 REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pu , pv, pw! 3 ocean transport components94 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: pt _lev1, pt_lev2! before and now tracer fields95 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pt _rhs! tracer trend92 REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pun, pvn, pwn ! 3 ocean transport components 93 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: ptb, ptn ! before and now tracer fields 94 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pta ! tracer trend 96 95 ! 97 96 INTEGER :: ji, jj, jk, jn ! dummy loop indices … … 127 126 DO jj = 1, jpjm1 ! First derivative (masked gradient) 128 127 DO ji = 1, fs_jpim1 ! vector opt. 129 zeeu = e2_e1u(ji,jj) * e3u (ji,jj,jk,ktlev) * umask(ji,jj,jk)130 zeev = e1_e2v(ji,jj) * e3v (ji,jj,jk,ktlev) * vmask(ji,jj,jk)131 ztu(ji,jj,jk) = zeeu * ( pt _lev1(ji+1,jj ,jk,jn) - pt_lev1(ji,jj,jk,jn) )132 ztv(ji,jj,jk) = zeev * ( pt _lev1(ji ,jj+1,jk,jn) - pt_lev1(ji,jj,jk,jn) )128 zeeu = e2_e1u(ji,jj) * e3u_n(ji,jj,jk) * umask(ji,jj,jk) 129 zeev = e1_e2v(ji,jj) * e3v_n(ji,jj,jk) * vmask(ji,jj,jk) 130 ztu(ji,jj,jk) = zeeu * ( ptb(ji+1,jj ,jk,jn) - ptb(ji,jj,jk,jn) ) 131 ztv(ji,jj,jk) = zeev * ( ptb(ji ,jj+1,jk,jn) - ptb(ji,jj,jk,jn) ) 133 132 END DO 134 133 END DO 135 134 DO jj = 2, jpjm1 ! Second derivative (divergence) 136 135 DO ji = fs_2, fs_jpim1 ! vector opt. 137 zcoef = 1._wp / ( 6._wp * e3t (ji,jj,jk,ktlev) )136 zcoef = 1._wp / ( 6._wp * e3t_n(ji,jj,jk) ) 138 137 zltu(ji,jj,jk) = ( ztu(ji,jj,jk) - ztu(ji-1,jj,jk) ) * zcoef 139 138 zltv(ji,jj,jk) = ( ztv(ji,jj,jk) - ztv(ji,jj-1,jk) ) * zcoef … … 147 146 DO jj = 1, jpjm1 148 147 DO ji = 1, fs_jpim1 ! vector opt. 149 zfp_ui = pu (ji,jj,jk) + ABS( pu(ji,jj,jk) ) ! upstream transport (x2)150 zfm_ui = pu (ji,jj,jk) - ABS( pu(ji,jj,jk) )151 zfp_vj = pv (ji,jj,jk) + ABS( pv(ji,jj,jk) )152 zfm_vj = pv (ji,jj,jk) - ABS( pv(ji,jj,jk) )148 zfp_ui = pun(ji,jj,jk) + ABS( pun(ji,jj,jk) ) ! upstream transport (x2) 149 zfm_ui = pun(ji,jj,jk) - ABS( pun(ji,jj,jk) ) 150 zfp_vj = pvn(ji,jj,jk) + ABS( pvn(ji,jj,jk) ) 151 zfm_vj = pvn(ji,jj,jk) - ABS( pvn(ji,jj,jk) ) 153 152 ! ! 2nd order centered advective fluxes (x2) 154 zcenut = pu (ji,jj,jk) * ( pt_lev2(ji,jj,jk,jn) + pt_lev2(ji+1,jj ,jk,jn) )155 zcenvt = pv (ji,jj,jk) * ( pt_lev2(ji,jj,jk,jn) + pt_lev2(ji ,jj+1,jk,jn) )153 zcenut = pun(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji+1,jj ,jk,jn) ) 154 zcenvt = pvn(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji ,jj+1,jk,jn) ) 156 155 ! ! UBS advective fluxes 157 156 ztu(ji,jj,jk) = 0.5 * ( zcenut - zfp_ui * zltu(ji,jj,jk) - zfm_ui * zltu(ji+1,jj,jk) ) … … 161 160 END DO 162 161 ! 163 zltu(:,:,:) = pt _rhs(:,:,:,jn) ! store the initial trends before its update162 zltu(:,:,:) = pta(:,:,:,jn) ! store the initial trends before its update 164 163 ! 165 164 DO jk = 1, jpkm1 !== add the horizontal advective trend ==! 166 165 DO jj = 2, jpjm1 167 166 DO ji = fs_2, fs_jpim1 ! vector opt. 168 pt _rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn) &167 pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) & 169 168 & - ( ztu(ji,jj,jk) - ztu(ji-1,jj ,jk) & 170 & + ztv(ji,jj,jk) - ztv(ji ,jj-1,jk) ) * r1_e1e2t(ji,jj) / e3t (ji,jj,jk,ktlev)169 & + ztv(ji,jj,jk) - ztv(ji ,jj-1,jk) ) * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) 171 170 END DO 172 171 END DO … … 174 173 END DO 175 174 ! 176 zltu(:,:,:) = pt _rhs(:,:,:,jn) - zltu(:,:,:) ! Horizontal advective trend used in vertical 2nd order FCT case175 zltu(:,:,:) = pta(:,:,:,jn) - zltu(:,:,:) ! Horizontal advective trend used in vertical 2nd order FCT case 177 176 ! ! and/or in trend diagnostic (l_trd=T) 178 177 ! 179 178 IF( l_trd ) THEN ! trend diagnostics 180 CALL trd_tra( kt, cdtype, jn, jptra_xad, ztu, pu , pt_lev2(:,:,:,jn) )181 CALL trd_tra( kt, cdtype, jn, jptra_yad, ztv, pv , pt_lev2(:,:,:,jn) )179 CALL trd_tra( kt, cdtype, jn, jptra_xad, ztu, pun, ptn(:,:,:,jn) ) 180 CALL trd_tra( kt, cdtype, jn, jptra_yad, ztv, pvn, ptn(:,:,:,jn) ) 182 181 END IF 183 182 ! … … 194 193 CASE( 2 ) ! 2nd order FCT 195 194 ! 196 IF( l_trd ) zltv(:,:,:) = pt _rhs(:,:,:,jn) ! store pt_rhsif trend diag.195 IF( l_trd ) zltv(:,:,:) = pta(:,:,:,jn) ! store pta if trend diag. 197 196 ! 198 197 ! !* upstream advection with initial mass fluxes & intermediate update ==! … … 200 199 DO jj = 1, jpj 201 200 DO ji = 1, jpi 202 zfp_wk = pw (ji,jj,jk) + ABS( pw(ji,jj,jk) )203 zfm_wk = pw (ji,jj,jk) - ABS( pw(ji,jj,jk) )204 ztw(ji,jj,jk) = 0.5_wp * ( zfp_wk * pt _lev1(ji,jj,jk,jn) + zfm_wk * pt_lev1(ji,jj,jk-1,jn) ) * wmask(ji,jj,jk)201 zfp_wk = pwn(ji,jj,jk) + ABS( pwn(ji,jj,jk) ) 202 zfm_wk = pwn(ji,jj,jk) - ABS( pwn(ji,jj,jk) ) 203 ztw(ji,jj,jk) = 0.5_wp * ( zfp_wk * ptb(ji,jj,jk,jn) + zfm_wk * ptb(ji,jj,jk-1,jn) ) * wmask(ji,jj,jk) 205 204 END DO 206 205 END DO … … 210 209 DO jj = 1, jpj 211 210 DO ji = 1, jpi 212 ztw(ji,jj, mikt(ji,jj) ) = pw (ji,jj,mikt(ji,jj)) * pt_lev1(ji,jj,mikt(ji,jj),jn) ! linear free surface211 ztw(ji,jj, mikt(ji,jj) ) = pwn(ji,jj,mikt(ji,jj)) * ptb(ji,jj,mikt(ji,jj),jn) ! linear free surface 213 212 END DO 214 213 END DO 215 214 ELSE ! no cavities: only at the ocean surface 216 ztw(:,:,1) = pw (:,:,1) * pt_lev1(:,:,1,jn)215 ztw(:,:,1) = pwn(:,:,1) * ptb(:,:,1,jn) 217 216 ENDIF 218 217 ENDIF … … 221 220 DO jj = 2, jpjm1 222 221 DO ji = fs_2, fs_jpim1 ! vector opt. 223 ztak = - ( ztw(ji,jj,jk) - ztw(ji,jj,jk+1) ) * r1_e1e2t(ji,jj) / e3t (ji,jj,jk,ktlev)224 pt _rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn) + ztak225 zti(ji,jj,jk) = ( pt _lev1(ji,jj,jk,jn) + p2dt * ( ztak + zltu(ji,jj,jk) ) ) * tmask(ji,jj,jk)222 ztak = - ( ztw(ji,jj,jk) - ztw(ji,jj,jk+1) ) * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) 223 pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztak 224 zti(ji,jj,jk) = ( ptb(ji,jj,jk,jn) + p2dt * ( ztak + zltu(ji,jj,jk) ) ) * tmask(ji,jj,jk) 226 225 END DO 227 226 END DO … … 233 232 DO jj = 1, jpj 234 233 DO ji = 1, jpi 235 ztw(ji,jj,jk) = ( 0.5_wp * pw (ji,jj,jk) * ( pt_lev2(ji,jj,jk,jn) + pt_lev2(ji,jj,jk-1,jn) ) &234 ztw(ji,jj,jk) = ( 0.5_wp * pwn(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji,jj,jk-1,jn) ) & 236 235 & - ztw(ji,jj,jk) ) * wmask(ji,jj,jk) 237 236 END DO … … 241 240 IF( ln_linssh ) ztw(:,:, 1 ) = 0._wp ! only ocean surface as interior zwz values have been w-masked 242 241 ! 243 CALL nonosc_z( pt _lev1(:,:,:,jn), ztw, zti, p2dt, e3t(:,:,:,ktlev)) ! monotonicity algorithm242 CALL nonosc_z( ptb(:,:,:,jn), ztw, zti, p2dt ) ! monotonicity algorithm 244 243 ! 245 244 CASE( 4 ) ! 4th order COMPACT 246 CALL interp_4th_cpt( pt _lev2(:,:,:,jn) , ztw ) ! 4th order compact interpolation of T at w-point245 CALL interp_4th_cpt( ptn(:,:,:,jn) , ztw ) ! 4th order compact interpolation of T at w-point 247 246 DO jk = 2, jpkm1 248 247 DO jj = 2, jpjm1 249 248 DO ji = fs_2, fs_jpim1 250 ztw(ji,jj,jk) = pw (ji,jj,jk) * ztw(ji,jj,jk) * wmask(ji,jj,jk)251 END DO 252 END DO 253 END DO 254 IF( ln_linssh ) ztw(:,:, 1 ) = pw (:,:,1) * pt_lev2(:,:,1,jn) !!gm ISF & 4th COMPACT doesn't work249 ztw(ji,jj,jk) = pwn(ji,jj,jk) * ztw(ji,jj,jk) * wmask(ji,jj,jk) 250 END DO 251 END DO 252 END DO 253 IF( ln_linssh ) ztw(:,:, 1 ) = pwn(:,:,1) * ptn(:,:,1,jn) !!gm ISF & 4th COMPACT doesn't work 255 254 ! 256 255 END SELECT … … 259 258 DO jj = 2, jpjm1 260 259 DO ji = fs_2, fs_jpim1 ! vector opt. 261 pt _rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn) - ( ztw(ji,jj,jk) - ztw(ji,jj,jk+1) ) * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,ktlev)260 pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) - ( ztw(ji,jj,jk) - ztw(ji,jj,jk+1) ) * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) 262 261 END DO 263 262 END DO … … 265 264 ! 266 265 IF( l_trd ) THEN ! vertical advective trend diagnostics 267 DO jk = 1, jpkm1 ! (compute -w.dk[ptn]= -dk[w.ptn] + pt _lev2.dk[w])266 DO jk = 1, jpkm1 ! (compute -w.dk[ptn]= -dk[w.ptn] + ptn.dk[w]) 268 267 DO jj = 2, jpjm1 269 268 DO ji = fs_2, fs_jpim1 ! vector opt. 270 zltv(ji,jj,jk) = pt _rhs(ji,jj,jk,jn) - zltv(ji,jj,jk) &271 & + pt _lev2(ji,jj,jk,jn) * ( pw(ji,jj,jk) - pw(ji,jj,jk+1) ) &272 & * r1_e1e2t(ji,jj) / e3t (ji,jj,jk,ktlev)269 zltv(ji,jj,jk) = pta(ji,jj,jk,jn) - zltv(ji,jj,jk) & 270 & + ptn(ji,jj,jk,jn) * ( pwn(ji,jj,jk) - pwn(ji,jj,jk+1) ) & 271 & * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) 273 272 END DO 274 273 END DO … … 282 281 283 282 284 SUBROUTINE nonosc_z( pbef, pcc, paft, p2dt , pe3t)283 SUBROUTINE nonosc_z( pbef, pcc, paft, p2dt ) 285 284 !!--------------------------------------------------------------------- 286 285 !! *** ROUTINE nonosc_z *** … … 297 296 REAL(wp), INTENT(in ) :: p2dt ! tracer time-step 298 297 REAL(wp), DIMENSION (jpi,jpj,jpk) :: pbef ! before field 299 REAL(wp), INTENT(in ), DIMENSION (jpi,jpj,jpk) :: pe3t ! now cell thickness field300 298 REAL(wp), INTENT(inout), DIMENSION (jpi,jpj,jpk) :: paft ! after field 301 299 REAL(wp), INTENT(inout), DIMENSION (jpi,jpj,jpk) :: pcc ! monotonic flux in the k direction … … 354 352 zneg = MAX( 0., pcc(ji ,jj ,jk ) ) - MIN( 0., pcc(ji ,jj ,jk+1) ) 355 353 ! up & down beta terms 356 zbt = e1e2t(ji,jj) * pe3t(ji,jj,jk) / p2dt354 zbt = e1e2t(ji,jj) * e3t_n(ji,jj,jk) / p2dt 357 355 zbetup(ji,jj,jk) = ( zbetup(ji,jj,jk) - paft(ji,jj,jk) ) / (zpos+zrtrn) * zbt 358 356 zbetdo(ji,jj,jk) = ( paft(ji,jj,jk) - zbetdo(ji,jj,jk) ) / (zneg+zrtrn) * zbt -
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/trabbc.F90
r10806 r10874 51 51 CONTAINS 52 52 53 SUBROUTINE tra_bbc( kt , ktlev, pts_rhs)53 SUBROUTINE tra_bbc( kt ) 54 54 !!---------------------------------------------------------------------- 55 55 !! *** ROUTINE tra_bbc *** … … 74 74 !!---------------------------------------------------------------------- 75 75 INTEGER, INTENT(in) :: kt ! ocean time-step index 76 INTEGER, INTENT(in) :: ktlev ! time level index for source terms77 REAL(wp), INTENT( inout), DIMENSION(jpi,jpj,jpk,jpts) :: pts_rhs ! temperature and salinity trends78 76 ! 79 77 INTEGER :: ji, jj ! dummy loop indices … … 85 83 IF( l_trdtra ) THEN ! Save the input temperature trend 86 84 ALLOCATE( ztrdt(jpi,jpj,jpk) ) 87 ztrdt(:,:,:) = pts_rhs(:,:,:,jp_tem)85 ztrdt(:,:,:) = tsa(:,:,:,jp_tem) 88 86 ENDIF 89 87 ! ! Add the geothermal trend on temperature 90 88 DO jj = 2, jpjm1 91 89 DO ji = 2, jpim1 92 pts_rhs(ji,jj,mbkt(ji,jj),jp_tem) = pts_rhs(ji,jj,mbkt(ji,jj),jp_tem) + qgh_trd0(ji,jj) / e3t(ji,jj,mbkt(ji,jj),ktlev)90 tsa(ji,jj,mbkt(ji,jj),jp_tem) = tsa(ji,jj,mbkt(ji,jj),jp_tem) + qgh_trd0(ji,jj) / e3t_n(ji,jj,mbkt(ji,jj)) 93 91 END DO 94 92 END DO 95 93 ! 96 CALL lbc_lnk( 'trabbc', pts_rhs(:,:,:,jp_tem) , 'T', 1. )94 CALL lbc_lnk( 'trabbc', tsa(:,:,:,jp_tem) , 'T', 1. ) 97 95 ! 98 96 IF( l_trdtra ) THEN ! Send the trend for diagnostics 99 ztrdt(:,:,:) = pts_rhs(:,:,:,jp_tem) - ztrdt(:,:,:)97 ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:) 100 98 CALL trd_tra( kt, 'TRA', jp_tem, jptra_bbc, ztrdt ) 101 99 DEALLOCATE( ztrdt ) -
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/trabbl.F90
r10806 r10874 89 89 90 90 91 SUBROUTINE tra_bbl( kt , ktlev1, ktlev2, kt2lev, pts_rhs)91 SUBROUTINE tra_bbl( kt ) 92 92 !!---------------------------------------------------------------------- 93 93 !! *** ROUTINE bbl *** … … 101 101 !! is added to the general tracer trend 102 102 !!---------------------------------------------------------------------- 103 INTEGER, INTENT( in ) :: kt ! ocean time-step 104 INTEGER, INTENT( in ) :: ktlev1, ktlev2 ! time level indices for 3-time-level source terms 105 INTEGER, INTENT( in ) :: kt2lev ! time level index for 2-time-level source terms 106 REAL(wp), INTENT( inout), DIMENSION(jpi,jpj,jpk,jpts) :: pts_rhs ! temperature and salinity trends 103 INTEGER, INTENT( in ) :: kt ! ocean time-step 107 104 ! 108 105 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: ztrdt, ztrds … … 113 110 IF( l_trdtra ) THEN !* Save the T-S input trends 114 111 ALLOCATE( ztrdt(jpi,jpj,jpk) , ztrds(jpi,jpj,jpk) ) 115 ztrdt(:,:,:) = pts_rhs(:,:,:,jp_tem)116 ztrds(:,:,:) = pts_rhs(:,:,:,jp_sal)117 ENDIF 118 119 IF( l_bbl ) CALL bbl( kt, nit000, ktlev1, ktlev2, kt2lev,'TRA' ) !* bbl coef. and transport (only if not already done in trcbbl)112 ztrdt(:,:,:) = tsa(:,:,:,jp_tem) 113 ztrds(:,:,:) = tsa(:,:,:,jp_sal) 114 ENDIF 115 116 IF( l_bbl ) CALL bbl( kt, nit000, 'TRA' ) !* bbl coef. and transport (only if not already done in trcbbl) 120 117 121 118 IF( nn_bbl_ldf == 1 ) THEN !* Diffusive bbl 122 119 ! 123 CALL tra_bbl_dif( ts (:,:,:,:,ktlev1), e3t(:,:,:,ktlev2), pts_rhs, jpts )120 CALL tra_bbl_dif( tsb, tsa, jpts ) 124 121 IF( ln_ctl ) & 125 122 CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' bbl_ldf - Ta: ', mask1=tmask, & … … 134 131 IF( nn_bbl_adv /= 0 ) THEN !* Advective bbl 135 132 ! 136 CALL tra_bbl_adv( ts (:,:,:,:,ktlev1), e3t(:,:,:,ktlev2), pts_rhs, jpts )133 CALL tra_bbl_adv( tsb, tsa, jpts ) 137 134 IF(ln_ctl) & 138 135 CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' bbl_adv - Ta: ', mask1=tmask, & … … 146 143 147 144 IF( l_trdtra ) THEN ! send the trends for further diagnostics 148 ztrdt(:,:,:) = pts_rhs(:,:,:,jp_tem) - ztrdt(:,:,:)149 ztrds(:,:,:) = pts_rhs(:,:,:,jp_sal) - ztrds(:,:,:)145 ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:) 146 ztrds(:,:,:) = tsa(:,:,:,jp_sal) - ztrds(:,:,:) 150 147 CALL trd_tra( kt, 'TRA', jp_tem, jptra_bbl, ztrdt ) 151 148 CALL trd_tra( kt, 'TRA', jp_sal, jptra_bbl, ztrds ) … … 158 155 159 156 160 SUBROUTINE tra_bbl_dif( pt , pe3t, pt_rhs, kjpt )157 SUBROUTINE tra_bbl_dif( ptb, pta, kjpt ) 161 158 !!---------------------------------------------------------------------- 162 159 !! *** ROUTINE tra_bbl_dif *** … … 174 171 !! convection is satified) 175 172 !! 176 !! ** Action : pt _rhsincreased by the bbl diffusive trend173 !! ** Action : pta increased by the bbl diffusive trend 177 174 !! 178 175 !! References : Beckmann, A., and R. Doscher, 1997, J. Phys.Oceanogr., 581-591. … … 180 177 !!---------------------------------------------------------------------- 181 178 INTEGER , INTENT(in ) :: kjpt ! number of tracers 182 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: pt ! tracer fields 183 REAL(wp), DIMENSION(jpi,jpj,jpk) , INTENT(in ) :: pe3t ! thickness fields 184 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pt_rhs ! tracer trend 179 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: ptb ! before and now tracer fields 180 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pta ! tracer trend 185 181 ! 186 182 INTEGER :: ji, jj, jn ! dummy loop indices … … 195 191 DO ji = 1, jpi 196 192 ik = mbkt(ji,jj) ! bottom T-level index 197 zptb(ji,jj) = pt (ji,jj,ik,jn) ! bottom before T and S193 zptb(ji,jj) = ptb(ji,jj,ik,jn) ! bottom before T and S 198 194 END DO 199 195 END DO … … 202 198 DO ji = 2, jpim1 203 199 ik = mbkt(ji,jj) ! bottom T-level index 204 pt _rhs(ji,jj,ik,jn) = pt_rhs(ji,jj,ik,jn) &200 pta(ji,jj,ik,jn) = pta(ji,jj,ik,jn) & 205 201 & + ( ahu_bbl(ji ,jj ) * ( zptb(ji+1,jj ) - zptb(ji ,jj ) ) & 206 202 & - ahu_bbl(ji-1,jj ) * ( zptb(ji ,jj ) - zptb(ji-1,jj ) ) & 207 203 & + ahv_bbl(ji ,jj ) * ( zptb(ji ,jj+1) - zptb(ji ,jj ) ) & 208 204 & - ahv_bbl(ji ,jj-1) * ( zptb(ji ,jj ) - zptb(ji ,jj-1) ) ) & 209 & * r1_e1e2t(ji,jj) / pe3t(ji,jj,ik)205 & * r1_e1e2t(ji,jj) / e3t_n(ji,jj,ik) 210 206 END DO 211 207 END DO … … 216 212 217 213 218 SUBROUTINE tra_bbl_adv( pt , pe3t, pt_rhs, kjpt )214 SUBROUTINE tra_bbl_adv( ptb, pta, kjpt ) 219 215 !!---------------------------------------------------------------------- 220 216 !! *** ROUTINE trc_bbl *** … … 232 228 !!---------------------------------------------------------------------- 233 229 INTEGER , INTENT(in ) :: kjpt ! number of tracers 234 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: pt ! before and now tracer fields 235 REAL(wp), DIMENSION(jpi,jpj,jpk) , INTENT(in ) :: pe3t ! thickness fields 236 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pt_rhs ! tracer trend 230 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: ptb ! before and now tracer fields 231 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pta ! tracer trend 237 232 ! 238 233 INTEGER :: ji, jj, jk, jn ! dummy loop indices … … 255 250 ! 256 251 ! ! up -slope T-point (shelf bottom point) 257 zbtr = r1_e1e2t(iis,jj) / pe3t(iis,jj,ikus)258 ztra = zu_bbl * ( pt (iid,jj,ikus,jn) - pt(iis,jj,ikus,jn) ) * zbtr259 pt _rhs(iis,jj,ikus,jn) = pt_rhs(iis,jj,ikus,jn) + ztra252 zbtr = r1_e1e2t(iis,jj) / e3t_n(iis,jj,ikus) 253 ztra = zu_bbl * ( ptb(iid,jj,ikus,jn) - ptb(iis,jj,ikus,jn) ) * zbtr 254 pta(iis,jj,ikus,jn) = pta(iis,jj,ikus,jn) + ztra 260 255 ! 261 256 DO jk = ikus, ikud-1 ! down-slope upper to down T-point (deep column) 262 zbtr = r1_e1e2t(iid,jj) / pe3t(iid,jj,jk)263 ztra = zu_bbl * ( pt (iid,jj,jk+1,jn) - pt(iid,jj,jk,jn) ) * zbtr264 pt _rhs(iid,jj,jk,jn) = pt_rhs(iid,jj,jk,jn) + ztra257 zbtr = r1_e1e2t(iid,jj) / e3t_n(iid,jj,jk) 258 ztra = zu_bbl * ( ptb(iid,jj,jk+1,jn) - ptb(iid,jj,jk,jn) ) * zbtr 259 pta(iid,jj,jk,jn) = pta(iid,jj,jk,jn) + ztra 265 260 END DO 266 261 ! 267 zbtr = r1_e1e2t(iid,jj) / pe3t(iid,jj,ikud)268 ztra = zu_bbl * ( pt (iis,jj,ikus,jn) - pt(iid,jj,ikud,jn) ) * zbtr269 pt _rhs(iid,jj,ikud,jn) = pt_rhs(iid,jj,ikud,jn) + ztra262 zbtr = r1_e1e2t(iid,jj) / e3t_n(iid,jj,ikud) 263 ztra = zu_bbl * ( ptb(iis,jj,ikus,jn) - ptb(iid,jj,ikud,jn) ) * zbtr 264 pta(iid,jj,ikud,jn) = pta(iid,jj,ikud,jn) + ztra 270 265 ENDIF 271 266 ! … … 277 272 ! 278 273 ! up -slope T-point (shelf bottom point) 279 zbtr = r1_e1e2t(ji,ijs) / pe3t(ji,ijs,ikvs)280 ztra = zv_bbl * ( pt (ji,ijd,ikvs,jn) - pt(ji,ijs,ikvs,jn) ) * zbtr281 pt _rhs(ji,ijs,ikvs,jn) = pt_rhs(ji,ijs,ikvs,jn) + ztra274 zbtr = r1_e1e2t(ji,ijs) / e3t_n(ji,ijs,ikvs) 275 ztra = zv_bbl * ( ptb(ji,ijd,ikvs,jn) - ptb(ji,ijs,ikvs,jn) ) * zbtr 276 pta(ji,ijs,ikvs,jn) = pta(ji,ijs,ikvs,jn) + ztra 282 277 ! 283 278 DO jk = ikvs, ikvd-1 ! down-slope upper to down T-point (deep column) 284 zbtr = r1_e1e2t(ji,ijd) / pe3t(ji,ijd,jk)285 ztra = zv_bbl * ( pt (ji,ijd,jk+1,jn) - pt(ji,ijd,jk,jn) ) * zbtr286 pt _rhs(ji,ijd,jk,jn) = pt_rhs(ji,ijd,jk,jn) + ztra279 zbtr = r1_e1e2t(ji,ijd) / e3t_n(ji,ijd,jk) 280 ztra = zv_bbl * ( ptb(ji,ijd,jk+1,jn) - ptb(ji,ijd,jk,jn) ) * zbtr 281 pta(ji,ijd,jk,jn) = pta(ji,ijd,jk,jn) + ztra 287 282 END DO 288 283 ! ! down-slope T-point (deep bottom point) 289 zbtr = r1_e1e2t(ji,ijd) / pe3t(ji,ijd,ikvd)290 ztra = zv_bbl * ( pt (ji,ijs,ikvs,jn) - pt(ji,ijd,ikvd,jn) ) * zbtr291 pt _rhs(ji,ijd,ikvd,jn) = pt_rhs(ji,ijd,ikvd,jn) + ztra284 zbtr = r1_e1e2t(ji,ijd) / e3t_n(ji,ijd,ikvd) 285 ztra = zv_bbl * ( ptb(ji,ijs,ikvs,jn) - ptb(ji,ijd,ikvd,jn) ) * zbtr 286 pta(ji,ijd,ikvd,jn) = pta(ji,ijd,ikvd,jn) + ztra 292 287 ENDIF 293 288 END DO … … 300 295 301 296 302 SUBROUTINE bbl( kt, kit000, ktlev1, ktlev2, kt2lev,cdtype )297 SUBROUTINE bbl( kt, kit000, cdtype ) 303 298 !!---------------------------------------------------------------------- 304 299 !! *** ROUTINE bbl *** … … 328 323 INTEGER , INTENT(in ) :: kt ! ocean time-step index 329 324 INTEGER , INTENT(in ) :: kit000 ! first time step index 330 INTEGER , INTENT(in ) :: ktlev1, ktlev2 ! time level indices for 3-time-levelsource terms331 INTEGER , INTENT(in ) :: kt2lev ! time level index for 2-time-level source terms332 325 CHARACTER(len=3), INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator) 333 326 ! … … 351 344 DO ji = 1, jpi 352 345 ik = mbkt(ji,jj) ! bottom T-level index 353 zts (ji,jj,jp_tem) = ts (ji,jj,ik,jp_tem,ktlev1) ! bottom before T and S354 zts (ji,jj,jp_sal) = ts (ji,jj,ik,jp_sal,ktlev1)346 zts (ji,jj,jp_tem) = tsb(ji,jj,ik,jp_tem) ! bottom before T and S 347 zts (ji,jj,jp_sal) = tsb(ji,jj,ik,jp_sal) 355 348 ! 356 zdep(ji,jj) = gdept (ji,jj,ik,kt2lev) ! bottom T-level reference depth357 zub (ji,jj) = u u(ji,jj,mbku(ji,jj),ktlev2) ! bottom velocity358 zvb (ji,jj) = v v(ji,jj,mbkv(ji,jj),ktlev2)349 zdep(ji,jj) = gdept_n(ji,jj,ik) ! bottom T-level reference depth 350 zub (ji,jj) = un(ji,jj,mbku(ji,jj)) ! bottom velocity 351 zvb (ji,jj) = vn(ji,jj,mbkv(ji,jj)) 359 352 END DO 360 353 END DO -
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/tradmp.F90
r10806 r10874 72 72 73 73 74 SUBROUTINE tra_dmp( kt , ktlev, kt2lev, pts_rhs)74 SUBROUTINE tra_dmp( kt ) 75 75 !!---------------------------------------------------------------------- 76 76 !! *** ROUTINE tra_dmp *** … … 90 90 !! ** Action : - tsa: tracer trends updated with the damping trend 91 91 !!---------------------------------------------------------------------- 92 INTEGER, INTENT(in) :: kt ! ocean time-step index 93 INTEGER, INTENT(in) :: ktlev ! time level index for 3-time-level source terms 94 INTEGER, INTENT(in) :: kt2lev ! time level index for 2-time-level source terms 95 REAL(wp), INTENT( inout), DIMENSION(jpi,jpj,jpk,jpts) :: pts_rhs ! temperature and salinity trends 92 INTEGER, INTENT(in) :: kt ! ocean time-step index 96 93 ! 97 94 INTEGER :: ji, jj, jk, jn ! dummy loop indices … … 104 101 IF( l_trdtra ) THEN !* Save ta and sa trends 105 102 ALLOCATE( ztrdts(jpi,jpj,jpk,jpts) ) 106 ztrdts(:,:,:,:) = pts_rhs(:,:,:,:)103 ztrdts(:,:,:,:) = tsa(:,:,:,:) 107 104 ENDIF 108 105 ! !== input T-S data at kt ==! … … 116 113 DO jj = 2, jpjm1 117 114 DO ji = fs_2, fs_jpim1 ! vector opt. 118 pts_rhs(ji,jj,jk,jn) = pts_rhs(ji,jj,jk,jn) + resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jn) - ts(ji,jj,jk,jn,ktlev) )115 tsa(ji,jj,jk,jn) = tsa(ji,jj,jk,jn) + resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jn) - tsb(ji,jj,jk,jn) ) 119 116 END DO 120 117 END DO … … 127 124 DO ji = fs_2, fs_jpim1 ! vector opt. 128 125 IF( avt(ji,jj,jk) <= avt_c ) THEN 129 pts_rhs(ji,jj,jk,jp_tem) = pts_rhs(ji,jj,jk,jp_tem) &130 & + resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jp_tem) - ts (ji,jj,jk,jp_tem,ktlev) )131 pts_rhs(ji,jj,jk,jp_sal) = pts_rhs(ji,jj,jk,jp_sal) &132 & + resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jp_sal) - ts (ji,jj,jk,jp_sal,ktlev) )126 tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem) & 127 & + resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jp_tem) - tsb(ji,jj,jk,jp_tem) ) 128 tsa(ji,jj,jk,jp_sal) = tsa(ji,jj,jk,jp_sal) & 129 & + resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jp_sal) - tsb(ji,jj,jk,jp_sal) ) 133 130 ENDIF 134 131 END DO … … 140 137 DO jj = 2, jpjm1 141 138 DO ji = fs_2, fs_jpim1 ! vector opt. 142 IF( gdept (ji,jj,jk,kt2lev) >= hmlp (ji,jj) ) THEN143 pts_rhs(ji,jj,jk,jp_tem) = pts_rhs(ji,jj,jk,jp_tem) &144 & + resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jp_tem) - ts (ji,jj,jk,jp_tem,ktlev) )145 pts_rhs(ji,jj,jk,jp_sal) = pts_rhs(ji,jj,jk,jp_sal) &146 & + resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jp_sal) - ts (ji,jj,jk,jp_sal,ktlev) )139 IF( gdept_n(ji,jj,jk) >= hmlp (ji,jj) ) THEN 140 tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem) & 141 & + resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jp_tem) - tsb(ji,jj,jk,jp_tem) ) 142 tsa(ji,jj,jk,jp_sal) = tsa(ji,jj,jk,jp_sal) & 143 & + resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jp_sal) - tsb(ji,jj,jk,jp_sal) ) 147 144 ENDIF 148 145 END DO … … 153 150 ! 154 151 IF( l_trdtra ) THEN ! trend diagnostic 155 ztrdts(:,:,:,:) = pts_rhs(:,:,:,:) - ztrdts(:,:,:,:)152 ztrdts(:,:,:,:) = tsa(:,:,:,:) - ztrdts(:,:,:,:) 156 153 CALL trd_tra( kt, 'TRA', jp_tem, jptra_dmp, ztrdts(:,:,:,jp_tem) ) 157 154 CALL trd_tra( kt, 'TRA', jp_sal, jptra_dmp, ztrdts(:,:,:,jp_sal) ) -
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traldf.F90
r10806 r10874 47 47 CONTAINS 48 48 49 SUBROUTINE tra_ldf( kt , ktlev1, ktlev2, kt2lev, pts_rhs)49 SUBROUTINE tra_ldf( kt ) 50 50 !!---------------------------------------------------------------------- 51 51 !! *** ROUTINE tra_ldf *** … … 53 53 !! ** Purpose : compute the lateral ocean tracer physics. 54 54 !!---------------------------------------------------------------------- 55 INTEGER, INTENT( in ) :: kt ! ocean time-step index 56 INTEGER, INTENT( in ) :: ktlev1, ktlev2 ! time level indices for 3-time-level source terms 57 INTEGER, INTENT( in ) :: kt2lev ! time level index for 2-time-level source terms 58 REAL(wp), INTENT( inout), DIMENSION(jpi,jpj,jpk,jpts) :: pts_rhs ! temperature and salinity trends 55 INTEGER, INTENT( in ) :: kt ! ocean time-step index 59 56 !! 60 57 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: ztrdt, ztrds … … 65 62 IF( l_trdtra ) THEN !* Save ta and sa trends 66 63 ALLOCATE( ztrdt(jpi,jpj,jpk) , ztrds(jpi,jpj,jpk) ) 67 ztrdt(:,:,:) = pts_rhs(:,:,:,jp_tem)68 ztrds(:,:,:) = pts_rhs(:,:,:,jp_sal)64 ztrdt(:,:,:) = tsa(:,:,:,jp_tem) 65 ztrds(:,:,:) = tsa(:,:,:,jp_sal) 69 66 ENDIF 70 67 ! 71 68 SELECT CASE ( nldf_tra ) !* compute lateral mixing trend and add it to the general trend 72 69 CASE ( np_lap ) ! laplacian: iso-level operator 73 CALL tra_ldf_lap ( kt, nit000, ktlev2, kt2lev, 'TRA', ahtu, ahtv, gtsu, gtsv, gtui, gtvi, ts(:,:,:,:,ktlev1), pts_rhs, jpts, 1 )70 CALL tra_ldf_lap ( kt, nit000,'TRA', ahtu, ahtv, gtsu, gtsv, gtui, gtvi, tsb, tsa, jpts, 1 ) 74 71 CASE ( np_lap_i ) ! laplacian: standard iso-neutral operator (Madec) 75 CALL tra_ldf_iso ( kt, nit000, ktlev2, kt2lev, 'TRA', ahtu, ahtv, gtsu, gtsv, gtui, gtvi, ts(:,:,:,:,ktlev1), ts(:,:,:,:,ktlev1), pts_rhs, jpts, 1 )72 CALL tra_ldf_iso ( kt, nit000,'TRA', ahtu, ahtv, gtsu, gtsv, gtui, gtvi, tsb, tsb, tsa, jpts, 1 ) 76 73 CASE ( np_lap_it ) ! laplacian: triad iso-neutral operator (griffies) 77 CALL tra_ldf_triad( kt, nit000, ktlev2, kt2lev, 'TRA', ahtu, ahtv, gtsu, gtsv, gtui, gtvi, ts(:,:,:,:,ktlev1), ts(:,:,:,:,ktlev1), pts_rhs, jpts, 1 )74 CALL tra_ldf_triad( kt, nit000,'TRA', ahtu, ahtv, gtsu, gtsv, gtui, gtvi, tsb, tsb, tsa, jpts, 1 ) 78 75 CASE ( np_blp , np_blp_i , np_blp_it ) ! bilaplacian: iso-level & iso-neutral operators 79 CALL tra_ldf_blp ( kt, nit000, ktlev2, kt2lev, 'TRA', ahtu, ahtv, gtsu, gtsv, gtui, gtvi, ts(:,:,:,:,ktlev1) , pts_rhs, jpts, nldf_tra )76 CALL tra_ldf_blp ( kt, nit000,'TRA', ahtu, ahtv, gtsu, gtsv, gtui, gtvi, tsb , tsa, jpts, nldf_tra ) 80 77 END SELECT 81 78 ! 82 79 IF( l_trdtra ) THEN !* save the horizontal diffusive trends for further diagnostics 83 ztrdt(:,:,:) = pts_rhs(:,:,:,jp_tem) - ztrdt(:,:,:)84 ztrds(:,:,:) = pts_rhs(:,:,:,jp_sal) - ztrds(:,:,:)80 ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:) 81 ztrds(:,:,:) = tsa(:,:,:,jp_sal) - ztrds(:,:,:) 85 82 CALL trd_tra( kt, 'TRA', jp_tem, jptra_ldf, ztrdt ) 86 83 CALL trd_tra( kt, 'TRA', jp_sal, jptra_ldf, ztrds ) -
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traldf_iso.F90
r10806 r10874 48 48 CONTAINS 49 49 50 SUBROUTINE tra_ldf_iso( kt, kit000, ktlev, kt2lev,cdtype, pahu, pahv, pgu , pgv , &50 SUBROUTINE tra_ldf_iso( kt, kit000, cdtype, pahu, pahv, pgu , pgv , & 51 51 & pgui, pgvi, & 52 & pt , pt_lev0, pt_rhs, kjpt, kpass )52 & ptb , ptbb, pta , kjpt, kpass ) 53 53 !!---------------------------------------------------------------------- 54 54 !! *** ROUTINE tra_ldf_iso *** … … 87 87 !! difft = 1/(e1e2t*e3t) dk[ zftw ] 88 88 !! Add this trend to the general trend (ta,sa): 89 !! pt _rhs = pt_rhs+ difft90 !! 91 !! ** Action : Update pt _rhsarrays with the before rotated diffusion89 !! pta = pta + difft 90 !! 91 !! ** Action : Update pta arrays with the before rotated diffusion 92 92 !!---------------------------------------------------------------------- 93 93 INTEGER , INTENT(in ) :: kt ! ocean time-step index 94 94 INTEGER , INTENT(in ) :: kit000 ! first time step index 95 INTEGER , INTENT(in ) :: ktlev ! time level index for e3t96 INTEGER , INTENT(in ) :: kt2lev ! time level index for 2-time-level thicknesses97 95 CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator) 98 96 INTEGER , INTENT(in ) :: kjpt ! number of tracers … … 101 99 REAL(wp), DIMENSION(jpi,jpj ,kjpt), INTENT(in ) :: pgu, pgv ! tracer gradient at pstep levels 102 100 REAL(wp), DIMENSION(jpi,jpj, kjpt), INTENT(in ) :: pgui, pgvi ! tracer gradient at top levels 103 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: pt ! tracer (kpass=1) or laplacian of tracer (kpass=2)104 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: pt _lev0! tracer (only used in kpass=2)105 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pt _rhs! tracer trend101 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: ptb ! tracer (kpass=1) or laplacian of tracer (kpass=2) 102 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: ptbb ! tracer (only used in kpass=2) 103 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pta ! tracer trend 106 104 ! 107 105 INTEGER :: ji, jj, jk, jn ! dummy loop indices … … 184 182 DO ji = 1, fs_jpim1 185 183 akz(ji,jj,jk) = 16._wp * ah_wslp2(ji,jj,jk) & 186 & * ( akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ( e3w (ji,jj,jk,kt2lev) * e3w(ji,jj,jk,kt2lev) ) )184 & * ( akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ( e3w_n(ji,jj,jk) * e3w_n(ji,jj,jk) ) ) 187 185 END DO 188 186 END DO … … 192 190 DO jj = 1, jpjm1 193 191 DO ji = 1, fs_jpim1 194 ze3w_2 = e3w (ji,jj,jk,kt2lev) * e3w(ji,jj,jk,kt2lev)192 ze3w_2 = e3w_n(ji,jj,jk) * e3w_n(ji,jj,jk) 195 193 zcoef0 = z2dt * ( akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ze3w_2 ) 196 194 akz(ji,jj,jk) = MAX( zcoef0 - 0.5_wp , 0._wp ) * ze3w_2 * z1_2dt … … 221 219 DO jj = 1, jpjm1 222 220 DO ji = 1, fs_jpim1 ! vector opt. 223 zdit(ji,jj,jk) = ( pt (ji+1,jj ,jk,jn) - pt(ji,jj,jk,jn) ) * umask(ji,jj,jk)224 zdjt(ji,jj,jk) = ( pt (ji ,jj+1,jk,jn) - pt(ji,jj,jk,jn) ) * vmask(ji,jj,jk)221 zdit(ji,jj,jk) = ( ptb(ji+1,jj ,jk,jn) - ptb(ji,jj,jk,jn) ) * umask(ji,jj,jk) 222 zdjt(ji,jj,jk) = ( ptb(ji ,jj+1,jk,jn) - ptb(ji,jj,jk,jn) ) * vmask(ji,jj,jk) 225 223 END DO 226 224 END DO … … 250 248 ! 251 249 ! !== Vertical tracer gradient 252 zdk1t(:,:) = ( pt (:,:,jk,jn) - pt(:,:,jk+1,jn) ) * wmask(:,:,jk+1) ! level jk+1250 zdk1t(:,:) = ( ptb(:,:,jk,jn) - ptb(:,:,jk+1,jn) ) * wmask(:,:,jk+1) ! level jk+1 253 251 ! 254 252 IF( jk == 1 ) THEN ; zdkt(:,:) = zdk1t(:,:) ! surface: zdkt(jk=1)=zdkt(jk=2) 255 ELSE ; zdkt(:,:) = ( pt (:,:,jk-1,jn) - pt(:,:,jk,jn) ) * wmask(:,:,jk)253 ELSE ; zdkt(:,:) = ( ptb(:,:,jk-1,jn) - ptb(:,:,jk,jn) ) * wmask(:,:,jk) 256 254 ENDIF 257 255 DO jj = 1 , jpjm1 !== Horizontal fluxes 258 256 DO ji = 1, fs_jpim1 ! vector opt. 259 zabe1 = pahu(ji,jj,jk) * e2_e1u(ji,jj) * e3u (ji,jj,jk,ktlev)260 zabe2 = pahv(ji,jj,jk) * e1_e2v(ji,jj) * e3v (ji,jj,jk,ktlev)257 zabe1 = pahu(ji,jj,jk) * e2_e1u(ji,jj) * e3u_n(ji,jj,jk) 258 zabe2 = pahv(ji,jj,jk) * e1_e2v(ji,jj) * e3v_n(ji,jj,jk) 261 259 ! 262 260 zmsku = 1. / MAX( wmask(ji+1,jj,jk ) + wmask(ji,jj,jk+1) & … … 278 276 END DO 279 277 ! 280 DO jj = 2 , jpjm1 !== horizontal divergence and add to pt _rhs278 DO jj = 2 , jpjm1 !== horizontal divergence and add to pta 281 279 DO ji = fs_2, fs_jpim1 ! vector opt. 282 pt _rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn) + zsign * ( zftu(ji,jj,jk) - zftu(ji-1,jj,jk) &280 pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + zsign * ( zftu(ji,jj,jk) - zftu(ji-1,jj,jk) & 283 281 & + zftv(ji,jj,jk) - zftv(ji,jj-1,jk) ) & 284 & * r1_e1e2t(ji,jj) / e3t (ji,jj,jk,ktlev)282 & * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) 285 283 END DO 286 284 END DO … … 327 325 DO jj = 1, jpjm1 328 326 DO ji = fs_2, fs_jpim1 329 ztfw(ji,jj,jk) = ztfw(ji,jj,jk) + e1e2t(ji,jj) / e3w (ji,jj,jk,kt2lev) * wmask(ji,jj,jk) &327 ztfw(ji,jj,jk) = ztfw(ji,jj,jk) + e1e2t(ji,jj) / e3w_n(ji,jj,jk) * wmask(ji,jj,jk) & 330 328 & * ( ah_wslp2(ji,jj,jk) - akz(ji,jj,jk) ) & 331 & * ( pt (ji,jj,jk-1,jn) - pt(ji,jj,jk,jn) )329 & * ( ptb(ji,jj,jk-1,jn) - ptb(ji,jj,jk,jn) ) 332 330 END DO 333 331 END DO … … 342 340 ztfw(ji,jj,jk) = ztfw(ji,jj,jk) & 343 341 & + ah_wslp2(ji,jj,jk) * e1e2t(ji,jj) & 344 & * ( pt (ji,jj,jk-1,jn) - pt(ji,jj,jk,jn) ) / e3w(ji,jj,jk,kt2lev) * wmask(ji,jj,jk)342 & * ( ptb(ji,jj,jk-1,jn) - ptb(ji,jj,jk,jn) ) / e3w_n(ji,jj,jk) * wmask(ji,jj,jk) 345 343 END DO 346 344 END DO 347 345 END DO 348 CASE( 2 ) ! 2nd pass : eddy flux = ah_wslp2 and akz applied on pt and pt_lev0gradients, resp.346 CASE( 2 ) ! 2nd pass : eddy flux = ah_wslp2 and akz applied on ptb and ptbb gradients, resp. 349 347 DO jk = 2, jpkm1 350 348 DO jj = 1, jpjm1 351 349 DO ji = fs_2, fs_jpim1 352 ztfw(ji,jj,jk) = ztfw(ji,jj,jk) + e1e2t(ji,jj) / e3w (ji,jj,jk,kt2lev) * wmask(ji,jj,jk) &353 & * ( ah_wslp2(ji,jj,jk) * ( pt (ji,jj,jk-1,jn) - pt(ji,jj,jk,jn) ) &354 & + akz (ji,jj,jk) * ( pt _lev0(ji,jj,jk-1,jn) - pt_lev0(ji,jj,jk,jn) ) )350 ztfw(ji,jj,jk) = ztfw(ji,jj,jk) + e1e2t(ji,jj) / e3w_n(ji,jj,jk) * wmask(ji,jj,jk) & 351 & * ( ah_wslp2(ji,jj,jk) * ( ptb (ji,jj,jk-1,jn) - ptb (ji,jj,jk,jn) ) & 352 & + akz (ji,jj,jk) * ( ptbb(ji,jj,jk-1,jn) - ptbb(ji,jj,jk,jn) ) ) 355 353 END DO 356 354 END DO … … 359 357 ENDIF 360 358 ! 361 DO jk = 1, jpkm1 !== Divergence of vertical fluxes added to pt _rhs==!359 DO jk = 1, jpkm1 !== Divergence of vertical fluxes added to pta ==! 362 360 DO jj = 2, jpjm1 363 361 DO ji = fs_2, fs_jpim1 ! vector opt. 364 pt _rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn) + zsign * ( ztfw (ji,jj,jk) - ztfw(ji,jj,jk+1) ) &365 & * r1_e1e2t(ji,jj) / e3t (ji,jj,jk,ktlev)362 pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + zsign * ( ztfw (ji,jj,jk) - ztfw(ji,jj,jk+1) ) & 363 & * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) 366 364 END DO 367 365 END DO -
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traldf_lap_blp.F90
r10806 r10874 45 45 CONTAINS 46 46 47 SUBROUTINE tra_ldf_lap( kt, kit000, ktlev, kt2lev,cdtype, pahu, pahv, pgu , pgv , &47 SUBROUTINE tra_ldf_lap( kt, kit000, cdtype, pahu, pahv, pgu , pgv , & 48 48 & pgui, pgvi, & 49 & pt , pt_rhs, kjpt, kpass )49 & ptb , pta , kjpt, kpass ) 50 50 !!---------------------------------------------------------------------- 51 51 !! *** ROUTINE tra_ldf_lap *** … … 59 59 !! difft = 1/(e1e2t*e3t) { di-1[ pahu e2u*e3u/e1u di(tb) ] 60 60 !! + dj-1[ pahv e1v*e3v/e2v dj(tb) ] } 61 !! Add this trend to the general tracer trend pt _rhs:62 !! pt _rhs = pt_rhs+ difft63 !! 64 !! ** Action : - Update pt _rhsarrays with the before iso-level61 !! Add this trend to the general tracer trend pta : 62 !! pta = pta + difft 63 !! 64 !! ** Action : - Update pta arrays with the before iso-level 65 65 !! harmonic mixing trend. 66 66 !!---------------------------------------------------------------------- 67 67 INTEGER , INTENT(in ) :: kt ! ocean time-step index 68 68 INTEGER , INTENT(in ) :: kit000 ! first time step index 69 INTEGER , INTENT(in ) :: ktlev ! time level index for e3t70 INTEGER , INTENT(in ) :: kt2lev ! time level index for 2-time-level thicknesses71 69 CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator) 72 70 INTEGER , INTENT(in ) :: kjpt ! number of tracers … … 75 73 REAL(wp), DIMENSION(jpi,jpj ,kjpt), INTENT(in ) :: pgu, pgv ! tracer gradient at pstep levels 76 74 REAL(wp), DIMENSION(jpi,jpj, kjpt), INTENT(in ) :: pgui, pgvi ! tracer gradient at top levels 77 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: pt ! before and now tracer fields78 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pt _rhs! tracer trend75 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: ptb ! before and now tracer fields 76 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pta ! tracer trend 79 77 ! 80 78 INTEGER :: ji, jj, jk, jn ! dummy loop indices … … 102 100 DO jj = 1, jpjm1 103 101 DO ji = 1, fs_jpim1 ! vector opt. 104 zaheeu(ji,jj,jk) = zsign * pahu(ji,jj,jk) * e2_e1u(ji,jj) * e3u (ji,jj,jk,ktlev) !!gm * umask(ji,jj,jk) pah masked!105 zaheev(ji,jj,jk) = zsign * pahv(ji,jj,jk) * e1_e2v(ji,jj) * e3v (ji,jj,jk,ktlev) !!gm * vmask(ji,jj,jk)102 zaheeu(ji,jj,jk) = zsign * pahu(ji,jj,jk) * e2_e1u(ji,jj) * e3u_n(ji,jj,jk) !!gm * umask(ji,jj,jk) pah masked! 103 zaheev(ji,jj,jk) = zsign * pahv(ji,jj,jk) * e1_e2v(ji,jj) * e3v_n(ji,jj,jk) !!gm * vmask(ji,jj,jk) 106 104 END DO 107 105 END DO … … 115 113 DO jj = 1, jpjm1 116 114 DO ji = 1, fs_jpim1 117 ztu(ji,jj,jk) = zaheeu(ji,jj,jk) * ( pt (ji+1,jj ,jk,jn) - pt(ji,jj,jk,jn) )118 ztv(ji,jj,jk) = zaheev(ji,jj,jk) * ( pt (ji ,jj+1,jk,jn) - pt(ji,jj,jk,jn) )115 ztu(ji,jj,jk) = zaheeu(ji,jj,jk) * ( ptb(ji+1,jj ,jk,jn) - ptb(ji,jj,jk,jn) ) 116 ztv(ji,jj,jk) = zaheev(ji,jj,jk) * ( ptb(ji ,jj+1,jk,jn) - ptb(ji,jj,jk,jn) ) 119 117 END DO 120 118 END DO … … 140 138 DO jj = 2, jpjm1 141 139 DO ji = fs_2, fs_jpim1 142 pt _rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn) + ( ztu(ji,jj,jk) - ztu(ji-1,jj,jk) &140 pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ( ztu(ji,jj,jk) - ztu(ji-1,jj,jk) & 143 141 & + ztv(ji,jj,jk) - ztv(ji,jj-1,jk) ) & 144 & / ( e1e2t(ji,jj) * e3t (ji,jj,jk,ktlev) )142 & / ( e1e2t(ji,jj) * e3t_n(ji,jj,jk) ) 145 143 END DO 146 144 END DO … … 161 159 162 160 163 SUBROUTINE tra_ldf_blp( kt, kit000, ktlev, kt2lev,cdtype, pahu, pahv, pgu , pgv , &161 SUBROUTINE tra_ldf_blp( kt, kit000, cdtype, pahu, pahv, pgu , pgv , & 164 162 & pgui, pgvi, & 165 & pt , pt_rhs, kjpt, kldf )163 & ptb , pta , kjpt, kldf ) 166 164 !!---------------------------------------------------------------------- 167 165 !! *** ROUTINE tra_ldf_blp *** … … 174 172 !! It is computed by two successive calls to laplacian routine 175 173 !! 176 !! ** Action : pt _rhsupdated with the before rotated bilaplacian diffusion174 !! ** Action : pta updated with the before rotated bilaplacian diffusion 177 175 !!---------------------------------------------------------------------- 178 176 INTEGER , INTENT(in ) :: kt ! ocean time-step index 179 177 INTEGER , INTENT(in ) :: kit000 ! first time step index 180 INTEGER , INTENT(in ) :: ktlev ! time level index for e3t181 INTEGER , INTENT(in ) :: kt2lev ! time level index for 2-time-level thicknesses182 178 CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator) 183 179 INTEGER , INTENT(in ) :: kjpt ! number of tracers … … 186 182 REAL(wp), DIMENSION(jpi,jpj ,kjpt), INTENT(in ) :: pgu, pgv ! tracer gradient at pstep levels 187 183 REAL(wp), DIMENSION(jpi,jpj, kjpt), INTENT(in ) :: pgui, pgvi ! tracer gradient at top levels 188 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: pt ! before and now tracer fields189 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pt _rhs! tracer trend184 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: ptb ! before and now tracer fields 185 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pta ! tracer trend 190 186 ! 191 187 INTEGER :: ji, jj, jk, jn ! dummy loop indices … … 207 203 zlap(:,:,:,:) = 0._wp 208 204 ! 209 SELECT CASE ( kldf ) !== 1st laplacian applied to pt (output in zlap) ==!205 SELECT CASE ( kldf ) !== 1st laplacian applied to ptb (output in zlap) ==! 210 206 ! 211 207 CASE ( np_blp ) ! iso-level bilaplacian 212 CALL tra_ldf_lap ( kt, kit000, ktlev, kt2lev, cdtype, pahu, pahv, pgu, pgv, pgui, pgvi, pt, zlap, kjpt, 1 )208 CALL tra_ldf_lap ( kt, kit000, cdtype, pahu, pahv, pgu, pgv, pgui, pgvi, ptb, zlap, kjpt, 1 ) 213 209 CASE ( np_blp_i ) ! rotated bilaplacian : standard operator (Madec) 214 CALL tra_ldf_iso ( kt, kit000, ktlev, kt2lev, cdtype, pahu, pahv, pgu, pgv, pgui, pgvi, pt, pt, zlap, kjpt, 1 )210 CALL tra_ldf_iso ( kt, kit000, cdtype, pahu, pahv, pgu, pgv, pgui, pgvi, ptb, ptb, zlap, kjpt, 1 ) 215 211 CASE ( np_blp_it ) ! rotated bilaplacian : triad operator (griffies) 216 CALL tra_ldf_triad( kt, kit000, ktlev, kt2lev, cdtype, pahu, pahv, pgu, pgv, pgui, pgvi, pt, pt, zlap, kjpt, 1 )212 CALL tra_ldf_triad( kt, kit000, cdtype, pahu, pahv, pgu, pgv, pgui, pgvi, ptb, ptb, zlap, kjpt, 1 ) 217 213 END SELECT 218 214 ! … … 223 219 ENDIF 224 220 ! 225 SELECT CASE ( kldf ) !== 2nd laplacian applied to zlap (output in pt _rhs) ==!221 SELECT CASE ( kldf ) !== 2nd laplacian applied to zlap (output in pta) ==! 226 222 ! 227 223 CASE ( np_blp ) ! iso-level bilaplacian 228 CALL tra_ldf_lap ( kt, kit000, ktlev, kt2lev, cdtype, pahu, pahv, zglu, zglv, zgui, zgvi, zlap, pt_rhs, kjpt, 2 )224 CALL tra_ldf_lap ( kt, kit000, cdtype, pahu, pahv, zglu, zglv, zgui, zgvi, zlap, pta, kjpt, 2 ) 229 225 CASE ( np_blp_i ) ! rotated bilaplacian : standard operator (Madec) 230 CALL tra_ldf_iso ( kt, kit000, ktlev, kt2lev, cdtype, pahu, pahv, zglu, zglv, zgui, zgvi, zlap, pt, pt_rhs, kjpt, 2 )226 CALL tra_ldf_iso ( kt, kit000, cdtype, pahu, pahv, zglu, zglv, zgui, zgvi, zlap, ptb, pta, kjpt, 2 ) 231 227 CASE ( np_blp_it ) ! rotated bilaplacian : triad operator (griffies) 232 CALL tra_ldf_triad( kt, kit000, ktlev, kt2lev, cdtype, pahu, pahv, zglu, zglv, zgui, zgvi, zlap, pt, pt_rhs, kjpt, 2 )228 CALL tra_ldf_triad( kt, kit000, cdtype, pahu, pahv, zglu, zglv, zgui, zgvi, zlap, ptb, pta, kjpt, 2 ) 233 229 END SELECT 234 230 ! -
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traldf_triad.F90
r10806 r10874 48 48 CONTAINS 49 49 50 SUBROUTINE tra_ldf_triad( kt, kit000, ktlev, kt2lev,cdtype, pahu, pahv, pgu , pgv , &50 SUBROUTINE tra_ldf_triad( kt, kit000, cdtype, pahu, pahv, pgu , pgv , & 51 51 & pgui, pgvi, & 52 & pt , pt_lev0, pt_rhs, kjpt, kpass )52 & ptb , ptbb, pta , kjpt, kpass ) 53 53 !!---------------------------------------------------------------------- 54 54 !! *** ROUTINE tra_ldf_triad *** … … 66 66 !! see documentation for the desciption 67 67 !! 68 !! ** Action : pt _rhsupdated with the before rotated diffusion68 !! ** Action : pta updated with the before rotated diffusion 69 69 !! ah_wslp2 .... 70 70 !! akz stabilizing vertical diffusivity coefficient (used in trazdf_imp) … … 72 72 INTEGER , INTENT(in ) :: kt ! ocean time-step index 73 73 INTEGER , INTENT(in ) :: kit000 ! first time step index 74 INTEGER , INTENT(in ) :: ktlev ! time level index for e3t75 INTEGER , INTENT(in ) :: kt2lev ! time level index for 2-time-level thicknesses76 74 CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator) 77 75 INTEGER , INTENT(in ) :: kjpt ! number of tracers … … 80 78 REAL(wp), DIMENSION(jpi,jpj ,kjpt), INTENT(in ) :: pgu , pgv ! tracer gradient at pstep levels 81 79 REAL(wp), DIMENSION(jpi,jpj, kjpt), INTENT(in ) :: pgui, pgvi ! tracer gradient at top levels 82 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: pt ! tracer (kpass=1) or laplacian of tracer (kpass=2)83 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: pt _lev0! tracer (only used in kpass=2)84 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pt _rhs! tracer trend80 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: ptb ! tracer (kpass=1) or laplacian of tracer (kpass=2) 81 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: ptbb ! tracer (only used in kpass=2) 82 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pta ! tracer trend 85 83 ! 86 84 INTEGER :: ji, jj, jk, jn ! dummy loop indices … … 144 142 DO jj = 1, jpjm1 145 143 DO ji = 1, fs_jpim1 146 ze3wr = 1._wp / e3w (ji+ip,jj,jk+kp,kt2lev)147 zbu = e1e2u(ji,jj) * e3u (ji,jj,jk,ktlev)144 ze3wr = 1._wp / e3w_n(ji+ip,jj,jk+kp) 145 zbu = e1e2u(ji,jj) * e3u_n(ji,jj,jk) 148 146 zah = 0.25_wp * pahu(ji,jj,jk) 149 147 zslope_skew = triadi_g(ji+ip,jj,jk,1-ip,kp) 150 148 ! Subtract s-coordinate slope at t-points to give slope rel to s-surfaces (do this by *adding* gradient of depth) 151 zslope2 = zslope_skew + ( gdept (ji+1,jj,jk,kt2lev) - gdept(ji,jj,jk,kt2lev) ) * r1_e1u(ji,jj) * umask(ji,jj,jk+kp)149 zslope2 = zslope_skew + ( gdept_n(ji+1,jj,jk) - gdept_n(ji,jj,jk) ) * r1_e1u(ji,jj) * umask(ji,jj,jk+kp) 152 150 zslope2 = zslope2 *zslope2 153 151 ah_wslp2(ji+ip,jj,jk+kp) = ah_wslp2(ji+ip,jj,jk+kp) + zah * zbu * ze3wr * r1_e1e2t(ji+ip,jj) * zslope2 … … 168 166 DO jj = 1, jpjm1 169 167 DO ji = 1, fs_jpim1 170 ze3wr = 1.0_wp / e3w (ji,jj+jp,jk+kp,kt2lev)171 zbv = e1e2v(ji,jj) * e3v (ji,jj,jk,ktlev)168 ze3wr = 1.0_wp / e3w_n(ji,jj+jp,jk+kp) 169 zbv = e1e2v(ji,jj) * e3v_n(ji,jj,jk) 172 170 zah = 0.25_wp * pahv(ji,jj,jk) 173 171 zslope_skew = triadj_g(ji,jj+jp,jk,1-jp,kp) 174 172 ! Subtract s-coordinate slope at t-points to give slope rel to s surfaces 175 173 ! (do this by *adding* gradient of depth) 176 zslope2 = zslope_skew + ( gdept (ji,jj+1,jk,kt2lev) - gdept(ji,jj,jk,kt2lev) ) * r1_e2v(ji,jj) * vmask(ji,jj,jk+kp)174 zslope2 = zslope_skew + ( gdept_n(ji,jj+1,jk) - gdept_n(ji,jj,jk) ) * r1_e2v(ji,jj) * vmask(ji,jj,jk+kp) 177 175 zslope2 = zslope2 * zslope2 178 176 ah_wslp2(ji,jj+jp,jk+kp) = ah_wslp2(ji,jj+jp,jk+kp) + zah * zbv * ze3wr * r1_e1e2t(ji,jj+jp) * zslope2 … … 195 193 DO ji = 1, fs_jpim1 196 194 akz(ji,jj,jk) = 16._wp * ah_wslp2(ji,jj,jk) & 197 & * ( akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ( e3w (ji,jj,jk,kt2lev) * e3w(ji,jj,jk,kt2lev) ) )195 & * ( akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ( e3w_n(ji,jj,jk) * e3w_n(ji,jj,jk) ) ) 198 196 END DO 199 197 END DO … … 203 201 DO jj = 1, jpjm1 204 202 DO ji = 1, fs_jpim1 205 ze3w_2 = e3w (ji,jj,jk,kt2lev) * e3w(ji,jj,jk,kt2lev)203 ze3w_2 = e3w_n(ji,jj,jk) * e3w_n(ji,jj,jk) 206 204 zcoef0 = z2dt * ( akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ze3w_2 ) 207 205 akz(ji,jj,jk) = MAX( zcoef0 - 0.5_wp , 0._wp ) * ze3w_2 * z1_2dt … … 231 229 DO jj = 1, jpjm1 232 230 DO ji = 1, fs_jpim1 ! vector opt. 233 zdit(ji,jj,jk) = ( pt (ji+1,jj ,jk,jn) - pt(ji,jj,jk,jn) ) * umask(ji,jj,jk)234 zdjt(ji,jj,jk) = ( pt (ji ,jj+1,jk,jn) - pt(ji,jj,jk,jn) ) * vmask(ji,jj,jk)231 zdit(ji,jj,jk) = ( ptb(ji+1,jj ,jk,jn) - ptb(ji,jj,jk,jn) ) * umask(ji,jj,jk) 232 zdjt(ji,jj,jk) = ( ptb(ji ,jj+1,jk,jn) - ptb(ji,jj,jk,jn) ) * vmask(ji,jj,jk) 235 233 END DO 236 234 END DO … … 259 257 DO jk = 1, jpkm1 260 258 ! !== Vertical tracer gradient at level jk and jk+1 261 zdkt3d(:,:,1) = ( pt (:,:,jk,jn) - pt(:,:,jk+1,jn) ) * tmask(:,:,jk+1)259 zdkt3d(:,:,1) = ( ptb(:,:,jk,jn) - ptb(:,:,jk+1,jn) ) * tmask(:,:,jk+1) 262 260 ! 263 261 ! ! surface boundary condition: zdkt3d(jk=0)=zdkt3d(jk=1) 264 262 IF( jk == 1 ) THEN ; zdkt3d(:,:,0) = zdkt3d(:,:,1) 265 ELSE ; zdkt3d(:,:,0) = ( pt (:,:,jk-1,jn) - pt(:,:,jk,jn) ) * tmask(:,:,jk)263 ELSE ; zdkt3d(:,:,0) = ( ptb(:,:,jk-1,jn) - ptb(:,:,jk,jn) ) * tmask(:,:,jk) 266 264 ENDIF 267 265 ! … … 275 273 ze1ur = r1_e1u(ji,jj) 276 274 zdxt = zdit(ji,jj,jk) * ze1ur 277 ze3wr = 1._wp / e3w (ji+ip,jj,jk+kp,kt2lev)275 ze3wr = 1._wp / e3w_n(ji+ip,jj,jk+kp) 278 276 zdzt = zdkt3d(ji+ip,jj,kp) * ze3wr 279 277 zslope_skew = triadi_g(ji+ip,jj,jk,1-ip,kp) 280 278 zslope_iso = triadi (ji+ip,jj,jk,1-ip,kp) 281 279 ! 282 zbu = 0.25_wp * e1e2u(ji,jj) * e3u (ji,jj,jk,ktlev)280 zbu = 0.25_wp * e1e2u(ji,jj) * e3u_n(ji,jj,jk) 283 281 ! ln_botmix_triad is .T. don't mask zah for bottom half cells !!gm ????? ahu is masked.... 284 282 zah = pahu(ji,jj,jk) … … 298 296 ze2vr = r1_e2v(ji,jj) 299 297 zdyt = zdjt(ji,jj,jk) * ze2vr 300 ze3wr = 1._wp / e3w (ji,jj+jp,jk+kp,kt2lev)298 ze3wr = 1._wp / e3w_n(ji,jj+jp,jk+kp) 301 299 zdzt = zdkt3d(ji,jj+jp,kp) * ze3wr 302 300 zslope_skew = triadj_g(ji,jj+jp,jk,1-jp,kp) 303 301 zslope_iso = triadj(ji,jj+jp,jk,1-jp,kp) 304 zbv = 0.25_wp * e1e2v(ji,jj) * e3v (ji,jj,jk,ktlev)302 zbv = 0.25_wp * e1e2v(ji,jj) * e3v_n(ji,jj,jk) 305 303 ! ln_botmix_triad is .T. don't mask zah for bottom half cells !!gm ????? ahv is masked... 306 304 zah = pahv(ji,jj,jk) … … 322 320 ze1ur = r1_e1u(ji,jj) 323 321 zdxt = zdit(ji,jj,jk) * ze1ur 324 ze3wr = 1._wp / e3w (ji+ip,jj,jk+kp,kt2lev)322 ze3wr = 1._wp / e3w_n(ji+ip,jj,jk+kp) 325 323 zdzt = zdkt3d(ji+ip,jj,kp) * ze3wr 326 324 zslope_skew = triadi_g(ji+ip,jj,jk,1-ip,kp) 327 325 zslope_iso = triadi(ji+ip,jj,jk,1-ip,kp) 328 326 ! 329 zbu = 0.25_wp * e1e2u(ji,jj) * e3u (ji,jj,jk,ktlev)327 zbu = 0.25_wp * e1e2u(ji,jj) * e3u_n(ji,jj,jk) 330 328 ! ln_botmix_triad is .F. mask zah for bottom half cells 331 329 zah = pahu(ji,jj,jk) * umask(ji,jj,jk+kp) ! pahu(ji+ip,jj,jk) ===>> ???? … … 345 343 ze2vr = r1_e2v(ji,jj) 346 344 zdyt = zdjt(ji,jj,jk) * ze2vr 347 ze3wr = 1._wp / e3w (ji,jj+jp,jk+kp,kt2lev)345 ze3wr = 1._wp / e3w_n(ji,jj+jp,jk+kp) 348 346 zdzt = zdkt3d(ji,jj+jp,kp) * ze3wr 349 347 zslope_skew = triadj_g(ji,jj+jp,jk,1-jp,kp) 350 348 zslope_iso = triadj(ji,jj+jp,jk,1-jp,kp) 351 zbv = 0.25_wp * e1e2v(ji,jj) * e3v (ji,jj,jk,ktlev)349 zbv = 0.25_wp * e1e2v(ji,jj) * e3v_n(ji,jj,jk) 352 350 ! ln_botmix_triad is .F. mask zah for bottom half cells 353 351 zah = pahv(ji,jj,jk) * vmask(ji,jj,jk+kp) ! pahv(ji,jj+jp,jk) ???? … … 364 362 DO jj = 2 , jpjm1 365 363 DO ji = fs_2, fs_jpim1 ! vector opt. 366 pt _rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn) + zsign * ( zftu(ji-1,jj,jk) - zftu(ji,jj,jk) &364 pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + zsign * ( zftu(ji-1,jj,jk) - zftu(ji,jj,jk) & 367 365 & + zftv(ji,jj-1,jk) - zftv(ji,jj,jk) ) & 368 & / ( e1e2t(ji,jj) * e3t (ji,jj,jk,ktlev) )366 & / ( e1e2t(ji,jj) * e3t_n(ji,jj,jk) ) 369 367 END DO 370 368 END DO … … 377 375 DO jj = 1, jpjm1 378 376 DO ji = fs_2, fs_jpim1 379 ztfw(ji,jj,jk) = ztfw(ji,jj,jk) - e1e2t(ji,jj) / e3w (ji,jj,jk,kt2lev) * tmask(ji,jj,jk) &377 ztfw(ji,jj,jk) = ztfw(ji,jj,jk) - e1e2t(ji,jj) / e3w_n(ji,jj,jk) * tmask(ji,jj,jk) & 380 378 & * ( ah_wslp2(ji,jj,jk) - akz(ji,jj,jk) ) & 381 & * ( pt (ji,jj,jk-1,jn) - pt(ji,jj,jk,jn) )379 & * ( ptb(ji,jj,jk-1,jn) - ptb(ji,jj,jk,jn) ) 382 380 END DO 383 381 END DO … … 389 387 DO jj = 1, jpjm1 390 388 DO ji = fs_2, fs_jpim1 391 ztfw(ji,jj,jk) = ztfw(ji,jj,jk) - e1e2t(ji,jj) / e3w (ji,jj,jk,kt2lev) * tmask(ji,jj,jk) &392 & * ah_wslp2(ji,jj,jk) * ( pt (ji,jj,jk-1,jn) - pt(ji,jj,jk,jn) )389 ztfw(ji,jj,jk) = ztfw(ji,jj,jk) - e1e2t(ji,jj) / e3w_n(ji,jj,jk) * tmask(ji,jj,jk) & 390 & * ah_wslp2(ji,jj,jk) * ( ptb(ji,jj,jk-1,jn) - ptb(ji,jj,jk,jn) ) 393 391 END DO 394 392 END DO 395 393 END DO 396 CASE( 2 ) ! 2nd pass : eddy flux = ah_wslp2 and akz applied on pt and pt_lev0gradients, resp.394 CASE( 2 ) ! 2nd pass : eddy flux = ah_wslp2 and akz applied on ptb and ptbb gradients, resp. 397 395 DO jk = 2, jpkm1 398 396 DO jj = 1, jpjm1 399 397 DO ji = fs_2, fs_jpim1 400 ztfw(ji,jj,jk) = ztfw(ji,jj,jk) - e1e2t(ji,jj) / e3w (ji,jj,jk,kt2lev) * tmask(ji,jj,jk) &401 & * ( ah_wslp2(ji,jj,jk) * ( pt (ji,jj,jk-1,jn) - pt(ji,jj,jk,jn) ) &402 & + akz (ji,jj,jk) * ( pt _lev0(ji,jj,jk-1,jn) - pt_lev0(ji,jj,jk,jn) ) )398 ztfw(ji,jj,jk) = ztfw(ji,jj,jk) - e1e2t(ji,jj) / e3w_n(ji,jj,jk) * tmask(ji,jj,jk) & 399 & * ( ah_wslp2(ji,jj,jk) * ( ptb (ji,jj,jk-1,jn) - ptb (ji,jj,jk,jn) ) & 400 & + akz (ji,jj,jk) * ( ptbb(ji,jj,jk-1,jn) - ptbb(ji,jj,jk,jn) ) ) 403 401 END DO 404 402 END DO … … 407 405 ENDIF 408 406 ! 409 DO jk = 1, jpkm1 !== Divergence of vertical fluxes added to pt _rhs==!407 DO jk = 1, jpkm1 !== Divergence of vertical fluxes added to pta ==! 410 408 DO jj = 2, jpjm1 411 409 DO ji = fs_2, fs_jpim1 ! vector opt. 412 pt _rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn) + zsign * ( ztfw(ji,jj,jk+1) - ztfw(ji,jj,jk) ) &413 & / ( e1e2t(ji,jj) * e3t (ji,jj,jk,ktlev) )410 pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + zsign * ( ztfw(ji,jj,jk+1) - ztfw(ji,jj,jk) ) & 411 & / ( e1e2t(ji,jj) * e3t_n(ji,jj,jk) ) 414 412 END DO 415 413 END DO -
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traqsr.F90
r10806 r10874 75 75 CONTAINS 76 76 77 SUBROUTINE tra_qsr( kt , ktlev, kt2lev, pts_rhs)77 SUBROUTINE tra_qsr( kt ) 78 78 !!---------------------------------------------------------------------- 79 79 !! *** ROUTINE tra_qsr *** … … 102 102 !!---------------------------------------------------------------------- 103 103 INTEGER, INTENT(in) :: kt ! ocean time-step 104 INTEGER, INTENT(in) :: ktlev ! time level index for 3-time-level source terms105 INTEGER, INTENT(in) :: kt2lev ! time level index for 2-time-level source terms106 REAL(wp), INTENT( inout), DIMENSION(jpi,jpj,jpk,jpts) :: pts_rhs ! temperature and salinity trends107 104 ! 108 105 INTEGER :: ji, jj, jk ! dummy loop indices … … 129 126 IF( l_trdtra ) THEN ! trends diagnostic: save the input temperature trend 130 127 ALLOCATE( ztrdt(jpi,jpj,jpk) ) 131 ztrdt(:,:,:) = pts_rhs(:,:,:,jp_tem)128 ztrdt(:,:,:) = tsa(:,:,:,jp_tem) 132 129 ENDIF 133 130 ! … … 175 172 zze = 568.2 * zCtot**(-0.746) 176 173 IF( zze > 102. ) zze = 200.0 * zCtot**(-0.293) 177 zpsi = gdepw (ji,jj,jk,kt2lev) / zze174 zpsi = gdepw_n(ji,jj,jk) / zze 178 175 ! 179 176 zlogc = LOG( zchl ) … … 221 218 DO jj = 2, jpjm1 222 219 DO ji = fs_2, fs_jpim1 223 zc0 = ze0(ji,jj,jk-1) * EXP( - e3t (ji,jj,jk-1,ktlev) * xsi0r )224 zc1 = ze1(ji,jj,jk-1) * EXP( - e3t (ji,jj,jk-1,ktlev) * zekb(ji,jj) )225 zc2 = ze2(ji,jj,jk-1) * EXP( - e3t (ji,jj,jk-1,ktlev) * zekg(ji,jj) )226 zc3 = ze3(ji,jj,jk-1) * EXP( - e3t (ji,jj,jk-1,ktlev) * zekr(ji,jj) )220 zc0 = ze0(ji,jj,jk-1) * EXP( - e3t_n(ji,jj,jk-1) * xsi0r ) 221 zc1 = ze1(ji,jj,jk-1) * EXP( - e3t_n(ji,jj,jk-1) * zekb(ji,jj) ) 222 zc2 = ze2(ji,jj,jk-1) * EXP( - e3t_n(ji,jj,jk-1) * zekg(ji,jj) ) 223 zc3 = ze3(ji,jj,jk-1) * EXP( - e3t_n(ji,jj,jk-1) * zekr(ji,jj) ) 227 224 ze0(ji,jj,jk) = zc0 228 225 ze1(ji,jj,jk) = zc1 … … 251 248 DO jj = 2, jpjm1 252 249 DO ji = fs_2, fs_jpim1 253 zc0 = zz0 * EXP( -gdepw (ji,jj,jk ,kt2lev)*xsi0r ) + zz1 * EXP( -gdepw(ji,jj,jk ,kt2lev)*xsi1r )254 zc1 = zz0 * EXP( -gdepw (ji,jj,jk+1,kt2lev)*xsi0r ) + zz1 * EXP( -gdepw(ji,jj,jk+1,kt2lev)*xsi1r )250 zc0 = zz0 * EXP( -gdepw_n(ji,jj,jk )*xsi0r ) + zz1 * EXP( -gdepw_n(ji,jj,jk )*xsi1r ) 251 zc1 = zz0 * EXP( -gdepw_n(ji,jj,jk+1)*xsi0r ) + zz1 * EXP( -gdepw_n(ji,jj,jk+1)*xsi1r ) 255 252 qsr_hc(ji,jj,jk) = qsr(ji,jj) * ( zc0 * wmask(ji,jj,jk) - zc1 * wmask(ji,jj,jk+1) ) 256 253 END DO … … 264 261 DO jj = 2, jpjm1 !-----------------------------! 265 262 DO ji = fs_2, fs_jpim1 ! vector opt. 266 pts_rhs(ji,jj,jk,jp_tem) = pts_rhs(ji,jj,jk,jp_tem) &267 & + z1_2 * ( qsr_hc_b(ji,jj,jk) + qsr_hc(ji,jj,jk) ) / e3t (ji,jj,jk,ktlev)263 tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem) & 264 & + z1_2 * ( qsr_hc_b(ji,jj,jk) + qsr_hc(ji,jj,jk) ) / e3t_n(ji,jj,jk) 268 265 END DO 269 266 END DO … … 298 295 ! 299 296 IF( l_trdtra ) THEN ! qsr tracers trends saved for diagnostics 300 ztrdt(:,:,:) = pts_rhs(:,:,:,jp_tem) - ztrdt(:,:,:)297 ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:) 301 298 CALL trd_tra( kt, 'TRA', jp_tem, jptra_qsr, ztrdt ) 302 299 DEALLOCATE( ztrdt ) -
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/trasbc.F90
r10806 r10874 51 51 CONTAINS 52 52 53 SUBROUTINE tra_sbc ( kt , ktlev, pts_rhs)53 SUBROUTINE tra_sbc ( kt ) 54 54 !!---------------------------------------------------------------------- 55 55 !! *** ROUTINE tra_sbc *** … … 63 63 !! (2) Fwe , tracer carried with the water that is exchanged with air+ice. 64 64 !! The input forcing fields (emp, rnf, sfx, isf) contain Fext+Fwe, 65 !! they are simply added to the tracer trend ( pts_rhs).65 !! they are simply added to the tracer trend (tsa). 66 66 !! In linear free surface case (ln_linssh=T), the volume of the 67 67 !! ocean does not change with the water exchanges at the (air+ice)-sea … … 69 69 !! concentration/dilution effect associated with water exchanges. 70 70 !! 71 !! ** Action : - Update pts_rhswith the surface boundary condition trend71 !! ** Action : - Update tsa with the surface boundary condition trend 72 72 !! - send trends to trdtra module for further diagnostics(l_trdtra=T) 73 73 !!---------------------------------------------------------------------- 74 INTEGER, INTENT(in) :: kt ! ocean time-step index 75 INTEGER, INTENT(in) :: ktlev ! time level index for source terms 76 REAL(wp), INTENT( inout), DIMENSION(jpi,jpj,jpk,jpts) :: pts_rhs ! temperature and salinity trends 74 INTEGER, INTENT(in) :: kt ! ocean time-step index 77 75 ! 78 76 INTEGER :: ji, jj, jk, jn ! dummy loop indices … … 92 90 IF( l_trdtra ) THEN !* Save ta and sa trends 93 91 ALLOCATE( ztrdt(jpi,jpj,jpk) , ztrds(jpi,jpj,jpk) ) 94 ztrdt(:,:,:) = pts_rhs(:,:,:,jp_tem)95 ztrds(:,:,:) = pts_rhs(:,:,:,jp_sal)92 ztrdt(:,:,:) = tsa(:,:,:,jp_tem) 93 ztrds(:,:,:) = tsa(:,:,:,jp_sal) 96 94 ENDIF 97 95 ! … … 133 131 DO jj = 2, jpj !==>> add concentration/dilution effect due to constant volume cell 134 132 DO ji = fs_2, fs_jpim1 ! vector opt. 135 sbc_tsc(ji,jj,jp_tem) = sbc_tsc(ji,jj,jp_tem) + r1_rau0 * emp(ji,jj) * ts (ji,jj,1,jp_tem,ktlev)136 sbc_tsc(ji,jj,jp_sal) = sbc_tsc(ji,jj,jp_sal) + r1_rau0 * emp(ji,jj) * ts (ji,jj,1,jp_sal,ktlev)133 sbc_tsc(ji,jj,jp_tem) = sbc_tsc(ji,jj,jp_tem) + r1_rau0 * emp(ji,jj) * tsn(ji,jj,1,jp_tem) 134 sbc_tsc(ji,jj,jp_sal) = sbc_tsc(ji,jj,jp_sal) + r1_rau0 * emp(ji,jj) * tsn(ji,jj,1,jp_sal) 137 135 END DO 138 136 END DO !==>> output c./d. term 139 IF( iom_use('emp_x_sst') ) CALL iom_put( "emp_x_sst", emp (:,:) * ts (:,:,1,jp_tem,ktlev) )140 IF( iom_use('emp_x_sss') ) CALL iom_put( "emp_x_sss", emp (:,:) * ts (:,:,1,jp_sal,ktlev) )137 IF( iom_use('emp_x_sst') ) CALL iom_put( "emp_x_sst", emp (:,:) * tsn(:,:,1,jp_tem) ) 138 IF( iom_use('emp_x_sss') ) CALL iom_put( "emp_x_sss", emp (:,:) * tsn(:,:,1,jp_sal) ) 141 139 ENDIF 142 140 ! … … 144 142 DO jj = 2, jpj 145 143 DO ji = fs_2, fs_jpim1 ! vector opt. 146 pts_rhs(ji,jj,1,jn) = pts_rhs(ji,jj,1,jn) + zfact * ( sbc_tsc_b(ji,jj,jn) + sbc_tsc(ji,jj,jn) ) / e3t(ji,jj,1,ktlev)144 tsa(ji,jj,1,jn) = tsa(ji,jj,1,jn) + zfact * ( sbc_tsc_b(ji,jj,jn) + sbc_tsc(ji,jj,jn) ) / e3t_n(ji,jj,1) 147 145 END DO 148 146 END DO … … 175 173 DO jk = ikt, ikb - 1 176 174 ! compute trend 177 pts_rhs(ji,jj,jk,jp_tem) = pts_rhs(ji,jj,jk,jp_tem) &175 tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem) & 178 176 & + zfact * ( risf_tsc_b(ji,jj,jp_tem) + risf_tsc(ji,jj,jp_tem) ) & 179 177 & * r1_hisf_tbl(ji,jj) … … 182 180 ! level partially include in ice shelf boundary layer 183 181 ! compute trend 184 pts_rhs(ji,jj,ikb,jp_tem) = pts_rhs(ji,jj,ikb,jp_tem) &182 tsa(ji,jj,ikb,jp_tem) = tsa(ji,jj,ikb,jp_tem) & 185 183 & + zfact * ( risf_tsc_b(ji,jj,jp_tem) + risf_tsc(ji,jj,jp_tem) ) & 186 184 & * r1_hisf_tbl(ji,jj) * ralpha(ji,jj) … … 201 199 zdep = zfact / h_rnf(ji,jj) 202 200 DO jk = 1, nk_rnf(ji,jj) 203 pts_rhs(ji,jj,jk,jp_tem) = pts_rhs(ji,jj,jk,jp_tem) &201 tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem) & 204 202 & + ( rnf_tsc_b(ji,jj,jp_tem) + rnf_tsc(ji,jj,jp_tem) ) * zdep 205 IF( ln_rnf_sal ) pts_rhs(ji,jj,jk,jp_sal) = pts_rhs(ji,jj,jk,jp_sal) &203 IF( ln_rnf_sal ) tsa(ji,jj,jk,jp_sal) = tsa(ji,jj,jk,jp_sal) & 206 204 & + ( rnf_tsc_b(ji,jj,jp_sal) + rnf_tsc(ji,jj,jp_sal) ) * zdep 207 205 END DO … … 211 209 ENDIF 212 210 213 IF( iom_use('rnf_x_sst') ) CALL iom_put( "rnf_x_sst", rnf*ts (:,:,1,jp_tem,ktlev) ) ! runoff term on sst214 IF( iom_use('rnf_x_sss') ) CALL iom_put( "rnf_x_sss", rnf*ts (:,:,1,jp_sal,ktlev) ) ! runoff term on sss211 IF( iom_use('rnf_x_sst') ) CALL iom_put( "rnf_x_sst", rnf*tsn(:,:,1,jp_tem) ) ! runoff term on sst 212 IF( iom_use('rnf_x_sss') ) CALL iom_put( "rnf_x_sss", rnf*tsn(:,:,1,jp_sal) ) ! runoff term on sss 215 213 216 214 #if defined key_asminc … … 225 223 DO jj = 2, jpj 226 224 DO ji = fs_2, fs_jpim1 227 ztim = ssh_iau(ji,jj) / e3t (ji,jj,1,ktlev)228 pts_rhs(ji,jj,1,jp_tem) = pts_rhs(ji,jj,1,jp_tem) + ts(ji,jj,1,jp_tem,ktlev) * ztim229 pts_rhs(ji,jj,1,jp_sal) = pts_rhs(ji,jj,1,jp_sal) + ts(ji,jj,1,jp_sal,ktlev) * ztim225 ztim = ssh_iau(ji,jj) / e3t_n(ji,jj,1) 226 tsa(ji,jj,1,jp_tem) = tsa(ji,jj,1,jp_tem) + tsn(ji,jj,1,jp_tem) * ztim 227 tsa(ji,jj,1,jp_sal) = tsa(ji,jj,1,jp_sal) + tsn(ji,jj,1,jp_sal) * ztim 230 228 END DO 231 229 END DO … … 234 232 DO ji = fs_2, fs_jpim1 235 233 ztim = ssh_iau(ji,jj) / ( ht_n(ji,jj) + 1. - ssmask(ji, jj) ) 236 pts_rhs(ji,jj,:,jp_tem) = pts_rhs(ji,jj,:,jp_tem) + ts(ji,jj,:,jp_tem,ktlev) * ztim237 pts_rhs(ji,jj,:,jp_sal) = pts_rhs(ji,jj,:,jp_sal) + ts(ji,jj,:,jp_sal,ktlev) * ztim234 tsa(ji,jj,:,jp_tem) = tsa(ji,jj,:,jp_tem) + tsn(ji,jj,:,jp_tem) * ztim 235 tsa(ji,jj,:,jp_sal) = tsa(ji,jj,:,jp_sal) + tsn(ji,jj,:,jp_sal) * ztim 238 236 END DO 239 237 END DO … … 252 250 DO jj = 2, jpj 253 251 DO ji = fs_2, fs_jpim1 254 zdep = 1._wp / e3t (ji,jj,jk,ktlev)255 pts_rhs(ji,jj,jk,jp_tem) = pts_rhs(ji,jj,jk,jp_tem) - htsc_iscpl(ji,jj,jk,jp_tem) * zdep256 pts_rhs(ji,jj,jk,jp_sal) = pts_rhs(ji,jj,jk,jp_sal) - htsc_iscpl(ji,jj,jk,jp_sal) * zdep252 zdep = 1._wp / e3t_n(ji,jj,jk) 253 tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem) - htsc_iscpl(ji,jj,jk,jp_tem) * zdep 254 tsa(ji,jj,jk,jp_sal) = tsa(ji,jj,jk,jp_sal) - htsc_iscpl(ji,jj,jk,jp_sal) * zdep 257 255 END DO 258 256 END DO … … 261 259 262 260 IF( l_trdtra ) THEN ! save the horizontal diffusive trends for further diagnostics 263 ztrdt(:,:,:) = pts_rhs(:,:,:,jp_tem) - ztrdt(:,:,:)264 ztrds(:,:,:) = pts_rhs(:,:,:,jp_sal) - ztrds(:,:,:)261 ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:) 262 ztrds(:,:,:) = tsa(:,:,:,jp_sal) - ztrds(:,:,:) 265 263 CALL trd_tra( kt, 'TRA', jp_tem, jptra_nsr, ztrdt ) 266 264 CALL trd_tra( kt, 'TRA', jp_sal, jptra_nsr, ztrds ) -
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/trazdf.F90
r10825 r10874 44 44 CONTAINS 45 45 46 SUBROUTINE tra_zdf( kt , ktlev1, ktlev2, ktlev3, kt2lev, pts_rhs)46 SUBROUTINE tra_zdf( kt ) 47 47 !!---------------------------------------------------------------------- 48 48 !! *** ROUTINE tra_zdf *** … … 51 51 !!--------------------------------------------------------------------- 52 52 INTEGER, INTENT(in) :: kt ! ocean time-step index 53 INTEGER, INTENT(in) :: ktlev1, ktlev2, ktlev3 ! time level indices for 3-time-level source terms54 INTEGER, INTENT(in) :: kt2lev ! time level index for 2-time-level source terms55 REAL(wp), INTENT( inout), DIMENSION(jpi,jpj,jpk,jpts) :: pts_rhs ! temperature and salinity trends56 53 ! 57 54 INTEGER :: jk ! Dummy loop indices … … 73 70 IF( l_trdtra ) THEN !* Save ta and sa trends 74 71 ALLOCATE( ztrdt(jpi,jpj,jpk) , ztrds(jpi,jpj,jpk) ) 75 ztrdt(:,:,:) = pts_rhs(:,:,:,jp_tem)76 ztrds(:,:,:) = pts_rhs(:,:,:,jp_sal)72 ztrdt(:,:,:) = tsa(:,:,:,jp_tem) 73 ztrds(:,:,:) = tsa(:,:,:,jp_sal) 77 74 ENDIF 78 75 ! 79 76 ! !* compute lateral mixing trend and add it to the general trend 80 CALL tra_zdf_imp( kt, nit000, ktlev1, ktlev2, ktlev3, kt2lev, 'TRA', r2dt, ts(:,:,:,:,ktlev1), pts_rhs, jpts )77 CALL tra_zdf_imp( kt, nit000, 'TRA', r2dt, tsb, tsa, jpts ) 81 78 82 79 !!gm WHY here ! and I don't like that ! … … 84 81 ! JMM avoid negative salinities near river outlet ! Ugly fix 85 82 ! JMM : restore negative salinities to small salinities: 86 WHERE( pts_rhs(:,:,:,jp_sal) < 0._wp ) pts_rhs(:,:,:,jp_sal) = 0.1_wp83 WHERE( tsa(:,:,:,jp_sal) < 0._wp ) tsa(:,:,:,jp_sal) = 0.1_wp 87 84 !!gm 88 85 89 86 IF( l_trdtra ) THEN ! save the vertical diffusive trends for further diagnostics 90 87 DO jk = 1, jpkm1 91 ztrdt(:,:,jk) = ( ( pts_rhs(:,:,jk,jp_tem)*e3t(:,:,jk,ktlev3) - ts(:,:,jk,jp_tem,ktlev1)*e3t(:,:,jk,ktlev1) ) &92 & / (e3t (:,:,jk,ktlev2)*r2dt) ) - ztrdt(:,:,jk)93 ztrds(:,:,jk) = ( ( pts_rhs(:,:,jk,jp_sal)*e3t(:,:,jk,ktlev3) - ts(:,:,jk,jp_sal,ktlev1)*e3t(:,:,jk,ktlev1) ) &94 & / (e3t (:,:,jk,ktlev2)*r2dt) ) - ztrds(:,:,jk)88 ztrdt(:,:,jk) = ( ( tsa(:,:,jk,jp_tem)*e3t_a(:,:,jk) - tsb(:,:,jk,jp_tem)*e3t_b(:,:,jk) ) & 89 & / (e3t_n(:,:,jk)*r2dt) ) - ztrdt(:,:,jk) 90 ztrds(:,:,jk) = ( ( tsa(:,:,jk,jp_sal)*e3t_a(:,:,jk) - tsb(:,:,jk,jp_sal)*e3t_b(:,:,jk) ) & 91 & / (e3t_n(:,:,jk)*r2dt) ) - ztrds(:,:,jk) 95 92 END DO 96 93 !!gm this should be moved in trdtra.F90 and done on all trends … … 110 107 111 108 112 SUBROUTINE tra_zdf_imp( kt, kit000, ktlev1, ktlev2, ktlev3, kt2lev, cdtype, p2dt, pt, pt_rhs, kjpt )109 SUBROUTINE tra_zdf_imp( kt, kit000, cdtype, p2dt, ptb, pta, kjpt ) 113 110 !!---------------------------------------------------------------------- 114 111 !! *** ROUTINE tra_zdf_imp *** … … 128 125 !! If iso-neutral mixing, add to avt the contribution due to lateral mixing. 129 126 !! 130 !! ** Action : - pt _rhsbecomes the after tracer127 !! ** Action : - pta becomes the after tracer 131 128 !!--------------------------------------------------------------------- 132 129 INTEGER , INTENT(in ) :: kt ! ocean time-step index 133 130 INTEGER , INTENT(in ) :: kit000 ! first time step index 134 INTEGER , INTENT(in ) :: ktlev1, ktlev2, ktlev3 ! time level indices for 3-time-level source terms135 INTEGER , INTENT(in ) :: kt2lev ! time level index for 2-time-level source terms136 131 CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator) 137 132 INTEGER , INTENT(in ) :: kjpt ! number of tracers 138 133 REAL(wp) , INTENT(in ) :: p2dt ! tracer time-step 139 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: pt ! before and now tracer fields140 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pt _rhs! in: tracer trend ; out: after tracer field134 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: ptb ! before and now tracer fields 135 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pta ! in: tracer trend ; out: after tracer field 141 136 ! 142 137 INTEGER :: ji, jj, jk, jn ! dummy loop indices … … 186 181 DO jj = 2, jpjm1 187 182 DO ji = fs_2, fs_jpim1 ! vector opt. (ensure same order of calculation as below if wi=0.) 188 zzwi = - p2dt * zwt(ji,jj,jk ) / e3w (ji,jj,jk ,kt2lev)189 zzws = - p2dt * zwt(ji,jj,jk+1) / e3w (ji,jj,jk+1,kt2lev)190 zwd(ji,jj,jk) = e3t (ji,jj,jk,ktlev3) - zzwi - zzws &183 zzwi = - p2dt * zwt(ji,jj,jk ) / e3w_n(ji,jj,jk ) 184 zzws = - p2dt * zwt(ji,jj,jk+1) / e3w_n(ji,jj,jk+1) 185 zwd(ji,jj,jk) = e3t_a(ji,jj,jk) - zzwi - zzws & 191 186 & + p2dt * ( MAX( wi(ji,jj,jk ) , 0._wp ) - MIN( wi(ji,jj,jk+1) , 0._wp ) ) 192 187 zwi(ji,jj,jk) = zzwi + p2dt * MIN( wi(ji,jj,jk ) , 0._wp ) … … 199 194 DO jj = 2, jpjm1 200 195 DO ji = fs_2, fs_jpim1 ! vector opt. 201 zwi(ji,jj,jk) = - p2dt * zwt(ji,jj,jk ) / e3w (ji,jj,jk,kt2lev)202 zws(ji,jj,jk) = - p2dt * zwt(ji,jj,jk+1) / e3w (ji,jj,jk+1,kt2lev)203 zwd(ji,jj,jk) = e3t (ji,jj,jk,ktlev3) - zwi(ji,jj,jk) - zws(ji,jj,jk)196 zwi(ji,jj,jk) = - p2dt * zwt(ji,jj,jk ) / e3w_n(ji,jj,jk) 197 zws(ji,jj,jk) = - p2dt * zwt(ji,jj,jk+1) / e3w_n(ji,jj,jk+1) 198 zwd(ji,jj,jk) = e3t_a(ji,jj,jk) - zwi(ji,jj,jk) - zws(ji,jj,jk) 204 199 END DO 205 200 END DO … … 221 216 ! Suffices i,s and d indicate "inferior" (below diagonal), diagonal 222 217 ! and "superior" (above diagonal) components of the tridiagonal system. 223 ! The solution will be in the 4d array pt _rhs.218 ! The solution will be in the 4d array pta. 224 219 ! The 3d array zwt is used as a work space array. 225 ! En route to the solution pt _rhsis used a to evaluate the rhs and then220 ! En route to the solution pta is used a to evaluate the rhs and then 226 221 ! used as a work space array: its value is modified. 227 222 ! … … 243 238 DO jj = 2, jpjm1 !* 2nd recurrence: Zk = Yk - Ik / Tk-1 Zk-1 244 239 DO ji = fs_2, fs_jpim1 245 pt _rhs(ji,jj,1,jn) = e3t(ji,jj,1,ktlev1) * pt(ji,jj,1,jn) + p2dt * e3t(ji,jj,1,ktlev2) * pt_rhs(ji,jj,1,jn)240 pta(ji,jj,1,jn) = e3t_b(ji,jj,1) * ptb(ji,jj,1,jn) + p2dt * e3t_n(ji,jj,1) * pta(ji,jj,1,jn) 246 241 END DO 247 242 END DO … … 249 244 DO jj = 2, jpjm1 250 245 DO ji = fs_2, fs_jpim1 251 zrhs = e3t (ji,jj,jk,ktlev1) * pt(ji,jj,jk,jn) + p2dt * e3t(ji,jj,jk,ktlev2) * pt_rhs(ji,jj,jk,jn) ! zrhs=right hand side252 pt _rhs(ji,jj,jk,jn) = zrhs - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) * pt_rhs(ji,jj,jk-1,jn)246 zrhs = e3t_b(ji,jj,jk) * ptb(ji,jj,jk,jn) + p2dt * e3t_n(ji,jj,jk) * pta(ji,jj,jk,jn) ! zrhs=right hand side 247 pta(ji,jj,jk,jn) = zrhs - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) * pta(ji,jj,jk-1,jn) 253 248 END DO 254 249 END DO … … 257 252 DO jj = 2, jpjm1 !* 3d recurrence: Xk = (Zk - Sk Xk+1 ) / Tk (result is the after tracer) 258 253 DO ji = fs_2, fs_jpim1 259 pt _rhs(ji,jj,jpkm1,jn) = pt_rhs(ji,jj,jpkm1,jn) / zwt(ji,jj,jpkm1) * tmask(ji,jj,jpkm1)254 pta(ji,jj,jpkm1,jn) = pta(ji,jj,jpkm1,jn) / zwt(ji,jj,jpkm1) * tmask(ji,jj,jpkm1) 260 255 END DO 261 256 END DO … … 263 258 DO jj = 2, jpjm1 264 259 DO ji = fs_2, fs_jpim1 265 pt _rhs(ji,jj,jk,jn) = ( pt_rhs(ji,jj,jk,jn) - zws(ji,jj,jk) * pt_rhs(ji,jj,jk+1,jn) ) &260 pta(ji,jj,jk,jn) = ( pta(ji,jj,jk,jn) - zws(ji,jj,jk) * pta(ji,jj,jk+1,jn) ) & 266 261 & / zwt(ji,jj,jk) * tmask(ji,jj,jk) 267 262 END DO
Note: See TracChangeset
for help on using the changeset viewer.
