Changeset 6140 for trunk/NEMOGCM/NEMO/OPA_SRC/TRA
- Timestamp:
- 2015-12-21T12:35:23+01:00 (8 years ago)
- Location:
- trunk/NEMOGCM/NEMO/OPA_SRC/TRA
- Files:
-
- 2 deleted
- 22 edited
- 1 copied
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMOGCM/NEMO/OPA_SRC/TRA/eosbn2.F90
r5541 r6140 2 2 !!============================================================================== 3 3 !! *** MODULE eosbn2 *** 4 !! Ocean diagnostic variable : equation of state - in situ and potential density 5 !! - Brunt-Vaisala frequency 4 !! Equation Of Seawater : in situ density - Brunt-Vaisala frequency 6 5 !!============================================================================== 7 6 !! History : OPA ! 1989-03 (O. Marti) Original code … … 26 25 27 26 !!---------------------------------------------------------------------- 28 !! eos 29 !! eos_insitu 30 !! eos_insitu_pot 31 !! eos_insitu_2d 32 !! bn2 33 !! eos_rab 34 !! eos_rab_3d 35 !! eos_rab_2d 36 !! eos_fzp_2d 37 !! eos_fzp_0d 38 !! eos_init 27 !! eos : generic interface of the equation of state 28 !! eos_insitu : Compute the in situ density 29 !! eos_insitu_pot: Compute the insitu and surface referenced potential volumic mass 30 !! eos_insitu_2d : Compute the in situ density for 2d fields 31 !! bn2 : Compute the Brunt-Vaisala frequency 32 !! eos_rab : generic interface of in situ thermal/haline expansion ratio 33 !! eos_rab_3d : compute in situ thermal/haline expansion ratio 34 !! eos_rab_2d : compute in situ thermal/haline expansion ratio for 2d fields 35 !! eos_fzp_2d : freezing temperature for 2d fields 36 !! eos_fzp_0d : freezing temperature for scalar 37 !! eos_init : set eos parameters (namelist) 39 38 !!---------------------------------------------------------------------- 40 USE dom_oce ! ocean space and time domain 41 USE phycst ! physical constants 39 USE dom_oce ! ocean space and time domain 40 USE phycst ! physical constants 41 USE stopar ! Stochastic T/S fluctuations 42 USE stopts ! Stochastic T/S fluctuations 42 43 ! 43 USE in_out_manager 44 USE lib_mpp 45 USE lib_fortran 46 USE prtctl 47 USE wrk_nemo 44 USE in_out_manager ! I/O manager 45 USE lib_mpp ! MPP library 46 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 47 USE prtctl ! Print control 48 USE wrk_nemo ! Memory Allocation 48 49 USE lbclnk ! ocean lateral boundary conditions 49 USE timing ! Timing 50 USE stopar ! Stochastic T/S fluctuations 51 USE stopts ! Stochastic T/S fluctuations 50 USE timing ! Timing 52 51 53 52 IMPLICIT NONE 54 53 PRIVATE 55 54 56 ! 55 ! !! * Interface 57 56 INTERFACE eos 58 57 MODULE PROCEDURE eos_insitu, eos_insitu_pot, eos_insitu_2d … … 75 74 PUBLIC eos_init ! called by istate module 76 75 77 ! !!* Namelist (nameos)*76 ! !!** Namelist nameos ** 78 77 INTEGER , PUBLIC :: nn_eos ! = 0/1/2 type of eq. of state and Brunt-Vaisala frequ. 79 78 LOGICAL , PUBLIC :: ln_useCT ! determine if eos_pt_from_ct is used to compute sst_m 80 79 81 ! !!! simplified eos coefficients 82 ! default value: Vallis 2006 80 ! !!! simplified eos coefficients (default value: Vallis 2006) 83 81 REAL(wp) :: rn_a0 = 1.6550e-1_wp ! thermal expansion coeff. 84 82 REAL(wp) :: rn_b0 = 7.6554e-1_wp ! saline expansion coeff. … … 172 170 173 171 !! * Substitutions 174 # include "domzgr_substitute.h90"175 172 # include "vectopt_loop_substitute.h90" 176 173 !!---------------------------------------------------------------------- … … 587 584 DO ji = 1, jpi 588 585 ! 589 zh = fsdept(ji,jj,jk) * r1_Z0 ! depth586 zh = gdept_n(ji,jj,jk) * r1_Z0 ! depth 590 587 zt = pts (ji,jj,jk,jp_tem) * r1_T0 ! temperature 591 588 zs = SQRT( ABS( pts(ji,jj,jk,jp_sal) + rdeltaS ) * r1_S0 ) ! square root salinity … … 645 642 zt = pts (ji,jj,jk,jp_tem) - 10._wp ! pot. temperature anomaly (t-T0) 646 643 zs = pts (ji,jj,jk,jp_sal) - 35._wp ! abs. salinity anomaly (s-S0) 647 zh = fsdept(ji,jj,jk)! depth in meters at t-point644 zh = gdept_n(ji,jj,jk) ! depth in meters at t-point 648 645 ztm = tmask(ji,jj,jk) ! land/sea bottom mask = surf. mask 649 646 ! … … 913 910 DO jj = 1, jpj ! surface and bottom value set to zero one for all in istate.F90 914 911 DO ji = 1, jpi 915 zrw = ( fsdepw(ji,jj,jk ) - fsdept(ji,jj,jk) ) &916 & / ( fsdept(ji,jj,jk-1) - fsdept(ji,jj,jk) )912 zrw = ( gdepw_n(ji,jj,jk ) - gdept_n(ji,jj,jk) ) & 913 & / ( gdept_n(ji,jj,jk-1) - gdept_n(ji,jj,jk) ) 917 914 ! 918 915 zaw = pab(ji,jj,jk,jp_tem) * (1. - zrw) + pab(ji,jj,jk-1,jp_tem) * zrw … … 921 918 pn2(ji,jj,jk) = grav * ( zaw * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) ) & 922 919 & - zbw * ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) ) ) & 923 & / fse3w(ji,jj,jk) * tmask(ji,jj,jk)920 & / e3w_n(ji,jj,jk) * tmask(ji,jj,jk) 924 921 END DO 925 922 END DO … … 1129 1126 DO ji = 1, jpi 1130 1127 ! 1131 zh = fsdept(ji,jj,jk) * r1_Z0 ! depth1128 zh = gdept_n(ji,jj,jk) * r1_Z0 ! depth 1132 1129 zt = pts (ji,jj,jk,jp_tem) * r1_T0 ! temperature 1133 1130 zs = SQRT( ABS( pts(ji,jj,jk,jp_sal) + rdeltaS ) * r1_S0 ) ! square root salinity … … 1193 1190 zt = pts(ji,jj,jk,jp_tem) - 10._wp ! temperature anomaly (t-T0) 1194 1191 zs = pts (ji,jj,jk,jp_sal) - 35._wp ! abs. salinity anomaly (s-S0) 1195 zh = fsdept(ji,jj,jk)! depth in meters at t-point1192 zh = gdept_n(ji,jj,jk) ! depth in meters at t-point 1196 1193 ztm = tmask(ji,jj,jk) ! tmask 1197 1194 zn = 0.5_wp * zh * r1_rau0 * ztm -
trunk/NEMOGCM/NEMO/OPA_SRC/TRA/traadv.F90
r5930 r6140 26 26 USE ldftra ! lateral diffusion: eddy diffusivity & EIV coeff. 27 27 USE ldfslp ! Lateral diffusion: slopes of neutral surfaces 28 USE c1d ! 1D vertical configuration29 28 ! 30 29 USE in_out_manager ! I/O manager … … 67 66 68 67 !! * Substitutions 69 # include "domzgr_substitute.h90"70 68 # include "vectopt_loop_substitute.h90" 71 69 !!---------------------------------------------------------------------- … … 96 94 ! ! set time step 97 95 IF( neuler == 0 .AND. kt == nit000 ) THEN ! at nit000 98 r2dt ra(:) = rdttra(:) ! = rdtra(restarting with Euler time stepping)96 r2dt = rdt ! = rdt (restarting with Euler time stepping) 99 97 ELSEIF( kt <= nit000 + 1) THEN ! at nit000 or nit000+1 100 r2dt ra(:) = 2._wp * rdttra(:) ! = 2 rdttra(leapfrog)98 r2dt = 2._wp * rdt ! = 2 rdt (leapfrog) 101 99 ENDIF 102 100 ! 103 101 ! !== effective transport ==! 104 102 DO jk = 1, jpkm1 105 zun(:,:,jk) = e2u (:,:) * fse3u(:,:,jk) * un(:,:,jk) ! eulerian transport only106 zvn(:,:,jk) = e1v (:,:) * fse3v(:,:,jk) * vn(:,:,jk)103 zun(:,:,jk) = e2u (:,:) * e3u_n(:,:,jk) * un(:,:,jk) ! eulerian transport only 104 zvn(:,:,jk) = e1v (:,:) * e3v_n(:,:,jk) * vn(:,:,jk) 107 105 zwn(:,:,jk) = e1e2t(:,:) * wn(:,:,jk) 108 106 END DO … … 135 133 CALL tra_adv_cen ( kt, nit000, 'TRA', zun, zvn, zwn , tsn, tsa, jpts, nn_cen_h, nn_cen_v ) 136 134 CASE ( np_FCT ) ! FCT scheme : 2nd / 4th order 137 CALL tra_adv_fct ( kt, nit000, 'TRA', r2dt ra, zun, zvn, zwn, tsb, tsn, tsa, jpts, nn_fct_h, nn_fct_v )135 CALL tra_adv_fct ( kt, nit000, 'TRA', r2dt, zun, zvn, zwn, tsb, tsn, tsa, jpts, nn_fct_h, nn_fct_v ) 138 136 CASE ( np_FCT_zts ) ! 2nd order FCT with vertical time-splitting 139 CALL tra_adv_fct_zts( kt, nit000, 'TRA', r2dt ra, zun, zvn, zwn, tsb, tsn, tsa, jpts , nn_fct_zts )137 CALL tra_adv_fct_zts( kt, nit000, 'TRA', r2dt, zun, zvn, zwn, tsb, tsn, tsa, jpts , nn_fct_zts ) 140 138 CASE ( np_MUS ) ! MUSCL 141 CALL tra_adv_mus ( kt, nit000, 'TRA', r2dt ra, zun, zvn, zwn, tsb, tsa, jpts , ln_mus_ups )139 CALL tra_adv_mus ( kt, nit000, 'TRA', r2dt, zun, zvn, zwn, tsb, tsa, jpts , ln_mus_ups ) 142 140 CASE ( np_UBS ) ! UBS 143 CALL tra_adv_ubs ( kt, nit000, 'TRA', r2dt ra, zun, zvn, zwn, tsb, tsn, tsa, jpts , nn_ubs_v )141 CALL tra_adv_ubs ( kt, nit000, 'TRA', r2dt, zun, zvn, zwn, tsb, tsn, tsa, jpts , nn_ubs_v ) 144 142 CASE ( np_QCK ) ! QUICKEST 145 CALL tra_adv_qck ( kt, nit000, 'TRA', r2dt ra, zun, zvn, zwn, tsb, tsn, tsa, jpts )143 CALL tra_adv_qck ( kt, nit000, 'TRA', r2dt, zun, zvn, zwn, tsb, tsn, tsa, jpts ) 146 144 ! 147 145 END SELECT 148 146 ! 149 ! 147 ! ! print mean trends (used for debugging) 150 148 IF(ln_ctl) CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' adv - Ta: ', mask1=tmask, & 151 149 & tab3d_2=tsa(:,:,:,jp_sal), clinfo2= ' Sa: ', mask2=tmask, clinfo3='tra' ) … … 175 173 ! 176 174 ! !== Namelist ==! 177 !178 175 REWIND( numnam_ref ) ! Namelist namtra_adv in reference namelist : Tracer advection scheme 179 176 READ ( numnam_ref, namtra_adv, IOSTAT = ios, ERR = 901) 180 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtra_adv in reference namelist', lwp )177 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtra_adv in reference namelist', lwp ) 181 178 ! 182 179 REWIND( numnam_cfg ) ! Namelist namtra_adv in configuration namelist : Tracer advection scheme 183 180 READ ( numnam_cfg, namtra_adv, IOSTAT = ios, ERR = 902 ) 184 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtra_adv in configuration namelist', lwp )185 IF(lwm) WRITE 186 181 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtra_adv in configuration namelist', lwp ) 182 IF(lwm) WRITE( numond, namtra_adv ) 183 ! 187 184 IF(lwp) THEN ! Namelist print 188 185 WRITE(numout,*) … … 203 200 WRITE(numout,*) ' QUICKEST scheme ln_traadv_qck = ', ln_traadv_qck 204 201 ENDIF 205 202 ! 206 203 ioptio = 0 !== Parameter control ==! 207 204 IF( ln_traadv_cen ) ioptio = ioptio + 1 … … 215 212 CALL ctl_warn( 'tra_adv_init: You are running without tracer advection.' ) 216 213 ENDIF 217 IF( (ioptio /= 1).AND. (.NOT. lk_c1d ) ) & 218 CALL ctl_stop( 'tra_adv_init: Choose ONE advection scheme in namelist namtra_adv' ) 214 IF( ioptio /= 1 ) CALL ctl_stop( 'tra_adv_init: Choose ONE advection scheme in namelist namtra_adv' ) 219 215 ! 220 216 IF( ln_traadv_cen .AND. ( nn_cen_h /= 2 .AND. nn_cen_h /= 4 ) & ! Centered … … 231 227 CALL ctl_stop( 'tra_adv_init: force 2nd order FCT scheme, 4th order does not exist with sub-timestepping' ) 232 228 ENDIF 233 IF( lk_vvl) THEN229 IF( .NOT.ln_linssh ) THEN 234 230 CALL ctl_stop( 'tra_adv_init: vertical sub-timestepping not allow in non-linear free surface' ) 235 231 ENDIF … … 243 239 ENDIF 244 240 IF( ln_isfcav ) THEN ! ice-shelf cavities 245 IF( ln_traadv_cen .AND. nn_cen_v /= 4 .OR. & ! NO 4th order with ISF246 & ln_traadv_fct .AND. nn_fct_v /= 4 ) CALL ctl_stop( 'tra_adv_init: 4th order COMPACT scheme not allowed with ISF' )241 IF( ln_traadv_cen .AND. nn_cen_v == 4 .OR. & ! NO 4th order with ISF 242 & ln_traadv_fct .AND. nn_fct_v == 4 ) CALL ctl_stop( 'tra_adv_init: 4th order COMPACT scheme not allowed with ISF' ) 247 243 ENDIF 248 244 ! … … 255 251 IF( ln_traadv_ubs ) nadv = np_UBS 256 252 IF( ln_traadv_qck ) nadv = np_QCK 257 253 ! 258 254 IF(lwp) THEN ! Print the choice 259 255 WRITE(numout,*) 260 IF( nadv == np_NO_adv ) WRITE(numout,*) ' NO T-S advection' 261 IF( nadv == np_CEN ) WRITE(numout,*) ' CEN scheme is used. Horizontal order: ', nn_cen_h, & 262 & ' Vertical order: ', nn_cen_v 263 IF( nadv == np_FCT ) WRITE(numout,*) ' FCT scheme is used. Horizontal order: ', nn_fct_h, & 264 & ' Vertical order: ', nn_fct_v 265 IF( nadv == np_FCT_zts ) WRITE(numout,*) ' use 2nd order FCT with ', nn_fct_zts,'vertical sub-timestepping' 266 IF( nadv == np_MUS ) WRITE(numout,*) ' MUSCL scheme is used' 267 IF( nadv == np_UBS ) WRITE(numout,*) ' UBS scheme is used' 268 IF( nadv == np_QCK ) WRITE(numout,*) ' QUICKEST scheme is used' 256 SELECT CASE ( nadv ) 257 CASE( np_NO_adv ) ; WRITE(numout,*) ' NO T-S advection' 258 CASE( np_CEN ) ; WRITE(numout,*) ' CEN scheme is used. Horizontal order: ', nn_cen_h, & 259 & ' Vertical order: ', nn_cen_v 260 CASE( np_FCT ) ; WRITE(numout,*) ' FCT scheme is used. Horizontal order: ', nn_fct_h, & 261 & ' Vertical order: ', nn_fct_v 262 CASE( np_FCT_zts ) ; WRITE(numout,*) ' use 2nd order FCT with ', nn_fct_zts,'vertical sub-timestepping' 263 CASE( np_MUS ) ; WRITE(numout,*) ' MUSCL scheme is used' 264 CASE( np_UBS ) ; WRITE(numout,*) ' UBS scheme is used' 265 CASE( np_QCK ) ; WRITE(numout,*) ' QUICKEST scheme is used' 266 END SELECT 269 267 ENDIF 270 268 ! -
trunk/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_cen.F90
r5836 r6140 2 2 !!====================================================================== 3 3 !! *** MODULE traadv_cen *** 4 !! Ocean tracers: horizontal & verticaladvective trend (2nd/4th order centered)4 !! Ocean tracers: advective trend (2nd/4th order centered) 5 5 !!====================================================================== 6 6 !! History : 3.7 ! 2014-05 (G. Madec) original code … … 8 8 9 9 !!---------------------------------------------------------------------- 10 !! tra_adv_cen : update the tracer trend with the advection trends using a centered or scheme (2nd or 4th order)11 !! NB: on the vertical it is actually a 4th order COMPACT scheme which is used12 !!---------------------------------------------------------------------- 13 USE oce , ONLY: tsn! now ocean temperature and salinity14 USE dom_oce 15 USE eosbn2 16 USE traadv_fct 17 USE trd_oce 18 USE trdtra 19 USE diaptr 10 !! tra_adv_cen : update the tracer trend with the advection trends using a centered or scheme (2nd or 4th order) 11 !! NB: on the vertical it is actually a 4th order COMPACT scheme which is used 12 !!---------------------------------------------------------------------- 13 USE oce , ONLY: tsn ! now ocean temperature and salinity 14 USE dom_oce ! ocean space and time domain 15 USE eosbn2 ! equation of state 16 USE traadv_fct ! acces to routine interp_4th_cpt 17 USE trd_oce ! trends: ocean variables 18 USE trdtra ! trends manager: tracers 19 USE diaptr ! poleward transport diagnostics 20 20 ! 21 USE in_out_manager 22 USE iom 23 USE trc_oce 24 USE lib_mpp 25 USE wrk_nemo 26 USE timing 21 USE in_out_manager ! I/O manager 22 USE iom ! IOM library 23 USE trc_oce ! share passive tracers/Ocean variables 24 USE lib_mpp ! MPP library 25 USE wrk_nemo ! Memory Allocation 26 USE timing ! Timing 27 27 28 28 IMPLICIT NONE … … 34 34 35 35 !! * Substitutions 36 # include "domzgr_substitute.h90"37 36 # include "vectopt_loop_substitute.h90" 38 37 !!---------------------------------------------------------------------- … … 53 52 !! ** Method : The advection is evaluated by a 2nd or 4th order scheme 54 53 !! using now fields (leap-frog scheme). 55 !!56 54 !! kn_cen_h = 2 ==>> 2nd order centered scheme on the horizontal 57 55 !! = 4 ==>> 4th order - - - - 58 !!59 56 !! kn_cen_v = 2 ==>> 2nd order centered scheme on the vertical 60 57 !! = 4 ==>> 4th order COMPACT scheme - - 61 58 !! 62 !! ** Action : - update pta with the now advective tracer trends 63 !! - send trends to trdtra module for further diagnostcs 59 !! ** Action : - update pta with the now advective tracer trends 60 !! - send trends to trdtra module for further diagnostcs (l_trdtra=T) 61 !! - htr_adv, str_adv : poleward advective heat and salt transport (ln_diaptr=T) 64 62 !!---------------------------------------------------------------------- 65 63 INTEGER , INTENT(in ) :: kt ! ocean time-step index … … 90 88 ENDIF 91 89 ! 92 ! ! surface & bottom values93 IF( lk_vvl ) zwz(:,:, 1 ) = 0._wp ! set to zero one for all94 zwz(:,:,jpk) = 0._wp ! except at the surface in linear free surface90 ! 91 zwz(:,:, 1 ) = 0._wp ! surface & bottom vertical flux set to zero for all tracers 92 zwz(:,:,jpk) = 0._wp 95 93 ! 96 94 DO jn = 1, kjpt !== loop over the tracers ==! 97 95 ! 98 SELECT CASE( kn_cen_h ) 99 ! 100 CASE( 2 ) !2nd order centered96 SELECT CASE( kn_cen_h ) !-- Horizontal fluxes --! 97 ! 98 CASE( 2 ) !* 2nd order centered 101 99 DO jk = 1, jpkm1 102 100 DO jj = 1, jpjm1 … … 108 106 END DO 109 107 ! 110 CASE( 4 ) !4th order centered111 ztu(:,:,jpk) = 0._wp 108 CASE( 4 ) !* 4th order centered 109 ztu(:,:,jpk) = 0._wp ! Bottom value : flux set to zero 112 110 ztv(:,:,jpk) = 0._wp 113 DO jk = 1, jpkm1 !gradient114 DO jj = 2, jpjm1 ! masked derivative111 DO jk = 1, jpkm1 ! masked gradient 112 DO jj = 2, jpjm1 115 113 DO ji = fs_2, fs_jpim1 ! vector opt. 116 114 ztu(ji,jj,jk) = ( ptn(ji+1,jj ,jk,jn) - ptn(ji,jj,jk,jn) ) * umask(ji,jj,jk) … … 121 119 CALL lbc_lnk( ztu, 'U', -1. ) ; CALL lbc_lnk( ztv, 'V', -1. ) ! Lateral boundary cond. (unchanged sgn) 122 120 ! 123 DO jk = 1, jpkm1 121 DO jk = 1, jpkm1 ! Horizontal advective fluxes 124 122 DO jj = 2, jpjm1 125 123 DO ji = 1, fs_jpim1 ! vector opt. … … 140 138 END SELECT 141 139 ! 142 ! !== Vertical fluxes ==! 143 ! 144 SELECT CASE( kn_cen_v ) !* interior fluxes 145 ! 146 CASE( 2 ) ! 2nd order centered 140 SELECT CASE( kn_cen_v ) !-- Vertical fluxes --! (interior) 141 ! 142 CASE( 2 ) !* 2nd order centered 147 143 DO jk = 2, jpk 148 144 DO jj = 2, jpjm1 … … 153 149 END DO 154 150 ! 155 CASE( 4 ) ! 4th order centered156 CALL interp_4th_cpt( ptn(:,:,:,jn) , ztw ) ! 4th order compact interpolationof T at w-point151 CASE( 4 ) !* 4th order compact 152 CALL interp_4th_cpt( ptn(:,:,:,jn) , ztw ) ! ztw = interpolated value of T at w-point 157 153 DO jk = 2, jpkm1 158 154 DO jj = 2, jpjm1 … … 165 161 END SELECT 166 162 ! 167 IF( .NOT.lk_vvl ) THEN !* top value (only in linear free surf.as zwz is multiplied by wmask)163 IF( ln_linssh ) THEN !* top value (linear free surf. only as zwz is multiplied by wmask) 168 164 IF( ln_isfcav ) THEN ! ice-shelf cavities (top of the ocean) 169 165 DO jj = 1, jpj 170 166 DO ji = 1, jpi 171 zwz(ji,jj, mikt(ji,jj) ) = pwn(ji,jj,mikt(ji,jj)) * ptn(ji,jj,mikt(ji,jj),jn) ! linear free surface167 zwz(ji,jj, mikt(ji,jj) ) = pwn(ji,jj,mikt(ji,jj)) * ptn(ji,jj,mikt(ji,jj),jn) 172 168 END DO 173 169 END DO … … 183 179 & - ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) & 184 180 & + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) & 185 & + zwz(ji,jj,jk) - zwz(ji ,jj ,jk+1) ) / ( e1e2t(ji,jj) * fse3t_n(ji,jj,jk))181 & + zwz(ji,jj,jk) - zwz(ji ,jj ,jk+1) ) * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) 186 182 END DO 187 183 END DO 188 184 END DO 189 ! 185 ! ! trend diagnostics 190 186 IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) ) THEN 191 187 CALL trd_tra( kt, cdtype, jn, jptra_xad, zwx, pun, ptn(:,:,:,jn) ) … … 193 189 CALL trd_tra( kt, cdtype, jn, jptra_zad, zwz, pwn, ptn(:,:,:,jn) ) 194 190 END IF 195 ! 191 ! ! "Poleward" heat and salt transports (contribution of upstream fluxes) 196 192 IF( cdtype == 'TRA' .AND. ln_diaptr ) THEN 197 IF( jn == jp_tem ) htr_adv(:) = ptr_sj( zwy(:,:,:) )198 IF( jn == jp_sal ) str_adv(:) = ptr_sj( zwy(:,:,:) )193 IF( jn == jp_tem ) htr_adv(:) = ptr_sj( zwy(:,:,:) ) 194 IF( jn == jp_sal ) str_adv(:) = ptr_sj( zwy(:,:,:) ) 199 195 ENDIF 200 196 ! -
trunk/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_fct.F90
r5930 r6140 39 39 40 40 !! * Substitutions 41 # include "domzgr_substitute.h90"42 41 # include "vectopt_loop_substitute.h90" 43 42 !!---------------------------------------------------------------------- … … 53 52 !! *** ROUTINE tra_adv_fct *** 54 53 !! 55 !! ** Purpose : Compute the now trend due to total advection of 56 !! tracersand add it to the general trend of tracer equations54 !! ** Purpose : Compute the now trend due to total advection of tracers 55 !! and add it to the general trend of tracer equations 57 56 !! 58 57 !! ** Method : - 2nd or 4th FCT scheme on the horizontal direction 59 58 !! (choice through the value of kn_fct) 60 !! - 4th order compact scheme on the vertical59 !! - on the vertical the 4th order is a compact scheme 61 60 !! - corrected flux (monotonic correction) 62 61 !! 63 !! ** Action : - update (pta) with the now advective tracer trends 64 !! - send the trends for further diagnostics 62 !! ** Action : - update pta with the now advective tracer trends 63 !! - send trends to trdtra module for further diagnostcs (l_trdtra=T) 64 !! - htr_adv, str_adv : poleward advective heat and salt transport (ln_diaptr=T) 65 65 !!---------------------------------------------------------------------- 66 66 INTEGER , INTENT(in ) :: kt ! ocean time-step index … … 70 70 INTEGER , INTENT(in ) :: kn_fct_h ! order of the FCT scheme (=2 or 4) 71 71 INTEGER , INTENT(in ) :: kn_fct_v ! order of the FCT scheme (=2 or 4) 72 REAL(wp) , DIMENSION( jpk ), INTENT(in ) :: p2dt ! vertical profile oftracer time-step72 REAL(wp) , INTENT(in ) :: p2dt ! tracer time-step 73 73 REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pun, pvn, pwn ! 3 ocean velocity components 74 74 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: ptb, ptn ! before and now tracer fields … … 76 76 ! 77 77 INTEGER :: ji, jj, jk, jn ! dummy loop indices 78 REAL(wp) :: z 2dtt, ztra! local scalar78 REAL(wp) :: ztra ! local scalar 79 79 REAL(wp) :: zfp_ui, zfp_vj, zfp_wk, zC2t_u, zC4t_u ! - - 80 80 REAL(wp) :: zfm_ui, zfm_vj, zfm_wk, zC2t_v, zC4t_v ! - - … … 101 101 ENDIF 102 102 ! 103 ! 104 IF( lk_vvl ) zwz(:,:, 1 ) = 0._wp ! except at the surface in linear free surface case103 ! ! surface & bottom value : flux set to zero one for all 104 zwz(:,:, 1 ) = 0._wp 105 105 zwx(:,:,jpk) = 0._wp ; zwy(:,:,jpk) = 0._wp ; zwz(:,:,jpk) = 0._wp 106 106 ! 107 107 zwi(:,:,:) = 0._wp 108 ! ! =========== 109 DO jn = 1, kjpt ! tracer loop 110 ! ! =========== 108 ! 109 DO jn = 1, kjpt !== loop over the tracers ==! 111 110 ! 112 111 ! !== upstream advection with initial mass fluxes & intermediate update ==! … … 126 125 END DO 127 126 ! !* upstream tracer flux in the k direction *! 128 DO jk = 2, jpkm1 127 DO jk = 2, jpkm1 ! Interior value ( multiplied by wmask) 129 128 DO jj = 1, jpj 130 129 DO ji = 1, jpi … … 135 134 END DO 136 135 END DO 137 ! 138 IF(.NOT.lk_vvl ) THEN ! top ocean value (only in linear free surface as zwz has been w-masked) 136 IF( ln_linssh ) THEN ! top ocean value (only in linear free surface as zwz has been w-masked) 139 137 IF( ln_isfcav ) THEN ! top of the ice-shelf cavities and at the ocean surface 140 138 DO jj = 1, jpj … … 149 147 ! 150 148 DO jk = 1, jpkm1 !* trend and after field with monotonic scheme 151 z2dtt = p2dt(jk)152 149 DO jj = 2, jpjm1 153 150 DO ji = fs_2, fs_jpim1 ! vector opt. … … 155 152 ztra = - ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) & 156 153 & + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) & 157 & + zwz(ji,jj,jk) - zwz(ji ,jj ,jk+1) ) / ( e1e2t(ji,jj) * fse3t_n(ji,jj,jk))154 & + zwz(ji,jj,jk) - zwz(ji ,jj ,jk+1) ) * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) 158 155 ! update and guess with monotonic sheme 159 156 !!gm why tmask added in the two following lines ??? the mask is done in tranxt ! 160 157 pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra * tmask(ji,jj,jk) 161 zwi(ji,jj,jk) = ( ptb(ji,jj,jk,jn) + z2dtt * ztra ) * tmask(ji,jj,jk)158 zwi(ji,jj,jk) = ( ptb(ji,jj,jk,jn) + p2dt * ztra ) * tmask(ji,jj,jk) 162 159 END DO 163 160 END DO … … 174 171 ENDIF 175 172 ! 176 !177 173 ! !== anti-diffusive flux : high order minus low order ==! 178 174 ! 179 SELECT CASE( kn_fct_h ) 180 ! 181 CASE( 2 ) !2nd order centered175 SELECT CASE( kn_fct_h ) !* horizontal anti-diffusive fluxes 176 ! 177 CASE( 2 ) !- 2nd order centered 182 178 DO jk = 1, jpkm1 183 179 DO jj = 1, jpjm1 … … 189 185 END DO 190 186 ! 191 CASE( 4 ) !4th order centered192 zltu(:,:,jpk) = 0._wp 187 CASE( 4 ) !- 4th order centered 188 zltu(:,:,jpk) = 0._wp ! Bottom value : flux set to zero 193 189 zltv(:,:,jpk) = 0._wp 194 DO jk = 1, jpkm1 195 DO jj = 1, jpjm1 ! First derivative (gradient)190 DO jk = 1, jpkm1 ! Laplacian 191 DO jj = 1, jpjm1 ! 1st derivative (gradient) 196 192 DO ji = 1, fs_jpim1 ! vector opt. 197 193 ztu(ji,jj,jk) = ( ptn(ji+1,jj ,jk,jn) - ptn(ji,jj,jk,jn) ) * umask(ji,jj,jk) … … 199 195 END DO 200 196 END DO 201 DO jj = 2, jpjm1 !197 DO jj = 2, jpjm1 ! 2nd derivative * 1/ 6 202 198 DO ji = fs_2, fs_jpim1 ! vector opt. 203 199 zltu(ji,jj,jk) = ( ztu(ji,jj,jk) + ztu(ji-1,jj,jk) ) * r1_6 … … 206 202 END DO 207 203 END DO 208 !209 204 CALL lbc_lnk( zltu, 'T', 1. ) ; CALL lbc_lnk( zltv, 'T', 1. ) ! Lateral boundary cond. (unchanged sgn) 210 205 ! 211 DO jk = 1, jpkm1 206 DO jk = 1, jpkm1 ! Horizontal advective fluxes 212 207 DO jj = 1, jpjm1 213 208 DO ji = 1, fs_jpim1 ! vector opt. … … 221 216 END DO 222 217 ! 223 CASE( 41 ) !4th order centered ==>> !!gm coding attempt need to be tested224 ztu(:,:,jpk) = 0._wp 218 CASE( 41 ) !- 4th order centered ==>> !!gm coding attempt need to be tested 219 ztu(:,:,jpk) = 0._wp ! Bottom value : flux set to zero 225 220 ztv(:,:,jpk) = 0._wp 226 DO jk = 1, jpkm1 ! gradient227 DO jj = 1, jpjm1 ! First derivative (gradient)221 DO jk = 1, jpkm1 ! 1st derivative (gradient) 222 DO jj = 1, jpjm1 228 223 DO ji = 1, fs_jpim1 ! vector opt. 229 224 ztu(ji,jj,jk) = ( ptn(ji+1,jj ,jk,jn) - ptn(ji,jj,jk,jn) ) * umask(ji,jj,jk) … … 234 229 CALL lbc_lnk( ztu, 'U', -1. ) ; CALL lbc_lnk( ztv, 'V', -1. ) ! Lateral boundary cond. (unchanged sgn) 235 230 ! 236 DO jk = 1, jpkm1 231 DO jk = 1, jpkm1 ! Horizontal advective fluxes 237 232 DO jj = 2, jpjm1 238 233 DO ji = 2, fs_jpim1 ! vector opt. … … 250 245 ! 251 246 END SELECT 252 ! !* vertical anti-diffusive fluxes253 SELECT CASE( kn_fct_v ) ! Interior values (w-masked)254 ! 255 CASE( 2 ) !2nd order centered247 ! 248 SELECT CASE( kn_fct_v ) !* vertical anti-diffusive fluxes (w-masked interior values) 249 ! 250 CASE( 2 ) !- 2nd order centered 256 251 DO jk = 2, jpkm1 257 252 DO jj = 2, jpjm1 258 253 DO ji = fs_2, fs_jpim1 259 zwz(ji,jj,jk) = ( 0.5_wp * pwn(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji,jj,jk-1,jn) ) & 260 - zwz(ji,jj,jk) ) * wmask(ji,jj,jk) 261 END DO 262 END DO 263 END DO 264 ! 265 CASE( 4 ) ! 4th order COMPACT 266 ! 267 CALL interp_4th_cpt( ptn(:,:,:,jn) , ztw ) ! COMPACT interpolation of T at w-point 268 ! 254 zwz(ji,jj,jk) = ( pwn(ji,jj,jk) * 0.5_wp * ( ptn(ji,jj,jk,jn) + ptn(ji,jj,jk-1,jn) ) & 255 & - zwz(ji,jj,jk) ) * wmask(ji,jj,jk) 256 END DO 257 END DO 258 END DO 259 ! 260 CASE( 4 ) !- 4th order COMPACT 261 CALL interp_4th_cpt( ptn(:,:,:,jn) , ztw ) ! zwt = COMPACT interpolation of T at w-point 269 262 DO jk = 2, jpkm1 270 263 DO jj = 2, jpjm1 … … 276 269 ! 277 270 END SELECT 278 ! ! top ocean value: high order = upstream ==>> zwz=0 279 zwz(:,:, 1 ) = 0._wp ! only ocean surface as interior zwz values have been w-masked 271 IF( ln_linssh ) THEN ! top ocean value: high order = upstream ==>> zwz=0 272 zwz(:,:,1) = 0._wp ! only ocean surface as interior zwz values have been w-masked 273 ENDIF 280 274 ! 281 275 CALL lbc_lnk( zwx, 'U', -1. ) ; CALL lbc_lnk( zwy, 'V', -1. ) ! Lateral bondary conditions 282 276 CALL lbc_lnk( zwz, 'W', 1. ) 283 277 ! 284 278 ! !== monotonicity algorithm ==! 285 279 ! 286 280 CALL nonosc( ptb(:,:,:,jn), zwx, zwy, zwz, zwi, p2dt ) 287 288 281 ! 289 282 ! !== final trend with corrected fluxes ==! 290 283 ! … … 295 288 & + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) & 296 289 & + zwz(ji,jj,jk) - zwz(ji ,jj ,jk+1) ) & 297 & / ( e1e2t(ji,jj) * fse3t(ji,jj,jk))298 END DO 299 END DO 300 END DO 301 ! 302 IF( l_trd ) THEN 290 & * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) 291 END DO 292 END DO 293 END DO 294 ! 295 IF( l_trd ) THEN ! trend diagnostics (contribution of upstream fluxes) 303 296 ztrdx(:,:,:) = ztrdx(:,:,:) + zwx(:,:,:) ! <<< Add to previously computed 304 297 ztrdy(:,:,:) = ztrdy(:,:,:) + zwy(:,:,:) ! <<< Add to previously computed … … 311 304 CALL wrk_dealloc( jpi,jpj,jpk, ztrdx, ztrdy, ztrdz ) 312 305 END IF 313 ! 306 ! ! "Poleward" heat and salt transports (contribution of upstream fluxes) 314 307 IF( cdtype == 'TRA' .AND. ln_diaptr ) THEN 315 IF( jn == jp_tem ) htr_adv(:) = ptr_sj( zwy(:,:,:) ) + htr_adv(:)316 IF( jn == jp_sal ) str_adv(:) = ptr_sj( zwy(:,:,:) ) + str_adv(:)308 IF( jn == jp_tem ) htr_adv(:) = htr_adv(:) + ptr_sj( zwy(:,:,:) ) 309 IF( jn == jp_sal ) str_adv(:) = str_adv(:) + ptr_sj( zwy(:,:,:) ) 317 310 ENDIF 318 311 ! 319 END DO 312 END DO ! end of tracer loop 320 313 ! 321 314 CALL wrk_dealloc( jpi,jpj,jpk, zwi, zwx, zwy, zwz, ztu, ztv, zltu, zltv, ztw ) … … 348 341 INTEGER , INTENT(in ) :: kjpt ! number of tracers 349 342 INTEGER , INTENT(in ) :: kn_fct_zts ! number of number of vertical sub-timesteps 350 REAL(wp) , DIMENSION( jpk ), INTENT(in ) :: p2dt ! vertical profile oftracer time-step343 REAL(wp) , INTENT(in ) :: p2dt ! tracer time-step 351 344 REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pun, pvn, pwn ! 3 ocean velocity components 352 345 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: ptb, ptn ! before and now tracer fields … … 354 347 ! 355 348 REAL(wp), DIMENSION( jpk ) :: zts ! length of sub-timestep for vertical advection 356 REAL(wp) , DIMENSION( jpk ):: zr_p2dt ! reciprocal of tracer timestep349 REAL(wp) :: zr_p2dt ! reciprocal of tracer timestep 357 350 INTEGER :: ji, jj, jk, jl, jn ! dummy loop indices 358 351 INTEGER :: jtb, jtn, jta ! sub timestep pointers for leap-frog/euler forward steps 359 352 INTEGER :: jtaken ! toggle for collecting appropriate fluxes from sub timesteps 360 353 REAL(wp) :: z_rzts ! Fractional length of Euler forward sub-timestep for vertical advection 361 REAL(wp) :: z 2dtt, ztra! local scalar354 REAL(wp) :: ztra ! local scalar 362 355 REAL(wp) :: zfp_ui, zfp_vj, zfp_wk ! - - 363 356 REAL(wp) :: zfm_ui, zfm_vj, zfm_wk ! - - … … 390 383 zwi(:,:,:) = 0._wp 391 384 z_rzts = 1._wp / REAL( kn_fct_zts, wp ) 392 zr_p2dt(:) = 1._wp / p2dt(:) 385 zr_p2dt = 1._wp / p2dt 386 ! 387 ! surface & Bottom value : flux set to zero for all tracers 388 zwz(:,:, 1 ) = 0._wp 389 zwx(:,:,jpk) = 0._wp ; zwz(:,:,jpk) = 0._wp 390 zwy(:,:,jpk) = 0._wp ; zwi(:,:,jpk) = 0._wp 393 391 ! 394 392 ! ! =========== 395 393 DO jn = 1, kjpt ! tracer loop 396 394 ! ! =========== 397 ! 1. Bottom value : flux set to zero 398 ! ---------------------------------- 399 zwx(:,:,jpk) = 0._wp ; zwz(:,:,jpk) = 0._wp 400 zwy(:,:,jpk) = 0._wp ; zwi(:,:,jpk) = 0._wp 401 402 ! 2. upstream advection with initial mass fluxes & intermediate update 403 ! -------------------------------------------------------------------- 404 ! upstream tracer flux in the i and j direction 405 DO jk = 1, jpkm1 395 ! 396 ! Upstream advection with initial mass fluxes & intermediate update 397 DO jk = 1, jpkm1 ! upstream tracer flux in the i and j direction 406 398 DO jj = 1, jpjm1 407 399 DO ji = 1, fs_jpim1 ! vector opt. … … 416 408 END DO 417 409 END DO 418 419 ! upstream tracer flux in the k direction 420 DO jk = 2, jpkm1 ! Interior value 410 ! ! upstream tracer flux in the k direction 411 DO jk = 2, jpkm1 ! Interior value 421 412 DO jj = 1, jpj 422 413 DO ji = 1, jpi … … 427 418 END DO 428 419 END DO 429 ! ! top value 430 IF( lk_vvl ) THEN ! variable volume: only k=1 as zwz is multiplied by wmask 431 zwz(:,:, 1 ) = 0._wp 432 ELSE ! linear free surface 433 IF( ln_isfcav ) THEN ! ice-shelf cavities 420 IF( ln_linssh ) THEN ! top value : linear free surface case only (as zwz is multiplied by wmask) 421 IF( ln_isfcav ) THEN ! ice-shelf cavities: top value 434 422 DO jj = 1, jpj 435 423 DO ji = 1, jpi … … 437 425 END DO 438 426 END DO 439 ELSE ! standard case427 ELSE ! no cavities, surface value 440 428 zwz(:,:,1) = pwn(:,:,1) * ptb(:,:,1,jn) 441 429 ENDIF … … 443 431 ! 444 432 DO jk = 1, jpkm1 ! total advective trend 445 z2dtt = p2dt(jk)446 433 DO jj = 2, jpjm1 447 434 DO ji = fs_2, fs_jpim1 ! vector opt. 448 ! total intermediate advective trends435 ! ! total intermediate advective trends 449 436 ztra = - ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) & 450 437 & + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) & 451 & + zwz(ji,jj,jk) - zwz(ji ,jj ,jk+1) ) / ( e1e2t(ji,jj) * fse3t_n(ji,jj,jk))452 ! update and guess with monotonic sheme438 & + zwz(ji,jj,jk) - zwz(ji ,jj ,jk+1) ) * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) 439 ! ! update and guess with monotonic sheme 453 440 pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra 454 zwi(ji,jj,jk) = ( ptb(ji,jj,jk,jn) + z2dtt * ztra ) * tmask(ji,jj,jk)441 zwi(ji,jj,jk) = ( ptb(ji,jj,jk,jn) + p2dt * ztra ) * tmask(ji,jj,jk) 455 442 END DO 456 443 END DO … … 497 484 END DO 498 485 END DO 499 486 ! 500 487 ! !* vertical anti-diffusive flux 501 488 zwz_sav(:,:,:) = zwz(:,:,:) 502 489 ztrs (:,:,:,1) = ptb(:,:,:,jn) 503 490 zwzts (:,:,:) = 0._wp 504 IF( lk_vvl ) zwz(:,:, 1 ) = 0._wp ! surface value set to zero in vvl case505 491 ! 506 492 DO jl = 1, kn_fct_zts ! Start of sub timestepping loop … … 508 494 IF( jl == 1 ) THEN ! Euler forward to kick things off 509 495 jtb = 1 ; jtn = 1 ; jta = 2 510 zts(:) = p2dt (:)* z_rzts496 zts(:) = p2dt * z_rzts 511 497 jtaken = MOD( kn_fct_zts + 1 , 2) ! Toggle to collect every second flux 512 498 ! ! starting at jl =1 if kn_fct_zts is odd; … … 514 500 ELSEIF( jl == 2 ) THEN ! First leapfrog step 515 501 jtb = 1 ; jtn = 2 ; jta = 3 516 zts(:) = 2._wp * p2dt (:)* z_rzts502 zts(:) = 2._wp * p2dt * z_rzts 517 503 ELSE ! Shuffle pointers for subsequent leapfrog steps 518 504 jtb = MOD(jtb,3) + 1 … … 528 514 END DO 529 515 END DO 530 IF( .NOT.lk_vvl) THEN ! top value (only in linear free surface case)516 IF( ln_linssh ) THEN ! top value (only in linear free surface case) 531 517 IF( ln_isfcav ) THEN ! ice-shelf cavities 532 518 DO jj = 1, jpj … … 535 521 END DO 536 522 END DO 537 ELSE ! standard case523 ELSE ! no ocean cavities 538 524 zwz(:,:,1) = pwn(:,:,1) * ptb(:,:,1,jn) 539 525 ENDIF … … 547 533 ztrs(ji,jj,jk,jta) = ztrs(ji,jj,jk,jtb) & 548 534 & - zts(jk) * ( zhdiv(ji,jj,jk) + zwz(ji,jj,jk) - zwz(ji,jj,jk+1) ) & 549 & / ( e1e2t(ji,jj) * fse3t(ji,jj,jk))535 & * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) 550 536 END DO 551 537 END DO … … 557 543 DO jj = 2, jpjm1 558 544 DO ji = fs_2, fs_jpim1 559 zwz(ji,jj,jk) = ( zwzts(ji,jj,jk) * zr_p2dt (jk)- zwz_sav(ji,jj,jk) ) * wmask(ji,jj,jk)545 zwz(ji,jj,jk) = ( zwzts(ji,jj,jk) * zr_p2dt - zwz_sav(ji,jj,jk) ) * wmask(ji,jj,jk) 560 546 END DO 561 547 END DO … … 576 562 pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ( zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) & 577 563 & + zwz(ji,jj,jk) - zwz(ji ,jj ,jk+1) ) & 578 & / ( e1e2t(ji,jj) * fse3t_n(ji,jj,jk))564 & * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) 579 565 END DO 580 566 END DO … … 623 609 !! in-space based differencing for fluid 624 610 !!---------------------------------------------------------------------- 625 REAL(wp) , DIMENSION(jpk) , INTENT(in ) :: p2dt ! vertical profile oftracer time-step611 REAL(wp) , INTENT(in ) :: p2dt ! tracer time-step 626 612 REAL(wp), DIMENSION (jpi,jpj,jpk), INTENT(in ) :: pbef, paft ! before & after field 627 613 REAL(wp), DIMENSION (jpi,jpj,jpk), INTENT(inout) :: paa, pbb, pcc ! monotonic fluxes in the 3 directions … … 629 615 INTEGER :: ji, jj, jk ! dummy loop indices 630 616 INTEGER :: ikm1 ! local integer 631 REAL(wp) :: zpos, zneg, zbt, za, zb, zc, zbig, zrtrn , z2dtt! local scalars617 REAL(wp) :: zpos, zneg, zbt, za, zb, zc, zbig, zrtrn ! local scalars 632 618 REAL(wp) :: zau, zbu, zcu, zav, zbv, zcv, zup, zdo ! - - 633 619 REAL(wp), POINTER, DIMENSION(:,:,:) :: zbetup, zbetdo, zbup, zbdo … … 652 638 DO jk = 1, jpkm1 653 639 ikm1 = MAX(jk-1,1) 654 z2dtt = p2dt(jk)655 640 DO jj = 2, jpjm1 656 641 DO ji = fs_2, fs_jpim1 ! vector opt. … … 679 664 680 665 ! up & down beta terms 681 zbt = e1 t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) / z2dtt666 zbt = e1e2t(ji,jj) * e3t_n(ji,jj,jk) / p2dt 682 667 zbetup(ji,jj,jk) = ( zup - paft(ji,jj,jk) ) / ( zpos + zrtrn ) * zbt 683 668 zbetdo(ji,jj,jk) = ( paft(ji,jj,jk) - zdo ) / ( zneg + zrtrn ) * zbt -
trunk/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_mle.F90
r5836 r6140 49 49 50 50 !! * Substitutions 51 # include "domzgr_substitute.h90"52 51 # include "vectopt_loop_substitute.h90" 53 52 !!---------------------------------------------------------------------- … … 125 124 DO jj = 1, jpj 126 125 DO ji = 1, jpi 127 zc = fse3t(ji,jj,jk) * REAL( MIN( MAX( 0, inml_mle(ji,jj)-jk ) , 1 ) ) ! zc being 0 outside the ML t-points126 zc = e3t_n(ji,jj,jk) * REAL( MIN( MAX( 0, inml_mle(ji,jj)-jk ) , 1 ) ) ! zc being 0 outside the ML t-points 128 127 zmld(ji,jj) = zmld(ji,jj) + zc 129 128 zbm (ji,jj) = zbm (ji,jj) + zc * (rau0 - rhop(ji,jj,jk) ) * r1_rau0 … … 157 156 END SELECT 158 157 ! ! convert density into buoyancy 159 zbm(:,:) = + grav * zbm(:,:) / MAX( fse3t(:,:,1), zmld(:,:) )158 zbm(:,:) = + grav * zbm(:,:) / MAX( e3t_n(:,:,1), zmld(:,:) ) 160 159 ! 161 160 ! … … 215 214 DO jj = 1, jpjm1 216 215 DO ji = 1, fs_jpim1 ! vector opt. 217 zcuw = 1._wp - ( fsdepw(ji+1,jj,jk) + fsdepw(ji,jj,jk) ) * zhu(ji,jj)218 zcvw = 1._wp - ( fsdepw(ji,jj+1,jk) + fsdepw(ji,jj,jk) ) * zhv(ji,jj)216 zcuw = 1._wp - ( gdepw_n(ji+1,jj,jk) + gdepw_n(ji,jj,jk) ) * zhu(ji,jj) 217 zcvw = 1._wp - ( gdepw_n(ji,jj+1,jk) + gdepw_n(ji,jj,jk) ) * zhv(ji,jj) 219 218 zcuw = zcuw * zcuw 220 219 zcvw = zcvw * zcvw -
trunk/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_mus.F90
r5930 r6140 41 41 42 42 !! * Substitutions 43 # include "domzgr_substitute.h90"44 43 # include "vectopt_loop_substitute.h90" 45 44 !!---------------------------------------------------------------------- … … 62 61 !! ld_msc_ups=T : 63 62 !! 64 !! ** Action : - update (ta,sa) with the now advective tracer trends 65 !! - send trends to trdtra module for further diagnostcs 63 !! ** Action : - update pta with the now advective tracer trends 64 !! - send trends to trdtra module for further diagnostcs (l_trdtra=T) 65 !! - htr_adv, str_adv : poleward advective heat and salt transport (ln_diaptr=T) 66 66 !! 67 67 !! References : Estubier, A., and M. Levy, Notes Techn. Pole de Modelisation … … 73 73 INTEGER , INTENT(in ) :: kjpt ! number of tracers 74 74 LOGICAL , INTENT(in ) :: ld_msc_ups ! use upstream scheme within muscl 75 REAL(wp) , DIMENSION( jpk ), INTENT(in ) :: p2dt ! vertical profile oftracer time-step75 REAL(wp) , INTENT(in ) :: p2dt ! tracer time-step 76 76 REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pun, pvn, pwn ! 3 ocean velocity components 77 77 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: ptb ! before tracer field … … 82 82 REAL(wp) :: zu, z0u, zzwx, zw ! local scalars 83 83 REAL(wp) :: zv, z0v, zzwy, z0w ! - - 84 REAL(wp) :: z dt, zalpha! - -84 REAL(wp) :: zalpha ! - - 85 85 REAL(wp), POINTER, DIMENSION(:,:,:) :: zslpx, zslpy ! 3D workspace 86 86 REAL(wp), POINTER, DIMENSION(:,:,:) :: zwx , zwy ! - - … … 116 116 ENDIF 117 117 ! 118 ! ! =========== 119 DO jn = 1, kjpt ! tracer loop 120 ! ! =========== 121 ! I. Horizontal advective fluxes 122 ! ------------------------------ 123 ! first guess of the slopes 124 zwx(:,:,jpk) = 0.e0 ; zwy(:,:,jpk) = 0.e0 ! bottom values 125 ! interior values 126 DO jk = 1, jpkm1 118 DO jn = 1, kjpt !== loop over the tracers ==! 119 ! 120 ! !* Horizontal advective fluxes 121 ! 122 ! !-- first guess of the slopes 123 zwx(:,:,jpk) = 0._wp ! bottom values 124 zwy(:,:,jpk) = 0._wp 125 DO jk = 1, jpkm1 ! interior values 127 126 DO jj = 1, jpjm1 128 127 DO ji = 1, fs_jpim1 ! vector opt. … … 132 131 END DO 133 132 END DO 134 ! 135 CALL lbc_lnk( zwx, 'U', -1. ) ! lateral boundary conditions on zwx, zwy (changed sign) 133 CALL lbc_lnk( zwx, 'U', -1. ) ! lateral boundary conditions (changed sign) 136 134 CALL lbc_lnk( zwy, 'V', -1. ) 137 ! !-- Slopes of tracer 138 zslpx(:,:,jpk) = 0.e0 ; zslpy(:,:,jpk) = 0.e0 ! bottom values 139 DO jk = 1, jpkm1 ! interior values 135 ! !-- Slopes of tracer 136 zslpx(:,:,jpk) = 0._wp ! bottom values 137 zslpy(:,:,jpk) = 0._wp 138 DO jk = 1, jpkm1 ! interior values 140 139 DO jj = 2, jpj 141 140 DO ji = fs_2, jpi ! vector opt. … … 148 147 END DO 149 148 ! 150 DO jk = 1, jpkm1 !Slopes limitation149 DO jk = 1, jpkm1 !-- Slopes limitation 151 150 DO jj = 2, jpj 152 151 DO ji = fs_2, jpi ! vector opt. … … 161 160 END DO 162 161 ! 163 ! !-- MUSCL horizontal advective fluxes 164 DO jk = 1, jpkm1 ! interior values 165 zdt = p2dt(jk) 162 DO jk = 1, jpkm1 !-- MUSCL horizontal advective fluxes 166 163 DO jj = 2, jpjm1 167 164 DO ji = fs_2, fs_jpim1 ! vector opt. … … 169 166 z0u = SIGN( 0.5, pun(ji,jj,jk) ) 170 167 zalpha = 0.5 - z0u 171 zu = z0u - 0.5 * pun(ji,jj,jk) * zdt / ( e1e2u(ji,jj) * fse3u(ji,jj,jk))168 zu = z0u - 0.5 * pun(ji,jj,jk) * p2dt * r1_e1e2u(ji,jj) / e3u_n(ji,jj,jk) 172 169 zzwx = ptb(ji+1,jj,jk,jn) + xind(ji,jj,jk) * zu * zslpx(ji+1,jj,jk) 173 170 zzwy = ptb(ji ,jj,jk,jn) + xind(ji,jj,jk) * zu * zslpx(ji ,jj,jk) … … 176 173 z0v = SIGN( 0.5, pvn(ji,jj,jk) ) 177 174 zalpha = 0.5 - z0v 178 zv = z0v - 0.5 * pvn(ji,jj,jk) * zdt / ( e1e2v(ji,jj) * fse3v(ji,jj,jk))175 zv = z0v - 0.5 * pvn(ji,jj,jk) * p2dt * r1_e1e2v(ji,jj) / e3v_n(ji,jj,jk) 179 176 zzwx = ptb(ji,jj+1,jk,jn) + xind(ji,jj,jk) * zv * zslpy(ji,jj+1,jk) 180 177 zzwy = ptb(ji,jj ,jk,jn) + xind(ji,jj,jk) * zv * zslpy(ji,jj ,jk) … … 185 182 CALL lbc_lnk( zwx, 'U', -1. ) ; CALL lbc_lnk( zwy, 'V', -1. ) ! lateral boundary conditions (changed sign) 186 183 ! 187 DO jk = 1, jpkm1 ! Tracer flux divergence at t-point added to the generaltrend184 DO jk = 1, jpkm1 !-- Tracer advective trend 188 185 DO jj = 2, jpjm1 189 186 DO ji = fs_2, fs_jpim1 ! vector opt. 190 187 pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) - ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) & 191 188 & + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) ) & 192 & / ( e1e2t(ji,jj) * fse3t(ji,jj,jk))189 & * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) 193 190 END DO 194 191 END DO 195 192 END DO 196 ! ! trend diagnostics (contribution of upstream fluxes)193 ! ! trend diagnostics 197 194 IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. & 198 195 &( cdtype == 'TRC' .AND. l_trdtrc ) ) THEN … … 200 197 CALL trd_tra( kt, cdtype, jn, jptra_yad, zwy, pvn, ptb(:,:,:,jn) ) 201 198 END IF 202 ! ! "Poleward" heat and salt transports (contribution of upstream fluxes)199 ! ! "Poleward" heat and salt transports 203 200 IF( cdtype == 'TRA' .AND. ln_diaptr ) THEN 204 201 IF( jn == jp_tem ) htr_adv(:) = ptr_sj( zwy(:,:,:) ) 205 202 IF( jn == jp_sal ) str_adv(:) = ptr_sj( zwy(:,:,:) ) 206 203 ENDIF 207 208 ! II.Vertical advective fluxes209 ! -----------------------------204 ! 205 ! !* Vertical advective fluxes 206 ! 210 207 ! !-- first guess of the slopes 211 208 zwx(:,:, 1 ) = 0._wp ! surface & bottom boundary conditions 212 zwx(:,:,jpk) = 0._wp ! surface & bottom boundary conditions213 DO jk = 2, jpkm1 209 zwx(:,:,jpk) = 0._wp 210 DO jk = 2, jpkm1 ! interior values 214 211 zwx(:,:,jk) = tmask(:,:,jk) * ( ptb(:,:,jk-1,jn) - ptb(:,:,jk,jn) ) 215 212 END DO 216 217 213 ! !-- Slopes of tracer 218 214 zslpx(:,:,1) = 0._wp ! surface values … … 220 216 DO jj = 1, jpj 221 217 DO ji = 1, jpi 222 zslpx(ji,jj,jk) = ( zwx(ji,jj,jk) + zwx(ji,jj,jk+1) ) & 223 & * ( 0.25 + SIGN( 0.25, zwx(ji,jj,jk) * zwx(ji,jj,jk+1) ) ) 224 END DO 225 END DO 226 END DO 227 ! !-- Slopes limitation 228 DO jk = 2, jpkm1 ! interior values 229 DO jj = 1, jpj 218 zslpx(ji,jj,jk) = ( zwx(ji,jj,jk) + zwx(ji,jj,jk+1) ) & 219 & * ( 0.25 + SIGN( 0.25, zwx(ji,jj,jk) * zwx(ji,jj,jk+1) ) ) 220 END DO 221 END DO 222 END DO 223 DO jk = 2, jpkm1 !-- Slopes limitation 224 DO jj = 1, jpj ! interior values 230 225 DO ji = 1, jpi 231 226 zslpx(ji,jj,jk) = SIGN( 1., zslpx(ji,jj,jk) ) * MIN( ABS( zslpx(ji,jj,jk ) ), & … … 235 230 END DO 236 231 END DO 237 ! !-- vertical advective flux 238 DO jk = 1, jpkm1 ! interior values 239 zdt = p2dt(jk) 232 DO jk = 1, jpk-2 !-- vertical advective flux 240 233 DO jj = 2, jpjm1 241 234 DO ji = fs_2, fs_jpim1 ! vector opt. 242 235 z0w = SIGN( 0.5, pwn(ji,jj,jk+1) ) 243 236 zalpha = 0.5 + z0w 244 zw = z0w - 0.5 * pwn(ji,jj,jk+1) * zdt / ( e1e2t(ji,jj) * fse3w(ji,jj,jk+1))237 zw = z0w - 0.5 * pwn(ji,jj,jk+1) * p2dt * r1_e1e2t(ji,jj) / e3w_n(ji,jj,jk+1) 245 238 zzwx = ptb(ji,jj,jk+1,jn) + xind(ji,jj,jk) * zw * zslpx(ji,jj,jk+1) 246 239 zzwy = ptb(ji,jj,jk ,jn) + xind(ji,jj,jk) * zw * zslpx(ji,jj,jk ) … … 249 242 END DO 250 243 END DO 251 ! ! top values (bottom already set to zero) 252 IF( lk_vvl ) THEN ! variable volume 253 zwx(:,:, 1 ) = 0._wp ! k=1 only as zwx has been multiplied by wmask 254 ELSE ! linear free surface 255 IF( ln_isfcav ) THEN ! ice-shelf cavities (top of the ocean) 244 IF( ln_linssh ) THEN ! top values, linear free surface only 245 IF( ln_isfcav ) THEN ! ice-shelf cavities (top of the ocean) 256 246 DO jj = 1, jpj 257 247 DO ji = 1, jpi … … 259 249 END DO 260 250 END DO 261 ELSE 251 ELSE ! no cavities: only at the ocean surface 262 252 zwx(:,:,1) = pwn(:,:,1) * ptb(:,:,1,jn) 263 253 ENDIF 264 254 ENDIF 265 255 ! 266 DO jk = 1, jpkm1 ! Compute & add thevertical advective trend256 DO jk = 1, jpkm1 !-- vertical advective trend 267 257 DO jj = 2, jpjm1 268 258 DO ji = fs_2, fs_jpim1 ! vector opt. 269 pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) - ( zwx(ji,jj,jk) - zwx(ji,jj,jk+1) ) / ( e1e2t(ji,jj) * fse3t(ji,jj,jk))270 END DO 271 END DO 272 END DO 273 ! ! Save the vertical advectivetrends for diagnostic259 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) 260 END DO 261 END DO 262 END DO 263 ! ! send trends for diagnostic 274 264 IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. & 275 265 &( cdtype == 'TRC' .AND. l_trdtrc ) ) & 276 266 CALL trd_tra( kt, cdtype, jn, jptra_zad, zwx, pwn, ptb(:,:,:,jn) ) 277 267 ! 278 END DO 268 END DO ! end of tracer loop 279 269 ! 280 270 CALL wrk_dealloc( jpi,jpj,jpk, zslpx, zslpy, zwx, zwy ) -
trunk/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_qck.F90
r5930 r6140 38 38 39 39 !! * Substitutions 40 # include "domzgr_substitute.h90"41 40 # include "vectopt_loop_substitute.h90" 42 41 !!---------------------------------------------------------------------- … … 78 77 !! prevent the appearance of spurious numerical oscillations 79 78 !! 80 !! ** Action : - update (pta) with the now advective tracer trends 81 !! - save the trends 79 !! ** Action : - update pta with the now advective tracer trends 80 !! - send trends to trdtra module for further diagnostcs (l_trdtra=T) 81 !! - htr_adv, str_adv : poleward advective heat and salt transport (ln_diaptr=T) 82 82 !! 83 83 !! ** Reference : Leonard (1979, 1991) … … 87 87 CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator) 88 88 INTEGER , INTENT(in ) :: kjpt ! number of tracers 89 REAL(wp) , DIMENSION( jpk ), INTENT(in ) :: p2dt ! vertical profile oftracer time-step89 REAL(wp) , INTENT(in ) :: p2dt ! tracer time-step 90 90 REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pun, pvn, pwn ! 3 ocean velocity components 91 91 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: ptb, ptn ! before and now tracer fields … … 105 105 IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) ) l_trd = .TRUE. 106 106 ! 107 ! I. Thehorizontal fluxes are computed with the QUICKEST + ULTIMATE scheme107 ! ! horizontal fluxes are computed with the QUICKEST + ULTIMATE scheme 108 108 CALL tra_adv_qck_i( kt, cdtype, p2dt, pun, ptb, ptn, pta, kjpt ) 109 109 CALL tra_adv_qck_j( kt, cdtype, p2dt, pvn, ptb, ptn, pta, kjpt ) 110 110 111 ! II. Thevertical fluxes are computed with the 2nd order centered scheme111 ! ! vertical fluxes are computed with the 2nd order centered scheme 112 112 CALL tra_adv_cen2_k( kt, cdtype, pwn, ptn, pta, kjpt ) 113 113 ! … … 125 125 CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator) 126 126 INTEGER , INTENT(in ) :: kjpt ! number of tracers 127 REAL(wp) , DIMENSION( jpk ), INTENT(in ) :: p2dt ! vertical profile oftracer time-step127 REAL(wp) , INTENT(in ) :: p2dt ! tracer time-step 128 128 REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pun ! i-velocity components 129 129 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: ptb, ptn ! before and now tracer fields … … 131 131 !! 132 132 INTEGER :: ji, jj, jk, jn ! dummy loop indices 133 REAL(wp) :: ztra, zbtr, zdir, zdx, z dt, zmsk ! local scalars133 REAL(wp) :: ztra, zbtr, zdir, zdx, zmsk ! local scalars 134 134 REAL(wp), POINTER, DIMENSION(:,:,:) :: zwx, zfu, zfc, zfd 135 135 !---------------------------------------------------------------------- … … 166 166 ! 167 167 DO jk = 1, jpkm1 168 zdt = p2dt(jk)169 168 DO jj = 2, jpjm1 170 169 DO ji = fs_2, fs_jpim1 ! vector opt. 171 170 zdir = 0.5 + SIGN( 0.5, pun(ji,jj,jk) ) ! if pun > 0 : zdir = 1 otherwise zdir = 0 172 zdx = ( zdir * e1t(ji,jj) + ( 1. - zdir ) * e1t(ji+1,jj) ) * e2u(ji,jj) * fse3u(ji,jj,jk)173 zwx(ji,jj,jk) = ABS( pun(ji,jj,jk) ) * zdt / zdx ! (0<zc_cfl<1 : Courant number on x-direction)171 zdx = ( zdir * e1t(ji,jj) + ( 1. - zdir ) * e1t(ji+1,jj) ) * e2u(ji,jj) * e3u_n(ji,jj,jk) 172 zwx(ji,jj,jk) = ABS( pun(ji,jj,jk) ) * p2dt / zdx ! (0<zc_cfl<1 : Courant number on x-direction) 174 173 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 175 174 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 … … 216 215 DO jj = 2, jpjm1 217 216 DO ji = fs_2, fs_jpim1 ! vector opt. 218 zbtr = 1. / ( e1e2t(ji,jj) * fse3t(ji,jj,jk))217 zbtr = r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) 219 218 ! horizontal advective trends 220 219 ztra = - zbtr * ( zwx(ji,jj,jk) - zwx(ji-1,jj,jk) ) … … 224 223 END DO 225 224 END DO 226 ! ! trend diagnostics (contribution of upstream fluxes)225 ! ! trend diagnostics 227 226 IF( l_trd ) CALL trd_tra( kt, cdtype, jn, jptra_xad, zwx, pun, ptn(:,:,:,jn) ) 228 227 ! … … 242 241 CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator) 243 242 INTEGER , INTENT(in ) :: kjpt ! number of tracers 244 REAL(wp) , DIMENSION( jpk ), INTENT(in ) :: p2dt ! vertical profile oftracer time-step243 REAL(wp) , INTENT(in ) :: p2dt ! tracer time-step 245 244 REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pvn ! j-velocity components 246 245 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: ptb, ptn ! before and now tracer fields … … 248 247 !! 249 248 INTEGER :: ji, jj, jk, jn ! dummy loop indices 250 REAL(wp) :: ztra, zbtr, zdir, zdx, z dt, zmsk ! local scalars249 REAL(wp) :: ztra, zbtr, zdir, zdx, zmsk ! local scalars 251 250 REAL(wp), POINTER, DIMENSION(:,:,:) :: zwy, zfu, zfc, zfd 252 251 !---------------------------------------------------------------------- … … 289 288 ! 290 289 DO jk = 1, jpkm1 291 zdt = p2dt(jk)292 290 DO jj = 2, jpjm1 293 291 DO ji = fs_2, fs_jpim1 ! vector opt. 294 292 zdir = 0.5 + SIGN( 0.5, pvn(ji,jj,jk) ) ! if pun > 0 : zdir = 1 otherwise zdir = 0 295 zdx = ( zdir * e2t(ji,jj) + ( 1. - zdir ) * e2t(ji,jj+1) ) * e1v(ji,jj) * fse3v(ji,jj,jk)296 zwy(ji,jj,jk) = ABS( pvn(ji,jj,jk) ) * zdt / zdx ! (0<zc_cfl<1 : Courant number on x-direction)293 zdx = ( zdir * e2t(ji,jj) + ( 1. - zdir ) * e2t(ji,jj+1) ) * e1v(ji,jj) * e3v_n(ji,jj,jk) 294 zwy(ji,jj,jk) = ABS( pvn(ji,jj,jk) ) * p2dt / zdx ! (0<zc_cfl<1 : Courant number on x-direction) 297 295 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 298 296 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 … … 340 338 DO jj = 2, jpjm1 341 339 DO ji = fs_2, fs_jpim1 ! vector opt. 342 zbtr = 1. / ( e1e2t(ji,jj) * fse3t(ji,jj,jk))340 zbtr = r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) 343 341 ! horizontal advective trends 344 342 ztra = - zbtr * ( zwy(ji,jj,jk) - zwy(ji,jj-1,jk) ) … … 348 346 END DO 349 347 END DO 350 ! ! trend diagnostics (contribution of upstream fluxes)348 ! ! trend diagnostics 351 349 IF( l_trd ) CALL trd_tra( kt, cdtype, jn, jptra_yad, zwy, pvn, ptn(:,:,:,jn) ) 352 350 ! ! "Poleward" heat and salt transports (contribution of upstream fluxes) … … 381 379 CALL wrk_alloc( jpi,jpj,jpk, zwz ) 382 380 ! 383 ! ! surface & bottom values 384 IF( lk_vvl ) zwz(:,:, 1 ) = 0._wp ! set to zero one for all 385 zwz(:,:,jpk) = 0._wp ! except at the surface in linear free surface 381 zwz(:,:, 1 ) = 0._wp ! surface & bottom values set to zero for all tracers 382 zwz(:,:,jpk) = 0._wp 386 383 ! 387 384 ! ! =========== … … 396 393 END DO 397 394 END DO 398 IF( .NOT.lk_vvl ) THEN!* top value (only in linear free surf. as zwz is multiplied by wmask)395 IF( ln_linssh ) THEN !* top value (only in linear free surf. as zwz is multiplied by wmask) 399 396 IF( ln_isfcav ) THEN ! ice-shelf cavities (top of the ocean) 400 397 DO jj = 1, jpj … … 403 400 END DO 404 401 END DO 405 ELSE ! no ice-shelfcavities (only ocean surface)402 ELSE ! no ocean cavities (only ocean surface) 406 403 zwz(:,:,1) = pwn(:,:,1) * ptn(:,:,1,jn) 407 404 ENDIF … … 412 409 DO ji = fs_2, fs_jpim1 ! vector opt. 413 410 pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) - ( zwz(ji,jj,jk) - zwz(ji,jj,jk+1) ) & 414 & / ( e1e2t(ji,jj) * fse3t(ji,jj,jk))415 END DO 416 END DO 417 END DO 418 ! ! S ave the vertical advectivetrends for diagnostic411 & * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) 412 END DO 413 END DO 414 END DO 415 ! ! Send trends for diagnostic 419 416 IF( l_trd ) CALL trd_tra( kt, cdtype, jn, jptra_zad, zwz, pwn, ptn(:,:,:,jn) ) 420 417 ! 421 418 END DO 422 419 ! 423 CALL wrk_dealloc( jpi, jpj, jpk,zwz )420 CALL wrk_dealloc( jpi,jpj,jpk, zwz ) 424 421 ! 425 422 END SUBROUTINE tra_adv_cen2_k -
trunk/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_ubs.F90
r5930 r6140 35 35 36 36 !! * Substitutions 37 # include "domzgr_substitute.h90"38 37 # include "vectopt_loop_substitute.h90" 39 38 !!---------------------------------------------------------------------- … … 71 70 !! On the vertical, the advection is evaluated using a FCT scheme, 72 71 !! as the UBS have been found to be too diffusive. 73 !!gm !! kn_ubs_v argument (not coded for the moment) 74 !! controles whether the FCT is based on a 2nd order centrered scheme (kn_ubs_v=2)75 !! or on a 4th order compactscheme (kn_ubs_v=4).72 !! kn_ubs_v argument controles whether the FCT is based on 73 !! a 2nd order centrered scheme (kn_ubs_v=2) or on a 4th order compact 74 !! scheme (kn_ubs_v=4). 76 75 !! 77 !! ** Action : - update (pta) with the now advective tracer trends 76 !! ** Action : - update pta with the now advective tracer trends 77 !! - send trends to trdtra module for further diagnostcs (l_trdtra=T) 78 !! - htr_adv, str_adv : poleward advective heat and salt transport (ln_diaptr=T) 78 79 !! 79 80 !! Reference : Shchepetkin, A. F., J. C. McWilliams, 2005, Ocean Modelling, 9, 347-404. … … 85 86 INTEGER , INTENT(in ) :: kjpt ! number of tracers 86 87 INTEGER , INTENT(in ) :: kn_ubs_v ! number of tracers 87 REAL(wp) , DIMENSION( jpk ), INTENT(in ) :: p2dt ! vertical profile oftracer time-step88 REAL(wp) , INTENT(in ) :: p2dt ! tracer time-step 88 89 REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pun, pvn, pwn ! 3 ocean transport components 89 90 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: ptb, ptn ! before and now tracer fields … … 91 92 ! 92 93 INTEGER :: ji, jj, jk, jn ! dummy loop indices 93 REAL(wp) :: ztra, zbtr, zcoef , z2dtt! local scalars94 REAL(wp) :: ztra, zbtr, zcoef ! local scalars 94 95 REAL(wp) :: zfp_ui, zfm_ui, zcenut, ztak, zfp_wk, zfm_wk ! - - 95 96 REAL(wp) :: zfp_vj, zfm_vj, zcenvt, zeeu, zeev, z_hdivn ! - - … … 110 111 IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) ) l_trd = .TRUE. 111 112 ! 112 zltu(:,:,jpk) = 0._wp ; zltv(:,:,jpk) = 0._wp ! Bottom value : set to zero one for all 113 ztw (:,:, 1 ) = 0._wp ! surface & bottom value : set to zero for all tracers 114 zltu(:,:,jpk) = 0._wp ; zltv(:,:,jpk) = 0._wp 113 115 ztw (:,:,jpk) = 0._wp ; zti (:,:,jpk) = 0._wp 114 IF( lk_vvl ) ztw(:,:, 1 ) = 0._wp ! surface value: set to zero only in vvl case115 116 ! 116 117 ! ! =========== … … 121 122 DO jj = 1, jpjm1 ! First derivative (masked gradient) 122 123 DO ji = 1, fs_jpim1 ! vector opt. 123 zeeu = e2_e1u(ji,jj) * fse3u(ji,jj,jk) * umask(ji,jj,jk)124 zeev = e1_e2v(ji,jj) * fse3v(ji,jj,jk) * vmask(ji,jj,jk)124 zeeu = e2_e1u(ji,jj) * e3u_n(ji,jj,jk) * umask(ji,jj,jk) 125 zeev = e1_e2v(ji,jj) * e3v_n(ji,jj,jk) * vmask(ji,jj,jk) 125 126 ztu(ji,jj,jk) = zeeu * ( ptb(ji+1,jj ,jk,jn) - ptb(ji,jj,jk,jn) ) 126 127 ztv(ji,jj,jk) = zeev * ( ptb(ji ,jj+1,jk,jn) - ptb(ji,jj,jk,jn) ) … … 129 130 DO jj = 2, jpjm1 ! Second derivative (divergence) 130 131 DO ji = fs_2, fs_jpim1 ! vector opt. 131 zcoef = 1._wp / ( 6._wp * fse3t(ji,jj,jk) )132 zcoef = 1._wp / ( 6._wp * e3t_n(ji,jj,jk) ) 132 133 zltu(ji,jj,jk) = ( ztu(ji,jj,jk) - ztu(ji-1,jj,jk) ) * zcoef 133 134 zltv(ji,jj,jk) = ( ztv(ji,jj,jk) - ztv(ji,jj-1,jk) ) * zcoef … … 162 163 pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) & 163 164 & - ( ztu(ji,jj,jk) - ztu(ji-1,jj ,jk) & 164 & + ztv(ji,jj,jk) - ztv(ji ,jj-1,jk) ) / ( e1e2t(ji,jj) * fse3t(ji,jj,jk))165 & + ztv(ji,jj,jk) - ztv(ji ,jj-1,jk) ) * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) 165 166 END DO 166 167 END DO … … 199 200 END DO 200 201 END DO 201 IF( .NOT.lk_vvl ) THEN! top ocean value (only in linear free surface as ztw has been w-masked)202 IF( ln_linssh ) THEN ! top ocean value (only in linear free surface as ztw has been w-masked) 202 203 IF( ln_isfcav ) THEN ! top of the ice-shelf cavities and at the ocean surface 203 204 DO jj = 1, jpj … … 212 213 ! 213 214 DO jk = 1, jpkm1 !* trend and after field with monotonic scheme 214 z2dtt = p2dt(jk)215 215 DO jj = 2, jpjm1 216 216 DO ji = fs_2, fs_jpim1 ! vector opt. 217 ztak = - ( ztw(ji,jj,jk) - ztw(ji,jj,jk+1) ) / ( e1e2t(ji,jj) * fse3t_n(ji,jj,jk))217 ztak = - ( ztw(ji,jj,jk) - ztw(ji,jj,jk+1) ) * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) 218 218 pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztak 219 zti(ji,jj,jk) = ( ptb(ji,jj,jk,jn) + z2dtt * ( ztak + zltu(ji,jj,jk) ) ) * tmask(ji,jj,jk)219 zti(ji,jj,jk) = ( ptb(ji,jj,jk,jn) + p2dt * ( ztak + zltu(ji,jj,jk) ) ) * tmask(ji,jj,jk) 220 220 END DO 221 221 END DO … … 233 233 END DO 234 234 ! ! top ocean value: high order == upstream ==>> zwz=0 235 IF( .NOT.lk_vvl ) ztw(:,:, 1 ) = 0._wp! only ocean surface as interior zwz values have been w-masked235 IF( ln_linssh ) ztw(:,:, 1 ) = 0._wp ! only ocean surface as interior zwz values have been w-masked 236 236 ! 237 237 CALL nonosc_z( ptb(:,:,:,jn), ztw, zti, p2dt ) ! monotonicity algorithm … … 246 246 END DO 247 247 END DO 248 IF( .NOT.lk_vvl) ztw(:,:, 1 ) = pwn(:,:,1) * ptn(:,:,1,jn) !!gm ISF & 4th COMPACT doesn't work248 IF( ln_linssh ) ztw(:,:, 1 ) = pwn(:,:,1) * ptn(:,:,1,jn) !!gm ISF & 4th COMPACT doesn't work 249 249 ! 250 250 END SELECT … … 253 253 DO jj = 2, jpjm1 254 254 DO ji = fs_2, fs_jpim1 ! vector opt. 255 pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) - ( ztw(ji,jj,jk) - ztw(ji,jj,jk+1) ) / ( e1e2t(ji,jj) * fse3t_n(ji,jj,jk))255 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) 256 256 END DO 257 257 END DO … … 264 264 zltv(ji,jj,jk) = pta(ji,jj,jk,jn) - zltv(ji,jj,jk) & 265 265 & + ptn(ji,jj,jk,jn) * ( pwn(ji,jj,jk) - pwn(ji,jj,jk+1) ) & 266 & / ( e1e2t(ji,jj) * fse3t_n(ji,jj,jk))266 & * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) 267 267 END DO 268 268 END DO … … 293 293 !! in-space based differencing for fluid 294 294 !!---------------------------------------------------------------------- 295 REAL(wp), INTENT(in ) , DIMENSION(jpk) :: p2dt ! vertical profile oftracer time-step295 REAL(wp), INTENT(in ) :: p2dt ! tracer time-step 296 296 REAL(wp), DIMENSION (jpi,jpj,jpk) :: pbef ! before field 297 297 REAL(wp), INTENT(inout), DIMENSION (jpi,jpj,jpk) :: paft ! after field … … 300 300 INTEGER :: ji, jj, jk ! dummy loop indices 301 301 INTEGER :: ikm1 ! local integer 302 REAL(wp) :: zpos, zneg, zbt, za, zb, zc, zbig, zrtrn , z2dtt! local scalars302 REAL(wp) :: zpos, zneg, zbt, za, zb, zc, zbig, zrtrn ! local scalars 303 303 REAL(wp), POINTER, DIMENSION(:,:,:) :: zbetup, zbetdo 304 304 !!---------------------------------------------------------------------- … … 349 349 ! --------------------------------------------------- 350 350 DO jk = 1, jpkm1 351 z2dtt = p2dt(jk)352 351 DO jj = 2, jpjm1 353 352 DO ji = fs_2, fs_jpim1 ! vector opt. … … 356 355 zneg = MAX( 0., pcc(ji ,jj ,jk ) ) - MIN( 0., pcc(ji ,jj ,jk+1) ) 357 356 ! up & down beta terms 358 zbt = e1e2t(ji,jj) * fse3t(ji,jj,jk) / z2dtt357 zbt = e1e2t(ji,jj) * e3t_n(ji,jj,jk) / p2dt 359 358 zbetup(ji,jj,jk) = ( zbetup(ji,jj,jk) - paft(ji,jj,jk) ) / (zpos+zrtrn) * zbt 360 359 zbetdo(ji,jj,jk) = ( paft(ji,jj,jk) - zbetdo(ji,jj,jk) ) / (zneg+zrtrn) * zbt -
trunk/NEMOGCM/NEMO/OPA_SRC/TRA/trabbc.F90
r5397 r6140 12 12 13 13 !!---------------------------------------------------------------------- 14 !! tra_bbc : update the tracer trend at ocean bottom15 !! tra_bbc_init : initialization of geothermal heat flux trend14 !! tra_bbc : update the tracer trend at ocean bottom 15 !! tra_bbc_init : initialization of geothermal heat flux trend 16 16 !!---------------------------------------------------------------------- 17 USE oce ! ocean variables 18 USE dom_oce ! domain: ocean 19 USE phycst ! physical constants 20 USE trd_oce ! trends: ocean variables 21 USE trdtra ! trends manager: tracers 22 USE in_out_manager ! I/O manager 23 USE iom ! I/O manager 24 USE fldread ! read input fields 25 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 26 USE lib_mpp ! distributed memory computing library 27 USE prtctl ! Print control 28 USE wrk_nemo ! Memory Allocation 29 USE timing ! Timing 17 USE oce ! ocean variables 18 USE dom_oce ! domain: ocean 19 USE phycst ! physical constants 20 USE trd_oce ! trends: ocean variables 21 USE trdtra ! trends manager: tracers 22 ! 23 USE in_out_manager ! I/O manager 24 USE iom ! xIOS 25 USE fldread ! read input fields 26 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 27 USE lib_mpp ! distributed memory computing library 28 USE prtctl ! Print control 29 USE wrk_nemo ! Memory Allocation 30 USE timing ! Timing 30 31 31 32 IMPLICIT NONE … … 40 41 REAL(wp) :: rn_geoflx_cst ! Constant value of geothermal heat flux 41 42 42 REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: qgh_trd0 ! geothermal heating trend 43 TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_qgh ! structure of input qgh (file informations, fields read) 43 REAL(wp), PUBLIC , ALLOCATABLE, DIMENSION(:,:) :: qgh_trd0 ! geothermal heating trend 44 45 TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_qgh ! structure of input qgh (file informations, fields read) 44 46 45 !! * Substitutions46 # include "domzgr_substitute.h90"47 47 !!---------------------------------------------------------------------- 48 48 !! NEMO/OPA 3.3 , NEMO Consortium (2010) … … 68 68 !! Where Qsf is the geothermal heat flux. 69 69 !! 70 !! ** Action : - update the temperature trends (ta) with the trend of71 !! the ocean bottom boundary condition70 !! ** Action : - update the temperature trends with geothermal heating trend 71 !! - send the trend for further diagnostics (ln_trdtra=T) 72 72 !! 73 73 !! References : Stein, C. A., and S. Stein, 1992, Nature, 359, 123-129. … … 75 75 !!---------------------------------------------------------------------- 76 76 INTEGER, INTENT(in) :: kt ! ocean time-step index 77 !! 78 INTEGER :: ji, jj, ik ! dummy loop indices 79 REAL(wp) :: zqgh_trd ! geothermal heat flux trend 77 ! 78 INTEGER :: ji, jj ! dummy loop indices 80 79 REAL(wp), POINTER, DIMENSION(:,:,:) :: ztrdt 81 80 !!---------------------------------------------------------------------- … … 83 82 IF( nn_timing == 1 ) CALL timing_start('tra_bbc') 84 83 ! 85 IF( l_trdtra ) THEN ! Save t a and sa trends86 CALL wrk_alloc( jpi, jpj, jpk,ztrdt )84 IF( l_trdtra ) THEN ! Save the input temperature trend 85 CALL wrk_alloc( jpi,jpj,jpk, ztrdt ) 87 86 ztrdt(:,:,:) = tsa(:,:,:,jp_tem) 88 87 ENDIF 89 ! 90 ! ! Add the geothermal heat flux trend on temperature 88 ! ! Add the geothermal trend on temperature 91 89 DO jj = 2, jpjm1 92 90 DO ji = 2, jpim1 93 ik = mbkt(ji,jj) 94 zqgh_trd = qgh_trd0(ji,jj) / fse3t(ji,jj,ik) 95 tsa(ji,jj,ik,jp_tem) = tsa(ji,jj,ik,jp_tem) + zqgh_trd 91 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)) 96 92 END DO 97 93 END DO … … 99 95 CALL lbc_lnk( tsa(:,:,:,jp_tem) , 'T', 1. ) 100 96 ! 101 IF( l_trdtra ) THEN ! S ave the geothermal heat fluxtrend for diagnostics97 IF( l_trdtra ) THEN ! Send the trend for diagnostics 102 98 ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:) 103 99 CALL trd_tra( kt, 'TRA', jp_tem, jptra_bbc, ztrdt ) 104 CALL wrk_dealloc( jpi, jpj, jpk,ztrdt )100 CALL wrk_dealloc( jpi,jpj,jpk, ztrdt ) 105 101 ENDIF 106 102 ! … … 127 123 !! ** Action : - read/fix the geothermal heat qgh_trd0 128 124 !!---------------------------------------------------------------------- 129 USE iom130 !!131 125 INTEGER :: ji, jj ! dummy loop indices 132 126 INTEGER :: inum ! temporary logical unit … … 139 133 NAMELIST/nambbc/ln_trabbc, nn_geoflx, rn_geoflx_cst, sn_qgh, cn_dir 140 134 !!---------------------------------------------------------------------- 141 135 ! 142 136 REWIND( numnam_ref ) ! Namelist nambbc in reference namelist : Bottom momentum boundary condition 143 137 READ ( numnam_ref, nambbc, IOSTAT = ios, ERR = 901) 144 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nambbc in reference namelist', lwp )145 138 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nambbc in reference namelist', lwp ) 139 ! 146 140 REWIND( numnam_cfg ) ! Namelist nambbc in configuration namelist : Bottom momentum boundary condition 147 141 READ ( numnam_cfg, nambbc, IOSTAT = ios, ERR = 902 ) 148 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nambbc in configuration namelist', lwp )142 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nambbc in configuration namelist', lwp ) 149 143 IF(lwm) WRITE ( numond, nambbc ) 150 144 ! 151 145 IF(lwp) THEN ! Control print 152 146 WRITE(numout,*) … … 159 153 WRITE(numout,*) 160 154 ENDIF 161 155 ! 162 156 IF( ln_trabbc ) THEN !== geothermal heating ==! 163 157 ! … … 190 184 WRITE(ctmp1,*) ' bad flag value for nn_geoflx = ', nn_geoflx 191 185 CALL ctl_stop( ctmp1 ) 192 !193 186 END SELECT 194 187 ! -
trunk/NEMOGCM/NEMO/OPA_SRC/TRA/trabbl.F90
r5836 r6140 70 70 71 71 !! * Substitutions 72 # include "domzgr_substitute.h90"73 72 # include "vectopt_loop_substitute.h90" 74 73 !!---------------------------------------------------------------------- … … 112 111 IF( nn_timing == 1 ) CALL timing_start( 'tra_bbl') 113 112 ! 114 IF( l_trdtra ) THEN !* Save t a and satrends113 IF( l_trdtra ) THEN !* Save the input trends 115 114 CALL wrk_alloc( jpi, jpj, jpk, ztrdt, ztrds ) 116 115 ztrdt(:,:,:) = tsa(:,:,:,jp_tem) … … 132 131 ! 133 132 END IF 134 133 ! 135 134 IF( nn_bbl_adv /= 0 ) THEN !* Advective bbl 136 135 ! … … 146 145 END IF 147 146 148 IF( l_trdtra ) THEN ! s ave the horizontal diffusive trends for further diagnostics147 IF( l_trdtra ) THEN ! send the trends for further diagnostics 149 148 ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:) 150 149 ztrds(:,:,:) = tsa(:,:,:,jp_sal) - ztrds(:,:,:) … … 211 210 & + ahv_bbl(ji ,jj ) * ( zptb(ji ,jj+1) - zptb(ji ,jj ) ) & 212 211 & - ahv_bbl(ji ,jj-1) * ( zptb(ji ,jj ) - zptb(ji ,jj-1) ) ) & 213 & / ( e1e2t(ji,jj) * fse3t(ji,jj,ik))212 & * r1_e1e2t(ji,jj) / e3t_n(ji,jj,ik) 214 213 END DO 215 214 END DO … … 263 262 ! 264 263 ! ! up -slope T-point (shelf bottom point) 265 zbtr = r1_e1e2t(iis,jj) / fse3t(iis,jj,ikus)264 zbtr = r1_e1e2t(iis,jj) / e3t_n(iis,jj,ikus) 266 265 ztra = zu_bbl * ( ptb(iid,jj,ikus,jn) - ptb(iis,jj,ikus,jn) ) * zbtr 267 266 pta(iis,jj,ikus,jn) = pta(iis,jj,ikus,jn) + ztra 268 267 ! 269 268 DO jk = ikus, ikud-1 ! down-slope upper to down T-point (deep column) 270 zbtr = r1_e1e2t(iid,jj) / fse3t(iid,jj,jk)269 zbtr = r1_e1e2t(iid,jj) / e3t_n(iid,jj,jk) 271 270 ztra = zu_bbl * ( ptb(iid,jj,jk+1,jn) - ptb(iid,jj,jk,jn) ) * zbtr 272 271 pta(iid,jj,jk,jn) = pta(iid,jj,jk,jn) + ztra 273 272 END DO 274 273 ! 275 zbtr = r1_e1e2t(iid,jj) / fse3t(iid,jj,ikud)274 zbtr = r1_e1e2t(iid,jj) / e3t_n(iid,jj,ikud) 276 275 ztra = zu_bbl * ( ptb(iis,jj,ikus,jn) - ptb(iid,jj,ikud,jn) ) * zbtr 277 276 pta(iid,jj,ikud,jn) = pta(iid,jj,ikud,jn) + ztra … … 285 284 ! 286 285 ! up -slope T-point (shelf bottom point) 287 zbtr = r1_e1e2t(ji,ijs) / fse3t(ji,ijs,ikvs)286 zbtr = r1_e1e2t(ji,ijs) / e3t_n(ji,ijs,ikvs) 288 287 ztra = zv_bbl * ( ptb(ji,ijd,ikvs,jn) - ptb(ji,ijs,ikvs,jn) ) * zbtr 289 288 pta(ji,ijs,ikvs,jn) = pta(ji,ijs,ikvs,jn) + ztra 290 289 ! 291 290 DO jk = ikvs, ikvd-1 ! down-slope upper to down T-point (deep column) 292 zbtr = r1_e1e2t(ji,ijd) / fse3t(ji,ijd,jk)291 zbtr = r1_e1e2t(ji,ijd) / e3t_n(ji,ijd,jk) 293 292 ztra = zv_bbl * ( ptb(ji,ijd,jk+1,jn) - ptb(ji,ijd,jk,jn) ) * zbtr 294 293 pta(ji,ijd,jk,jn) = pta(ji,ijd,jk,jn) + ztra 295 294 END DO 296 295 ! ! down-slope T-point (deep bottom point) 297 zbtr = r1_e1e2t(ji,ijd) / fse3t(ji,ijd,ikvd)296 zbtr = r1_e1e2t(ji,ijd) / e3t_n(ji,ijd,ikvd) 298 297 ztra = zv_bbl * ( ptb(ji,ijs,ikvs,jn) - ptb(ji,ijd,ikvd,jn) ) * zbtr 299 298 pta(ji,ijd,ikvd,jn) = pta(ji,ijd,ikvd,jn) + ztra … … 302 301 ! 303 302 END DO 304 ! ! =========== 305 END DO ! end tracer 306 ! ! =========== 307 ! 303 ! ! =========== 304 END DO ! end tracer 305 ! ! =========== 308 306 IF( nn_timing == 1 ) CALL timing_stop( 'tra_bbl_adv') 309 307 ! … … 340 338 INTEGER , INTENT(in ) :: kit000 ! first time step index 341 339 CHARACTER(len=3), INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator) 342 ! !340 ! 343 341 INTEGER :: ji, jj ! dummy loop indices 344 342 INTEGER :: ik ! local integers … … 365 363 zts (ji,jj,jp_sal) = tsb(ji,jj,ik,jp_sal) 366 364 ! 367 zdep(ji,jj) = fsdept(ji,jj,ik)! bottom T-level reference depth365 zdep(ji,jj) = gdept_n(ji,jj,ik) ! bottom T-level reference depth 368 366 zub (ji,jj) = un(ji,jj,mbku(ji,jj)) ! bottom velocity 369 367 zvb (ji,jj) = vn(ji,jj,mbkv(ji,jj)) … … 401 399 ! 402 400 ENDIF 403 401 ! 404 402 ! !-------------------! 405 403 IF( nn_bbl_adv /= 0 ) THEN ! advective bbl ! … … 500 498 INTEGER :: ios ! - - 501 499 REAL(wp), POINTER, DIMENSION(:,:) :: zmbk 502 ! !500 ! 503 501 NAMELIST/nambbl/ nn_bbl_ldf, nn_bbl_adv, rn_ahtbbl, rn_gambbl 504 502 !!---------------------------------------------------------------------- … … 506 504 IF( nn_timing == 1 ) CALL timing_start( 'tra_bbl_init') 507 505 ! 508 CALL wrk_alloc( jpi, jpj, zmbk )509 !510 511 506 REWIND( numnam_ref ) ! Namelist nambbl in reference namelist : Bottom boundary layer scheme 512 507 READ ( numnam_ref, nambbl, IOSTAT = ios, ERR = 901) 513 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nambbl in reference namelist', lwp )514 508 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nambbl in reference namelist', lwp ) 509 ! 515 510 REWIND( numnam_cfg ) ! Namelist nambbl in configuration namelist : Bottom boundary layer scheme 516 511 READ ( numnam_cfg, nambbl, IOSTAT = ios, ERR = 902 ) 517 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nambbl in configuration namelist', lwp )512 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nambbl in configuration namelist', lwp ) 518 513 IF(lwm) WRITE ( numond, nambbl ) 519 514 ! … … 545 540 END DO 546 541 ! converte into REAL to use lbc_lnk ; impose a min value of 1 as a zero can be set in lbclnk 542 CALL wrk_alloc( jpi, jpj, zmbk ) 547 543 zmbk(:,:) = REAL( mbku_d(:,:), wp ) ; CALL lbc_lnk(zmbk,'U',1.) ; mbku_d(:,:) = MAX( INT( zmbk(:,:) ), 1 ) 548 544 zmbk(:,:) = REAL( mbkv_d(:,:), wp ) ; CALL lbc_lnk(zmbk,'V',1.) ; mbkv_d(:,:) = MAX( INT( zmbk(:,:) ), 1 ) 545 CALL wrk_dealloc( jpi, jpj, zmbk ) 549 546 550 547 !* sign of grad(H) at u- and v-points … … 593 590 ENDIF 594 591 ! 595 CALL wrk_dealloc( jpi, jpj, zmbk )596 !597 592 IF( nn_timing == 1 ) CALL timing_stop( 'tra_bbl_init') 598 593 ! -
trunk/NEMOGCM/NEMO/OPA_SRC/TRA/tradmp.F90
r5836 r6140 31 31 USE dtatsd ! data: temperature & salinity 32 32 USE zdfmxl ! vertical physics: mixed layer depth 33 ! 33 34 USE in_out_manager ! I/O manager 34 35 USE lib_mpp ! MPP library … … 41 42 PRIVATE 42 43 43 PUBLIC tra_dmp ! routinecalled by step.F9044 PUBLIC tra_dmp_init ! routinecalled by nemogcm.F9044 PUBLIC tra_dmp ! called by step.F90 45 PUBLIC tra_dmp_init ! called by nemogcm.F90 45 46 46 47 ! !!* Namelist namtra_dmp : T & S newtonian damping * … … 52 53 53 54 !! * Substitutions 54 # include "domzgr_substitute.h90"55 55 # include "vectopt_loop_substitute.h90" 56 56 !!---------------------------------------------------------------------- … … 89 89 !! below the well mixed layer (nlmdmp=2) 90 90 !! 91 !! ** Action : - (ta,sa)tracer trends updated with the damping trend91 !! ** Action : - tsa: tracer trends updated with the damping trend 92 92 !!---------------------------------------------------------------------- 93 93 INTEGER, INTENT(in) :: kt ! ocean time-step index … … 100 100 ! 101 101 CALL wrk_alloc( jpi,jpj,jpk,jpts, zts_dta ) 102 !103 102 IF( l_trdtra ) THEN !* Save ta and sa trends 104 103 CALL wrk_alloc( jpi,jpj,jpk,jpts, ztrdts ) … … 139 138 DO jj = 2, jpjm1 140 139 DO ji = fs_2, fs_jpim1 ! vector opt. 141 IF( fsdept(ji,jj,jk) >= hmlp (ji,jj) ) THEN140 IF( gdept_n(ji,jj,jk) >= hmlp (ji,jj) ) THEN 142 141 tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem) & 143 142 & + resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jp_tem) - tsb(ji,jj,jk,jp_tem) ) … … 177 176 !!---------------------------------------------------------------------- 178 177 INTEGER :: ios, imask ! local integers 179 ! !178 ! 180 179 NAMELIST/namtra_dmp/ ln_tradmp, nn_zdmp, cn_resto 181 180 !!---------------------------------------------------------------------- … … 229 228 END SUBROUTINE tra_dmp_init 230 229 230 !!====================================================================== 231 231 END MODULE tradmp -
trunk/NEMOGCM/NEMO/OPA_SRC/TRA/traldf.F90
r5836 r6140 12 12 13 13 !!---------------------------------------------------------------------- 14 !! tra_ldf : update the tracer trend with the lateral diffusion 15 !! tra_ldf_init : initialization, namelist read, and parameters control 16 !!---------------------------------------------------------------------- 17 USE oce ! ocean dynamics and tracers 18 USE dom_oce ! ocean space and time domain 19 USE phycst ! physical constants 20 USE ldftra ! lateral diffusion: eddy diffusivity & EIV coeff. 21 USE ldfslp ! lateral diffusion: iso-neutral slope 22 USE traldf_lap ! lateral diffusion: laplacian iso-level operator (tra_ldf_lap routine) 23 USE traldf_iso ! lateral diffusion: laplacian iso-neutral standard operator (tra_ldf_iso routine) 24 USE traldf_triad ! lateral diffusion: laplacian iso-neutral triad operator (tra_ldf_triad routine) 25 USE traldf_blp ! lateral diffusion (iso-level lap/blp) (tra_ldf_lap routine) 26 USE trd_oce ! trends: ocean variables 27 USE trdtra ! ocean active tracers trends 14 !! tra_ldf : update the tracer trend with the lateral diffusion trend 15 !! tra_ldf_init : initialization, namelist read, and parameters control 16 !!---------------------------------------------------------------------- 17 USE oce ! ocean dynamics and tracers 18 USE dom_oce ! ocean space and time domain 19 USE phycst ! physical constants 20 USE ldftra ! lateral diffusion: eddy diffusivity & EIV coeff. 21 USE ldfslp ! lateral diffusion: iso-neutral slope 22 USE traldf_lap_blp ! lateral diffusion: laplacian iso-level operator (tra_ldf_lap/_blp routines) 23 USE traldf_iso ! lateral diffusion: laplacian iso-neutral standard operator (tra_ldf_iso routine ) 24 USE traldf_triad ! lateral diffusion: laplacian iso-neutral triad operator (tra_ldf_triad routine ) 25 USE trd_oce ! trends: ocean variables 26 USE trdtra ! ocean active tracers trends 28 27 ! 29 28 USE prtctl ! Print control … … 43 42 44 43 !! * Substitutions 45 # include "domzgr_substitute.h90"46 44 # include "vectopt_loop_substitute.h90" 47 45 !!---------------------------------------------------------------------- … … 72 70 ! 73 71 SELECT CASE ( nldf ) !* compute lateral mixing trend and add it to the general trend 74 !75 72 CASE ( np_lap ) ! laplacian: iso-level operator 76 73 CALL tra_ldf_lap ( kt, nit000,'TRA', ahtu, ahtv, gtsu, gtsv, gtui, gtvi, tsb, tsa, jpts, 1 ) … … 82 79 CALL tra_ldf_blp ( kt, nit000,'TRA', ahtu, ahtv, gtsu, gtsv, gtui, gtvi, tsb , tsa, jpts, nldf ) 83 80 END SELECT 84 81 ! 85 82 IF( l_trdtra ) THEN !* save the horizontal diffusive trends for further diagnostics 86 83 ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:) … … 114 111 WRITE(numout,*) 'tra_ldf_init : lateral tracer diffusive operator' 115 112 WRITE(numout,*) '~~~~~~~~~~~' 116 WRITE(numout,*) ' Namelist namtra_ldf already read in ldftra module'117 WRITE(numout,*) ' see ldf_tra_init report for lateral mixing parameters'113 WRITE(numout,*) ' Namelist namtra_ldf: already read in ldftra module' 114 WRITE(numout,*) ' see ldf_tra_init report for lateral mixing parameters' 118 115 WRITE(numout,*) 119 116 ENDIF … … 178 175 ENDIF 179 176 ! 180 IF( ierr == 1 ) CALL ctl_stop( ' iso-level in z-coordinate -partial step, not allowed' )177 IF( ierr == 1 ) CALL ctl_stop( 'iso-level in z-partial step, not allowed' ) 181 178 IF( ln_ldfeiv .AND. .NOT.( ln_traldf_iso .OR. ln_traldf_triad ) ) & 182 & CALL ctl_stop( ' eddy induced velocity on tracers requires isopycnal', &183 & ' laplacian diffusion' )179 & CALL ctl_stop( 'eddy induced velocity on tracers requires iso-neutral laplacian diffusion' ) 180 ! 184 181 IF( nldf == np_lap_i .OR. nldf == np_lap_it .OR. & 185 182 & nldf == np_blp_i .OR. nldf == np_blp_it ) l_ldfslp = .TRUE. ! slope of neutral surfaces required … … 187 184 IF(lwp) THEN 188 185 WRITE(numout,*) 189 IF( nldf == np_no_ldf ) WRITE(numout,*) ' NO lateral diffusion' 190 IF( nldf == np_lap ) WRITE(numout,*) ' laplacian iso-level operator' 191 IF( nldf == np_lap_i ) WRITE(numout,*) ' Rotated laplacian operator (standard)' 192 IF( nldf == np_lap_it ) WRITE(numout,*) ' Rotated laplacian operator (triad)' 193 IF( nldf == np_blp ) WRITE(numout,*) ' bilaplacian iso-level operator' 194 IF( nldf == np_blp_i ) WRITE(numout,*) ' Rotated bilaplacian operator (standard)' 195 IF( nldf == np_blp_it ) WRITE(numout,*) ' Rotated bilaplacian operator (triad)' 186 SELECT CASE( nldf ) 187 CASE( np_no_ldf ) ; WRITE(numout,*) ' NO lateral diffusion' 188 CASE( np_lap ) ; WRITE(numout,*) ' laplacian iso-level operator' 189 CASE( np_lap_i ) ; WRITE(numout,*) ' Rotated laplacian operator (standard)' 190 CASE( np_lap_it ) ; WRITE(numout,*) ' Rotated laplacian operator (triad)' 191 CASE( np_blp ) ; WRITE(numout,*) ' bilaplacian iso-level operator' 192 CASE( np_blp_i ) ; WRITE(numout,*) ' Rotated bilaplacian operator (standard)' 193 CASE( np_blp_it ) ; WRITE(numout,*) ' Rotated bilaplacian operator (triad)' 194 END SELECT 196 195 ENDIF 197 196 ! -
trunk/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_iso.F90
r5836 r6140 14 14 15 15 !!---------------------------------------------------------------------- 16 !! tra_ldf_iso : update the tracer trend with the horizontal component of a iso-neutral laplacian operator17 !! and with the vertical part of the isopycnal or geopotential s-coord. operator16 !! tra_ldf_iso : update the tracer trend with the horizontal component of a iso-neutral laplacian operator 17 !! and with the vertical part of the isopycnal or geopotential s-coord. operator 18 18 !!---------------------------------------------------------------------- 19 USE oce 20 USE dom_oce 21 USE trc_oce 22 USE zdf_oce 23 USE ldftra 24 USE ldfslp 25 USE diaptr 19 USE oce ! ocean dynamics and active tracers 20 USE dom_oce ! ocean space and time domain 21 USE trc_oce ! share passive tracers/Ocean variables 22 USE zdf_oce ! ocean vertical physics 23 USE ldftra ! lateral diffusion: tracer eddy coefficients 24 USE ldfslp ! iso-neutral slopes 25 USE diaptr ! poleward transport diagnostics 26 26 ! 27 USE in_out_manager 28 USE iom 29 USE phycst 30 USE lbclnk 31 USE wrk_nemo 32 USE timing 27 USE in_out_manager ! I/O manager 28 USE iom ! I/O library 29 USE phycst ! physical constants 30 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 31 USE wrk_nemo ! Memory Allocation 32 USE timing ! Timing 33 33 34 34 IMPLICIT NONE … … 38 38 39 39 !! * Substitutions 40 # include "domzgr_substitute.h90"41 40 # include "vectopt_loop_substitute.h90" 42 41 !!---------------------------------------------------------------------- … … 103 102 ! 104 103 INTEGER :: ji, jj, jk, jn ! dummy loop indices 104 INTEGER :: ikt 105 105 INTEGER :: ierr ! local integer 106 106 REAL(wp) :: zmsku, zahu_w, zabe1, zcof1, zcoef3 ! local scalars … … 127 127 ah_wslp2(:,:,:) = 0._wp 128 128 ENDIF 129 !130 129 ! ! set time step size (Euler/Leapfrog) 131 IF( neuler == 0 .AND. kt == nit000 ) THEN ; z2dt = rdt tra(1)! at nit000 (Euler)132 ELSE ; z2dt = 2.* rdt tra(1)! (Leapfrog)130 IF( neuler == 0 .AND. kt == nit000 ) THEN ; z2dt = rdt ! at nit000 (Euler) 131 ELSE ; z2dt = 2.* rdt ! (Leapfrog) 133 132 ENDIF 134 133 z1_2dt = 1._wp / z2dt … … 138 137 ENDIF 139 138 140 141 139 !!---------------------------------------------------------------------- 142 140 !! 0 - calculate ah_wslp2 and akz … … 149 147 DO ji = fs_2, fs_jpim1 ! vector opt. 150 148 ! 151 zmsku = tmask(ji,jj,jk) / MAX( umask(ji ,jj,jk-1) + umask(ji-1,jj,jk) &149 zmsku = wmask(ji,jj,jk) / MAX( umask(ji ,jj,jk-1) + umask(ji-1,jj,jk) & 152 150 & + umask(ji-1,jj,jk-1) + umask(ji ,jj,jk) , 1._wp ) 153 zmskv = tmask(ji,jj,jk) / MAX( vmask(ji,jj ,jk-1) + vmask(ji,jj-1,jk) &151 zmskv = wmask(ji,jj,jk) / MAX( vmask(ji,jj ,jk-1) + vmask(ji,jj-1,jk) & 154 152 & + vmask(ji,jj-1,jk-1) + vmask(ji,jj ,jk) , 1._wp ) 155 153 ! … … 183 181 DO ji = 1, fs_jpim1 184 182 akz(ji,jj,jk) = 16._wp * ah_wslp2(ji,jj,jk) & 185 & * ( akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ( fse3w(ji,jj,jk) * fse3w(ji,jj,jk) ) )183 & * ( akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ( e3w_n(ji,jj,jk) * e3w_n(ji,jj,jk) ) ) 186 184 END DO 187 185 END DO … … 191 189 DO jj = 1, jpjm1 192 190 DO ji = 1, fs_jpim1 193 ze3w_2 = fse3w(ji,jj,jk) * fse3w(ji,jj,jk)191 ze3w_2 = e3w_n(ji,jj,jk) * e3w_n(ji,jj,jk) 194 192 zcoef0 = z2dt * ( akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ze3w_2 ) 195 193 akz(ji,jj,jk) = MAX( zcoef0 - 0.5_wp , 0._wp ) * ze3w_2 * z1_2dt … … 228 226 DO jj = 1, jpjm1 ! bottom correction (partial bottom cell) 229 227 DO ji = 1, fs_jpim1 ! vector opt. 230 !!gm the following anonymous remark is to considered: ! IF useless if zpshde defines pgu everywhere231 228 zdit(ji,jj,mbku(ji,jj)) = pgu(ji,jj,jn) 232 229 zdjt(ji,jj,mbkv(ji,jj)) = pgv(ji,jj,jn) … … 242 239 ENDIF 243 240 ENDIF 244 241 ! 245 242 !!---------------------------------------------------------------------- 246 243 !! II - horizontal trend (full) … … 255 252 ELSE ; zdkt(:,:) = ( ptb(:,:,jk-1,jn) - ptb(:,:,jk,jn) ) * wmask(:,:,jk) 256 253 ENDIF 257 !!gm I don't understand why we should need this.... since wmask is used instead of tmask258 ! IF ( ln_isfcav ) THEN259 ! DO jj = 1, jpj260 ! DO ji = 1, jpi ! vector opt.261 ! ikt = mikt(ji,jj) ! surface level262 ! zdk1t(ji,jj,ikt) = ( ptb(ji,jj,ikt,jn ) - ptb(ji,jj,ikt+1,jn) ) * wmask(ji,jj,ikt+1)263 ! zdkt (ji,jj,ikt) = zdk1t(ji,jj,ikt)264 ! END DO265 ! END DO266 ! END IF267 !!gm268 269 254 DO jj = 1 , jpjm1 !== Horizontal fluxes 270 255 DO ji = 1, fs_jpim1 ! vector opt. 271 zabe1 = pahu(ji,jj,jk) * e2_e1u(ji,jj) * fse3u_n(ji,jj,jk)272 zabe2 = pahv(ji,jj,jk) * e1_e2v(ji,jj) * fse3v_n(ji,jj,jk)273 ! 274 zmsku = 1. / MAX( tmask(ji+1,jj,jk ) + tmask(ji,jj,jk+1) &275 & + tmask(ji+1,jj,jk+1) + tmask(ji,jj,jk ), 1. )276 ! 277 zmskv = 1. / MAX( tmask(ji,jj+1,jk ) + tmask(ji,jj,jk+1) &278 & + tmask(ji,jj+1,jk+1) + tmask(ji,jj,jk ), 1. )256 zabe1 = pahu(ji,jj,jk) * e2_e1u(ji,jj) * e3u_n(ji,jj,jk) 257 zabe2 = pahv(ji,jj,jk) * e1_e2v(ji,jj) * e3v_n(ji,jj,jk) 258 ! 259 zmsku = 1. / MAX( wmask(ji+1,jj,jk ) + wmask(ji,jj,jk+1) & 260 & + wmask(ji+1,jj,jk+1) + wmask(ji,jj,jk ), 1. ) 261 ! 262 zmskv = 1. / MAX( wmask(ji,jj+1,jk ) + wmask(ji,jj,jk+1) & 263 & + wmask(ji,jj+1,jk+1) + wmask(ji,jj,jk ), 1. ) 279 264 ! 280 265 zcof1 = - pahu(ji,jj,jk) * e2u(ji,jj) * uslp(ji,jj,jk) * zmsku … … 294 279 pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + zsign * ( zftu(ji,jj,jk) - zftu(ji-1,jj,jk) & 295 280 & + zftv(ji,jj,jk) - zftv(ji,jj-1,jk) ) & 296 & / ( e1e2t(ji,jj) * fse3t(ji,jj,jk))281 & * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) 297 282 END DO 298 283 END DO 299 284 END DO ! End of slab 300 285 301 302 286 !!---------------------------------------------------------------------- 303 287 !! III - vertical trend (full) 304 288 !!---------------------------------------------------------------------- 305 289 ! 306 290 ztfw(1,:,:) = 0._wp ; ztfw(jpi,:,:) = 0._wp 307 291 ! 308 292 ! Vertical fluxes 309 293 ! --------------- 310 311 ! Surface and bottom vertical fluxes set to zero 294 ! ! Surface and bottom vertical fluxes set to zero 312 295 ztfw(:,:, 1 ) = 0._wp ; ztfw(:,:,jpk) = 0._wp 313 296 314 ! interior (2=<jk=<jpk-1) 315 DO jk = 2, jpkm1 297 DO jk = 2, jpkm1 ! interior (2=<jk=<jpk-1) 316 298 DO jj = 2, jpjm1 317 299 DO ji = fs_2, fs_jpim1 ! vector opt. 318 300 ! 319 zmsku = tmask(ji,jj,jk) / MAX( umask(ji ,jj,jk-1) + umask(ji-1,jj,jk) &301 zmsku = wmask(ji,jj,jk) / MAX( umask(ji ,jj,jk-1) + umask(ji-1,jj,jk) & 320 302 & + umask(ji-1,jj,jk-1) + umask(ji ,jj,jk) , 1._wp ) 321 zmskv = tmask(ji,jj,jk) / MAX( vmask(ji,jj ,jk-1) + vmask(ji,jj-1,jk) &303 zmskv = wmask(ji,jj,jk) / MAX( vmask(ji,jj ,jk-1) + vmask(ji,jj-1,jk) & 322 304 & + vmask(ji,jj-1,jk-1) + vmask(ji,jj ,jk) , 1._wp ) 323 305 ! … … 337 319 END DO 338 320 END DO 339 !340 321 ! !== add the vertical 33 flux ==! 341 322 IF( ln_traldf_lap ) THEN ! laplacian case: eddy coef = ah_wslp2 - akz … … 343 324 DO jj = 1, jpjm1 344 325 DO ji = fs_2, fs_jpim1 345 ztfw(ji,jj,jk) = ztfw(ji,jj,jk) + e1e2t(ji,jj) / fse3w(ji,jj,jk) * tmask(ji,jj,jk) &326 ztfw(ji,jj,jk) = ztfw(ji,jj,jk) + e1e2t(ji,jj) / e3w_n(ji,jj,jk) * wmask(ji,jj,jk) & 346 327 & * ( ah_wslp2(ji,jj,jk) - akz(ji,jj,jk) ) & 347 328 & * ( ptb(ji,jj,jk-1,jn) - ptb(ji,jj,jk,jn) ) … … 358 339 ztfw(ji,jj,jk) = ztfw(ji,jj,jk) & 359 340 & + ah_wslp2(ji,jj,jk) * e1e2t(ji,jj) & 360 & * ( ptb(ji,jj,jk-1,jn) - ptb(ji,jj,jk,jn) ) * tmask(ji,jj,jk) / fse3w(ji,jj,jk)341 & * ( ptb(ji,jj,jk-1,jn) - ptb(ji,jj,jk,jn) ) / e3w_n(ji,jj,jk) * wmask(ji,jj,jk) 361 342 END DO 362 343 END DO … … 366 347 DO jj = 1, jpjm1 367 348 DO ji = fs_2, fs_jpim1 368 ztfw(ji,jj,jk) = ztfw(ji,jj,jk) + e1e2t(ji,jj) / fse3w(ji,jj,jk) * tmask(ji,jj,jk) &349 ztfw(ji,jj,jk) = ztfw(ji,jj,jk) + e1e2t(ji,jj) / e3w_n(ji,jj,jk) * wmask(ji,jj,jk) & 369 350 & * ( ah_wslp2(ji,jj,jk) * ( ptb (ji,jj,jk-1,jn) - ptb (ji,jj,jk,jn) ) & 370 351 & + akz (ji,jj,jk) * ( ptbb(ji,jj,jk-1,jn) - ptbb(ji,jj,jk,jn) ) ) … … 379 360 DO ji = fs_2, fs_jpim1 ! vector opt. 380 361 pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + zsign * ( ztfw (ji,jj,jk) - ztfw(ji,jj,jk+1) ) & 381 & / ( e1e2t(ji,jj) * fse3t_n(ji,jj,jk))362 & * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) 382 363 END DO 383 364 END DO -
trunk/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_triad.F90
r5836 r6140 11 11 !! tra_ldf_triad : update the tracer trend with the iso-neutral laplacian triad-operator 12 12 !!---------------------------------------------------------------------- 13 USE oce 14 USE dom_oce 15 USE phycst 16 USE trc_oce 17 USE zdf_oce 18 USE ldftra 19 USE ldfslp 20 USE traldf_iso 21 USE diaptr 22 USE zpshde 13 USE oce ! ocean dynamics and active tracers 14 USE dom_oce ! ocean space and time domain 15 USE phycst ! physical constants 16 USE trc_oce ! share passive tracers/Ocean variables 17 USE zdf_oce ! ocean vertical physics 18 USE ldftra ! lateral physics: eddy diffusivity 19 USE ldfslp ! lateral physics: iso-neutral slopes 20 USE traldf_iso ! lateral diffusion (Madec operator) (tra_ldf_iso routine) 21 USE diaptr ! poleward transport diagnostics 22 USE zpshde ! partial step: hor. derivative (zps_hde routine) 23 23 ! 24 USE in_out_manager 25 USE iom 26 USE lbclnk 27 USE lib_mpp 28 USE wrk_nemo 29 USE timing 24 USE in_out_manager ! I/O manager 25 USE iom ! I/O library 26 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 27 USE lib_mpp ! MPP library 28 USE wrk_nemo ! Memory Allocation 29 USE timing ! Timing 30 30 31 31 IMPLICIT NONE … … 37 37 38 38 !! * Substitutions 39 # include "domzgr_substitute.h90"40 39 # include "vectopt_loop_substitute.h90" 41 40 !!---------------------------------------------------------------------- … … 113 112 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~' 114 113 ENDIF 115 !116 114 ! ! set time step size (Euler/Leapfrog) 117 IF( neuler == 0 .AND. kt == kit000 ) THEN ; z2dt = rdt tra(1)! at nit000 (Euler)118 ELSE ; z2dt = 2.* rdt tra(1)! (Leapfrog)115 IF( neuler == 0 .AND. kt == kit000 ) THEN ; z2dt = rdt ! at nit000 (Euler) 116 ELSE ; z2dt = 2.* rdt ! (Leapfrog) 119 117 ENDIF 120 118 z1_2dt = 1._wp / z2dt … … 123 121 ELSE ; zsign = -1._wp 124 122 ENDIF 125 123 ! 126 124 !!---------------------------------------------------------------------- 127 125 !! 0 - calculate ah_wslp2, akz, and optionally zpsi_uw, zpsi_vw … … 142 140 DO jj = 1, jpjm1 143 141 DO ji = 1, fs_jpim1 144 ze3wr = 1._wp / fse3w(ji+ip,jj,jk+kp)145 zbu = e1e2u(ji,jj) * fse3u(ji,jj,jk)142 ze3wr = 1._wp / e3w_n(ji+ip,jj,jk+kp) 143 zbu = e1e2u(ji,jj) * e3u_n(ji,jj,jk) 146 144 zah = 0.25_wp * pahu(ji,jj,jk) 147 145 zslope_skew = triadi_g(ji+ip,jj,jk,1-ip,kp) 148 146 ! Subtract s-coordinate slope at t-points to give slope rel to s-surfaces (do this by *adding* gradient of depth) 149 zslope2 = zslope_skew + ( fsdept(ji+1,jj,jk) - fsdept(ji,jj,jk) ) * r1_e1u(ji,jj) * umask(ji,jj,jk+kp)147 zslope2 = zslope_skew + ( gdept_n(ji+1,jj,jk) - gdept_n(ji,jj,jk) ) * r1_e1u(ji,jj) * umask(ji,jj,jk+kp) 150 148 zslope2 = zslope2 *zslope2 151 149 ah_wslp2(ji+ip,jj,jk+kp) = ah_wslp2(ji+ip,jj,jk+kp) + zah * zbu * ze3wr * r1_e1e2t(ji+ip,jj) * zslope2 152 150 akz (ji+ip,jj,jk+kp) = akz (ji+ip,jj,jk+kp) + zah * r1_e1u(ji,jj) & 153 151 & * r1_e1u(ji,jj) * umask(ji,jj,jk+kp) 154 !152 ! 155 153 IF( ln_ldfeiv_dia ) zpsi_uw(ji,jj,jk+kp) = zpsi_uw(ji,jj,jk+kp) & 156 154 & + 0.25_wp * aeiu(ji,jj,jk) * e2u(ji,jj) * zslope_skew … … 166 164 DO jj = 1, jpjm1 167 165 DO ji = 1, fs_jpim1 168 ze3wr = 1.0_wp / fse3w(ji,jj+jp,jk+kp)169 zbv = e1e2v(ji,jj) * fse3v(ji,jj,jk)166 ze3wr = 1.0_wp / e3w_n(ji,jj+jp,jk+kp) 167 zbv = e1e2v(ji,jj) * e3v_n(ji,jj,jk) 170 168 zah = 0.25_wp * pahv(ji,jj,jk) 171 169 zslope_skew = triadj_g(ji,jj+jp,jk,1-jp,kp) 172 170 ! Subtract s-coordinate slope at t-points to give slope rel to s surfaces 173 171 ! (do this by *adding* gradient of depth) 174 zslope2 = zslope_skew + ( fsdept(ji,jj+1,jk) - fsdept(ji,jj,jk) ) * r1_e2v(ji,jj) * vmask(ji,jj,jk+kp)172 zslope2 = zslope_skew + ( gdept_n(ji,jj+1,jk) - gdept_n(ji,jj,jk) ) * r1_e2v(ji,jj) * vmask(ji,jj,jk+kp) 175 173 zslope2 = zslope2 * zslope2 176 174 ah_wslp2(ji,jj+jp,jk+kp) = ah_wslp2(ji,jj+jp,jk+kp) + zah * zbv * ze3wr * r1_e1e2t(ji,jj+jp) * zslope2 … … 193 191 DO ji = 1, fs_jpim1 194 192 akz(ji,jj,jk) = 16._wp * ah_wslp2(ji,jj,jk) & 195 & * ( akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ( fse3w(ji,jj,jk) * fse3w(ji,jj,jk) ) )193 & * ( akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ( e3w_n(ji,jj,jk) * e3w_n(ji,jj,jk) ) ) 196 194 END DO 197 195 END DO … … 201 199 DO jj = 1, jpjm1 202 200 DO ji = 1, fs_jpim1 203 ze3w_2 = fse3w(ji,jj,jk) * fse3w(ji,jj,jk)201 ze3w_2 = e3w_n(ji,jj,jk) * e3w_n(ji,jj,jk) 204 202 zcoef0 = z2dt * ( akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ze3w_2 ) 205 203 akz(ji,jj,jk) = MAX( zcoef0 - 0.5_wp , 0._wp ) * ze3w_2 * z1_2dt … … 250 248 ENDIF 251 249 ENDIF 252 250 ! 253 251 !!---------------------------------------------------------------------- 254 252 !! II - horizontal trend (full) … … 256 254 ! 257 255 DO jk = 1, jpkm1 258 !259 256 ! !== Vertical tracer gradient at level jk and jk+1 260 257 zdkt3d(:,:,1) = ( ptb(:,:,jk,jn) - ptb(:,:,jk+1,jn) ) * tmask(:,:,jk+1) … … 274 271 ze1ur = r1_e1u(ji,jj) 275 272 zdxt = zdit(ji,jj,jk) * ze1ur 276 ze3wr = 1._wp / fse3w(ji+ip,jj,jk+kp)273 ze3wr = 1._wp / e3w_n(ji+ip,jj,jk+kp) 277 274 zdzt = zdkt3d(ji+ip,jj,kp) * ze3wr 278 275 zslope_skew = triadi_g(ji+ip,jj,jk,1-ip,kp) 279 276 zslope_iso = triadi (ji+ip,jj,jk,1-ip,kp) 280 281 zbu = 0.25_wp * e1e2u(ji,jj) * fse3u(ji,jj,jk)277 ! 278 zbu = 0.25_wp * e1e2u(ji,jj) * e3u_n(ji,jj,jk) 282 279 ! ln_botmix_triad is .T. don't mask zah for bottom half cells !!gm ????? ahu is masked.... 283 280 zah = pahu(ji,jj,jk) … … 290 287 END DO 291 288 END DO 292 289 ! 293 290 DO jp = 0, 1 294 291 DO kp = 0, 1 … … 297 294 ze2vr = r1_e2v(ji,jj) 298 295 zdyt = zdjt(ji,jj,jk) * ze2vr 299 ze3wr = 1._wp / fse3w(ji,jj+jp,jk+kp)296 ze3wr = 1._wp / e3w_n(ji,jj+jp,jk+kp) 300 297 zdzt = zdkt3d(ji,jj+jp,kp) * ze3wr 301 298 zslope_skew = triadj_g(ji,jj+jp,jk,1-jp,kp) 302 299 zslope_iso = triadj(ji,jj+jp,jk,1-jp,kp) 303 zbv = 0.25_wp * e1e2v(ji,jj) * fse3v(ji,jj,jk)300 zbv = 0.25_wp * e1e2v(ji,jj) * e3v_n(ji,jj,jk) 304 301 ! ln_botmix_triad is .T. don't mask zah for bottom half cells !!gm ????? ahv is masked... 305 302 zah = pahv(ji,jj,jk) … … 312 309 END DO 313 310 END DO 314 311 ! 315 312 ELSE 316 313 ! 317 314 DO ip = 0, 1 !== Horizontal & vertical fluxes 318 315 DO kp = 0, 1 … … 321 318 ze1ur = r1_e1u(ji,jj) 322 319 zdxt = zdit(ji,jj,jk) * ze1ur 323 ze3wr = 1._wp / fse3w(ji+ip,jj,jk+kp)320 ze3wr = 1._wp / e3w_n(ji+ip,jj,jk+kp) 324 321 zdzt = zdkt3d(ji+ip,jj,kp) * ze3wr 325 322 zslope_skew = triadi_g(ji+ip,jj,jk,1-ip,kp) 326 323 zslope_iso = triadi(ji+ip,jj,jk,1-ip,kp) 327 328 zbu = 0.25_wp * e1e2u(ji,jj) * fse3u(ji,jj,jk)324 ! 325 zbu = 0.25_wp * e1e2u(ji,jj) * e3u_n(ji,jj,jk) 329 326 ! ln_botmix_triad is .F. mask zah for bottom half cells 330 327 zah = pahu(ji,jj,jk) * umask(ji,jj,jk+kp) ! pahu(ji+ip,jj,jk) ===>> ???? 331 328 zah_slp = zah * zslope_iso 332 IF( ln_ldfeiv ) zaei_slp = aeiu(ji,jj,jk) * zslope_skew ! fsaeit(ji+ip,jj,jk)*zslope_skew329 IF( ln_ldfeiv ) zaei_slp = aeiu(ji,jj,jk) * zslope_skew ! aeit(ji+ip,jj,jk)*zslope_skew 333 330 zftu(ji ,jj,jk ) = zftu(ji ,jj,jk ) - ( zah * zdxt + (zah_slp - zaei_slp) * zdzt ) * zbu * ze1ur 334 331 ztfw(ji+ip,jj,jk+kp) = ztfw(ji+ip,jj,jk+kp) - (zah_slp + zaei_slp) * zdxt * zbu * ze3wr … … 337 334 END DO 338 335 END DO 339 336 ! 340 337 DO jp = 0, 1 341 338 DO kp = 0, 1 … … 344 341 ze2vr = r1_e2v(ji,jj) 345 342 zdyt = zdjt(ji,jj,jk) * ze2vr 346 ze3wr = 1._wp / fse3w(ji,jj+jp,jk+kp)343 ze3wr = 1._wp / e3w_n(ji,jj+jp,jk+kp) 347 344 zdzt = zdkt3d(ji,jj+jp,kp) * ze3wr 348 345 zslope_skew = triadj_g(ji,jj+jp,jk,1-jp,kp) 349 346 zslope_iso = triadj(ji,jj+jp,jk,1-jp,kp) 350 zbv = 0.25_wp * e1e2v(ji,jj) * fse3v(ji,jj,jk)347 zbv = 0.25_wp * e1e2v(ji,jj) * e3v_n(ji,jj,jk) 351 348 ! ln_botmix_triad is .F. mask zah for bottom half cells 352 349 zah = pahv(ji,jj,jk) * vmask(ji,jj,jk+kp) ! pahv(ji,jj+jp,jk) ???? 353 350 zah_slp = zah * zslope_iso 354 IF( ln_ldfeiv ) zaei_slp = aeiv(ji,jj,jk) * zslope_skew ! fsaeit(ji,jj+jp,jk)*zslope_skew351 IF( ln_ldfeiv ) zaei_slp = aeiv(ji,jj,jk) * zslope_skew ! aeit(ji,jj+jp,jk)*zslope_skew 355 352 zftv(ji,jj,jk) = zftv(ji,jj,jk) - ( zah * zdyt + (zah_slp - zaei_slp) * zdzt ) * zbv * ze2vr 356 353 ztfw(ji,jj+jp,jk+kp) = ztfw(ji,jj+jp,jk+kp) - (zah_slp + zaei_slp) * zdyt * zbv * ze3wr … … 365 362 pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + zsign * ( zftu(ji-1,jj,jk) - zftu(ji,jj,jk) & 366 363 & + zftv(ji,jj-1,jk) - zftv(ji,jj,jk) ) & 367 & / ( e1e2t(ji,jj) * fse3t(ji,jj,jk) )364 & / ( e1e2t(ji,jj) * e3t_n(ji,jj,jk) ) 368 365 END DO 369 366 END DO … … 376 373 DO jj = 1, jpjm1 377 374 DO ji = fs_2, fs_jpim1 378 ztfw(ji,jj,jk) = ztfw(ji,jj,jk) - e1e2t(ji,jj) / fse3w(ji,jj,jk) * tmask(ji,jj,jk) &375 ztfw(ji,jj,jk) = ztfw(ji,jj,jk) - e1e2t(ji,jj) / e3w_n(ji,jj,jk) * tmask(ji,jj,jk) & 379 376 & * ( ah_wslp2(ji,jj,jk) - akz(ji,jj,jk) ) & 380 377 & * ( ptb(ji,jj,jk-1,jn) - ptb(ji,jj,jk,jn) ) … … 388 385 DO jj = 1, jpjm1 389 386 DO ji = fs_2, fs_jpim1 390 ztfw(ji,jj,jk) = ztfw(ji,jj,jk) - e1e2t(ji,jj) / fse3w(ji,jj,jk) * tmask(ji,jj,jk) &387 ztfw(ji,jj,jk) = ztfw(ji,jj,jk) - e1e2t(ji,jj) / e3w_n(ji,jj,jk) * tmask(ji,jj,jk) & 391 388 & * ah_wslp2(ji,jj,jk) * ( ptb(ji,jj,jk-1,jn) - ptb(ji,jj,jk,jn) ) 392 389 END DO … … 397 394 DO jj = 1, jpjm1 398 395 DO ji = fs_2, fs_jpim1 399 ztfw(ji,jj,jk) = ztfw(ji,jj,jk) - e1e2t(ji,jj) / fse3w(ji,jj,jk) * tmask(ji,jj,jk) &396 ztfw(ji,jj,jk) = ztfw(ji,jj,jk) - e1e2t(ji,jj) / e3w_n(ji,jj,jk) * tmask(ji,jj,jk) & 400 397 & * ( ah_wslp2(ji,jj,jk) * ( ptb (ji,jj,jk-1,jn) - ptb (ji,jj,jk,jn) ) & 401 398 & + akz (ji,jj,jk) * ( ptbb(ji,jj,jk-1,jn) - ptbb(ji,jj,jk,jn) ) ) … … 410 407 DO ji = fs_2, fs_jpim1 ! vector opt. 411 408 pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + zsign * ( ztfw(ji,jj,jk+1) - ztfw(ji,jj,jk) ) & 412 & / ( e1e2t(ji,jj) * fse3t(ji,jj,jk) )409 & / ( e1e2t(ji,jj) * e3t_n(ji,jj,jk) ) 413 410 END DO 414 411 END DO -
trunk/NEMOGCM/NEMO/OPA_SRC/TRA/tranpc.F90
r5386 r6140 13 13 14 14 !!---------------------------------------------------------------------- 15 !! tra_npc : apply the non penetrative convection scheme16 !!---------------------------------------------------------------------- 17 USE oce 18 USE dom_oce 19 USE phycst 20 USE zdf_oce 21 USE trd_oce 22 USE trdtra 23 USE eosbn2 15 !! tra_npc : apply the non penetrative convection scheme 16 !!---------------------------------------------------------------------- 17 USE oce ! ocean dynamics and active tracers 18 USE dom_oce ! ocean space and time domain 19 USE phycst ! physical constants 20 USE zdf_oce ! ocean vertical physics 21 USE trd_oce ! ocean active tracer trends 22 USE trdtra ! ocean active tracer trends 23 USE eosbn2 ! equation of state (eos routine) 24 24 ! 25 USE lbclnk 26 USE in_out_manager 27 USE lib_mpp 28 USE wrk_nemo 29 USE timing 25 USE lbclnk ! lateral boundary conditions (or mpp link) 26 USE in_out_manager ! I/O manager 27 USE lib_mpp ! MPP library 28 USE wrk_nemo ! Memory Allocation 29 USE timing ! Timing 30 30 31 31 IMPLICIT NONE … … 35 35 36 36 !! * Substitutions 37 # include "domzgr_substitute.h90"38 37 # include "vectopt_loop_substitute.h90" 39 38 !!---------------------------------------------------------------------- … … 55 54 !! (i.e. static stability computed locally) 56 55 !! 57 !! ** Action : - (ta,sa) after the application odthe npc scheme56 !! ** Action : - tsa: after tracers with the application of the npc scheme 58 57 !! - send the associated trends for on-line diagnostics (l_trdtra=T) 59 58 !! … … 115 114 zvts(:,jp_tem) = tsa(ji,jj,:,jp_tem) ! temperature 116 115 zvts(:,jp_sal) = tsa(ji,jj,:,jp_sal) ! salinity 117 116 ! 118 117 zvab(:,jp_tem) = zab(ji,jj,:,jp_tem) ! Alpha 119 118 zvab(:,jp_sal) = zab(ji,jj,:,jp_sal) ! Beta 120 119 zvn2(:) = zn2(ji,jj,:) ! N^2 121 120 ! 122 121 IF( l_LB_debug ) THEN !LB debug: 123 122 lp_monitor_point = .FALSE. … … 126 125 lp_monitor_point = (narea == nncpu).AND.lp_monitor_point 127 126 ENDIF !LB debug end 128 127 ! 129 128 ikbot = mbkt(ji,jj) ! ikbot: ocean bottom T-level 130 129 ikp = 1 ! because N2 is irrelevant at the surface level (will start at ikp=2) … … 132 131 jiter = 0 133 132 l_column_treated = .FALSE. 134 133 ! 135 134 DO WHILE ( .NOT. l_column_treated ) 136 135 ! 137 136 jiter = jiter + 1 138 137 ! 139 138 IF( jiter >= 400 ) EXIT 140 139 ! 141 140 l_bottom_reached = .FALSE. 142 141 ! 143 142 DO WHILE ( .NOT. l_bottom_reached ) 144 143 ! 145 144 ikp = ikp + 1 146 145 ! 147 146 !! Testing level ikp for instability 148 147 !! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 149 148 IF( zvn2(ikp) < -zn2_zero ) THEN ! Instability found! 150 149 ! 151 150 ilayer = ilayer + 1 ! yet another instable portion of the water column found.... 152 151 ! 153 152 IF( lp_monitor_point ) THEN 154 153 WRITE(numout,*) … … 165 164 WRITE(numout,*) 166 165 ENDIF 167 168 166 ! 169 167 IF( jiter == 1 ) inpcc = inpcc + 1 170 168 ! 171 169 IF( lp_monitor_point ) WRITE(numout, *) 'Negative N2 at ikp =',ikp,' for layer #', ilayer 172 170 ! 173 171 !! ikup is the uppermost point where mixing will start: 174 172 ikup = ikp - 1 ! ikup is always "at most at ikp-1", less if neutral levels overlying 175 173 ! 176 174 !! If the points above ikp-1 have N2 == 0 they must also be mixed: 177 175 IF( ikp > 2 ) THEN … … 184 182 END DO 185 183 ENDIF 186 184 ! 187 185 IF( ikup < 1 ) CALL ctl_stop( 'tra_npc : PROBLEM #1') 188 186 ! 189 187 zsum_temp = 0._wp 190 188 zsum_sali = 0._wp … … 195 193 DO jk = ikup, ikbot ! Inside the instable (and overlying neutral) portion of the column 196 194 ! 197 zdz = fse3t(ji,jj,jk)195 zdz = e3t_n(ji,jj,jk) 198 196 zsum_temp = zsum_temp + zvts(jk,jp_tem)*zdz 199 197 zsum_sali = zsum_sali + zvts(jk,jp_sal)*zdz … … 244 242 245 243 !! Interpolating alfa and beta at W point: 246 zrw = ( fsdepw(ji,jj,jk ) - fsdept(ji,jj,jk)) &247 & / ( fsdept(ji,jj,jk-1) - fsdept(ji,jj,jk))244 zrw = (gdepw_n(ji,jj,jk ) - gdept_n(ji,jj,jk)) & 245 & / (gdept_n(ji,jj,jk-1) - gdept_n(ji,jj,jk)) 248 246 zaw = zvab(jk,jp_tem) * (1._wp - zrw) + zvab(jk-1,jp_tem) * zrw 249 247 zbw = zvab(jk,jp_sal) * (1._wp - zrw) + zvab(jk-1,jp_sal) * zrw … … 252 250 zvn2(jk) = grav*( zaw * ( zvts(jk-1,jp_tem) - zvts(jk,jp_tem) ) & 253 251 & - zbw * ( zvts(jk-1,jp_sal) - zvts(jk,jp_sal) ) ) & 254 & / fse3w(ji,jj,jk) * tmask(ji,jj,jk)252 & / e3w_n(ji,jj,jk) * tmask(ji,jj,jk) 255 253 256 254 !! OR, faster => just considering the vertical gradient of density -
trunk/NEMOGCM/NEMO/OPA_SRC/TRA/tranxt.F90
r5930 r6140 28 28 USE sbc_oce ! surface boundary condition: ocean 29 29 USE sbcrnf ! river runoffs 30 USE sbcisf ! ice shelf melting /freezing30 USE sbcisf ! ice shelf melting 31 31 USE zdf_oce ! ocean vertical mixing 32 32 USE domvvl ! variable volume … … 56 56 PUBLIC tra_nxt_vvl ! to be used in trcnxt 57 57 58 59 58 !! * Substitutions 60 # include " domzgr_substitute.h90"59 # include "vectopt_loop_substitute.h90" 61 60 !!---------------------------------------------------------------------- 62 61 !! NEMO/OPA 3.3 , NEMO-Consortium (2010) … … 86 85 !! domains (lk_agrif=T) 87 86 !! 88 !! ** Action : - (tb,sb) and (tn,sn) ready for the next time step 89 !! 87 !! ** Action : - tsb & tsn ready for the next time step 90 88 !!---------------------------------------------------------------------- 91 89 INTEGER, INTENT(in) :: kt ! ocean time-step index 92 90 !! 93 INTEGER :: j k, jn! dummy loop indices94 REAL(wp) :: zfact ! local scalars91 INTEGER :: ji, jj, jk, jn ! dummy loop indices 92 REAL(wp) :: zfact ! local scalars 95 93 REAL(wp), POINTER, DIMENSION(:,:,:) :: ztrdt, ztrds 96 94 !!---------------------------------------------------------------------- … … 106 104 ! Update after tracer on domain lateral boundaries 107 105 ! 108 !109 106 #if defined key_agrif 110 107 CALL Agrif_tra ! AGRIF zoom boundaries … … 119 116 120 117 ! set time step size (Euler/Leapfrog) 121 IF( neuler == 0 .AND. kt == nit000 ) THEN ; r2dt ra(:) = rdttra(:)! at nit000 (Euler)122 ELSEIF( kt <= nit000 + 1 ) THEN ; r2dt ra(:) = 2._wp* rdttra(:)! at nit000 or nit000+1 (Leapfrog)118 IF( neuler == 0 .AND. kt == nit000 ) THEN ; r2dt = rdt ! at nit000 (Euler) 119 ELSEIF( kt <= nit000 + 1 ) THEN ; r2dt = 2._wp* rdt ! at nit000 or nit000+1 (Leapfrog) 123 120 ENDIF 124 121 … … 140 137 END DO 141 138 END DO 139 ! 142 140 ELSE ! Leap-Frog + Asselin filter time stepping 143 141 ! 144 IF( l k_vvl ) THEN ; CALL tra_nxt_vvl( kt, nit000, rdttra, 'TRA', tsb, tsn, tsa, &145 & sbc_tsc, sbc_tsc_b, jpts ) ! variable volume level (vvl)146 ELSE ; CALL tra_nxt_fix( kt, nit000, 'TRA', tsb, tsn, tsa, jpts ) ! fixed volume level142 IF( ln_linssh ) THEN ; CALL tra_nxt_fix( kt, nit000, 'TRA', tsb, tsn, tsa, jpts ) ! linear free surface 143 ELSE ; CALL tra_nxt_vvl( kt, nit000, rdt, 'TRA', tsb, tsn, tsa, & 144 & sbc_tsc, sbc_tsc_b, jpts ) ! non-linear free surface 147 145 ENDIF 146 ! 147 DO jn = 1, jpts 148 CALL lbc_lnk( tsb(:,:,:,jn), 'T', 1._wp ) 149 CALL lbc_lnk( tsn(:,:,:,jn), 'T', 1._wp ) 150 CALL lbc_lnk( tsa(:,:,:,jn), 'T', 1._wp ) 151 END DO 148 152 ENDIF 149 153 ! 150 ! trends computation151 154 IF( l_trdtra ) THEN ! trend of the Asselin filter (tb filtered - tb)/dt 152 155 DO jk = 1, jpkm1 153 zfact = 1._wp / r2dt ra(jk)156 zfact = 1._wp / r2dt 154 157 ztrdt(:,:,jk) = ( tsb(:,:,jk,jp_tem) - ztrdt(:,:,jk) ) * zfact 155 158 ztrds(:,:,jk) = ( tsb(:,:,jk,jp_sal) - ztrds(:,:,jk) ) * zfact … … 179 182 !! - swap tracer fields to prepare the next time_step. 180 183 !! 181 !! ** Action : - (tb,sb) and (tn,sn) ready for the next time step 182 !! 183 !!---------------------------------------------------------------------- 184 INTEGER , INTENT(in ) :: kt ! ocean time-step index 185 INTEGER , INTENT(in ) :: kit000 ! first time step index 186 CHARACTER(len=3), INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator) 187 INTEGER , INTENT(in ) :: kjpt ! number of tracers 188 REAL(wp) , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt) :: ptb ! before tracer fields 189 REAL(wp) , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt) :: ptn ! now tracer fields 190 REAL(wp) , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt) :: pta ! tracer trend 184 !! ** Action : - tsb & tsn ready for the next time step 185 !!---------------------------------------------------------------------- 186 INTEGER , INTENT(in ) :: kt ! ocean time-step index 187 INTEGER , INTENT(in ) :: kit000 ! first time step index 188 CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator) 189 INTEGER , INTENT(in ) :: kjpt ! number of tracers 190 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: ptb ! before tracer fields 191 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: ptn ! now tracer fields 192 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pta ! tracer trend 191 193 ! 192 194 INTEGER :: ji, jj, jk, jn ! dummy loop indices 193 195 REAL(wp) :: ztn, ztd ! local scalars 194 196 !!---------------------------------------------------------------------- 195 197 ! 196 198 IF( kt == kit000 ) THEN 197 199 IF(lwp) WRITE(numout,*) … … 200 202 ENDIF 201 203 ! 202 !203 204 DO jn = 1, kjpt 204 205 ! 205 206 DO jk = 1, jpkm1 206 DO jj = 1, jpj207 DO ji = 1, jpi207 DO jj = 2, jpjm1 208 DO ji = fs_2, fs_jpim1 208 209 ztn = ptn(ji,jj,jk,jn) 209 ztd = pta(ji,jj,jk,jn) - 2. * ztn + ptb(ji,jj,jk,jn) ! time laplacian on tracers 210 ! 211 ptb(ji,jj,jk,jn) = ztn + atfp * ztd ! ptb <-- filtered ptn 212 ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn) ! ptn <-- pta 213 ! 210 ztd = pta(ji,jj,jk,jn) - 2._wp * ztn + ptb(ji,jj,jk,jn) ! time laplacian on tracers 211 ! 212 ptb(ji,jj,jk,jn) = ztn + atfp * ztd ! ptb <-- filtered ptn 213 ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn) ! ptn <-- pta 214 214 END DO 215 215 END DO … … 230 230 !! ** Method : - Apply a thickness weighted Asselin time filter on now fields. 231 231 !! - swap tracer fields to prepare the next time_step. 232 !! This can be summurized for tempearture as:233 !! 234 !! ** Action : - (tb,sb) and (tn,sn) ready for the next time step235 !! 236 !! ----------------------------------------------------------------------237 INTEGER , INTENT(in ) :: kt ! ocean time-step index238 INTEGER , INTENT(in ) :: kit000 ! first timestep index239 REAL(wp) , INTENT(in ), DIMENSION(jpk) :: p2dt ! time-step240 CHARACTER(len=3), INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator)241 INTEGER , INTENT(in ) :: kjpt ! number of tracers242 REAL(wp) , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt) :: ptb ! before tracer fields243 REAL(wp) , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt) :: ptn ! nowtracer fields244 REAL(wp) , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt) :: pta ! tracer trend245 REAL(wp) , INTENT(in ), DIMENSION(jpi,jpj,kjpt) :: psbc_tc ! surface tracer content246 REAL(wp) , INTENT(in ), DIMENSION(jpi,jpj,kjpt) :: psbc_tc_b ! beforesurface tracer content247 248 ! !232 !! tb = ( e3t_n*tn + atfp*[ e3t_b*tb - 2 e3t_n*tn + e3t_a*ta ] ) 233 !! /( e3t_n + atfp*[ e3t_b - 2 e3t_n + e3t_a ] ) 234 !! tn = ta 235 !! 236 !! ** Action : - tsb & tsn ready for the next time step 237 !!---------------------------------------------------------------------- 238 INTEGER , INTENT(in ) :: kt ! ocean time-step index 239 INTEGER , INTENT(in ) :: kit000 ! first time step index 240 REAL(wp) , INTENT(in ) :: p2dt ! time-step 241 CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator) 242 INTEGER , INTENT(in ) :: kjpt ! number of tracers 243 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: ptb ! before tracer fields 244 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: ptn ! now tracer fields 245 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pta ! tracer trend 246 REAL(wp), DIMENSION(jpi,jpj ,kjpt), INTENT(in ) :: psbc_tc ! surface tracer content 247 REAL(wp), DIMENSION(jpi,jpj ,kjpt), INTENT(in ) :: psbc_tc_b ! before surface tracer content 248 ! 249 249 LOGICAL :: ll_traqsr, ll_rnf, ll_isf ! local logical 250 250 INTEGER :: ji, jj, jk, jn ! dummy loop indices … … 262 262 ll_traqsr = ln_traqsr ! active tracers case and solar penetration 263 263 ll_rnf = ln_rnf ! active tracers case and river runoffs 264 IF (nn_isf .GE. 1) THEN 265 ll_isf = .TRUE. ! active tracers case and ice shelf melting/freezing 266 ELSE 267 ll_isf = .FALSE. 268 END IF 269 ELSE 270 ll_traqsr = .FALSE. ! active tracers case and NO solar penetration 271 ll_rnf = .FALSE. ! passive tracers or NO river runoffs 272 ll_isf = .FALSE. ! passive tracers or NO ice shelf melting/freezing 264 ll_isf = ln_isf ! active tracers case and ice shelf melting 265 ELSE ! passive tracers case 266 ll_traqsr = .FALSE. ! NO solar penetration 267 ll_rnf = .FALSE. ! NO river runoffs ???? !!gm BUG ? 268 ll_isf = .FALSE. ! NO ice shelf melting/freezing !!gm BUG ?? 273 269 ENDIF 274 270 ! 275 271 DO jn = 1, kjpt 276 272 DO jk = 1, jpkm1 277 zfact1 = atfp * p2dt (jk)278 zfact2 = zfact1 /rau0279 DO jj = 1, jpj280 DO ji = 1, jpi281 ze3t_b = fse3t_b(ji,jj,jk)282 ze3t_n = fse3t_n(ji,jj,jk)283 ze3t_a = fse3t_a(ji,jj,jk)273 zfact1 = atfp * p2dt 274 zfact2 = zfact1 * r1_rau0 275 DO jj = 2, jpjm1 276 DO ji = fs_2, fs_jpim1 277 ze3t_b = e3t_b(ji,jj,jk) 278 ze3t_n = e3t_n(ji,jj,jk) 279 ze3t_a = e3t_a(ji,jj,jk) 284 280 ! ! tracer content at Before, now and after 285 281 ztc_b = ptb(ji,jj,jk,jn) * ze3t_b … … 299 295 ztc_f = ztc_f - zfact1 * ( psbc_tc(ji,jj,jn) - psbc_tc_b(ji,jj,jn) ) 300 296 ENDIF 301 297 ! 302 298 ! solar penetration (temperature only) 303 299 IF( ll_traqsr .AND. jn == jp_tem .AND. jk <= nksr ) & 304 300 & ztc_f = ztc_f - zfact1 * ( qsr_hc(ji,jj,jk) - qsr_hc_b(ji,jj,jk) ) 305 301 ! 306 302 ! river runoff 307 303 IF( ll_rnf .AND. jk <= nk_rnf(ji,jj) ) & 308 304 & ztc_f = ztc_f - zfact1 * ( rnf_tsc(ji,jj,jn) - rnf_tsc_b(ji,jj,jn) ) & 309 & * fse3t_n(ji,jj,jk) / h_rnf(ji,jj)310 305 & * e3t_n(ji,jj,jk) / h_rnf(ji,jj) 306 ! 311 307 ! ice shelf 312 308 IF( ll_isf ) THEN … … 314 310 IF ( jk >= misfkt(ji,jj) .AND. jk < misfkb(ji,jj) ) & 315 311 ztc_f = ztc_f - zfact1 * ( risf_tsc(ji,jj,jn) - risf_tsc_b(ji,jj,jn) ) & 316 & * fse3t_n(ji,jj,jk) * r1_hisf_tbl (ji,jj)312 & * e3t_n(ji,jj,jk) * r1_hisf_tbl (ji,jj) 317 313 ! level partially include in Losch_2008 ice shelf boundary layer 318 314 IF ( jk == misfkb(ji,jj) ) & 319 315 ztc_f = ztc_f - zfact1 * ( risf_tsc(ji,jj,jn) - risf_tsc_b(ji,jj,jn) ) & 320 & * fse3t_n(ji,jj,jk) * r1_hisf_tbl (ji,jj) * ralpha(ji,jj)316 & * e3t_n(ji,jj,jk) * r1_hisf_tbl (ji,jj) * ralpha(ji,jj) 321 317 END IF 322 318 ! 323 319 ze3t_f = 1.e0 / ze3t_f 324 320 ptb(ji,jj,jk,jn) = ztc_f * ze3t_f ! ptb <-- ptn filtered -
trunk/NEMOGCM/NEMO/OPA_SRC/TRA/traqsr.F90
r5836 r6140 2 2 !!====================================================================== 3 3 !! *** MODULE traqsr *** 4 !! Ocean physics: solar radiation penetration in the top ocean levels4 !! Ocean physics: solar radiation penetration in the top ocean levels 5 5 !!====================================================================== 6 6 !! History : OPA ! 1990-10 (B. Blanke) Original code … … 10 10 !! - ! 2005-11 (G. Madec) zco, zps, sco coordinate 11 11 !! 3.2 ! 2009-04 (G. Madec & NEMO team) 12 !! 4.0 ! 2012-05 (C. Rousset) store attenuation coef for use in ice model 12 !! 3.6 ! 2012-05 (C. Rousset) store attenuation coef for use in ice model 13 !! 3.7 ! 2015-11 (G. Madec, A. Coward) remove optimisation for fix volume 13 14 !!---------------------------------------------------------------------- 14 15 15 16 !!---------------------------------------------------------------------- 16 !! tra_qsr : trend due to the solar radiation penetration17 !! tra_qsr_init : solar radiation penetration initialization17 !! tra_qsr : temperature trend due to the penetration of solar radiation 18 !! tra_qsr_init : initialization of the qsr penetration 18 19 !!---------------------------------------------------------------------- 19 USE oce ! ocean dynamics and active tracers 20 USE dom_oce ! ocean space and time domain 21 USE sbc_oce ! surface boundary condition: ocean 22 USE trc_oce ! share SMS/Ocean variables 20 USE oce ! ocean dynamics and active tracers 21 USE phycst ! physical constants 22 USE dom_oce ! ocean space and time domain 23 USE sbc_oce ! surface boundary condition: ocean 24 USE trc_oce ! share SMS/Ocean variables 23 25 USE trd_oce ! trends: ocean variables 24 26 USE trdtra ! trends manager: tracers 25 USE in_out_manager ! I/O manager 26 USE phycst ! physical constants 27 USE prtctl ! Print control 28 USE iom ! I/O manager 29 USE fldread ! read input fields 30 USE restart ! ocean restart 31 USE lib_mpp ! MPP library 27 ! 28 USE in_out_manager ! I/O manager 29 USE prtctl ! Print control 30 USE iom ! I/O manager 31 USE fldread ! read input fields 32 USE restart ! ocean restart 33 USE lib_mpp ! MPP library 34 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 32 35 USE wrk_nemo ! Memory Allocation 33 36 USE timing ! Timing … … 49 52 REAL(wp), PUBLIC :: rn_si0 !: very near surface depth of extinction (RGB & 2 bands) 50 53 REAL(wp), PUBLIC :: rn_si1 !: deepest depth of extinction (water type I) (2 bands) 54 ! 55 INTEGER , PUBLIC :: nksr !: levels below which the light cannot penetrate (depth larger than 391 m) 51 56 52 ! Module variables 53 REAL(wp) :: xsi0r !: inverse of rn_si0 54 REAL(wp) :: xsi1r !: inverse of rn_si1 57 INTEGER, PARAMETER :: np_RGB = 1 ! R-G-B light penetration with constant Chlorophyll 58 INTEGER, PARAMETER :: np_RGBc = 2 ! R-G-B light penetration with Chlorophyll data 59 INTEGER, PARAMETER :: np_2BD = 3 ! 2 bands light penetration 60 INTEGER, PARAMETER :: np_BIO = 4 ! bio-model light penetration 61 ! 62 INTEGER :: nqsr ! user choice of the type of light penetration 63 REAL(wp) :: xsi0r ! inverse of rn_si0 64 REAL(wp) :: xsi1r ! inverse of rn_si1 65 ! 66 REAL(wp) , DIMENSION(3,61) :: rkrgb ! tabulated attenuation coefficients for RGB absorption 55 67 TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_chl ! structure of input Chl (file informations, fields read) 56 INTEGER, PUBLIC :: nksr ! levels below which the light cannot penetrate ( depth larger than 391 m)57 REAL(wp), DIMENSION(3,61) :: rkrgb !: tabulated attenuation coefficients for RGB absorption58 68 59 69 !! * Substitutions 60 # include "domzgr_substitute.h90"61 70 # include "vectopt_loop_substitute.h90" 62 71 !!---------------------------------------------------------------------- … … 72 81 !! 73 82 !! ** Purpose : Compute the temperature trend due to the solar radiation 74 !! penetration and add it to the general temperature trend.83 !! penetration and add it to the general temperature trend. 75 84 !! 76 85 !! ** Method : The profile of the solar radiation within the ocean is defined … … 83 92 !! all heat which has not been absorbed in the above levels is put 84 93 !! in the last ocean level. 85 !! In z-coordinate case, the computation is only done down to the 86 !! level where I(k) < 1.e-15 W/m2. In addition, the coefficients 87 !! used for the computation are calculated one for once as they 88 !! depends on k only. 94 !! The computation is only done down to the level where 95 !! I(k) < 1.e-15 W/m2 (i.e. over the top nksr levels) . 89 96 !! 90 97 !! ** Action : - update ta with the penetrative solar radiation trend 91 !! - s ave the trend in ttrd ('key_trdtra')98 !! - send trend for further diagnostics (l_trdtra=T) 92 99 !! 93 100 !! Reference : Jerlov, N. G., 1968 Optical Oceanography, Elsevier, 194pp. 94 101 !! Lengaigne et al. 2007, Clim. Dyn., V28, 5, 503-516. 95 102 !!---------------------------------------------------------------------- 96 !97 103 INTEGER, INTENT(in) :: kt ! ocean time-step 98 104 ! 99 INTEGER :: ji, jj, jk ! dummy loop indices100 INTEGER :: irgb ! local integers101 REAL(wp) :: zchl, zcoef, z fact! local scalars102 REAL(wp) :: zc0 , zc1, zc2, zc3! - -105 INTEGER :: ji, jj, jk ! dummy loop indices 106 INTEGER :: irgb ! local integers 107 REAL(wp) :: zchl, zcoef, z1_2 ! local scalars 108 REAL(wp) :: zc0 , zc1 , zc2 , zc3 ! - - 103 109 REAL(wp) :: zzc0, zzc1, zzc2, zzc3 ! - - 104 REAL(wp) :: zz0 , zz1, z1_e3t! - -105 REAL(wp), POINTER, DIMENSION(:,: ):: zekb, zekg, zekr110 REAL(wp) :: zz0 , zz1 ! - - 111 REAL(wp), POINTER, DIMENSION(:,:) :: zekb, zekg, zekr 106 112 REAL(wp), POINTER, DIMENSION(:,:,:) :: ze0, ze1, ze2, ze3, zea, ztrdt 113 REAL(wp), POINTER, DIMENSION(:,:,:) :: zetot 107 114 !!---------------------------------------------------------------------- 108 115 ! 109 116 IF( nn_timing == 1 ) CALL timing_start('tra_qsr') 110 !111 CALL wrk_alloc( jpi, jpj, zekb, zekg, zekr )112 CALL wrk_alloc( jpi, jpj, jpk, ze0, ze1, ze2, ze3, zea )113 117 ! 114 118 IF( kt == nit000 ) THEN … … 116 120 IF(lwp) WRITE(numout,*) 'tra_qsr : penetration of the surface solar radiation' 117 121 IF(lwp) WRITE(numout,*) '~~~~~~~' 118 IF( .NOT.ln_traqsr ) RETURN 119 ENDIF 120 121 IF( l_trdtra ) THEN ! Save ta and sa trends 122 CALL wrk_alloc( jpi, jpj, jpk, ztrdt ) 122 ENDIF 123 ! 124 IF( l_trdtra ) THEN ! trends diagnostic: save the input temperature trend 125 CALL wrk_alloc( jpi,jpj,jpk, ztrdt ) 123 126 ztrdt(:,:,:) = tsa(:,:,:,jp_tem) 124 127 ENDIF 125 126 ! Set before qsr tracer content field 127 ! *********************************** 128 IF( kt == nit000 ) THEN ! Set the forcing field at nit000 - 1 129 ! ! ----------------------------------- 130 qsr_hc(:,:,:) = 0.e0 131 ! 132 IF( ln_rstart .AND. & ! Restart: read in restart file 133 & iom_varid( numror, 'qsr_hc_b', ldstop = .FALSE. ) > 0 ) THEN 134 IF(lwp) WRITE(numout,*) ' nit000-1 qsr tracer content forcing field red in the restart file' 135 zfact = 0.5e0 128 ! 129 ! !-----------------------------------! 130 ! ! before qsr induced heat content ! 131 ! !-----------------------------------! 132 IF( kt == nit000 ) THEN !== 1st time step ==! 133 !!gm case neuler not taken into account.... 134 IF( ln_rstart .AND. iom_varid( numror, 'qsr_hc_b', ldstop = .FALSE. ) > 0 ) THEN ! read in restart 135 IF(lwp) WRITE(numout,*) ' nit000-1 qsr tracer content forcing field read in the restart file' 136 z1_2 = 0.5_wp 136 137 CALL iom_get( numror, jpdom_autoglo, 'qsr_hc_b', qsr_hc_b ) ! before heat content trend due to Qsr flux 137 138 ELSE ! No restart or restart not found: Euler forward time stepping 138 z fact = 1.e0139 qsr_hc_b(:,:,:) = 0. e0139 z1_2 = 1._wp 140 qsr_hc_b(:,:,:) = 0._wp 140 141 ENDIF 141 ELSE ! Swap of forcing field 142 ! ! --------------------- 143 zfact = 0.5e0 142 ELSE !== Swap of qsr heat content ==! 143 z1_2 = 0.5_wp 144 144 qsr_hc_b(:,:,:) = qsr_hc(:,:,:) 145 145 ENDIF 146 ! Compute now qsr tracer content field 147 ! ************************************ 148 149 ! ! ============================================== ! 150 IF( lk_qsr_bio .AND. ln_qsr_bio ) THEN ! bio-model fluxes : all vertical coordinates ! 151 ! ! ============================================== ! 152 DO jk = 1, jpkm1 146 ! 147 ! !--------------------------------! 148 SELECT CASE( nqsr ) ! now qsr induced heat content ! 149 ! !--------------------------------! 150 ! 151 CASE( np_BIO ) !== bio-model fluxes ==! 152 ! 153 DO jk = 1, nksr 153 154 qsr_hc(:,:,jk) = r1_rau0_rcp * ( etot3(:,:,jk) - etot3(:,:,jk+1) ) 154 155 END DO 155 ! Add to the general trend 156 DO jk = 1, jpkm1 157 DO jj = 2, jpjm1 158 DO ji = fs_2, fs_jpim1 ! vector opt. 159 z1_e3t = zfact / fse3t(ji,jj,jk) 160 tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem) + ( qsr_hc_b(ji,jj,jk) + qsr_hc(ji,jj,jk) ) * z1_e3t 156 ! 157 CASE( np_RGB , np_RGBc ) !== R-G-B fluxes ==! 158 ! 159 CALL wrk_alloc( jpi,jpj, zekb, zekg, zekr ) 160 CALL wrk_alloc( jpi,jpj,jpk, ze0, ze1, ze2, ze3, zea ) 161 ! 162 IF( nqsr == np_RGBc ) THEN !* Variable Chlorophyll 163 CALL fld_read( kt, 1, sf_chl ) ! Read Chl data and provides it at the current time step 164 DO jj = 2, jpjm1 ! Separation in R-G-B depending of the surface Chl 165 DO ji = fs_2, fs_jpim1 166 zchl = MIN( 10. , MAX( 0.03, sf_chl(1)%fnow(ji,jj,1) ) ) 167 irgb = NINT( 41 + 20.*LOG10(zchl) + 1.e-15 ) 168 zekb(ji,jj) = rkrgb(1,irgb) 169 zekg(ji,jj) = rkrgb(2,irgb) 170 zekr(ji,jj) = rkrgb(3,irgb) 161 171 END DO 162 172 END DO 163 END DO 164 CALL iom_put( 'qsr3d', etot3 ) ! Shortwave Radiation 3D distribution 165 ! clem: store attenuation coefficient of the first ocean level 166 IF ( ln_qsr_ice ) THEN 167 DO jj = 1, jpj 168 DO ji = 1, jpi 169 IF ( qsr(ji,jj) /= 0._wp ) THEN 170 fraqsr_1lev(ji,jj) = ( qsr_hc(ji,jj,1) / ( r1_rau0_rcp * qsr(ji,jj) ) ) 171 ELSE 172 fraqsr_1lev(ji,jj) = 1. 173 ENDIF 173 ELSE !* constant chrlorophyll 174 zchl = 0.05 ! constant chlorophyll 175 ! ! Separation in R-G-B depending of the chlorophyll 176 irgb = NINT( 41 + 20.*LOG10( zchl ) + 1.e-15 ) 177 DO jj = 2, jpjm1 178 DO ji = fs_2, fs_jpim1 179 zekb(ji,jj) = rkrgb(1,irgb) 180 zekg(ji,jj) = rkrgb(2,irgb) 181 zekr(ji,jj) = rkrgb(3,irgb) 174 182 END DO 175 183 END DO 176 184 ENDIF 177 ! ! ============================================== ! 178 ELSE ! Ocean alone : 179 ! ! ============================================== ! 180 ! 181 ! ! ------------------------- ! 182 IF( ln_qsr_rgb) THEN ! R-G-B light penetration ! 183 ! ! ------------------------- ! 184 ! Set chlorophyl concentration 185 IF( nn_chldta == 1 .OR. lk_vvl ) THEN !* Variable Chlorophyll or ocean volume 186 ! 187 IF( nn_chldta == 1 ) THEN !* Variable Chlorophyll 188 ! 189 CALL fld_read( kt, 1, sf_chl ) ! Read Chl data and provides it at the current time step 190 ! 191 DO jj = 1, jpj ! Separation in R-G-B depending of the surface Chl 192 DO ji = 1, jpi 193 zchl = MIN( 10. , MAX( 0.03, sf_chl(1)%fnow(ji,jj,1) ) ) 194 irgb = NINT( 41 + 20.*LOG10(zchl) + 1.e-15 ) 195 zekb(ji,jj) = rkrgb(1,irgb) 196 zekg(ji,jj) = rkrgb(2,irgb) 197 zekr(ji,jj) = rkrgb(3,irgb) 198 END DO 199 END DO 200 ELSE ! Variable ocean volume but constant chrlorophyll 201 zchl = 0.05 ! constant chlorophyll 202 irgb = NINT( 41 + 20.*LOG10( zchl ) + 1.e-15 ) 203 zekb(:,:) = rkrgb(1,irgb) ! Separation in R-G-B depending of the chlorophyll 204 zekg(:,:) = rkrgb(2,irgb) 205 zekr(:,:) = rkrgb(3,irgb) 185 ! 186 zcoef = ( 1. - rn_abs ) / 3._wp !* surface equi-partition in R-G-B 187 DO jj = 2, jpjm1 188 DO ji = fs_2, fs_jpim1 189 ze0(ji,jj,1) = rn_abs * qsr(ji,jj) 190 ze1(ji,jj,1) = zcoef * qsr(ji,jj) 191 ze2(ji,jj,1) = zcoef * qsr(ji,jj) 192 ze3(ji,jj,1) = zcoef * qsr(ji,jj) 193 zea(ji,jj,1) = qsr(ji,jj) 194 END DO 195 END DO 196 ! 197 DO jk = 2, nksr+1 !* interior equi-partition in R-G-B 198 DO jj = 2, jpjm1 199 DO ji = fs_2, fs_jpim1 200 zc0 = ze0(ji,jj,jk-1) * EXP( - e3t_n(ji,jj,jk-1) * xsi0r ) 201 zc1 = ze1(ji,jj,jk-1) * EXP( - e3t_n(ji,jj,jk-1) * zekb(ji,jj) ) 202 zc2 = ze2(ji,jj,jk-1) * EXP( - e3t_n(ji,jj,jk-1) * zekg(ji,jj) ) 203 zc3 = ze3(ji,jj,jk-1) * EXP( - e3t_n(ji,jj,jk-1) * zekr(ji,jj) ) 204 ze0(ji,jj,jk) = zc0 205 ze1(ji,jj,jk) = zc1 206 ze2(ji,jj,jk) = zc2 207 ze3(ji,jj,jk) = zc3 208 zea(ji,jj,jk) = ( zc0 + zc1 + zc2 + zc3 ) * wmask(ji,jj,jk) 209 END DO 210 END DO 211 END DO 212 ! 213 DO jk = 1, nksr !* now qsr induced heat content 214 DO jj = 2, jpjm1 215 DO ji = fs_2, fs_jpim1 216 qsr_hc(ji,jj,jk) = r1_rau0_rcp * ( zea(ji,jj,jk) - zea(ji,jj,jk+1) ) 217 END DO 218 END DO 219 END DO 220 ! 221 CALL wrk_dealloc( jpi,jpj, zekb, zekg, zekr ) 222 CALL wrk_dealloc( jpi,jpj,jpk, ze0, ze1, ze2, ze3, zea ) 223 ! 224 CASE( np_2BD ) !== 2-bands fluxes ==! 225 ! 226 zz0 = rn_abs * r1_rau0_rcp ! surface equi-partition in 2-bands 227 zz1 = ( 1. - rn_abs ) * r1_rau0_rcp 228 DO jk = 1, nksr ! solar heat absorbed at T-point in the top 400m 229 DO jj = 2, jpjm1 230 DO ji = fs_2, fs_jpim1 231 zc0 = zz0 * EXP( -gdepw_n(ji,jj,jk )*xsi0r ) + zz1 * EXP( -gdepw_n(ji,jj,jk )*xsi1r ) 232 zc1 = zz0 * EXP( -gdepw_n(ji,jj,jk+1)*xsi0r ) + zz1 * EXP( -gdepw_n(ji,jj,jk+1)*xsi1r ) 233 qsr_hc(ji,jj,jk) = qsr(ji,jj) * ( zc0 * wmask(ji,jj,jk) - zc1 * wmask(ji,jj,jk+1) ) 234 END DO 235 END DO 236 END DO 237 ! 238 END SELECT 239 ! 240 ! !-----------------------------! 241 DO jk = 1, nksr ! update to the temp. trend ! 242 DO jj = 2, jpjm1 !-----------------------------! 243 DO ji = fs_2, fs_jpim1 ! vector opt. 244 tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem) & 245 & + z1_2 * ( qsr_hc_b(ji,jj,jk) + qsr_hc(ji,jj,jk) ) / e3t_n(ji,jj,jk) 246 END DO 247 END DO 248 END DO 249 ! 250 IF( ln_qsr_ice ) THEN ! sea-ice: store the 1st ocean level attenuation coefficient 251 DO jj = 2, jpjm1 252 DO ji = fs_2, fs_jpim1 ! vector opt. 253 IF( qsr(ji,jj) /= 0._wp ) THEN ; fraqsr_1lev(ji,jj) = qsr_hc(ji,jj,1) / ( r1_rau0_rcp * qsr(ji,jj) ) 254 ELSE ; fraqsr_1lev(ji,jj) = 1._wp 206 255 ENDIF 207 ! 208 zcoef = ( 1. - rn_abs ) / 3.e0 ! equi-partition in R-G-B 209 ze0(:,:,1) = rn_abs * qsr(:,:) 210 ze1(:,:,1) = zcoef * qsr(:,:) 211 ze2(:,:,1) = zcoef * qsr(:,:) 212 ze3(:,:,1) = zcoef * qsr(:,:) 213 zea(:,:,1) = qsr(:,:) 214 ! 215 DO jk = 2, nksr+1 216 DO jj = 1, jpj 217 DO ji = 1, jpi 218 zc0 = ze0(ji,jj,jk-1) * EXP( - fse3t(ji,jj,jk-1) * xsi0r ) 219 zc1 = ze1(ji,jj,jk-1) * EXP( - fse3t(ji,jj,jk-1) * zekb(ji,jj) ) 220 zc2 = ze2(ji,jj,jk-1) * EXP( - fse3t(ji,jj,jk-1) * zekg(ji,jj) ) 221 zc3 = ze3(ji,jj,jk-1) * EXP( - fse3t(ji,jj,jk-1) * zekr(ji,jj) ) 222 ze0(ji,jj,jk) = zc0 223 ze1(ji,jj,jk) = zc1 224 ze2(ji,jj,jk) = zc2 225 ze3(ji,jj,jk) = zc3 226 zea(ji,jj,jk) = ( zc0 + zc1 + zc2 + zc3 ) * tmask(ji,jj,jk) 227 END DO 228 END DO 229 END DO 230 ! clem: store attenuation coefficient of the first ocean level 231 IF ( ln_qsr_ice ) THEN 232 DO jj = 1, jpj 233 DO ji = 1, jpi 234 zzc0 = rn_abs * EXP( - fse3t(ji,jj,1) * xsi0r ) 235 zzc1 = zcoef * EXP( - fse3t(ji,jj,1) * zekb(ji,jj) ) 236 zzc2 = zcoef * EXP( - fse3t(ji,jj,1) * zekg(ji,jj) ) 237 zzc3 = zcoef * EXP( - fse3t(ji,jj,1) * zekr(ji,jj) ) 238 fraqsr_1lev(ji,jj) = 1.0 - ( zzc0 + zzc1 + zzc2 + zzc3 ) * tmask(ji,jj,2) 239 END DO 240 END DO 241 ENDIF 242 ! 243 DO jk = 1, nksr ! compute and add qsr trend to ta 244 qsr_hc(:,:,jk) = r1_rau0_rcp * ( zea(:,:,jk) - zea(:,:,jk+1) ) 245 END DO 246 zea(:,:,nksr+1:jpk) = 0.e0 ! below 400m set to zero 247 CALL iom_put( 'qsr3d', zea ) ! Shortwave Radiation 3D distribution 248 ! 249 ELSE !* Constant Chlorophyll 250 DO jk = 1, nksr 251 qsr_hc(:,:,jk) = etot3(:,:,jk) * qsr(:,:) 252 END DO 253 ! clem: store attenuation coefficient of the first ocean level 254 IF ( ln_qsr_ice ) THEN 255 fraqsr_1lev(:,:) = etot3(:,:,1) / r1_rau0_rcp 256 ENDIF 257 ENDIF 258 259 ENDIF 260 ! ! ------------------------- ! 261 IF( ln_qsr_2bd ) THEN ! 2 band light penetration ! 262 ! ! ------------------------- ! 263 ! 264 IF( lk_vvl ) THEN !* variable volume 265 zz0 = rn_abs * r1_rau0_rcp 266 zz1 = ( 1. - rn_abs ) * r1_rau0_rcp 267 DO jk = 1, nksr ! solar heat absorbed at T-point in the top 400m 268 DO jj = 1, jpj 269 DO ji = 1, jpi 270 zc0 = zz0 * EXP( -fsdepw(ji,jj,jk )*xsi0r ) + zz1 * EXP( -fsdepw(ji,jj,jk )*xsi1r ) 271 zc1 = zz0 * EXP( -fsdepw(ji,jj,jk+1)*xsi0r ) + zz1 * EXP( -fsdepw(ji,jj,jk+1)*xsi1r ) 272 qsr_hc(ji,jj,jk) = qsr(ji,jj) * ( zc0*tmask(ji,jj,jk) - zc1*tmask(ji,jj,jk+1) ) 273 END DO 274 END DO 275 END DO 276 ! clem: store attenuation coefficient of the first ocean level 277 IF ( ln_qsr_ice ) THEN 278 DO jj = 1, jpj 279 DO ji = 1, jpi 280 zc0 = zz0 * EXP( -fsdepw(ji,jj,1)*xsi0r ) + zz1 * EXP( -fsdepw(ji,jj,1)*xsi1r ) 281 zc1 = zz0 * EXP( -fsdepw(ji,jj,2)*xsi0r ) + zz1 * EXP( -fsdepw(ji,jj,2)*xsi1r ) 282 fraqsr_1lev(ji,jj) = ( zc0*tmask(ji,jj,1) - zc1*tmask(ji,jj,2) ) / r1_rau0_rcp 283 END DO 284 END DO 285 ENDIF 286 ELSE !* constant volume: coef. computed one for all 287 DO jk = 1, nksr 288 DO jj = 2, jpjm1 289 DO ji = fs_2, fs_jpim1 ! vector opt. 290 ! (ISF) no light penetration below the ice shelves 291 qsr_hc(ji,jj,jk) = etot3(ji,jj,jk) * qsr(ji,jj) * tmask(ji,jj,1) 292 END DO 293 END DO 294 END DO 295 ! clem: store attenuation coefficient of the first ocean level 296 IF ( ln_qsr_ice ) THEN 297 fraqsr_1lev(:,:) = etot3(:,:,1) / r1_rau0_rcp 298 ENDIF 299 ! 300 ENDIF 301 ! 302 ENDIF 303 ! 304 ! Add to the general trend 305 DO jk = 1, nksr 306 DO jj = 2, jpjm1 307 DO ji = fs_2, fs_jpim1 ! vector opt. 308 z1_e3t = zfact / fse3t(ji,jj,jk) 309 tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem) + ( qsr_hc_b(ji,jj,jk) + qsr_hc(ji,jj,jk) ) * z1_e3t 310 END DO 311 END DO 312 END DO 313 ! 314 ENDIF 315 ! 316 IF( lrst_oce ) THEN ! Write in the ocean restart file 317 ! ******************************* 318 IF(lwp) WRITE(numout,*) 319 IF(lwp) WRITE(numout,*) 'qsr tracer content forcing field written in ocean restart file ', & 320 & 'at it= ', kt,' date= ', ndastp 321 IF(lwp) WRITE(numout,*) '~~~~' 256 END DO 257 END DO 258 ! Update haloes since lim_thd needs fraqsr_1lev to be defined everywhere 259 CALL lbc_lnk( fraqsr_1lev(:,:), 'T', 1._wp ) 260 ENDIF 261 ! 262 IF( iom_use('qsr3d') ) THEN ! output the shortwave Radiation distribution 263 CALL wrk_alloc( jpi,jpj,jpk, zetot ) 264 ! 265 zetot(:,:,nksr+1:jpk) = 0._wp ! below ~400m set to zero 266 DO jk = nksr, 1, -1 267 zetot(:,:,jk) = zetot(:,:,jk+1) + qsr_hc(:,:,jk) / r1_rau0_rcp 268 END DO 269 CALL iom_put( 'qsr3d', zetot ) ! 3D distribution of shortwave Radiation 270 ! 271 CALL wrk_dealloc( jpi,jpj,jpk, zetot ) 272 ENDIF 273 ! 274 IF( lrst_oce ) THEN ! write in the ocean restart file 322 275 CALL iom_rstput( kt, nitrst, numrow, 'qsr_hc_b' , qsr_hc ) 323 CALL iom_rstput( kt, nitrst, numrow, 'fraqsr_1lev', fraqsr_1lev ) ! default definition in sbcssm 324 ! 325 ENDIF 326 276 CALL iom_rstput( kt, nitrst, numrow, 'fraqsr_1lev', fraqsr_1lev ) 277 ENDIF 278 ! 327 279 IF( l_trdtra ) THEN ! qsr tracers trends saved for diagnostics 328 280 ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:) 329 281 CALL trd_tra( kt, 'TRA', jp_tem, jptra_qsr, ztrdt ) 330 CALL wrk_dealloc( jpi, jpj, jpk,ztrdt )282 CALL wrk_dealloc( jpi,jpj,jpk, ztrdt ) 331 283 ENDIF 332 284 ! ! print mean trends (used for debugging) 333 285 IF(ln_ctl) CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' qsr - Ta: ', mask1=tmask, clinfo3='tra-ta' ) 334 !335 CALL wrk_dealloc( jpi, jpj, zekb, zekg, zekr )336 CALL wrk_dealloc( jpi, jpj, jpk, ze0, ze1, ze2, ze3, zea )337 286 ! 338 287 IF( nn_timing == 1 ) CALL timing_stop('tra_qsr') … … 358 307 !! Reference : Jerlov, N. G., 1968 Optical Oceanography, Elsevier, 194pp. 359 308 !!---------------------------------------------------------------------- 360 ! 361 INTEGER :: ji, jj, jk ! dummy loop indices 362 INTEGER :: irgb, ierror, ioptio, nqsr ! local integer 363 INTEGER :: ios ! Local integer output status for namelist read 364 REAL(wp) :: zz0, zc0 , zc1, zcoef ! local scalars 365 REAL(wp) :: zz1, zc2 , zc3, zchl ! - - 366 REAL(wp), POINTER, DIMENSION(:,: ) :: zekb, zekg, zekr 367 REAL(wp), POINTER, DIMENSION(:,:,:) :: ze0, ze1, ze2, ze3, zea 309 INTEGER :: ji, jj, jk ! dummy loop indices 310 INTEGER :: ios, irgb, ierror, ioptio ! local integer 311 REAL(wp) :: zz0, zc0 , zc1, zcoef ! local scalars 312 REAL(wp) :: zz1, zc2 , zc3, zchl ! - - 368 313 ! 369 314 CHARACTER(len=100) :: cn_dir ! Root directory for location of ssr files 370 315 TYPE(FLD_N) :: sn_chl ! informations about the chlorofyl field to be read 371 316 !! 372 NAMELIST/namtra_qsr/ sn_chl, cn_dir, ln_ traqsr, ln_qsr_rgb, ln_qsr_2bd, ln_qsr_bio, ln_qsr_ice, &317 NAMELIST/namtra_qsr/ sn_chl, cn_dir, ln_qsr_rgb, ln_qsr_2bd, ln_qsr_bio, ln_qsr_ice, & 373 318 & nn_chldta, rn_abs, rn_si0, rn_si1 374 319 !!---------------------------------------------------------------------- 375 376 ! 377 IF( nn_timing == 1 ) CALL timing_start('tra_qsr_init') 378 ! 379 CALL wrk_alloc( jpi, jpj, zekb, zekg, zekr ) 380 CALL wrk_alloc( jpi, jpj, jpk, ze0, ze1, ze2, ze3, zea ) 381 ! 382 383 REWIND( numnam_ref ) ! Namelist namtra_qsr in reference namelist : Ratio and length of penetration 320 ! 321 IF( nn_timing == 1 ) CALL timing_start('tra_qsr_init') 322 ! 323 REWIND( numnam_ref ) ! Namelist namtra_qsr in reference namelist 384 324 READ ( numnam_ref, namtra_qsr, IOSTAT = ios, ERR = 901) 385 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtra_qsr in reference namelist', lwp )386 387 REWIND( numnam_cfg ) ! Namelist namtra_qsr in configuration namelist : Ratio and length of penetration325 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtra_qsr in reference namelist', lwp ) 326 ! 327 REWIND( numnam_cfg ) ! Namelist namtra_qsr in configuration namelist 388 328 READ ( numnam_cfg, namtra_qsr, IOSTAT = ios, ERR = 902 ) 389 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtra_qsr in configuration namelist', lwp )329 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtra_qsr in configuration namelist', lwp ) 390 330 IF(lwm) WRITE ( numond, namtra_qsr ) 391 331 ! … … 395 335 WRITE(numout,*) '~~~~~~~~~~~~' 396 336 WRITE(numout,*) ' Namelist namtra_qsr : set the parameter of penetration' 397 WRITE(numout,*) ' Light penetration (T) or not (F) ln_traqsr = ', ln_traqsr 398 WRITE(numout,*) ' RGB (Red-Green-Blue) light penetration ln_qsr_rgb = ', ln_qsr_rgb 399 WRITE(numout,*) ' 2 band light penetration ln_qsr_2bd = ', ln_qsr_2bd 400 WRITE(numout,*) ' bio-model light penetration ln_qsr_bio = ', ln_qsr_bio 401 WRITE(numout,*) ' light penetration for ice-model LIM3 ln_qsr_ice = ', ln_qsr_ice 402 WRITE(numout,*) ' RGB : Chl data (=1) or cst value (=0) nn_chldta = ', nn_chldta 403 WRITE(numout,*) ' RGB & 2 bands: fraction of light (rn_si1) rn_abs = ', rn_abs 404 WRITE(numout,*) ' RGB & 2 bands: shortess depth of extinction rn_si0 = ', rn_si0 405 WRITE(numout,*) ' 2 bands: longest depth of extinction rn_si1 = ', rn_si1 406 ENDIF 407 408 IF( ln_traqsr ) THEN ! control consistency 409 ! 410 IF( .NOT.lk_qsr_bio .AND. ln_qsr_bio ) THEN 411 CALL ctl_warn( 'No bio model : force ln_qsr_bio = FALSE ' ) 412 ln_qsr_bio = .FALSE. 337 WRITE(numout,*) ' RGB (Red-Green-Blue) light penetration ln_qsr_rgb = ', ln_qsr_rgb 338 WRITE(numout,*) ' 2 band light penetration ln_qsr_2bd = ', ln_qsr_2bd 339 WRITE(numout,*) ' bio-model light penetration ln_qsr_bio = ', ln_qsr_bio 340 WRITE(numout,*) ' light penetration for ice-model (LIM3) ln_qsr_ice = ', ln_qsr_ice 341 WRITE(numout,*) ' RGB : Chl data (=1) or cst value (=0) nn_chldta = ', nn_chldta 342 WRITE(numout,*) ' RGB & 2 bands: fraction of light (rn_si1) rn_abs = ', rn_abs 343 WRITE(numout,*) ' RGB & 2 bands: shortess depth of extinction rn_si0 = ', rn_si0 344 WRITE(numout,*) ' 2 bands: longest depth of extinction rn_si1 = ', rn_si1 345 WRITE(numout,*) 346 ENDIF 347 ! 348 ioptio = 0 ! Parameter control 349 IF( ln_qsr_rgb ) ioptio = ioptio + 1 350 IF( ln_qsr_2bd ) ioptio = ioptio + 1 351 IF( ln_qsr_bio ) ioptio = ioptio + 1 352 ! 353 IF( ioptio /= 1 ) CALL ctl_stop( 'Choose ONE type of light penetration in namelist namtra_qsr', & 354 & ' 2 bands, 3 RGB bands or bio-model light penetration' ) 355 ! 356 IF( ln_qsr_rgb .AND. nn_chldta == 0 ) nqsr = np_RGB 357 IF( ln_qsr_rgb .AND. nn_chldta == 1 ) nqsr = np_RGBc 358 IF( ln_qsr_2bd ) nqsr = np_2BD 359 IF( ln_qsr_bio ) nqsr = np_BIO 360 ! 361 ! ! Initialisation 362 xsi0r = 1._wp / rn_si0 363 xsi1r = 1._wp / rn_si1 364 ! 365 SELECT CASE( nqsr ) 366 ! 367 CASE( np_RGB , np_RGBc ) !== Red-Green-Blue light penetration ==! 368 ! 369 IF(lwp) WRITE(numout,*) ' R-G-B light penetration ' 370 ! 371 CALL trc_oce_rgb( rkrgb ) ! tabulated attenuation coef. 372 ! 373 nksr = trc_oce_ext_lev( r_si2, 33._wp ) ! level of light extinction 374 ! 375 IF(lwp) WRITE(numout,*) ' level of light extinction = ', nksr, ' ref depth = ', gdepw_1d(nksr+1), ' m' 376 ! 377 IF( nqsr == np_RGBc ) THEN ! Chl data : set sf_chl structure 378 IF(lwp) WRITE(numout,*) ' Chlorophyll read in a file' 379 ALLOCATE( sf_chl(1), STAT=ierror ) 380 IF( ierror > 0 ) THEN 381 CALL ctl_stop( 'tra_qsr_init: unable to allocate sf_chl structure' ) ; RETURN 382 ENDIF 383 ALLOCATE( sf_chl(1)%fnow(jpi,jpj,1) ) 384 IF( sn_chl%ln_tint ) ALLOCATE( sf_chl(1)%fdta(jpi,jpj,1,2) ) 385 ! ! fill sf_chl with sn_chl and control print 386 CALL fld_fill( sf_chl, (/ sn_chl /), cn_dir, 'tra_qsr_init', & 387 & 'Solar penetration function of read chlorophyll', 'namtra_qsr' ) 413 388 ENDIF 414 ! 415 ioptio = 0 ! Parameter control 416 IF( ln_qsr_rgb ) ioptio = ioptio + 1 417 IF( ln_qsr_2bd ) ioptio = ioptio + 1 418 IF( ln_qsr_bio ) ioptio = ioptio + 1 419 ! 420 IF( ioptio /= 1 ) & 421 CALL ctl_stop( ' Choose ONE type of light penetration in namelist namtra_qsr', & 422 & ' 2 bands, 3 RGB bands or bio-model light penetration' ) 423 ! 424 IF( ln_qsr_rgb .AND. nn_chldta == 0 ) nqsr = 1 425 IF( ln_qsr_rgb .AND. nn_chldta == 1 ) nqsr = 2 426 IF( ln_qsr_2bd ) nqsr = 3 427 IF( ln_qsr_bio ) nqsr = 4 428 ! 429 IF(lwp) THEN ! Print the choice 430 WRITE(numout,*) 431 IF( nqsr == 1 ) WRITE(numout,*) ' R-G-B light penetration - Constant Chlorophyll' 432 IF( nqsr == 2 ) WRITE(numout,*) ' R-G-B light penetration - Chl data ' 433 IF( nqsr == 3 ) WRITE(numout,*) ' 2 bands light penetration' 434 IF( nqsr == 4 ) WRITE(numout,*) ' bio-model light penetration' 389 IF( nqsr == np_RGB ) THEN ! constant Chl 390 IF(lwp) WRITE(numout,*) ' Constant Chlorophyll concentration = 0.05' 435 391 ENDIF 436 392 ! 437 ENDIF 438 ! ! ===================================== ! 439 IF( ln_traqsr ) THEN ! Initialisation of Light Penetration ! 440 ! ! ===================================== ! 441 ! 442 xsi0r = 1.e0 / rn_si0 443 xsi1r = 1.e0 / rn_si1 444 ! ! ---------------------------------- ! 445 IF( ln_qsr_rgb ) THEN ! Red-Green-Blue light penetration ! 446 ! ! ---------------------------------- ! 447 ! 448 CALL trc_oce_rgb( rkrgb ) !* tabulated attenuation coef. 449 ! 450 ! !* level of light extinction 451 IF( ln_sco ) THEN ; nksr = jpkm1 452 ELSE ; nksr = trc_oce_ext_lev( r_si2, 0.33e2 ) 453 ENDIF 454 455 IF(lwp) WRITE(numout,*) ' level of light extinction = ', nksr, ' ref depth = ', gdepw_1d(nksr+1), ' m' 456 ! 457 IF( nn_chldta == 1 ) THEN !* Chl data : set sf_chl structure 458 IF(lwp) WRITE(numout,*) 459 IF(lwp) WRITE(numout,*) ' Chlorophyll read in a file' 460 ALLOCATE( sf_chl(1), STAT=ierror ) 461 IF( ierror > 0 ) THEN 462 CALL ctl_stop( 'tra_qsr_init: unable to allocate sf_chl structure' ) ; RETURN 463 ENDIF 464 ALLOCATE( sf_chl(1)%fnow(jpi,jpj,1) ) 465 IF( sn_chl%ln_tint )ALLOCATE( sf_chl(1)%fdta(jpi,jpj,1,2) ) 466 ! ! fill sf_chl with sn_chl and control print 467 CALL fld_fill( sf_chl, (/ sn_chl /), cn_dir, 'tra_qsr_init', & 468 & 'Solar penetration function of read chlorophyll', 'namtra_qsr' ) 469 ! 470 ELSE !* constant Chl : compute once for all the distribution of light (etot3) 471 IF(lwp) WRITE(numout,*) 472 IF(lwp) WRITE(numout,*) ' Constant Chlorophyll concentration = 0.05' 473 IF( lk_vvl ) THEN ! variable volume 474 IF(lwp) WRITE(numout,*) ' key_vvl: light distribution will be computed at each time step' 475 ELSE ! constant volume: computes one for all 476 IF(lwp) WRITE(numout,*) ' fixed volume: light distribution computed one for all' 477 ! 478 zchl = 0.05 ! constant chlorophyll 479 irgb = NINT( 41 + 20.*LOG10(zchl) + 1.e-15 ) 480 zekb(:,:) = rkrgb(1,irgb) ! Separation in R-G-B depending of the chlorophyll 481 zekg(:,:) = rkrgb(2,irgb) 482 zekr(:,:) = rkrgb(3,irgb) 483 ! 484 zcoef = ( 1. - rn_abs ) / 3.e0 ! equi-partition in R-G-B 485 ze0(:,:,1) = rn_abs 486 ze1(:,:,1) = zcoef 487 ze2(:,:,1) = zcoef 488 ze3(:,:,1) = zcoef 489 zea(:,:,1) = tmask(:,:,1) ! = ( ze0+ze1+z2+ze3 ) * tmask 490 491 DO jk = 2, nksr+1 492 DO jj = 1, jpj 493 DO ji = 1, jpi 494 zc0 = ze0(ji,jj,jk-1) * EXP( - e3t_0(ji,jj,jk-1) * xsi0r ) 495 zc1 = ze1(ji,jj,jk-1) * EXP( - e3t_0(ji,jj,jk-1) * zekb(ji,jj) ) 496 zc2 = ze2(ji,jj,jk-1) * EXP( - e3t_0(ji,jj,jk-1) * zekg(ji,jj) ) 497 zc3 = ze3(ji,jj,jk-1) * EXP( - e3t_0(ji,jj,jk-1) * zekr(ji,jj) ) 498 ze0(ji,jj,jk) = zc0 499 ze1(ji,jj,jk) = zc1 500 ze2(ji,jj,jk) = zc2 501 ze3(ji,jj,jk) = zc3 502 zea(ji,jj,jk) = ( zc0 + zc1 + zc2 + zc3 ) * tmask(ji,jj,jk) 503 END DO 504 END DO 505 END DO 506 ! 507 DO jk = 1, nksr 508 ! (ISF) no light penetration below the ice shelves 509 etot3(:,:,jk) = r1_rau0_rcp * ( zea(:,:,jk) - zea(:,:,jk+1) ) * tmask(:,:,1) 510 END DO 511 etot3(:,:,nksr+1:jpk) = 0.e0 ! below 400m set to zero 512 ENDIF 513 ENDIF 514 ! 515 ENDIF 516 ! ! ---------------------------------- ! 517 IF( ln_qsr_2bd ) THEN ! 2 bands light penetration ! 518 ! ! ---------------------------------- ! 519 ! 520 ! ! level of light extinction 521 nksr = trc_oce_ext_lev( rn_si1, 1.e2 ) 522 IF(lwp) THEN 523 WRITE(numout,*) 524 IF(lwp) WRITE(numout,*) ' level of light extinction = ', nksr, ' ref depth = ', gdepw_1d(nksr+1), ' m' 525 ENDIF 526 ! 527 IF( lk_vvl ) THEN ! variable volume 528 IF(lwp) WRITE(numout,*) ' key_vvl: light distribution will be computed at each time step' 529 ELSE ! constant volume: computes one for all 530 zz0 = rn_abs * r1_rau0_rcp 531 zz1 = ( 1. - rn_abs ) * r1_rau0_rcp 532 DO jk = 1, nksr !* solar heat absorbed at T-point computed once for all 533 DO jj = 1, jpj ! top 400 meters 534 DO ji = 1, jpi 535 zc0 = zz0 * EXP( -fsdepw(ji,jj,jk )*xsi0r ) + zz1 * EXP( -fsdepw(ji,jj,jk )*xsi1r ) 536 zc1 = zz0 * EXP( -fsdepw(ji,jj,jk+1)*xsi0r ) + zz1 * EXP( -fsdepw(ji,jj,jk+1)*xsi1r ) 537 etot3(ji,jj,jk) = ( zc0 * tmask(ji,jj,jk) - zc1 * tmask(ji,jj,jk+1) ) * tmask(ji,jj,1) 538 END DO 539 END DO 540 END DO 541 etot3(:,:,nksr+1:jpk) = 0.e0 ! below 400m set to zero 542 ! 543 ENDIF 544 ENDIF 545 ! ! ===================================== ! 546 ELSE ! No light penetration ! 547 ! ! ===================================== ! 548 IF(lwp) THEN 549 WRITE(numout,*) 550 WRITE(numout,*) 'tra_qsr_init : NO solar flux penetration' 551 WRITE(numout,*) '~~~~~~~~~~~~' 552 ENDIF 553 ENDIF 554 ! 555 ! initialisation of fraqsr_1lev used in sbcssm 393 CASE( np_2BD ) !== 2 bands light penetration ==! 394 ! 395 IF(lwp) WRITE(numout,*) ' 2 bands light penetration' 396 ! 397 nksr = trc_oce_ext_lev( rn_si1, 100._wp ) ! level of light extinction 398 IF(lwp) WRITE(numout,*) ' level of light extinction = ', nksr, ' ref depth = ', gdepw_1d(nksr+1), ' m' 399 ! 400 CASE( np_BIO ) !== BIO light penetration ==! 401 ! 402 IF(lwp) WRITE(numout,*) ' bio-model light penetration' 403 IF( .NOT.lk_qsr_bio ) CALL ctl_stop( 'No bio model : ln_qsr_bio = true impossible ' ) 404 ! 405 END SELECT 406 ! 407 qsr_hc(:,:,:) = 0._wp ! now qsr heat content set to zero where it will not be computed 408 ! 409 ! 1st ocean level attenuation coefficient (used in sbcssm) 556 410 IF( iom_varid( numror, 'fraqsr_1lev', ldstop = .FALSE. ) > 0 ) THEN 557 411 CALL iom_get( numror, jpdom_autoglo, 'fraqsr_1lev' , fraqsr_1lev ) 558 412 ELSE 559 fraqsr_1lev(:,:) = 1._wp ! default definition 560 ENDIF 561 ! 562 CALL wrk_dealloc( jpi, jpj, zekb, zekg, zekr ) 563 CALL wrk_dealloc( jpi, jpj, jpk, ze0, ze1, ze2, ze3, zea ) 564 ! 565 IF( nn_timing == 1 ) CALL timing_stop('tra_qsr_init') 413 fraqsr_1lev(:,:) = 1._wp ! default : no penetration 414 ENDIF 415 ! 416 IF( nn_timing == 1 ) CALL timing_stop('tra_qsr_init') 566 417 ! 567 418 END SUBROUTINE tra_qsr_init -
trunk/NEMOGCM/NEMO/OPA_SRC/TRA/trasbc.F90
r5643 r6140 13 13 14 14 !!---------------------------------------------------------------------- 15 !! tra_sbc : update the tracer trend at ocean surface 16 !!---------------------------------------------------------------------- 17 USE oce ! ocean dynamics and active tracers 18 USE sbc_oce ! surface boundary condition: ocean 19 USE dom_oce ! ocean space domain variables 20 USE phycst ! physical constant 21 USE sbcmod ! ln_rnf 22 USE sbcrnf ! River runoff 23 USE sbcisf ! Ice shelf 24 USE traqsr ! solar radiation penetration 25 USE trd_oce ! trends: ocean variables 26 USE trdtra ! trends manager: tracers 15 !! tra_sbc : update the tracer trend at ocean surface 16 !!---------------------------------------------------------------------- 17 USE oce ! ocean dynamics and active tracers 18 USE sbc_oce ! surface boundary condition: ocean 19 USE dom_oce ! ocean space domain variables 20 USE phycst ! physical constant 21 USE eosbn2 ! Equation Of State 22 USE sbcmod ! ln_rnf 23 USE sbcrnf ! River runoff 24 USE sbcisf ! Ice shelf 25 USE iscplini ! Ice sheet coupling 26 USE traqsr ! solar radiation penetration 27 USE trd_oce ! trends: ocean variables 28 USE trdtra ! trends manager: tracers 27 29 ! 28 USE in_out_manager ! I/O manager 29 USE prtctl ! Print control 30 USE iom 31 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 32 USE wrk_nemo ! Memory Allocation 33 USE timing ! Timing 34 USE eosbn2 30 USE in_out_manager ! I/O manager 31 USE prtctl ! Print control 32 USE iom ! xIOS server 33 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 34 USE wrk_nemo ! Memory Allocation 35 USE timing ! Timing 35 36 36 37 IMPLICIT NONE 37 38 PRIVATE 38 39 39 PUBLIC tra_sbc 40 PUBLIC tra_sbc ! routine called by step.F90 40 41 41 42 !! * Substitutions 42 # include "domzgr_substitute.h90"43 43 # include "vectopt_loop_substitute.h90" 44 44 !!---------------------------------------------------------------------- … … 57 57 !! and add it to the general trend of tracer equations. 58 58 !! 59 !! ** Method : 60 !! Following Roullet and Madec (2000), the air-sea flux can be divided 61 !! into three effects: (1) Fext, external forcing; 62 !! (2) Fwi, concentration/dilution effect due to water exchanged 63 !! at the surface by evaporation, precipitations and runoff (E-P-R); 64 !! (3) Fwe, tracer carried with the water that is exchanged. 65 !! - salinity : salt flux only due to freezing/melting 66 !! sa = sa + sfx / rau0 / e3t for k=1 59 !! ** Method : The (air+ice)-sea flux has two components: 60 !! (1) Fext, external forcing (i.e. flux through the (air+ice)-sea interface); 61 !! (2) Fwe , tracer carried with the water that is exchanged with air+ice. 62 !! The input forcing fields (emp, rnf, sfx, isf) contain Fext+Fwe, 63 !! they are simply added to the tracer trend (tsa). 64 !! In linear free surface case (ln_linssh=T), the volume of the 65 !! ocean does not change with the water exchanges at the (air+ice)-sea 66 !! interface. Therefore another term has to be added, to mimic the 67 !! concentration/dilution effect associated with water exchanges. 67 68 !! 68 !! Fext, flux through the air-sea interface for temperature and salt: 69 !! - temperature : heat flux q (w/m2). If penetrative solar 70 !! radiation q is only the non solar part of the heat flux, the 71 !! solar part is added in traqsr.F routine. 72 !! ta = ta + q /(rau0 rcp e3t) for k=1 73 !! - salinity : no salt flux 74 !! 75 !! The formulation for Fwb and Fwi vary according to the free 76 !! surface formulation (linear or variable volume). 77 !! * Linear free surface 78 !! The surface freshwater flux modifies the ocean volume 79 !! and thus the concentration of a tracer and the temperature. 80 !! First order of the effect of surface freshwater exchange 81 !! for salinity, it can be neglected on temperature (especially 82 !! as the temperature of precipitations and runoffs is usually 83 !! unknown). 84 !! - temperature : we assume that the temperature of both 85 !! precipitations and runoffs is equal to the SST, thus there 86 !! is no additional flux since in this case, the concentration 87 !! dilution effect is balanced by the net heat flux associated 88 !! to the freshwater exchange (Fwe+Fwi=0): 89 !! (Tp P - Te E) + SST (P-E) = 0 when Tp=Te=SST 90 !! - salinity : evaporation, precipitation and runoff 91 !! water has a zero salinity but there is a salt flux due to 92 !! freezing/melting, thus: 93 !! sa = sa + emp * sn / rau0 / e3t for k=1 94 !! + sfx / rau0 / e3t 95 !! where emp, the surface freshwater budget (evaporation minus 96 !! precipitation minus runoff) given in kg/m2/s is divided 97 !! by rau0 (density of sea water) to obtain m/s. 98 !! Note: even though Fwe does not appear explicitly for 99 !! temperature in this routine, the heat carried by the water 100 !! exchanged through the surface is part of the total heat flux 101 !! forcing and must be taken into account in the global heat 102 !! balance). 103 !! * nonlinear free surface (variable volume, lk_vvl) 104 !! contrary to the linear free surface case, Fwi is properly 105 !! taken into account by using the true layer thicknesses to 106 !! calculate tracer content and advection. There is no need to 107 !! deal with it in this routine. 108 !! - temperature: Fwe=SST (P-E+R) is added to Fext. 109 !! - salinity: Fwe = 0, there is no surface flux of salt. 110 !! 111 !! ** Action : - Update the 1st level of (ta,sa) with the trend associated 112 !! with the tracer surface boundary condition 113 !! - send trends to trdtra module (l_trdtra=T) 69 !! ** Action : - Update tsa with the surface boundary condition trend 70 !! - send trends to trdtra module for further diagnostics(l_trdtra=T) 114 71 !!---------------------------------------------------------------------- 115 72 INTEGER, INTENT(in) :: kt ! ocean time-step index 116 !! 117 INTEGER :: ji, jj, jk, jn ! dummy loop indices 118 INTEGER :: ikt, ikb 119 INTEGER :: nk_isf 120 REAL(wp) :: zfact, z1_e3t, zdep 121 REAL(wp) :: zalpha, zhk 73 ! 74 INTEGER :: ji, jj, jk, jn ! dummy loop indices 75 INTEGER :: ikt, ikb ! local integers 76 REAL(wp) :: zfact, z1_e3t, zdep ! local scalar 122 77 REAL(wp), POINTER, DIMENSION(:,:,:) :: ztrdt, ztrds 123 78 !!---------------------------------------------------------------------- … … 130 85 IF(lwp) WRITE(numout,*) '~~~~~~~ ' 131 86 ENDIF 132 87 ! 133 88 IF( l_trdtra ) THEN !* Save ta and sa trends 134 89 CALL wrk_alloc( jpi, jpj, jpk, ztrdt, ztrds ) … … 136 91 ztrds(:,:,:) = tsa(:,:,:,jp_sal) 137 92 ENDIF 138 139 !!gm IF( .NOT.ln_traqsr ) qsr(:,:) = 0.e0 ! no solar radiation penetration93 ! 94 !!gm This should be moved into sbcmod.F90 module ? (especially now that ln_traqsr is read in namsbc namelist) 140 95 IF( .NOT.ln_traqsr ) THEN ! no solar radiation penetration 141 96 qns(:,:) = qns(:,:) + qsr(:,:) ! total heat flux in qns 142 qsr(:,:) = 0. e0! qsr set to zero97 qsr(:,:) = 0._wp ! qsr set to zero 143 98 ENDIF 144 99 … … 146 101 ! EMP, SFX and QNS effects 147 102 !---------------------------------------- 148 ! Set before sbc tracer content fields 149 ! ************************************ 150 IF( kt == nit000 ) THEN ! Set the forcing field at nit000 - 1 151 ! ! ----------------------------------- 152 IF( ln_rstart .AND. & ! Restart: read in restart file 103 ! !== Set before sbc tracer content fields ==! 104 IF( kt == nit000 ) THEN !* 1st time-step 105 IF( ln_rstart .AND. & ! Restart: read in restart file 153 106 & iom_varid( numror, 'sbc_hc_b', ldstop = .FALSE. ) > 0 ) THEN 154 IF(lwp) WRITE(numout,*) ' nit000-1 s urface tracer content forcing fields red in the restart file'107 IF(lwp) WRITE(numout,*) ' nit000-1 sbc tracer content field read in the restart file' 155 108 zfact = 0.5_wp 156 109 CALL iom_get( numror, jpdom_autoglo, 'sbc_hc_b', sbc_tsc_b(:,:,jp_tem) ) ! before heat content sbc trend 157 110 CALL iom_get( numror, jpdom_autoglo, 'sbc_sc_b', sbc_tsc_b(:,:,jp_sal) ) ! before salt content sbc trend 158 ELSE 111 ELSE ! No restart or restart not found: Euler forward time stepping 159 112 zfact = 1._wp 160 113 sbc_tsc_b(:,:,:) = 0._wp 161 114 ENDIF 162 ELSE ! Swap of forcing fields 163 ! ! ---------------------- 115 ELSE !* other time-steps: swap of forcing fields 164 116 zfact = 0.5_wp 165 117 sbc_tsc_b(:,:,:) = sbc_tsc(:,:,:) 166 118 ENDIF 167 ! Compute now sbc tracer content fields 168 ! ************************************* 169 170 ! Concentration dilution effect on (t,s) due to 171 ! evaporation, precipitation and qns, but not river runoff 172 173 IF( lk_vvl ) THEN ! Variable Volume case ==>> heat content of mass flux is in qns 174 DO jj = 1, jpj 175 DO ji = 1, jpi 176 sbc_tsc(ji,jj,jp_tem) = r1_rau0_rcp * qns(ji,jj) ! non solar heat flux 177 sbc_tsc(ji,jj,jp_sal) = r1_rau0 * sfx(ji,jj) ! salt flux due to freezing/melting 119 ! !== Now sbc tracer content fields ==! 120 DO jj = 2, jpj 121 DO ji = fs_2, fs_jpim1 ! vector opt. 122 sbc_tsc(ji,jj,jp_tem) = r1_rau0_rcp * qns(ji,jj) ! non solar heat flux 123 sbc_tsc(ji,jj,jp_sal) = r1_rau0 * sfx(ji,jj) ! salt flux due to freezing/melting 124 END DO 125 END DO 126 IF( ln_linssh ) THEN !* linear free surface 127 DO jj = 2, jpj !==>> add concentration/dilution effect due to constant volume cell 128 DO ji = fs_2, fs_jpim1 ! vector opt. 129 sbc_tsc(ji,jj,jp_tem) = sbc_tsc(ji,jj,jp_tem) + r1_rau0 * emp(ji,jj) * tsn(ji,jj,1,jp_tem) 130 sbc_tsc(ji,jj,jp_sal) = sbc_tsc(ji,jj,jp_sal) + r1_rau0 * emp(ji,jj) * tsn(ji,jj,1,jp_sal) 178 131 END DO 179 END DO 180 ELSE ! Constant Volume case ==>> Concentration dilution effect 132 END DO !==>> output c./d. term 133 IF( iom_use('emp_x_sst') ) CALL iom_put( "emp_x_sst", emp (:,:) * tsn(:,:,1,jp_tem) ) 134 IF( iom_use('emp_x_sss') ) CALL iom_put( "emp_x_sss", emp (:,:) * tsn(:,:,1,jp_sal) ) 135 ENDIF 136 ! 137 DO jn = 1, jpts !== update tracer trend ==! 181 138 DO jj = 2, jpj 182 DO ji = fs_2, fs_jpim1 ! vector opt. 183 ! temperature : heat flux 184 sbc_tsc(ji,jj,jp_tem) = r1_rau0_rcp * qns(ji,jj) & ! non solar heat flux 185 & + r1_rau0 * emp(ji,jj) * tsn(ji,jj,1,jp_tem) ! concent./dilut. effect 186 ! salinity : salt flux + concent./dilut. effect (both in sfx) 187 sbc_tsc(ji,jj,jp_sal) = r1_rau0 * ( sfx(ji,jj) & ! salt flux (freezing/melting) 188 & + emp(ji,jj) * tsn(ji,jj,1,jp_sal) ) ! concent./dilut. effect 139 DO ji = fs_2, fs_jpim1 ! vector opt. 140 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) 189 141 END DO 190 142 END DO 191 IF( iom_use('emp_x_sst') ) CALL iom_put( "emp_x_sst", emp (:,:) * tsn(:,:,1,jp_tem) ) ! c/d term on sst192 IF( iom_use('emp_x_sss') ) CALL iom_put( "emp_x_sss", emp (:,:) * tsn(:,:,1,jp_sal) ) ! c/d term on sss193 ENDIF194 ! Concentration dilution effect on (t,s) due to evapouration, precipitation and qns, but not river runoff195 DO jn = 1, jpts196 DO jj = 2, jpj197 DO ji = fs_2, fs_jpim1 ! vector opt.198 z1_e3t = zfact / fse3t(ji,jj,1)199 tsa(ji,jj,1,jn) = tsa(ji,jj,1,jn) + ( sbc_tsc_b(ji,jj,jn) + sbc_tsc(ji,jj,jn) ) * z1_e3t200 END DO201 END DO202 143 END DO 203 ! Write in the ocean restart file 204 ! ******************************* 205 IF( lrst_oce ) THEN 206 IF(lwp) WRITE(numout,*) 207 IF(lwp) WRITE(numout,*) 'sbc : ocean surface tracer content forcing fields written in ocean restart file ', & 208 & 'at it= ', kt,' date= ', ndastp 209 IF(lwp) WRITE(numout,*) '~~~~' 144 ! 145 IF( lrst_oce ) THEN !== write sbc_tsc in the ocean restart file ==! 210 146 CALL iom_rstput( kt, nitrst, numrow, 'sbc_hc_b', sbc_tsc(:,:,jp_tem) ) 211 147 CALL iom_rstput( kt, nitrst, numrow, 'sbc_sc_b', sbc_tsc(:,:,jp_sal) ) 212 148 ENDIF 213 149 ! 214 !215 150 !---------------------------------------- 216 151 ! Ice Shelf effects (ISF) … … 218 153 !---------------------------------------- 219 154 ! 220 IF( nn_isf > 0 ) THEN 221 zfact = 0.5e0 155 !!gm BUG ? Why no differences between non-linear and linear free surface ? 156 !!gm probably taken into account in r1_hisf_tbl : to be verified 157 IF( ln_isf ) THEN 158 zfact = 0.5_wp 222 159 DO jj = 2, jpj 223 160 DO ji = fs_2, fs_jpim1 224 161 ! 225 162 ikt = misfkt(ji,jj) 226 163 ikb = misfkb(ji,jj) 227 164 ! 228 165 ! level fully include in the ice shelf boundary layer 229 ! if isfdiv, we have to remove heat flux due to inflow at 0oC (as in rnf when you add rnf at sst)230 166 ! sign - because fwf sign of evapo (rnf sign of precip) 231 167 DO jk = ikt, ikb - 1 232 ! compute tfreez for the temperature correction (we add water at freezing temperature)233 168 ! compute trend 234 tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem) & 235 & + zfact * (risf_tsc_b(ji,jj,jp_tem) + risf_tsc(ji,jj,jp_tem)) * r1_hisf_tbl(ji,jj) 236 tsa(ji,jj,jk,jp_sal) = tsa(ji,jj,jk,jp_sal) & 237 & + zfact * (risf_tsc_b(ji,jj,jp_sal) + risf_tsc(ji,jj,jp_sal)) * r1_hisf_tbl(ji,jj) 169 tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem) & 170 & + zfact * ( risf_tsc_b(ji,jj,jp_tem) + risf_tsc(ji,jj,jp_tem) ) & 171 & * r1_hisf_tbl(ji,jj) 238 172 END DO 239 173 240 174 ! level partially include in ice shelf boundary layer 241 ! compute tfreez for the temperature correction (we add water at freezing temperature)242 175 ! compute trend 243 tsa(ji,jj,ikb,jp_tem) = tsa(ji,jj,ikb,jp_tem) &244 & + zfact * ( risf_tsc_b(ji,jj,jp_tem) + risf_tsc(ji,jj,jp_tem)) * r1_hisf_tbl(ji,jj) * ralpha(ji,jj)245 tsa(ji,jj,ikb,jp_sal) = tsa(ji,jj,ikb,jp_sal) &246 & + zfact * (risf_tsc_b(ji,jj,jp_sal) + risf_tsc(ji,jj,jp_sal)) * r1_hisf_tbl(ji,jj) * ralpha(ji,jj) 176 tsa(ji,jj,ikb,jp_tem) = tsa(ji,jj,ikb,jp_tem) & 177 & + zfact * ( risf_tsc_b(ji,jj,jp_tem) + risf_tsc(ji,jj,jp_tem) ) & 178 & * r1_hisf_tbl(ji,jj) * ralpha(ji,jj) 179 247 180 END DO 248 181 END DO 249 182 IF( lrst_oce ) THEN 250 IF(lwp) WRITE(numout,*) 251 IF(lwp) WRITE(numout,*) 'sbc : isf surface tracer content forcing fields written in ocean restart file ', & 252 & 'at it= ', kt,' date= ', ndastp 253 IF(lwp) WRITE(numout,*) '~~~~' 254 CALL iom_rstput( kt, nitrst, numrow, 'fwf_isf_b', fwfisf(:,:) ) 183 CALL iom_rstput( kt, nitrst, numrow, 'fwf_isf_b', fwfisf (:,:) ) 255 184 CALL iom_rstput( kt, nitrst, numrow, 'isf_hc_b' , risf_tsc(:,:,jp_tem) ) 256 185 CALL iom_rstput( kt, nitrst, numrow, 'isf_sc_b' , risf_tsc(:,:,jp_sal) ) … … 269 198 zdep = zfact / h_rnf(ji,jj) 270 199 DO jk = 1, nk_rnf(ji,jj) 271 tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem) &272 &+ ( rnf_tsc_b(ji,jj,jp_tem) + rnf_tsc(ji,jj,jp_tem) ) * zdep273 IF( ln_rnf_sal ) tsa(ji,jj,jk,jp_sal) = tsa(ji,jj,jk,jp_sal) &274 &+ ( rnf_tsc_b(ji,jj,jp_sal) + rnf_tsc(ji,jj,jp_sal) ) * zdep200 tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem) & 201 & + ( rnf_tsc_b(ji,jj,jp_tem) + rnf_tsc(ji,jj,jp_tem) ) * zdep 202 IF( ln_rnf_sal ) tsa(ji,jj,jk,jp_sal) = tsa(ji,jj,jk,jp_sal) & 203 & + ( rnf_tsc_b(ji,jj,jp_sal) + rnf_tsc(ji,jj,jp_sal) ) * zdep 275 204 END DO 276 205 ENDIF … … 278 207 END DO 279 208 ENDIF 280 281 IF( l_trdtra ) THEN ! send trends for further diagnostics 209 ! 210 !---------------------------------------- 211 ! Ice Sheet coupling imbalance correction to have conservation 212 !---------------------------------------- 213 ! 214 IF( ln_iscpl .AND. ln_hsb) THEN ! input of heat and salt due to river runoff 215 DO jk = 1,jpk 216 DO jj = 2, jpj 217 DO ji = fs_2, fs_jpim1 218 zdep = 1._wp / e3t_n(ji,jj,jk) 219 tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem) - htsc_iscpl(ji,jj,jk,jp_tem) & 220 & * zdep 221 tsa(ji,jj,jk,jp_sal) = tsa(ji,jj,jk,jp_sal) - htsc_iscpl(ji,jj,jk,jp_sal) & 222 & * zdep 223 END DO 224 END DO 225 END DO 226 ENDIF 227 228 IF( l_trdtra ) THEN ! save the horizontal diffusive trends for further diagnostics 282 229 ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:) 283 230 ztrds(:,:,:) = tsa(:,:,:,jp_sal) - ztrds(:,:,:) -
trunk/NEMOGCM/NEMO/OPA_SRC/TRA/trazdf.F90
r5930 r6140 9 9 10 10 !!---------------------------------------------------------------------- 11 !! tra_zdf : Update the tracer trend with the vertical diffusion12 !! tra_zdf_init : initialisation of the computation11 !! tra_zdf : Update the tracer trend with the vertical diffusion 12 !! tra_zdf_init : initialisation of the computation 13 13 !!---------------------------------------------------------------------- 14 USE oce ! ocean dynamics and tracers variables 15 USE dom_oce ! ocean space and time domain variables 16 USE domvvl ! variable volume 17 USE phycst ! physical constant 18 USE zdf_oce ! ocean vertical physics variables 19 USE sbc_oce ! surface boundary condition: ocean 14 USE oce ! ocean dynamics and tracers variables 15 USE dom_oce ! ocean space and time domain variables 16 USE domvvl ! variable volume 17 USE phycst ! physical constant 18 USE zdf_oce ! ocean vertical physics variables 19 USE sbc_oce ! surface boundary condition: ocean 20 USE ldftra ! lateral diffusion: eddy diffusivity 21 USE ldfslp ! lateral diffusion: iso-neutral slope 22 USE trazdf_exp ! vertical diffusion: explicit (tra_zdf_exp routine) 23 USE trazdf_imp ! vertical diffusion: implicit (tra_zdf_imp routine) 24 USE trd_oce ! trends: ocean variables 25 USE trdtra ! trends: tracer trend manager 20 26 ! 21 USE ldftra ! lateral diffusion: eddy diffusivity 22 USE ldfslp ! lateral diffusion: iso-neutral slope 23 USE trazdf_exp ! vertical diffusion: explicit (tra_zdf_exp routine) 24 USE trazdf_imp ! vertical diffusion: implicit (tra_zdf_imp routine) 25 ! 26 USE trd_oce ! trends: ocean variables 27 USE trdtra ! trends manager: tracers 28 ! 29 USE in_out_manager ! I/O manager 30 USE prtctl ! Print control 31 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 32 USE lib_mpp ! MPP library 33 USE wrk_nemo ! Memory allocation 34 USE timing ! Timing 27 USE in_out_manager ! I/O manager 28 USE prtctl ! Print control 29 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 30 USE lib_mpp ! MPP library 31 USE wrk_nemo ! Memory allocation 32 USE timing ! Timing 35 33 36 34 IMPLICIT NONE … … 43 41 44 42 !! * Substitutions 45 # include "domzgr_substitute.h90"46 43 # include "zdfddm_substitute.h90" 47 44 # include "vectopt_loop_substitute.h90" … … 60 57 !!--------------------------------------------------------------------- 61 58 INTEGER, INTENT( in ) :: kt ! ocean time-step index 62 ! !59 ! 63 60 INTEGER :: jk ! Dummy loop indices 64 61 REAL(wp), POINTER, DIMENSION(:,:,:) :: ztrdt, ztrds ! 3D workspace … … 68 65 ! 69 66 IF( neuler == 0 .AND. kt == nit000 ) THEN ! at nit000 70 r2dt ra(:) = rdttra(:) ! = rdtra(restarting with Euler time stepping)67 r2dt = rdt ! = rdt (restarting with Euler time stepping) 71 68 ELSEIF( kt <= nit000 + 1) THEN ! at nit000 or nit000+1 72 r2dt ra(:) = 2. * rdttra(:) ! = 2 rdttra(leapfrog)69 r2dt = 2. * rdt ! = 2 rdt (leapfrog) 73 70 ENDIF 74 71 ! 75 72 IF( l_trdtra ) THEN !* Save ta and sa trends 76 73 CALL wrk_alloc( jpi, jpj, jpk, ztrdt, ztrds ) … … 78 75 ztrds(:,:,:) = tsa(:,:,:,jp_sal) 79 76 ENDIF 80 77 ! 81 78 SELECT CASE ( nzdf ) ! compute lateral mixing trend and add it to the general trend 82 CASE ( 0 ) ; CALL tra_zdf_exp( kt, nit000, 'TRA', r2dt ra, nn_zdfexp, tsb, tsa, jpts ) ! explicit scheme83 CASE ( 1 ) ; CALL tra_zdf_imp( kt, nit000, 'TRA', r2dt ra, tsb, tsa, jpts ) ! implicit scheme79 CASE ( 0 ) ; CALL tra_zdf_exp( kt, nit000, 'TRA', r2dt, nn_zdfexp, tsb, tsa, jpts ) ! explicit scheme 80 CASE ( 1 ) ; CALL tra_zdf_imp( kt, nit000, 'TRA', r2dt, tsb, tsa, jpts ) ! implicit scheme 84 81 END SELECT 85 82 !!gm WHY here ! and I don't like that ! … … 87 84 ! JMM avoid negative salinities near river outlet ! Ugly fix 88 85 ! JMM : restore negative salinities to small salinities: 89 WHERE 86 WHERE( tsa(:,:,:,jp_sal) < 0._wp ) tsa(:,:,:,jp_sal) = 0.1_wp 90 87 !!gm 91 88 92 89 IF( l_trdtra ) THEN ! save the vertical diffusive trends for further diagnostics 93 90 DO jk = 1, jpkm1 94 ztrdt(:,:,jk) = ( ( tsa(:,:,jk,jp_tem) - tsb(:,:,jk,jp_tem) ) / r2dt ra(jk)) - ztrdt(:,:,jk)95 ztrds(:,:,jk) = ( ( tsa(:,:,jk,jp_sal) - tsb(:,:,jk,jp_sal) ) / r2dt ra(jk)) - ztrds(:,:,jk)91 ztrdt(:,:,jk) = ( ( tsa(:,:,jk,jp_tem) - tsb(:,:,jk,jp_tem) ) / r2dt ) - ztrdt(:,:,jk) 92 ztrds(:,:,jk) = ( ( tsa(:,:,jk,jp_sal) - tsb(:,:,jk,jp_sal) ) / r2dt ) - ztrds(:,:,jk) 96 93 END DO 97 94 !!gm this should be moved in trdtra.F90 and done on all trends … … 103 100 CALL wrk_dealloc( jpi, jpj, jpk, ztrdt, ztrds ) 104 101 ENDIF 105 106 102 ! ! print mean trends (used for debugging) 107 103 IF(ln_ctl) CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' zdf - Ta: ', mask1=tmask, & … … 128 124 USE zdfgls 129 125 !!---------------------------------------------------------------------- 130 126 ! 131 127 ! Choice from ln_zdfexp already read in namelist in zdfini module 132 128 IF( ln_zdfexp ) THEN ; nzdf = 0 ! use explicit scheme 133 129 ELSE ; nzdf = 1 ! use implicit scheme 134 130 ENDIF 135 131 ! 136 132 ! Force implicit schemes 137 133 IF( lk_zdftke .OR. lk_zdfgls ) nzdf = 1 ! TKE, or GLS physics … … 140 136 IF( ln_zdfexp .AND. nzdf == 1 ) CALL ctl_stop( 'tra_zdf : If using the rotation of lateral mixing operator', & 141 137 & ' GLS or TKE scheme, the implicit scheme is required, set ln_zdfexp = .false.' ) 142 138 ! 143 139 IF(lwp) THEN 144 140 WRITE(numout,*) -
trunk/NEMOGCM/NEMO/OPA_SRC/TRA/trazdf_exp.F90
r3294 r6140 20 20 21 21 !!---------------------------------------------------------------------- 22 !! tra_zdf_exp : compute the tracer the vertical diffusion trend using a23 !! split-explicit time stepping and provide the after tracer22 !! tra_zdf_exp : compute the tracer the vertical diffusion trend using a 23 !! split-explicit time stepping and provide the after tracer 24 24 !!---------------------------------------------------------------------- 25 USE oce ! ocean dynamics and active tracers 26 USE dom_oce ! ocean space and time domain 27 USE domvvl ! variable volume levels 28 USE zdf_oce ! ocean vertical physics 29 USE zdfddm ! ocean vertical physics: double diffusion 30 USE trc_oce ! share passive tracers/Ocean variables 31 USE in_out_manager ! I/O manager 32 USE lib_mpp ! MPP library 33 USE wrk_nemo ! Memory Allocation 34 USE timing ! Timing 25 USE oce ! ocean dynamics and active tracers 26 USE dom_oce ! ocean space and time domain 27 USE domvvl ! variable volume levels 28 USE zdf_oce ! ocean vertical physics 29 USE zdfddm ! ocean vertical physics: double diffusion 30 USE trc_oce ! share passive tracers/Ocean variables 31 ! 32 USE in_out_manager ! I/O manager 33 USE lib_mpp ! MPP library 34 USE wrk_nemo ! Memory Allocation 35 USE timing ! Timing 35 36 36 37 IMPLICIT NONE … … 40 41 41 42 !! * Substitutions 42 # include "domzgr_substitute.h90"43 43 # include "zdfddm_substitute.h90" 44 44 # include "vectopt_loop_substitute.h90" … … 50 50 CONTAINS 51 51 52 SUBROUTINE tra_zdf_exp( kt, kit000, cdtype, p2dt, k n_zdfexp, &53 & ptb , pta, kjpt )52 SUBROUTINE tra_zdf_exp( kt, kit000, cdtype, p2dt, ksts, & 53 & ptb , pta , kjpt ) 54 54 !!---------------------------------------------------------------------- 55 55 !! *** ROUTINE tra_zdf_exp *** … … 60 60 !! ** Method : - The after tracer fields due to the vertical diffusion 61 61 !! of tracers alone is given by: 62 !! z wx= ptb + p2dt difft62 !! ztb = ptb + p2dt difft 63 63 !! where difft = dz( avt dz(ptb) ) = 1/e3t dk+1( avt/e3w dk(ptb) ) 64 64 !! (if lk_zdfddm=T use avs on salinity and passive tracers instead of avt) … … 67 67 !! (N.B. bottom condition is applied through the masked field avt). 68 68 !! - the after tracer fields due to the whole trend is 69 !! obtained in leap-frog environment by : 70 !! pta = zwx + p2dt pta 71 !! - in case of variable level thickness (lk_vvl=T) the 72 !! the leap-frog is applied on thickness weighted tracer. That is: 73 !! pta = [ ptb*e3tb + e3tn*( zwx - ptb + p2dt pta ) ] / e3tn 69 !! obtained in leap-frog environment applied on thickness weighted tracer by : 70 !! pta = [ ptb*e3tb + e3tn*( ztb - ptb + p2dt pta ) ] / e3tn 74 71 !! 75 72 !! ** Action : - after tracer fields pta 76 73 !!--------------------------------------------------------------------- 74 INTEGER , INTENT(in ) :: kt ! ocean time-step index 75 INTEGER , INTENT(in ) :: kit000 ! first time step index 76 CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator) 77 INTEGER , INTENT(in ) :: kjpt ! number of tracers 78 INTEGER , INTENT(in ) :: ksts ! number of sub-time step 79 REAL(wp) , INTENT(in ) :: p2dt ! vertical profile of tracer time-step 80 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: ptb ! before and now tracer fields 81 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pta ! in: tracer trend ; out: after tracer field 77 82 ! 78 INTEGER , INTENT(in ) :: kt ! ocean time-step index 79 INTEGER , INTENT(in ) :: kit000 ! first time step index 80 CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator) 81 INTEGER , INTENT(in ) :: kjpt ! number of tracers 82 INTEGER , INTENT(in ) :: kn_zdfexp ! number of sub-time step 83 REAL(wp), DIMENSION( jpk ), INTENT(in ) :: p2dt ! vertical profile of tracer time-step 84 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: ptb ! before and now tracer fields 85 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pta ! tracer trend 86 ! 87 INTEGER :: ji, jj, jk, jn, jl ! dummy loop indices 88 REAL(wp) :: zlavmr, zave3r, ze3tr ! local scalars 89 REAL(wp) :: ztra, ze3tb ! - - 90 REAL(wp), POINTER, DIMENSION(:,:,:) :: zwx, zwy 83 INTEGER :: ji, jj, jk, jn, jl ! dummy loop indices 84 REAL(wp) :: z1_ksts, ze3tr ! local scalars 85 REAL(wp) :: ztra, ze3tb ! - - 86 REAL(wp), POINTER, DIMENSION(:,:,:) :: ztb, zwf 91 87 !!--------------------------------------------------------------------- 92 88 ! 93 89 IF( nn_timing == 1 ) CALL timing_start('tra_zdf_exp') 94 90 ! 95 CALL wrk_alloc( jpi, jpj, jpk, zwx, zwy)91 CALL wrk_alloc( jpi,jpj,jpk, ztb, zwf ) 96 92 ! 97 98 93 IF( kt == kit000 ) THEN 99 94 IF(lwp) WRITE(numout,*) … … 104 99 ! Initializations 105 100 ! --------------- 106 zlavmr = 1. / float( kn_zdfexp ) ! Local constant 101 z1_ksts = 1._wp / REAL( ksts, wp ) 102 zwf(:,:, 1 ) = 0._wp ! no flux at the surface and at bottom level 103 zwf(:,:,jpk) = 0._wp 107 104 ! 108 105 ! 109 DO jn = 1, kjpt ! loop over tracers106 DO jn = 1, kjpt !== loop over tracers ==! 110 107 ! 111 zwy(:,:, 1 ) = 0.e0 ! surface boundary conditions: no flux 112 zwy(:,:,jpk) = 0.e0 ! bottom boundary conditions: no flux 113 ! 114 zwx(:,:,:) = ptb(:,:,:,jn) ! zwx array set to before tracer values 115 116 ! Split-explicit loop (after tracer due to the vertical diffusion alone) 117 ! ------------------- 118 ! 119 DO jl = 1, kn_zdfexp 120 ! ! first vertical derivative 121 DO jk = 2, jpk 108 ztb(:,:,:) = ptb(:,:,:,jn) ! initial before value for tracer 109 ! 110 DO jl = 1, ksts !== Split-explicit loop ==! 111 ! 112 DO jk = 2, jpk ! 1st vertical derivative (w-flux) 122 113 DO jj = 2, jpjm1 123 114 DO ji = fs_2, fs_jpim1 ! vector opt. 124 zave3r = 1.e0 / fse3w_n(ji,jj,jk)125 115 IF( cdtype == 'TRA' .AND. jn == jp_tem ) THEN ! temperature : use of avt 126 zw y(ji,jj,jk) = avt(ji,jj,jk) * ( zwx(ji,jj,jk-1) - zwx(ji,jj,jk) ) * zave3r116 zwf(ji,jj,jk) = avt(ji,jj,jk) * ( ztb(ji,jj,jk-1) - ztb(ji,jj,jk) ) / e3w_b(ji,jj,jk) 127 117 ELSE ! salinity or pass. tracer : use of avs 128 zw y(ji,jj,jk) = fsavs(ji,jj,jk) * ( zwx(ji,jj,jk-1) - zwx(ji,jj,jk) ) * zave3r118 zwf(ji,jj,jk) = fsavs(ji,jj,jk) * ( ztb(ji,jj,jk-1) - ztb(ji,jj,jk) ) / e3w_b(ji,jj,jk) 129 119 END IF 130 120 END DO … … 132 122 END DO 133 123 ! 134 DO jk = 1, jpkm1 ! second vertical derivative ==> tracer at kt+l*2*rdt/nn_zdfexp124 DO jk = 1, jpkm1 ! 2nd vertical derivative ==> tracer at kt+l*2*rdt/nn_zdfexp 135 125 DO jj = 2, jpjm1 136 126 DO ji = fs_2, fs_jpim1 ! vector opt. 137 ze3tr = zlavmr / fse3t_n(ji,jj,jk) 138 zwx(ji,jj,jk) = zwx(ji,jj,jk) + p2dt(jk) * ( zwy(ji,jj,jk) - zwy(ji,jj,jk+1) ) * ze3tr 127 ztb(ji,jj,jk) = ztb(ji,jj,jk) + p2dt * ( zwf(ji,jj,jk) - zwf(ji,jj,jk+1) ) / e3t_n(ji,jj,jk) 139 128 END DO 140 129 END DO 141 130 END DO 142 131 ! 143 END DO 132 END DO ! end sub-time stepping 144 133 145 ! After tracer due to all trends 146 ! ------------------------------ 147 IF( lk_vvl ) THEN ! variable level thickness : leap-frog on tracer*e3t 148 DO jk = 1, jpkm1 149 DO jj = 2, jpjm1 150 DO ji = fs_2, fs_jpim1 ! vector opt. 151 ze3tb = fse3t_b(ji,jj,jk) / fse3t(ji,jj,jk) ! before e3t 152 ztra = zwx(ji,jj,jk) - ptb(ji,jj,jk,jn) + p2dt(jk) * pta(ji,jj,jk,jn) ! total trends * 2*rdt 153 pta(ji,jj,jk,jn) = ( ze3tb * ptb(ji,jj,jk,jn) + ztra ) * tmask(ji,jj,jk) 154 END DO 134 DO jk = 1, jpkm1 !== After tracer due to all trends 135 DO jj = 2, jpjm1 136 DO ji = fs_2, fs_jpim1 ! vector opt. 137 ze3tb = e3t_b(ji,jj,jk) / e3t_n(ji,jj,jk) 138 ztra = ( ztb(ji,jj,jk) - ptb(ji,jj,jk,jn) ) + p2dt * pta(ji,jj,jk,jn) ! total trend * 2dt 139 pta(ji,jj,jk,jn) = ( ze3tb * ptb(ji,jj,jk,jn) + ztra ) * tmask(ji,jj,jk) ! after tracer 155 140 END DO 156 141 END DO 157 ELSE ! fixed level thickness : leap-frog on tracers 158 DO jk = 1, jpkm1 159 DO jj = 2, jpjm1 160 DO ji = fs_2, fs_jpim1 ! vector opt. 161 pta(ji,jj,jk,jn) = ( zwx(ji,jj,jk) + p2dt(jk) * pta(ji,jj,jk,jn) ) * tmask(ji,jj,jk) 162 END DO 163 END DO 164 END DO 165 ENDIF 142 END DO 166 143 ! 167 END DO 144 END DO ! end of tracer loop 168 145 ! 169 CALL wrk_dealloc( jpi, jpj, jpk, zwx, zwy)146 CALL wrk_dealloc( jpi,jpj,jpk, ztb, zwf ) 170 147 ! 171 148 IF( nn_timing == 1 ) CALL timing_stop('tra_zdf_exp') -
trunk/NEMOGCM/NEMO/OPA_SRC/TRA/trazdf_imp.F90
r5836 r6140 16 16 !! 3.3 ! 2010-06 (C. Ethe, G. Madec) Merge TRA-TRC 17 17 !! - ! 2011-02 (A. Coward, C. Ethe, G. Madec) improvment of surface boundary condition 18 !! 3.7 ! 2015-11 (G. Madec, A. Coward) non linear free surface by default 18 19 !!---------------------------------------------------------------------- 19 20 20 21 !!---------------------------------------------------------------------- 21 !! tra_zdf_imp : Update the tracer trend with the diagonal vertical part of the mixing tensor.22 !! tra_zdf_imp : Update the tracer trend with vertical mixing, nad compute the after tracer field 22 23 !!---------------------------------------------------------------------- 23 24 USE oce ! ocean dynamics and tracers variables … … 42 43 PUBLIC tra_zdf_imp ! routine called by step.F90 43 44 44 REAL(wp) :: r_vvl ! variable volume indicator, =1 if lk_vvl=T, =0 otherwise45 46 45 !! * Substitutions 47 # include "domzgr_substitute.h90"48 46 # include "zdfddm_substitute.h90" 49 47 # include "vectopt_loop_substitute.h90" … … 64 62 !! it is already computed and add to the general trend in traldf) 65 63 !! 66 !! ** Method : The vertical diffusion of the tracer t is given by: 67 !! difft = dz( avt dz(t) ) = 1/e3t dk+1( avt/e3w dk(t) ) 68 !! It is computed using a backward time scheme (t=ta). 64 !! ** Method : The vertical diffusion of a tracer ,t , is given by: 65 !! difft = dz( avt dz(t) ) = 1/e3t dk+1( avt/e3w dk(t) ) 66 !! It is computed using a backward time scheme (t=after field) 67 !! which provide directly the after tracer field. 69 68 !! If lk_zdfddm=T, use avs for salinity or for passive tracers 70 69 !! Surface and bottom boundary conditions: no diffusive flux on … … 75 74 !!--------------------------------------------------------------------- 76 75 INTEGER , INTENT(in ) :: kt ! ocean time-step index 77 INTEGER , INTENT(in ) :: kit000 76 INTEGER , INTENT(in ) :: kit000 ! first time step index 78 77 CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator) 79 78 INTEGER , INTENT(in ) :: kjpt ! number of tracers 80 REAL(wp) , DIMENSION( jpk ), INTENT(in ) :: p2dt ! vertical profile oftracer time-step79 REAL(wp) , INTENT(in ) :: p2dt ! tracer time-step 81 80 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: ptb ! before and now tracer fields 82 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pta ! tracer trend81 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pta ! in: tracer trend ; out: after tracer field 83 82 ! 84 83 INTEGER :: ji, jj, jk, jn ! dummy loop indices 85 REAL(wp) :: zrhs , ze3tb, ze3tn, ze3ta! local scalars84 REAL(wp) :: zrhs ! local scalars 86 85 REAL(wp), POINTER, DIMENSION(:,:,:) :: zwi, zwt, zwd, zws 87 86 !!--------------------------------------------------------------------- … … 95 94 IF(lwp)WRITE(numout,*) 'tra_zdf_imp : implicit vertical mixing on ', cdtype 96 95 IF(lwp)WRITE(numout,*) '~~~~~~~~~~~ ' 97 !98 IF( lk_vvl ) THEN ; r_vvl = 1._wp ! Variable volume indicator99 ELSE ; r_vvl = 0._wp100 ENDIF101 96 ENDIF 102 !103 97 ! ! ============= ! 104 98 DO jn = 1, kjpt ! tracer loop ! 105 99 ! ! ============= ! 106 !107 100 ! Matrix construction 108 101 ! -------------------- … … 142 135 DO jj = 2, jpjm1 143 136 DO ji = fs_2, fs_jpim1 ! vector opt. 144 ze3ta = ( 1. - r_vvl ) + r_vvl * fse3t_a(ji,jj,jk) ! after scale factor at T-point 145 ze3tn = r_vvl + ( 1. - r_vvl ) * fse3t_n(ji,jj,jk) ! now scale factor at T-point 146 zwi(ji,jj,jk) = - p2dt(jk) * zwt(ji,jj,jk ) / ( ze3tn * fse3w(ji,jj,jk ) ) 147 zws(ji,jj,jk) = - p2dt(jk) * zwt(ji,jj,jk+1) / ( ze3tn * fse3w(ji,jj,jk+1) ) 148 zwd(ji,jj,jk) = ze3ta - zwi(ji,jj,jk) - zws(ji,jj,jk) 137 !!gm BUG I think, use e3w_a instead of e3w_n 138 zwi(ji,jj,jk) = - p2dt * zwt(ji,jj,jk ) / e3w_n(ji,jj,jk ) 139 zws(ji,jj,jk) = - p2dt * zwt(ji,jj,jk+1) / e3w_n(ji,jj,jk+1) 140 zwd(ji,jj,jk) = e3t_a(ji,jj,jk) - zwi(ji,jj,jk) - zws(ji,jj,jk) 149 141 END DO 150 142 END DO … … 170 162 ! used as a work space array: its value is modified. 171 163 ! 172 ! first recurrence: Tk = Dk - Ik Sk-1 / Tk-1 (increasing k) 173 ! done once for all passive tracers (so included in the IF instruction) 174 DO jj = 2, jpjm1 175 DO ji = fs_2, fs_jpim1 164 DO jj = 2, jpjm1 !* 1st recurrence: Tk = Dk - Ik Sk-1 / Tk-1 (increasing k) 165 DO ji = fs_2, fs_jpim1 ! done one for all passive tracers (so included in the IF instruction) 176 166 zwt(ji,jj,1) = zwd(ji,jj,1) 177 167 END DO … … 185 175 END DO 186 176 ! 187 END 177 ENDIF 188 178 ! 189 ! second recurrence: Zk = Yk - Ik / Tk-1 Zk-1 190 DO jj = 2, jpjm1 179 DO jj = 2, jpjm1 !* 2nd recurrence: Zk = Yk - Ik / Tk-1 Zk-1 191 180 DO ji = fs_2, fs_jpim1 192 ze3tb = ( 1. - r_vvl ) + r_vvl * fse3t_b(ji,jj,1) 193 ze3tn = ( 1. - r_vvl ) + r_vvl * fse3t(ji,jj,1) 194 pta(ji,jj,1,jn) = ze3tb * ptb(ji,jj,1,jn) + p2dt(1) * ze3tn * pta(ji,jj,1,jn) 181 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) 195 182 END DO 196 183 END DO … … 198 185 DO jj = 2, jpjm1 199 186 DO ji = fs_2, fs_jpim1 200 ze3tb = ( 1. - r_vvl ) + r_vvl * fse3t_b(ji,jj,jk) 201 ze3tn = ( 1. - r_vvl ) + r_vvl * fse3t (ji,jj,jk) 202 zrhs = ze3tb * ptb(ji,jj,jk,jn) + p2dt(jk) * ze3tn * pta(ji,jj,jk,jn) ! zrhs=right hand side 187 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 203 188 pta(ji,jj,jk,jn) = zrhs - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) * pta(ji,jj,jk-1,jn) 204 189 END DO 205 190 END DO 206 191 END DO 207 208 ! third recurrence: Xk = (Zk - Sk Xk+1 ) / Tk (result is the after tracer) 209 DO jj = 2, jpjm1 192 ! 193 DO jj = 2, jpjm1 !* 3d recurrence: Xk = (Zk - Sk Xk+1 ) / Tk (result is the after tracer) 210 194 DO ji = fs_2, fs_jpim1 211 195 pta(ji,jj,jpkm1,jn) = pta(ji,jj,jpkm1,jn) / zwt(ji,jj,jpkm1) * tmask(ji,jj,jpkm1) -
trunk/NEMOGCM/NEMO/OPA_SRC/TRA/zpshde.F90
r5836 r6140 32 32 33 33 !! * Substitutions 34 # include "domzgr_substitute.h90"35 34 # include "vectopt_loop_substitute.h90" 36 35 !!---------------------------------------------------------------------- … … 111 110 iku = mbku(ji,jj) ; ikum1 = MAX( iku - 1 , 1 ) ! last and before last ocean level at u- & v-points 112 111 ikv = mbkv(ji,jj) ; ikvm1 = MAX( ikv - 1 , 1 ) ! if level first is a p-step, ik.m1=1 113 ze3wu = fse3w(ji+1,jj ,iku) - fse3w(ji,jj,iku) 114 ze3wv = fse3w(ji ,jj+1,ikv) - fse3w(ji,jj,ikv) 112 !!gm BUG ? when applied to before fields, e3w_b should be used.... 113 ze3wu = e3w_n(ji+1,jj ,iku) - e3w_n(ji,jj,iku) 114 ze3wv = e3w_n(ji ,jj+1,ikv) - e3w_n(ji,jj,ikv) 115 115 ! 116 116 ! i- direction 117 117 IF( ze3wu >= 0._wp ) THEN ! case 1 118 zmaxu = ze3wu / fse3w(ji+1,jj,iku)118 zmaxu = ze3wu / e3w_n(ji+1,jj,iku) 119 119 ! interpolated values of tracers 120 120 zti (ji,jj,jn) = pta(ji+1,jj,iku,jn) + zmaxu * ( pta(ji+1,jj,ikum1,jn) - pta(ji+1,jj,iku,jn) ) … … 122 122 pgtu(ji,jj,jn) = umask(ji,jj,1) * ( zti(ji,jj,jn) - pta(ji,jj,iku,jn) ) 123 123 ELSE ! case 2 124 zmaxu = -ze3wu / fse3w(ji,jj,iku)124 zmaxu = -ze3wu / e3w_n(ji,jj,iku) 125 125 ! interpolated values of tracers 126 126 zti (ji,jj,jn) = pta(ji,jj,iku,jn) + zmaxu * ( pta(ji,jj,ikum1,jn) - pta(ji,jj,iku,jn) ) … … 131 131 ! j- direction 132 132 IF( ze3wv >= 0._wp ) THEN ! case 1 133 zmaxv = ze3wv / fse3w(ji,jj+1,ikv)133 zmaxv = ze3wv / e3w_n(ji,jj+1,ikv) 134 134 ! interpolated values of tracers 135 135 ztj (ji,jj,jn) = pta(ji,jj+1,ikv,jn) + zmaxv * ( pta(ji,jj+1,ikvm1,jn) - pta(ji,jj+1,ikv,jn) ) … … 137 137 pgtv(ji,jj,jn) = vmask(ji,jj,1) * ( ztj(ji,jj,jn) - pta(ji,jj,ikv,jn) ) 138 138 ELSE ! case 2 139 zmaxv = -ze3wv / fse3w(ji,jj,ikv)139 zmaxv = -ze3wv / e3w_n(ji,jj,ikv) 140 140 ! interpolated values of tracers 141 141 ztj (ji,jj,jn) = pta(ji,jj,ikv,jn) + zmaxv * ( pta(ji,jj,ikvm1,jn) - pta(ji,jj,ikv,jn) ) … … 156 156 iku = mbku(ji,jj) 157 157 ikv = mbkv(ji,jj) 158 ze3wu = fse3w(ji+1,jj ,iku) - fse3w(ji,jj,iku)159 ze3wv = fse3w(ji ,jj+1,ikv) - fse3w(ji,jj,ikv)160 IF( ze3wu >= 0._wp ) THEN ; zhi(ji,jj) = fsdept(ji ,jj,iku) ! i-direction: case 1161 ELSE ; zhi(ji,jj) = fsdept(ji+1,jj,iku) ! - - case 2162 ENDIF 163 IF( ze3wv >= 0._wp ) THEN ; zhj(ji,jj) = fsdept(ji,jj ,ikv) ! j-direction: case 1164 ELSE ; zhj(ji,jj) = fsdept(ji,jj+1,ikv) ! - - case 2158 ze3wu = e3w_n(ji+1,jj ,iku) - e3w_n(ji,jj,iku) 159 ze3wv = e3w_n(ji ,jj+1,ikv) - e3w_n(ji,jj,ikv) 160 IF( ze3wu >= 0._wp ) THEN ; zhi(ji,jj) = gdept_n(ji ,jj,iku) ! i-direction: case 1 161 ELSE ; zhi(ji,jj) = gdept_n(ji+1,jj,iku) ! - - case 2 162 ENDIF 163 IF( ze3wv >= 0._wp ) THEN ; zhj(ji,jj) = gdept_n(ji,jj ,ikv) ! j-direction: case 1 164 ELSE ; zhj(ji,jj) = gdept_n(ji,jj+1,ikv) ! - - case 2 165 165 ENDIF 166 166 END DO … … 174 174 iku = mbku(ji,jj) 175 175 ikv = mbkv(ji,jj) 176 ze3wu = fse3w(ji+1,jj ,iku) - fse3w(ji,jj,iku)177 ze3wv = fse3w(ji ,jj+1,ikv) - fse3w(ji,jj,ikv)176 ze3wu = e3w_n(ji+1,jj ,iku) - e3w_n(ji,jj,iku) 177 ze3wv = e3w_n(ji ,jj+1,ikv) - e3w_n(ji,jj,ikv) 178 178 IF( ze3wu >= 0._wp ) THEN ; pgru(ji,jj) = umask(ji,jj,1) * ( zri(ji ,jj ) - prd(ji,jj,iku) ) ! i: 1 179 179 ELSE ; pgru(ji,jj) = umask(ji,jj,1) * ( prd(ji+1,jj,iku) - zri(ji,jj ) ) ! i: 2 … … 191 191 ! 192 192 END SUBROUTINE zps_hde 193 194 195 SUBROUTINE zps_hde_isf( kt, kjpt, pta, pgtu , pgtv , pgtui, pgtvi, & 196 & prd, pgru , pgrv , pmru , pmrv , pgzu , pgzv , pge3ru , pge3rv , & 197 & pgrui, pgrvi, pmrui, pmrvi, pgzui, pgzvi, pge3rui, pge3rvi ) 198 !!---------------------------------------------------------------------- 199 !! *** ROUTINE zps_hde *** 193 ! 194 SUBROUTINE zps_hde_isf( kt, kjpt, pta, pgtu, pgtv, pgtui, pgtvi, & 195 & prd, pgru, pgrv, pgrui, pgrvi ) 196 !!---------------------------------------------------------------------- 197 !! *** ROUTINE zps_hde_isf *** 200 198 !! 201 199 !! ** Purpose : Compute the horizontal derivative of T, S and rho 202 200 !! at u- and v-points with a linear interpolation for z-coordinate 203 !! with partial steps .201 !! with partial steps for top (ice shelf) and bottom. 204 202 !! 205 203 !! ** Method : In z-coord with partial steps, scale factors on last 206 204 !! levels are different for each grid point, so that T, S and rd 207 205 !! points are not at the same depth as in z-coord. To have horizontal 208 !! gradients again, we interpolate T and S at the good depth : 206 !! gradients again, we interpolate T and S at the good depth : 207 !! For the bottom case: 209 208 !! Linear interpolation of T, S 210 209 !! Computation of di(tb) and dj(tb) by vertical interpolation: … … 235 234 !! di(rho) = rd~ - rd(i,j,k) or rd(i+1,j,k) - rd~ 236 235 !! 236 !! For the top case (ice shelf): As for the bottom case but upside down 237 !! 237 238 !! ** Action : compute for top and bottom interfaces 238 239 !! - pgtu, pgtv, pgtui, pgtvi: horizontal gradient of tracer at u- & v-points 239 240 !! - pgru, pgrv, pgrui, pgtvi: horizontal gradient of rho (if present) at u- & v-points 240 !! - pmru, pmrv, pmrui, pmrvi: horizontal sum of rho at u- & v- point (used in dynhpg with vvl) 241 !! - pgzu, pgzv, pgzui, pgzvi: horizontal gradient of z at u- and v- point (used in dynhpg with vvl) 242 !! - pge3ru, pge3rv, pge3rui, pge3rvi: horizontal gradient of rho weighted by local e3w at u- & v-points 243 !!---------------------------------------------------------------------- 244 INTEGER , INTENT(in ) :: kt ! ocean time-step index 245 INTEGER , INTENT(in ) :: kjpt ! number of tracers 246 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: pta ! 4D tracers fields 247 ! !! u-point ! v-point ! 248 REAL(wp), DIMENSION(jpi,jpj, kjpt), INTENT( out) :: pgtu , pgtv ! bottom GRADh( ptra ) 249 REAL(wp), DIMENSION(jpi,jpj, kjpt), INTENT( out) :: pgtui , pgtvi ! top GRADh( ptra ) 250 REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ), OPTIONAL :: prd ! 3D density anomaly fields 251 ! !! u-point ! v-point ! 252 REAL(wp), DIMENSION(jpi,jpj ), INTENT( out), OPTIONAL :: pgru , pgrv ! bottom GRADh( prd ) 253 REAL(wp), DIMENSION(jpi,jpj ), INTENT( out), OPTIONAL :: pmru , pmrv ! bottom SUM ( prd ) 254 REAL(wp), DIMENSION(jpi,jpj ), INTENT( out), OPTIONAL :: pgzu , pgzv ! bottom GRADh( z ) 255 REAL(wp), DIMENSION(jpi,jpj ), INTENT( out), OPTIONAL :: pge3ru , pge3rv ! bottom GRADh( prd ) weighted by e3w 256 REAL(wp), DIMENSION(jpi,jpj ), INTENT( out), OPTIONAL :: pgrui , pgrvi ! top GRADh( prd ) 257 REAL(wp), DIMENSION(jpi,jpj ), INTENT( out), OPTIONAL :: pmrui , pmrvi ! top SUM ( prd ) 258 REAL(wp), DIMENSION(jpi,jpj ), INTENT( out), OPTIONAL :: pgzui , pgzvi ! top GRADh( z ) 259 REAL(wp), DIMENSION(jpi,jpj ), INTENT( out), OPTIONAL :: pge3rui , pge3rvi ! top GRADh( prd ) weighted by e3w 241 !!---------------------------------------------------------------------- 242 INTEGER , INTENT(in ) :: kt ! ocean time-step index 243 INTEGER , INTENT(in ) :: kjpt ! number of tracers 244 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: pta ! 4D tracers fields 245 REAL(wp), DIMENSION(jpi,jpj, kjpt), INTENT( out) :: pgtu, pgtv ! hor. grad. of ptra at u- & v-pts 246 REAL(wp), DIMENSION(jpi,jpj, kjpt), INTENT( out) :: pgtui, pgtvi ! hor. grad. of stra at u- & v-pts (ISF) 247 REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ), OPTIONAL :: prd ! 3D density anomaly fields 248 REAL(wp), DIMENSION(jpi,jpj ), INTENT( out), OPTIONAL :: pgru, pgrv ! hor. grad of prd at u- & v-pts (bottom) 249 REAL(wp), DIMENSION(jpi,jpj ), INTENT( out), OPTIONAL :: pgrui, pgrvi ! hor. grad of prd at u- & v-pts (top) 260 250 ! 261 251 INTEGER :: ji, jj, jn ! Dummy loop indices 262 252 INTEGER :: iku, ikv, ikum1, ikvm1,ikup1, ikvp1 ! partial step level (ocean bottom level) at u- and v-points 263 REAL(wp) :: ze3wu, ze3wv, zmaxu, zmaxv , zdzwu, zdzwv, zdzwuip1, zdzwvjp1! temporary scalars253 REAL(wp) :: ze3wu, ze3wv, zmaxu, zmaxv ! temporary scalars 264 254 REAL(wp), DIMENSION(jpi,jpj) :: zri, zrj, zhi, zhj ! NB: 3rd dim=1 to use eos 265 255 REAL(wp), DIMENSION(jpi,jpj,kjpt) :: zti, ztj ! … … 277 267 DO jj = 1, jpjm1 278 268 DO ji = 1, jpim1 279 iku = mbku(ji,jj) ; ikum1 = MAX( iku - 1 , 1 ) ! last and before last ocean level at u- & v-points 280 ikv = mbkv(ji,jj) ; ikvm1 = MAX( ikv - 1 , 1 ) ! if level first is a p-step, ik.m1=1 269 270 iku = mbku(ji,jj); ikum1 = MAX( iku - 1 , 1 ) ! last and before last ocean level at u- & v-points 271 ikv = mbkv(ji,jj); ikvm1 = MAX( ikv - 1 , 1 ) ! if level first is a p-step, ik.m1=1 272 ze3wu = gdept_n(ji+1,jj,iku) - gdept_n(ji,jj,iku) 273 ze3wv = gdept_n(ji,jj+1,ikv) - gdept_n(ji,jj,ikv) 274 ! 275 ! i- direction 276 IF( ze3wu >= 0._wp ) THEN ! case 1 277 zmaxu = ze3wu / e3w_n(ji+1,jj,iku) 278 ! interpolated values of tracers 279 zti (ji,jj,jn) = pta(ji+1,jj,iku,jn) + zmaxu * ( pta(ji+1,jj,ikum1,jn) - pta(ji+1,jj,iku,jn) ) 280 ! gradient of tracers 281 pgtu(ji,jj,jn) = ssumask(ji,jj) * ( zti(ji,jj,jn) - pta(ji,jj,iku,jn) ) 282 ELSE ! case 2 283 zmaxu = -ze3wu / e3w_n(ji,jj,iku) 284 ! interpolated values of tracers 285 zti (ji,jj,jn) = pta(ji,jj,iku,jn) + zmaxu * ( pta(ji,jj,ikum1,jn) - pta(ji,jj,iku,jn) ) 286 ! gradient of tracers 287 pgtu(ji,jj,jn) = ssumask(ji,jj) * ( pta(ji+1,jj,iku,jn) - zti(ji,jj,jn) ) 288 ENDIF 289 ! 290 ! j- direction 291 IF( ze3wv >= 0._wp ) THEN ! case 1 292 zmaxv = ze3wv / e3w_n(ji,jj+1,ikv) 293 ! interpolated values of tracers 294 ztj (ji,jj,jn) = pta(ji,jj+1,ikv,jn) + zmaxv * ( pta(ji,jj+1,ikvm1,jn) - pta(ji,jj+1,ikv,jn) ) 295 ! gradient of tracers 296 pgtv(ji,jj,jn) = ssvmask(ji,jj) * ( ztj(ji,jj,jn) - pta(ji,jj,ikv,jn) ) 297 ELSE ! case 2 298 zmaxv = -ze3wv / e3w_n(ji,jj,ikv) 299 ! interpolated values of tracers 300 ztj (ji,jj,jn) = pta(ji,jj,ikv,jn) + zmaxv * ( pta(ji,jj,ikvm1,jn) - pta(ji,jj,ikv,jn) ) 301 ! gradient of tracers 302 pgtv(ji,jj,jn) = ssvmask(ji,jj) * ( pta(ji,jj+1,ikv,jn) - ztj(ji,jj,jn) ) 303 ENDIF 304 305 END DO 306 END DO 307 CALL lbc_lnk( pgtu(:,:,jn), 'U', -1. ) ; CALL lbc_lnk( pgtv(:,:,jn), 'V', -1. ) ! Lateral boundary cond. 308 ! 309 END DO 310 311 ! horizontal derivative of density anomalies (rd) 312 IF( PRESENT( prd ) ) THEN ! depth of the partial step level 313 pgru(:,:)=0.0_wp ; pgrv(:,:)=0.0_wp ; 314 ! 315 DO jj = 1, jpjm1 316 DO ji = 1, jpim1 317 318 iku = mbku(ji,jj) 319 ikv = mbkv(ji,jj) 320 ze3wu = gdept_n(ji+1,jj,iku) - gdept_n(ji,jj,iku) 321 ze3wv = gdept_n(ji,jj+1,ikv) - gdept_n(ji,jj,ikv) 322 ! 323 IF( ze3wu >= 0._wp ) THEN ; zhi(ji,jj) = gdept_n(ji ,jj,iku) ! i-direction: case 1 324 ELSE ; zhi(ji,jj) = gdept_n(ji+1,jj,iku) ! - - case 2 325 ENDIF 326 IF( ze3wv >= 0._wp ) THEN ; zhj(ji,jj) = gdept_n(ji,jj ,ikv) ! j-direction: case 1 327 ELSE ; zhj(ji,jj) = gdept_n(ji,jj+1,ikv) ! - - case 2 328 ENDIF 329 330 END DO 331 END DO 332 333 ! Compute interpolated rd from zti, ztj for the 2 cases at the depth of the partial 334 ! step and store it in zri, zrj for each case 335 CALL eos( zti, zhi, zri ) 336 CALL eos( ztj, zhj, zrj ) 337 338 DO jj = 1, jpjm1 ! Gradient of density at the last level 339 DO ji = 1, jpim1 340 iku = mbku(ji,jj) 341 ikv = mbkv(ji,jj) 342 ze3wu = gdept_n(ji+1,jj,iku) - gdept_n(ji,jj,iku) 343 ze3wv = gdept_n(ji,jj+1,ikv) - gdept_n(ji,jj,ikv) 344 345 IF( ze3wu >= 0._wp ) THEN ; pgru(ji,jj) = ssumask(ji,jj) * ( zri(ji ,jj ) - prd(ji,jj,iku) ) ! i: 1 346 ELSE ; pgru(ji,jj) = ssumask(ji,jj) * ( prd(ji+1,jj,iku) - zri(ji,jj ) ) ! i: 2 347 ENDIF 348 IF( ze3wv >= 0._wp ) THEN ; pgrv(ji,jj) = ssvmask(ji,jj) * ( zrj(ji,jj ) - prd(ji,jj,ikv) ) ! j: 1 349 ELSE ; pgrv(ji,jj) = ssvmask(ji,jj) * ( prd(ji,jj+1,ikv) - zrj(ji,jj ) ) ! j: 2 350 ENDIF 351 352 END DO 353 END DO 354 355 CALL lbc_lnk( pgru , 'U', -1. ) ; CALL lbc_lnk( pgrv , 'V', -1. ) ! Lateral boundary conditions 356 ! 357 END IF 358 ! 359 ! !== (ISH) compute grui and gruvi ==! 360 ! 361 DO jn = 1, kjpt !== Interpolation of tracers at the last ocean level ==! ! 362 DO jj = 1, jpjm1 363 DO ji = 1, jpim1 364 iku = miku(ji,jj); ikup1 = miku(ji,jj) + 1 365 ikv = mikv(ji,jj); ikvp1 = mikv(ji,jj) + 1 366 ! 281 367 ! (ISF) case partial step top and bottom in adjacent cell in vertical 282 368 ! cannot used e3w because if 2 cell water column, we have ps at top and bottom 283 369 ! in this case e3w(i,j) - e3w(i,j+1) is not the distance between Tj~ and Tj 284 370 ! the only common depth between cells (i,j) and (i,j+1) is gdepw_0 285 ze3wu = (gdept_0(ji+1,jj,iku) - gdepw_0(ji+1,jj,iku)) - (gdept_0(ji,jj,iku) - gdepw_0(ji,jj,iku))286 ze3wv = (gdept_0(ji,jj+1,ikv) - gdepw_0(ji,jj+1,ikv)) - (gdept_0(ji,jj,ikv) - gdepw_0(ji,jj,ikv))287 ! 371 ze3wu = gdept_n(ji,jj,iku) - gdept_n(ji+1,jj,iku) 372 ze3wv = gdept_n(ji,jj,ikv) - gdept_n(ji,jj+1,ikv) 373 288 374 ! i- direction 289 375 IF( ze3wu >= 0._wp ) THEN ! case 1 290 zmaxu = ze3wu / fse3w(ji+1,jj,iku) 291 ! interpolated values of tracers 292 zti (ji,jj,jn) = pta(ji+1,jj,iku,jn) + zmaxu * ( pta(ji+1,jj,ikum1,jn) - pta(ji+1,jj,iku,jn) ) 376 zmaxu = ze3wu / e3w_n(ji+1,jj,ikup1) 377 ! interpolated values of tracers 378 zti(ji,jj,jn) = pta(ji+1,jj,iku,jn) + zmaxu * ( pta(ji+1,jj,ikup1,jn) - pta(ji+1,jj,iku,jn) ) 379 ! gradient of tracers 380 pgtui(ji,jj,jn) = ssumask(ji,jj) * ( zti(ji,jj,jn) - pta(ji,jj,iku,jn) ) 381 ELSE ! case 2 382 zmaxu = - ze3wu / e3w_n(ji,jj,ikup1) 383 ! interpolated values of tracers 384 zti(ji,jj,jn) = pta(ji,jj,iku,jn) + zmaxu * ( pta(ji,jj,ikup1,jn) - pta(ji,jj,iku,jn) ) 293 385 ! gradient of tracers 294 pgtu(ji,jj,jn) = umask(ji,jj,iku) * ( zti(ji,jj,jn) - pta(ji,jj,iku,jn) ) 295 ELSE ! case 2 296 zmaxu = -ze3wu / fse3w(ji,jj,iku) 297 ! interpolated values of tracers 298 zti (ji,jj,jn) = pta(ji,jj,iku,jn) + zmaxu * ( pta(ji,jj,ikum1,jn) - pta(ji,jj,iku,jn) ) 299 ! gradient of tracers 300 pgtu(ji,jj,jn) = umask(ji,jj,iku) * ( pta(ji+1,jj,iku,jn) - zti(ji,jj,jn) ) 386 pgtui(ji,jj,jn) = ssumask(ji,jj) * ( pta(ji+1,jj,iku,jn) - zti(ji,jj,jn) ) 301 387 ENDIF 302 388 ! 303 389 ! j- direction 304 390 IF( ze3wv >= 0._wp ) THEN ! case 1 305 zmaxv = ze3wv / fse3w(ji,jj+1,ikv) 306 ! interpolated values of tracers 307 ztj (ji,jj,jn) = pta(ji,jj+1,ikv,jn) + zmaxv * ( pta(ji,jj+1,ikvm1,jn) - pta(ji,jj+1,ikv,jn) ) 308 ! gradient of tracers 309 pgtv(ji,jj,jn) = vmask(ji,jj,ikv) * ( ztj(ji,jj,jn) - pta(ji,jj,ikv,jn) ) 310 ELSE ! case 2 311 zmaxv = -ze3wv / fse3w(ji,jj,ikv) 312 ! interpolated values of tracers 313 ztj (ji,jj,jn) = pta(ji,jj,ikv,jn) + zmaxv * ( pta(ji,jj,ikvm1,jn) - pta(ji,jj,ikv,jn) ) 314 ! gradient of tracers 315 pgtv(ji,jj,jn) = vmask(ji,jj,ikv) * ( pta(ji,jj+1,ikv,jn) - ztj(ji,jj,jn) ) 316 ENDIF 317 END DO 318 END DO 319 CALL lbc_lnk( pgtu(:,:,jn), 'U', -1. ) ; CALL lbc_lnk( pgtv(:,:,jn), 'V', -1. ) ! Lateral boundary cond. 391 zmaxv = ze3wv / e3w_n(ji,jj+1,ikvp1) 392 ! interpolated values of tracers 393 ztj(ji,jj,jn) = pta(ji,jj+1,ikv,jn) + zmaxv * ( pta(ji,jj+1,ikvp1,jn) - pta(ji,jj+1,ikv,jn) ) 394 ! gradient of tracers 395 pgtvi(ji,jj,jn) = ssvmask(ji,jj) * ( ztj(ji,jj,jn) - pta(ji,jj,ikv,jn) ) 396 ELSE ! case 2 397 zmaxv = - ze3wv / e3w_n(ji,jj,ikvp1) 398 ! interpolated values of tracers 399 ztj(ji,jj,jn) = pta(ji,jj,ikv,jn) + zmaxv * ( pta(ji,jj,ikvp1,jn) - pta(ji,jj,ikv,jn) ) 400 ! gradient of tracers 401 pgtvi(ji,jj,jn) = ssvmask(ji,jj) * ( pta(ji,jj+1,ikv,jn) - ztj(ji,jj,jn) ) 402 ENDIF 403 404 END DO 405 END DO 406 CALL lbc_lnk( pgtui(:,:,jn), 'U', -1. ); CALL lbc_lnk( pgtvi(:,:,jn), 'V', -1. ) ! Lateral boundary cond. 320 407 ! 321 408 END DO … … 323 410 IF( PRESENT( prd ) ) THEN !== horizontal derivative of density anomalies (rd) ==! (optional part) 324 411 ! 325 pgru (:,:)=0._wp ; pgrv (:,:) = 0._wp 326 pgzu (:,:)=0._wp ; pgzv (:,:) = 0._wp 327 pmru (:,:)=0._wp ; pmru (:,:) = 0._wp 328 pge3ru(:,:)=0._wp ; pge3rv(:,:) = 0._wp 329 ! 330 DO jj = 1, jpjm1 ! depth of the partial step level 331 DO ji = 1, jpim1 332 iku = mbku(ji,jj) 333 ikv = mbkv(ji,jj) 334 ze3wu = (gdept_0(ji+1,jj,iku) - gdepw_0(ji+1,jj,iku)) - (gdept_0(ji,jj,iku) - gdepw_0(ji,jj,iku)) 335 ze3wv = (gdept_0(ji,jj+1,ikv) - gdepw_0(ji,jj+1,ikv)) - (gdept_0(ji,jj,ikv) - gdepw_0(ji,jj,ikv)) 336 ! 337 IF( ze3wu >= 0._wp ) THEN ; zhi(ji,jj) = fsdept(ji+1,jj,iku) - ze3wu ! i-direction: case 1 338 ELSE ; zhi(ji,jj) = fsdept(ji ,jj,iku) + ze3wu ! - - case 2 339 ENDIF 340 IF( ze3wv >= 0._wp ) THEN ; zhj(ji,jj) = fsdept(ji,jj+1,ikv) - ze3wv ! j-direction: case 1 341 ELSE ; zhj(ji,jj) = fsdept(ji,jj ,ikv) + ze3wv ! - - case 2 342 ENDIF 412 pgrui(:,:) =0.0_wp; pgrvi(:,:) =0.0_wp; 413 DO jj = 1, jpjm1 414 DO ji = 1, jpim1 415 416 iku = miku(ji,jj) 417 ikv = mikv(ji,jj) 418 ze3wu = gdept_n(ji,jj,iku) - gdept_n(ji+1,jj,iku) 419 ze3wv = gdept_n(ji,jj,ikv) - gdept_n(ji,jj+1,ikv) 420 ! 421 IF( ze3wu >= 0._wp ) THEN ; zhi(ji,jj) = gdept_n(ji ,jj,iku) ! i-direction: case 1 422 ELSE ; zhi(ji,jj) = gdept_n(ji+1,jj,iku) ! - - case 2 423 ENDIF 424 425 IF( ze3wv >= 0._wp ) THEN ; zhj(ji,jj) = gdept_n(ji,jj ,ikv) ! j-direction: case 1 426 ELSE ; zhj(ji,jj) = gdept_n(ji,jj+1,ikv) ! - - case 2 427 ENDIF 428 343 429 END DO 344 430 END DO … … 346 432 CALL eos( zti, zhi, zri ) ! interpolated density from zti, ztj 347 433 CALL eos( ztj, zhj, zrj ) ! at the partial step depth output in zri, zrj 348 434 ! 349 435 DO jj = 1, jpjm1 ! Gradient of density at the last level 350 436 DO ji = 1, jpim1 351 iku = mbku(ji,jj) ; ikum1 = MAX( iku - 1 , 1 ) ! last and before last ocean level at u- & v-points 352 ikv = mbkv(ji,jj) ; ikvm1 = MAX( ikv - 1 , 1 ) ! last and before last ocean level at u- & v-points 353 ze3wu = (gdept_0(ji+1,jj,iku) - gdepw_0(ji+1,jj,iku)) - (gdept_0(ji,jj,iku) - gdepw_0(ji,jj,iku)) 354 ze3wv = (gdept_0(ji,jj+1,ikv) - gdepw_0(ji,jj+1,ikv)) - (gdept_0(ji,jj,ikv) - gdepw_0(ji,jj,ikv)) 355 IF( ze3wu >= 0._wp ) THEN 356 pgzu(ji,jj) = (fsde3w(ji+1,jj,iku) - ze3wu) - fsde3w(ji,jj,iku) 357 pgru(ji,jj) = umask(ji,jj,iku) * ( zri(ji ,jj) - prd(ji,jj,iku) ) ! i: 1 358 pmru(ji,jj) = umask(ji,jj,iku) * ( zri(ji ,jj) + prd(ji,jj,iku) ) ! i: 1 359 pge3ru(ji,jj) = umask(ji,jj,iku) & 360 * ( (fse3w(ji+1,jj,iku) - ze3wu )* ( zri(ji ,jj ) + prd(ji+1,jj,ikum1) + 2._wp) & 361 - fse3w(ji ,jj,iku) * ( prd(ji ,jj,iku) + prd(ji ,jj,ikum1) + 2._wp) ) ! j: 2 362 ELSE 363 pgzu(ji,jj) = fsde3w(ji+1,jj,iku) - (fsde3w(ji,jj,iku) + ze3wu) 364 pgru(ji,jj) = umask(ji,jj,iku) * ( prd(ji+1,jj,iku) - zri(ji,jj) ) ! i: 2 365 pmru(ji,jj) = umask(ji,jj,iku) * ( prd(ji+1,jj,iku) + zri(ji,jj) ) ! i: 2 366 pge3ru(ji,jj) = umask(ji,jj,iku) & 367 * ( fse3w(ji+1,jj,iku) * ( prd(ji+1,jj,iku) + prd(ji+1,jj,ikum1) + 2._wp) & 368 -(fse3w(ji ,jj,iku) + ze3wu) * ( zri(ji ,jj ) + prd(ji ,jj,ikum1) + 2._wp) ) ! j: 2 369 ENDIF 370 IF( ze3wv >= 0._wp ) THEN 371 pgzv(ji,jj) = (fsde3w(ji,jj+1,ikv) - ze3wv) - fsde3w(ji,jj,ikv) 372 pgrv(ji,jj) = vmask(ji,jj,ikv) * ( zrj(ji,jj ) - prd(ji,jj,ikv) ) ! j: 1 373 pmrv(ji,jj) = vmask(ji,jj,ikv) * ( zrj(ji,jj ) + prd(ji,jj,ikv) ) ! j: 1 374 pge3rv(ji,jj) = vmask(ji,jj,ikv) & 375 * ( (fse3w(ji,jj+1,ikv) - ze3wv )* ( zrj(ji,jj ) + prd(ji,jj+1,ikvm1) + 2._wp) & 376 - fse3w(ji,jj ,ikv) * ( prd(ji,jj ,ikv) + prd(ji,jj ,ikvm1) + 2._wp) ) ! j: 2 377 ELSE 378 pgzv(ji,jj) = fsde3w(ji,jj+1,ikv) - (fsde3w(ji,jj,ikv) + ze3wv) 379 pgrv(ji,jj) = vmask(ji,jj,ikv) * ( prd(ji,jj+1,ikv) - zrj(ji,jj) ) ! j: 2 380 pmrv(ji,jj) = vmask(ji,jj,ikv) * ( prd(ji,jj+1,ikv) + zrj(ji,jj) ) ! j: 2 381 pge3rv(ji,jj) = vmask(ji,jj,ikv) & 382 * ( fse3w(ji,jj+1,ikv) * ( prd(ji,jj+1,ikv) + prd(ji,jj+1,ikvm1) + 2._wp) & 383 -(fse3w(ji,jj ,ikv) + ze3wv) * ( zrj(ji,jj ) + prd(ji,jj ,ikvm1) + 2._wp) ) ! j: 2 384 ENDIF 385 END DO 386 END DO 387 CALL lbc_lnk( pgru , 'U', -1. ) ; CALL lbc_lnk( pgrv , 'V', -1. ) ! Lateral boundary conditions 388 CALL lbc_lnk( pmru , 'U', 1. ) ; CALL lbc_lnk( pmrv , 'V', 1. ) ! Lateral boundary conditions 389 CALL lbc_lnk( pgzu , 'U', -1. ) ; CALL lbc_lnk( pgzv , 'V', -1. ) ! Lateral boundary conditions 390 CALL lbc_lnk( pge3ru , 'U', -1. ) ; CALL lbc_lnk( pge3rv , 'V', -1. ) ! Lateral boundary conditions 391 ! 392 END IF 393 ! 394 ! !== (ISH) compute grui and gruvi ==! 395 ! 396 DO jn = 1, kjpt !== Interpolation of tracers at the last ocean level ==! ! 397 DO jj = 1, jpjm1 398 DO ji = 1, jpim1 399 iku = miku(ji,jj) ; ikup1 = miku(ji,jj) + 1 400 ikv = mikv(ji,jj) ; ikvp1 = mikv(ji,jj) + 1 401 ! 402 ! (ISF) case partial step top and bottom in adjacent cell in vertical 403 ! cannot used e3w because if 2 cell water column, we have ps at top and bottom 404 ! in this case e3w(i,j) - e3w(i,j+1) is not the distance between Tj~ and Tj 405 ! the only common depth between cells (i,j) and (i,j+1) is gdepw_0 406 ze3wu = (gdepw_0(ji+1,jj,iku+1) - gdept_0(ji+1,jj,iku)) - (gdepw_0(ji,jj,iku+1) - gdept_0(ji,jj,iku)) 407 ze3wv = (gdepw_0(ji,jj+1,ikv+1) - gdept_0(ji,jj+1,ikv)) - (gdepw_0(ji,jj,ikv+1) - gdept_0(ji,jj,ikv)) 408 ! i- direction 409 IF( ze3wu >= 0._wp ) THEN ! case 1 410 zmaxu = ze3wu / fse3w(ji+1,jj,iku+1) 411 ! interpolated values of tracers 412 zti(ji,jj,jn) = pta(ji+1,jj,iku,jn) + zmaxu * ( pta(ji+1,jj,iku+1,jn) - pta(ji+1,jj,iku,jn) ) 413 ! gradient of tracers 414 pgtui(ji,jj,jn) = umask(ji,jj,iku) * ( zti(ji,jj,jn) - pta(ji,jj,iku,jn) ) 415 ELSE ! case 2 416 zmaxu = - ze3wu / fse3w(ji,jj,iku+1) 417 ! interpolated values of tracers 418 zti(ji,jj,jn) = pta(ji,jj,iku,jn) + zmaxu * ( pta(ji,jj,iku+1,jn) - pta(ji,jj,iku,jn) ) 419 ! gradient of tracers 420 pgtui(ji,jj,jn) = umask(ji,jj,iku) * ( pta(ji+1,jj,iku,jn) - zti(ji,jj,jn) ) 421 ENDIF 422 ! 423 ! j- direction 424 IF( ze3wv >= 0._wp ) THEN ! case 1 425 zmaxv = ze3wv / fse3w(ji,jj+1,ikv+1) 426 ! interpolated values of tracers 427 ztj(ji,jj,jn) = pta(ji,jj+1,ikv,jn) + zmaxv * ( pta(ji,jj+1,ikv+1,jn) - pta(ji,jj+1,ikv,jn) ) 428 ! gradient of tracers 429 pgtvi(ji,jj,jn) = vmask(ji,jj,ikv) * ( ztj(ji,jj,jn) - pta(ji,jj,ikv,jn) ) 430 ELSE ! case 2 431 zmaxv = - ze3wv / fse3w(ji,jj,ikv+1) 432 ! interpolated values of tracers 433 ztj(ji,jj,jn) = pta(ji,jj,ikv,jn) + zmaxv * ( pta(ji,jj,ikv+1,jn) - pta(ji,jj,ikv,jn) ) 434 ! gradient of tracers 435 pgtvi(ji,jj,jn) = vmask(ji,jj,ikv) * ( pta(ji,jj+1,ikv,jn) - ztj(ji,jj,jn) ) 436 ENDIF 437 END DO!! 438 END DO!! 439 CALL lbc_lnk( pgtui(:,:,jn), 'U', -1. ) ; CALL lbc_lnk( pgtvi(:,:,jn), 'V', -1. ) ! Lateral boundary cond. 440 ! 441 END DO 442 443 IF( PRESENT( prd ) ) THEN !== horizontal derivative of density anomalies (rd) ==! (optional part) 444 ! 445 pgrui(:,:) =0.0_wp ; pgrvi(:,:) =0.0_wp ; 446 pgzui(:,:) =0.0_wp ; pgzvi(:,:) =0.0_wp ; 447 pmrui(:,:) =0.0_wp ; pmrui(:,:) =0.0_wp ; 448 pge3rui(:,:)=0.0_wp ; pge3rvi(:,:)=0.0_wp ; 449 ! 450 DO jj = 1, jpjm1 ! depth of the partial step level 451 DO ji = 1, jpim1 452 iku = miku(ji,jj) 453 ikv = mikv(ji,jj) 454 ze3wu = (gdepw_0(ji+1,jj,iku+1) - gdept_0(ji+1,jj,iku)) - (gdepw_0(ji,jj,iku+1) - gdept_0(ji,jj,iku)) 455 ze3wv = (gdepw_0(ji,jj+1,ikv+1) - gdept_0(ji,jj+1,ikv)) - (gdepw_0(ji,jj,ikv+1) - gdept_0(ji,jj,ikv)) 456 ! 457 IF( ze3wu >= 0._wp ) THEN ; zhi(ji,jj) = fsdept(ji+1,jj,iku) + ze3wu ! i-direction: case 1 458 ELSE ; zhi(ji,jj) = fsdept(ji ,jj,iku) - ze3wu ! - - case 2 459 ENDIF 460 IF( ze3wv >= 0._wp ) THEN ; zhj(ji,jj) = fsdept(ji,jj+1,ikv) + ze3wv ! j-direction: case 1 461 ELSE ; zhj(ji,jj) = fsdept(ji,jj ,ikv) - ze3wv ! - - case 2 462 ENDIF 463 END DO 464 END DO 465 ! 466 CALL eos( zti, zhi, zri ) ! interpolated density from zti, ztj 467 CALL eos( ztj, zhj, zrj ) ! at the partial step depth output in zri, zrj 468 ! 469 DO jj = 1, jpjm1 ! Gradient of density at the last level 470 DO ji = 1, jpim1 471 iku = miku(ji,jj) ; ikup1 = miku(ji,jj) + 1 472 ikv = mikv(ji,jj) ; ikvp1 = mikv(ji,jj) + 1 473 ze3wu = (gdepw_0(ji+1,jj,iku+1) - gdept_0(ji+1,jj,iku)) - (gdepw_0(ji,jj,iku+1) - gdept_0(ji,jj,iku)) 474 ze3wv = (gdepw_0(ji,jj+1,ikv+1) - gdept_0(ji,jj+1,ikv)) - (gdepw_0(ji,jj,ikv+1) - gdept_0(ji,jj,ikv)) 475 IF( ze3wu >= 0._wp ) THEN 476 pgzui (ji,jj) = (fsde3w(ji+1,jj,iku) + ze3wu) - fsde3w(ji,jj,iku) 477 pgrui (ji,jj) = umask(ji,jj,iku) * ( zri(ji,jj) - prd(ji,jj,iku) ) ! i: 1 478 pmrui (ji,jj) = umask(ji,jj,iku) * ( zri(ji,jj) + prd(ji,jj,iku) ) ! i: 1 479 pge3rui(ji,jj) = umask(ji,jj,iku+1) & 480 & * ( (fse3w(ji+1,jj,iku+1) - ze3wu) * (zri(ji,jj ) + prd(ji+1,jj,iku+1) + 2._wp) & 481 & - fse3w(ji ,jj,iku+1) * (prd(ji,jj,iku) + prd(ji ,jj,iku+1) + 2._wp) ) ! i: 1 482 ELSE 483 pgzui (ji,jj) = fsde3w(ji+1,jj,iku) - (fsde3w(ji,jj,iku) - ze3wu) 484 pgrui (ji,jj) = umask(ji,jj,iku) * ( prd(ji+1,jj,iku) - zri(ji,jj) ) ! i: 2 485 pmrui (ji,jj) = umask(ji,jj,iku) * ( prd(ji+1,jj,iku) + zri(ji,jj) ) ! i: 2 486 pge3rui(ji,jj) = umask(ji,jj,iku+1) & 487 & * ( fse3w(ji+1,jj,iku+1) * (prd(ji+1,jj,iku) + prd(ji+1,jj,iku+1) + 2._wp) & 488 & -(fse3w(ji ,jj,iku+1) + ze3wu) * (zri(ji,jj ) + prd(ji ,jj,iku+1) + 2._wp) ) ! i: 2 489 ENDIF 490 IF( ze3wv >= 0._wp ) THEN 491 pgzvi (ji,jj) = (fsde3w(ji,jj+1,ikv) + ze3wv) - fsde3w(ji,jj,ikv) 492 pgrvi (ji,jj) = vmask(ji,jj,ikv) * ( zrj(ji,jj ) - prd(ji,jj,ikv) ) ! j: 1 493 pmrvi (ji,jj) = vmask(ji,jj,ikv) * ( zrj(ji,jj ) + prd(ji,jj,ikv) ) ! j: 1 494 pge3rvi(ji,jj) = vmask(ji,jj,ikv+1) & 495 & * ( (fse3w(ji,jj+1,ikv+1) - ze3wv) * ( zrj(ji,jj ) + prd(ji,jj+1,ikv+1) + 2._wp) & 496 & - fse3w(ji,jj ,ikv+1) * ( prd(ji,jj,ikv) + prd(ji,jj ,ikv+1) + 2._wp) ) ! j: 1 497 ! + 2 due to the formulation in density and not in anomalie in hpg sco 498 ELSE 499 pgzvi (ji,jj) = fsde3w(ji,jj+1,ikv) - (fsde3w(ji,jj,ikv) - ze3wv) 500 pgrvi (ji,jj) = vmask(ji,jj,ikv) * ( prd(ji,jj+1,ikv) - zrj(ji,jj) ) ! j: 2 501 pmrvi (ji,jj) = vmask(ji,jj,ikv) * ( prd(ji,jj+1,ikv) + zrj(ji,jj) ) ! j: 2 502 pge3rvi(ji,jj) = vmask(ji,jj,ikv+1) & 503 & * ( fse3w(ji,jj+1,ikv+1) * ( prd(ji,jj+1,ikv) + prd(ji,jj+1,ikv+1) + 2._wp) & 504 & -(fse3w(ji,jj ,ikv+1) + ze3wv) * ( zrj(ji,jj ) + prd(ji,jj ,ikv+1) + 2._wp) ) ! j: 2 505 ENDIF 506 END DO 507 END DO 508 CALL lbc_lnk( pgrui , 'U', -1. ) ; CALL lbc_lnk( pgrvi , 'V', -1. ) ! Lateral boundary conditions 509 CALL lbc_lnk( pmrui , 'U', 1. ) ; CALL lbc_lnk( pmrvi , 'V', 1. ) ! Lateral boundary conditions 510 CALL lbc_lnk( pgzui , 'U', -1. ) ; CALL lbc_lnk( pgzvi , 'V', -1. ) ! Lateral boundary conditions 511 CALL lbc_lnk( pge3rui , 'U', -1. ) ; CALL lbc_lnk( pge3rvi , 'V', -1. ) ! Lateral boundary conditions 437 iku = miku(ji,jj) 438 ikv = mikv(ji,jj) 439 ze3wu = gdept_n(ji,jj,iku) - gdept_n(ji+1,jj,iku) 440 ze3wv = gdept_n(ji,jj,ikv) - gdept_n(ji,jj+1,ikv) 441 442 IF( ze3wu >= 0._wp ) THEN ; pgrui(ji,jj) = ssumask(ji,jj) * ( zri(ji ,jj ) - prd(ji,jj,iku) ) ! i: 1 443 ELSE ; pgrui(ji,jj) = ssumask(ji,jj) * ( prd(ji+1,jj ,iku) - zri(ji,jj ) ) ! i: 2 444 ENDIF 445 IF( ze3wv >= 0._wp ) THEN ; pgrvi(ji,jj) = ssvmask(ji,jj) * ( zrj(ji ,jj ) - prd(ji,jj,ikv) ) ! j: 1 446 ELSE ; pgrvi(ji,jj) = ssvmask(ji,jj) * ( prd(ji ,jj+1,ikv) - zrj(ji,jj ) ) ! j: 2 447 ENDIF 448 449 END DO 450 END DO 451 CALL lbc_lnk( pgrui , 'U', -1. ); CALL lbc_lnk( pgrvi , 'V', -1. ) ! Lateral boundary conditions 512 452 ! 513 453 END IF
Note: See TracChangeset
for help on using the changeset viewer.