Changeset 2450 for branches/nemo_v3_3_beta
- Timestamp:
- 2010-12-04T16:20:50+01:00 (14 years ago)
- Location:
- branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC
- Files:
-
- 23 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/DIA/diaar5.F90
r2287 r2450 16 16 USE oce ! ocean dynamics and active tracers 17 17 USE dom_oce ! ocean space and time domain 18 USE eosbn2 18 USE eosbn2 ! equation of state (eos_bn2 routine) 19 19 USE lib_mpp ! distribued memory computing library 20 20 USE iom ! I/O manager library … … 28 28 LOGICAL, PUBLIC, PARAMETER :: lk_diaar5 = .TRUE. ! coupled flag 29 29 30 REAL(wp) :: vol0 31 REAL(wp) :: area_tot 32 REAL(wp), DIMENSION(jpi,jpj ) :: area 33 REAL(wp), DIMENSION(jpi,jpj ) :: thick0 34 REAL(wp), DIMENSION(jpi,jpj,jpk) :: sn0 30 REAL(wp) :: vol0 ! ocean volume (interior domain) 31 REAL(wp) :: area_tot ! total ocean surface (interior domain) 32 REAL(wp), DIMENSION(jpi,jpj ) :: area ! cell surface (interior domain) 33 REAL(wp), DIMENSION(jpi,jpj ) :: thick0 ! ocean thickness (interior domain) 34 REAL(wp), DIMENSION(jpi,jpj,jpk) :: sn0 ! initial salinity 35 35 36 36 !! * Substitutions … … 39 39 !! NEMO/OPA 3.3 , NEMO Consortium (2010) 40 40 !! $Id$ 41 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 42 !!---------------------------------------------------------------------- 43 41 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 42 !!---------------------------------------------------------------------- 44 43 CONTAINS 45 44 … … 77 76 CALL eos( ztsn, zrhd ) ! now in situ density using initial salinity 78 77 ! 79 zbotpres(:,:) = 0. e0! no atmospheric surface pressure, levitating sea-ice78 zbotpres(:,:) = 0._wp ! no atmospheric surface pressure, levitating sea-ice 80 79 DO jk = 1, jpkm1 81 80 zbotpres(:,:) = zbotpres(:,:) + fse3t(:,:,jk) * zrhd(:,:,jk) 82 81 END DO 83 IF( .NOT. 82 IF( .NOT.lk_vvl ) zbotpres(:,:) = zbotpres(:,:) + sshn(:,:) * zrhd(:,:,1) 84 83 ! 85 84 zarho = SUM( area(:,:) * zbotpres(:,:) ) … … 90 89 ! ! steric sea surface height 91 90 CALL eos( tsn, zrhd, zrhop ) ! now in situ and potential density 92 zrhop(:,:,jpk) = 0. e091 zrhop(:,:,jpk) = 0._wp 93 92 CALL iom_put( 'rhop', zrhop ) 94 93 ! 95 zbotpres(:,:) = 0. e0! no atmospheric surface pressure, levitating sea-ice94 zbotpres(:,:) = 0._wp ! no atmospheric surface pressure, levitating sea-ice 96 95 DO jk = 1, jpkm1 97 96 zbotpres(:,:) = zbotpres(:,:) + fse3t(:,:,jk) * zrhd(:,:,jk) 98 97 END DO 99 IF( .NOT. 98 IF( .NOT.lk_vvl ) zbotpres(:,:) = zbotpres(:,:) + sshn(:,:) * zrhd(:,:,1) 100 99 ! 101 100 zarho = SUM( area(:,:) * zbotpres(:,:) ) … … 105 104 106 105 ! ! ocean bottom pressure 107 zztmp = rau0 * grav * 1.e-4 106 zztmp = rau0 * grav * 1.e-4_wp ! recover pressure from pressure anomaly and cover to dbar = 1.e4 Pa 108 107 zbotpres(:,:) = zztmp * ( zbotpres(:,:) + sshn(:,:) + thick0(:,:) ) 109 108 CALL iom_put( 'botpres', zbotpres ) 110 109 111 110 ! ! Mean density anomalie, temperature and salinity 112 ztemp = 0. e0113 zsal = 0. e0111 ztemp = 0._wp 112 zsal = 0._wp 114 113 DO jk = 1, jpkm1 115 114 DO jj = 1, jpj 116 115 DO ji = 1, jpi 117 116 zztmp = area(ji,jj) * fse3t(ji,jj,jk) 118 ztemp = ztemp + zztmp * tn 119 zsal = zsal + zztmp * sn 117 ztemp = ztemp + zztmp * tn(ji,jj,jk) 118 zsal = zsal + zztmp * sn(ji,jj,jk) 120 119 END DO 121 120 END DO 122 121 END DO 123 IF( .NOT. 124 ztemp = ztemp + SUM( zarea_ssh(:,:) * tn 125 zsal = zsal + SUM( zarea_ssh(:,:) * sn 122 IF( .NOT.lk_vvl ) THEN 123 ztemp = ztemp + SUM( zarea_ssh(:,:) * tn(:,:,1) ) 124 zsal = zsal + SUM( zarea_ssh(:,:) * sn(:,:,1) ) 126 125 ENDIF 127 126 IF( lk_mpp ) THEN … … 146 145 !! 147 146 !! ** Purpose : initialization for AR5 diagnostic computation 148 !!149 147 !!---------------------------------------------------------------------- 150 148 INTEGER :: inum … … 159 157 area_tot = SUM( area(:,:) ) ; IF( lk_mpp ) CALL mpp_sum( area_tot ) 160 158 161 vol0 = 0. e0162 thick0(:,:) = 0. e0159 vol0 = 0._wp 160 thick0(:,:) = 0._wp 163 161 DO jk = 1, jpkm1 164 162 vol0 = vol0 + SUM( area (:,:) * tmask(:,:,jk) * fse3t_0(:,:,jk) ) … … 171 169 CALL iom_get ( inum, jpdom_data, 'vosaline', zsaldta(:,:,:,2), 12 ) 172 170 CALL iom_close( inum ) 173 sn0(:,:,:) = 0.5 * ( zsaldta(:,:,:,1) + zsaldta(:,:,:,2) )174 sn0(:,:,:) = sn0(:,:,:) *tmask(:,:,:)171 sn0(:,:,:) = 0.5_wp * ( zsaldta(:,:,:,1) + zsaldta(:,:,:,2) ) 172 sn0(:,:,:) = sn0(:,:,:) * tmask(:,:,:) 175 173 IF( ln_zps ) THEN ! z-coord. partial steps 176 174 DO jj = 1, jpj ! interpolation of salinity at the last ocean level (i.e. the partial step) 177 175 DO ji = 1, jpi 178 ik = mb athy(ji,jj) - 1179 IF( ik > 2) THEN176 ik = mbkt(ji,jj) 177 IF( ik > 1 ) THEN 180 178 zztmp = ( gdept_0(ik) - fsdept_0(ji,jj,ik) ) / ( gdept_0(ik) - gdept_0(ik-1) ) 181 sn0(ji,jj,ik) = ( 1.-zztmp) * sn0(ji,jj,ik) + zztmp * sn0(ji,jj,ik-1)179 sn0(ji,jj,ik) = ( 1._wp - zztmp ) * sn0(ji,jj,ik) + zztmp * sn0(ji,jj,ik-1) 182 180 ENDIF 183 181 END DO -
branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/DIA/diahth.F90
r2287 r2450 11 11 !! NEMO 3.2 ! 2009-07 (S. Masson) hc300 bugfix + cleaning + add new diag 12 12 !!---------------------------------------------------------------------- 13 14 13 #if defined key_diahth || defined key_esopa 15 14 !!---------------------------------------------------------------------- … … 18 17 !! dia_hth : Compute varius diagnostics associated with the mixed layer 19 18 !!---------------------------------------------------------------------- 20 !! * Modules used21 19 USE oce ! ocean dynamics and tracers 22 20 USE dom_oce ! ocean space and time domain 23 21 USE phycst ! physical constants 24 22 USE in_out_manager ! I/O manager 25 USE iom 23 USE iom ! I/O library 26 24 27 25 IMPLICIT NONE 28 26 PRIVATE 29 27 30 !! * Routine accessibility 31 PUBLIC dia_hth ! routine called by step.F90 32 33 !! * Shared module variables 28 PUBLIC dia_hth ! routine called by step.F90 29 34 30 LOGICAL , PUBLIC, PARAMETER :: lk_diahth = .TRUE. !: thermocline-20d depths flag 35 31 ! note: following variables should move to local variables once iom_put is always used … … 44 40 !! NEMO/OPA 3.3 , NEMO Consortium (2010) 45 41 !! $Id$ 46 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 47 !!---------------------------------------------------------------------- 48 42 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 43 !!---------------------------------------------------------------------- 49 44 CONTAINS 50 45 … … 68 63 !! 69 64 !! ** Method : 70 !!71 65 !!------------------------------------------------------------------- 72 66 INTEGER, INTENT( in ) :: kt ! ocean time-step index 73 67 !! 74 68 INTEGER :: ji, jj, jk ! dummy loop arguments 75 INTEGER :: iid, i if, ilevel! temporary integers69 INTEGER :: iid, ilevel ! temporary integers 76 70 INTEGER, DIMENSION(jpi,jpj) :: ik20, ik28 ! levels 77 71 REAL(wp) :: zavt5 = 5.e-4_wp ! Kz criterion for the turbocline depth … … 103 97 104 98 ! initialization 105 ztinv (:,:) = 0. e0_wp106 zdepinv(:,:) = 0. e0_wp107 zmaxdzT(:,:) = 0. e0_wp99 ztinv (:,:) = 0._wp 100 zdepinv(:,:) = 0._wp 101 zmaxdzT(:,:) = 0._wp 108 102 DO jj = 1, jpj 109 103 DO ji = 1, jpi … … 128 122 ! Preliminary computation 129 123 ! computation of zdelr = (dr/dT)(T,S,10m)*(-0.2 degC) 130 DO jj =1, jpj131 DO ji =1, jpi124 DO jj = 1, jpj 125 DO ji = 1, jpi 132 126 IF( tmask(ji,jj,nla10) == 1. ) THEN 133 127 zu = 1779.50 + 11.250*tn(ji,jj,nla10) - 3.80*sn(ji,jj,nla10) - 0.0745*tn(ji,jj,nla10)*tn(ji,jj,nla10) & … … 139 133 zdelr(ji,jj) = ztem2 * (1000.*(zut*zv - zvt*zu)/zw) 140 134 ELSE 141 zdelr(ji,jj) = 0. e0135 zdelr(ji,jj) = 0._wp 142 136 ENDIF 143 137 END DO … … 153 147 DO jj = 1, jpj 154 148 DO ji = 1, jpi 155 149 ! 156 150 zzdep = fsdepw(ji,jj,jk) 157 151 zztmp = ( tn(ji,jj,jk-1) - tn(ji,jj,jk) ) / zzdep * tmask(ji,jj,jk) ! vertical gradient of temperature (dT/dz) … … 189 183 DO jj = 1, jpj 190 184 DO ji = 1, jpi 191 185 ! 192 186 zzdep = fsdepw(ji,jj,jk) * tmask(ji,jj,1) 193 187 ! 194 188 zztmp = tn(ji,jj,nla10) - tn(ji,jj,jk) ! - delta T(10m) 195 189 IF( ABS(zztmp) > ztem2 ) zabs2 (ji,jj) = zzdep ! abs > 0.2 … … 203 197 IF( zztmp > zrho3 ) zrho10_3(ji,jj) = zzdep ! > 0.03 204 198 IF( zztmp > zdelr(ji,jj) ) zpycn (ji,jj) = zzdep ! > equi. delta T(10m) - 0.2 205 199 ! 206 200 END DO 207 201 END DO … … 237 231 DO jj = 1, jpj 238 232 DO ji = 1, jpi 239 240 iif = mbathy(ji,jj) 241 zzdep = fsdepw(ji,jj,iif) 242 233 ! 234 zzdep = fsdepw(ji,jj,mbkt(ji,jj)+1) ! depth of the oean bottom 235 ! 243 236 iid = ik20(ji,jj) 244 237 IF( iid /= 1 ) THEN 245 ! linear interpolation 246 zztmp = fsdept(ji,jj,iid ) & 238 zztmp = fsdept(ji,jj,iid ) & ! linear interpolation 247 239 & + ( fsdept(ji,jj,iid+1) - fsdept(ji,jj,iid) ) & 248 240 & * ( 20.*tmask(ji,jj,iid+1) - tn(ji,jj,iid) ) & 249 241 & / ( tn(ji,jj,iid+1) - tn(ji,jj,iid) + (1.-tmask(ji,jj,1)) ) 250 ! bound by the ocean depth, minimum value, first T-point depth 251 hd20(ji,jj) = MIN( zztmp*tmask(ji,jj,1), zzdep) 242 hd20(ji,jj) = MIN( zztmp , zzdep) * tmask(ji,jj,1) ! bound by the ocean depth 252 243 ELSE 253 hd20(ji,jj) =0.244 hd20(ji,jj) = 0._wp 254 245 ENDIF 255 246 ! 256 247 iid = ik28(ji,jj) 257 248 IF( iid /= 1 ) THEN 258 ! linear interpolation 259 zztmp = fsdept(ji,jj,iid ) & 249 zztmp = fsdept(ji,jj,iid ) & ! linear interpolation 260 250 & + ( fsdept(ji,jj,iid+1) - fsdept(ji,jj,iid) ) & 261 251 & * ( 28.*tmask(ji,jj,iid+1) - tn(ji,jj,iid) ) & 262 252 & / ( tn(ji,jj,iid+1) - tn(ji,jj,iid) + (1.-tmask(ji,jj,1)) ) 263 ! bound by the ocean depth, minimum value, first T-point depth 264 hd28(ji,jj) = MIN( zztmp*tmask(ji,jj,1), zzdep ) 253 hd28(ji,jj) = MIN( zztmp , zzdep ) * tmask(ji,jj,1) ! bound by the ocean depth 265 254 ELSE 266 hd28(ji,jj) = 0. 255 hd28(ji,jj) = 0._wp 267 256 ENDIF 268 257 … … 277 266 278 267 ! find ilevel with (ilevel+1) the deepest W-level above 300m (we assume we can use e3t_0 to do this search...) 279 ilevel = 0280 zthick_0 = 0. e0_wp268 ilevel = 0 269 zthick_0 = 0._wp 281 270 DO jk = 1, jpkm1 282 271 zthick_0 = zthick_0 + e3t_0(jk) … … 284 273 END DO 285 274 ! surface boundary condition 286 IF( lk_vvl ) THEN ; zthick(:,:) = 0. e0_wp ; htc3(:,:) = 0.e0_wp275 IF( lk_vvl ) THEN ; zthick(:,:) = 0._wp ; htc3(:,:) = 0._wp 287 276 ELSE ; zthick(:,:) = sshn(:,:) ; htc3(:,:) = tn(:,:,jk) * sshn(:,:) * tmask(:,:,jk) 288 277 ENDIF … … 303 292 htc3(:,:) = zcoef * htc3(:,:) 304 293 CALL iom_put( "hc300", htc3 ) ! first 300m heat content 305 306 294 ! 307 295 END SUBROUTINE dia_hth 308 296 -
branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/DOM/domwri.F90
r2287 r2450 4 4 !! Ocean initialization : write the ocean domain mesh ask file(s) 5 5 !!====================================================================== 6 !! History : OPA ! 1997-02 (G. Madec) Original code 7 !! 8.1 ! 1999-11 (M. Imbard) NetCDF FORMAT with IOIPSL 8 !! NEMO 1.0 ! 2002-08 (G. Madec) F90 and several file 9 !!---------------------------------------------------------------------- 6 10 7 11 !!---------------------------------------------------------------------- … … 11 15 !! = 3 : mesh_hgr, mesh_zgr and mask 12 16 !!---------------------------------------------------------------------- 13 !! * Modules used14 17 USE dom_oce ! ocean space and time domain 15 USE in_out_manager 16 USE iom 17 USE lbclnk ! ocean lateral boundary conditions (or mpp link)18 USE lib_mpp 18 USE in_out_manager ! I/O manager 19 USE iom ! I/O library 20 USE lbclnk ! lateral boundary conditions - mpp exchanges 21 USE lib_mpp ! MPP library 19 22 20 23 IMPLICIT NONE … … 28 31 !! NEMO/OPA 3.3 , NEMO Consortium (2010) 29 32 !! $Id$ 30 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 31 !!---------------------------------------------------------------------- 32 33 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 34 !!---------------------------------------------------------------------- 33 35 CONTAINS 34 36 … … 58 60 !! thickness of the bottom points hdep[tw] and e3[tw]_ps 59 61 !! 60 !! ** output file : 61 !! meshmask.nc : domain size, horizontal grid-point position, 62 !! masks, depth and vertical scale factors 63 !! 64 !! History : 65 !! ! 97-02 (G. Madec) Original code 66 !! ! 99-11 (M. Imbard) NetCDF FORMAT with IOIPSL 67 !! 9.0 ! 02-08 (G. Madec) F90 and several file 62 !! ** output file : meshmask.nc : domain size, horizontal grid-point position, 63 !! masks, depth and vertical scale factors 68 64 !!---------------------------------------------------------------------- 69 65 INTEGER :: inum0 ! temprary units for 'mesh_mask.nc' file … … 72 68 INTEGER :: inum3 ! temprary units for 'mesh_hgr.nc' file 73 69 INTEGER :: inum4 ! temprary units for 'mesh_zgr.nc' file 74 INTEGER :: ji, jj, jk , ik70 INTEGER :: ji, jj, jk ! dummy loop indices 75 71 REAL(wp), DIMENSION(jpi,jpj) :: zprt ! temporary array for bathymetry 76 72 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdepu ! 3D depth of U point … … 118 114 CALL iom_open( TRIM(clnam3), inum3, ldwrt = .TRUE., kiolib = jprstlib ) 119 115 CALL iom_open( TRIM(clnam4), inum4, ldwrt = .TRUE., kiolib = jprstlib ) 120 116 ! 121 117 END SELECT 122 118 … … 160 156 CALL iom_rstput( 0, 0, inum3, 'ff', ff, ktype = jp_r8 ) ! ! coriolis factor 161 157 162 ! note that mbathy has been modified in dommsk or in solver. 163 ! it is the number of non-zero "w" levels in the water, and the minimum 164 ! value (on land) is 2. We define zprt as the number of "T" points in the ocean 165 ! at any location, and zero on land. 166 ! 167 zprt = tmask(:,:,1)*(mbathy-1) 168 CALL iom_rstput( 0, 0, inum4, 'mbathy', zprt, ktype = jp_i2 ) 158 ! note that mbkt is set to 1 over land ==> use surface tmask 159 zprt(:,:) = tmask(:,:,1) * REAL( mbkt(:,:) , wp ) 160 CALL iom_rstput( 0, 0, inum4, 'mbathy', zprt, ktype = jp_i2 ) ! ! nb of ocean T-points 169 161 170 162 IF( ln_sco ) THEN ! s-coordinate … … 173 165 CALL iom_rstput( 0, 0, inum4, 'hbatv', hbatv ) 174 166 CALL iom_rstput( 0, 0, inum4, 'hbatf', hbatf ) 175 167 ! 176 168 CALL iom_rstput( 0, 0, inum4, 'gsigt', gsigt ) ! ! scaling coef. 177 169 CALL iom_rstput( 0, 0, inum4, 'gsigw', gsigw ) … … 179 171 CALL iom_rstput( 0, 0, inum4, 'esigt', esigt ) 180 172 CALL iom_rstput( 0, 0, inum4, 'esigw', esigw ) 181 173 ! 182 174 CALL iom_rstput( 0, 0, inum4, 'e3t', e3t ) ! ! scale factors 183 175 CALL iom_rstput( 0, 0, inum4, 'e3u', e3u ) 184 176 CALL iom_rstput( 0, 0, inum4, 'e3v', e3v ) 185 177 CALL iom_rstput( 0, 0, inum4, 'e3w', e3w ) 186 178 ! 187 179 CALL iom_rstput( 0, 0, inum4, 'gdept_0' , gdept_0 ) ! ! stretched system 188 180 CALL iom_rstput( 0, 0, inum4, 'gdepw_0' , gdepw_0 ) … … 190 182 191 183 IF( ln_zps ) THEN ! z-coordinate - partial steps 192 184 ! 193 185 IF( nmsh <= 6 ) THEN ! ! 3D vertical scale factors 194 186 CALL iom_rstput( 0, 0, inum4, 'e3t', e3t ) … … 196 188 CALL iom_rstput( 0, 0, inum4, 'e3v', e3v ) 197 189 CALL iom_rstput( 0, 0, inum4, 'e3w', e3w ) 198 ELSE ! ! 2D bottomscale factors199 DO jj = 1,jpj ; DO ji = 1,jpi200 ik = NINT( zprt(ji,jj) ) ! take care that mbathy is not what you think it is here !201 IF ( ik /= 0 ) THEN ; e3tp(ji,jj) = e3t(ji,jj,ik) ; e3wp(ji,jj) = e3w(ji,jj,ik)202 ELSE ; e3tp(ji,jj) = 0. ; e3wp(ji,jj) = 0.203 END IF204 END DO ; END DO190 ELSE ! ! 2D masked bottom ocean scale factors 191 DO jj = 1,jpj 192 DO ji = 1,jpi 193 e3tp(ji,jj) = e3t(ji,jj,mbkt(ji,jj)) * tmask(ji,jj,1) 194 e3wp(ji,jj) = e3w(ji,jj,mbkt(ji,jj)) * tmask(ji,jj,1) 195 END DO 196 END DO 205 197 CALL iom_rstput( 0, 0, inum4, 'e3t_ps', e3tp ) 206 198 CALL iom_rstput( 0, 0, inum4, 'e3w_ps', e3wp ) 207 199 END IF 208 200 ! 209 201 IF( nmsh <= 3 ) THEN ! ! 3D depth 210 202 CALL iom_rstput( 0, 0, inum4, 'gdept', gdept, ktype = jp_r4 ) 211 DO jk = 1,jpk ; DO jj = 1, jpjm1 ; DO ji = 1, fs_jpim1 ! vector opt. 212 zdepu(ji,jj,jk) = MIN( gdept(ji,jj,jk), gdept(ji+1,jj ,jk) ) 213 zdepv(ji,jj,jk) = MIN( gdept(ji,jj,jk), gdept(ji ,jj+1,jk) ) 214 END DO ; END DO ; END DO 203 DO jk = 1,jpk 204 DO jj = 1, jpjm1 205 DO ji = 1, fs_jpim1 ! vector opt. 206 zdepu(ji,jj,jk) = MIN( gdept(ji,jj,jk) , gdept(ji+1,jj ,jk) ) 207 zdepv(ji,jj,jk) = MIN( gdept(ji,jj,jk) , gdept(ji ,jj+1,jk) ) 208 END DO 209 END DO 210 END DO 215 211 CALL lbc_lnk( zdepu, 'U', 1. ) ; CALL lbc_lnk( zdepv, 'V', 1. ) 216 212 CALL iom_rstput( 0, 0, inum4, 'gdepu', zdepu, ktype = jp_r4 ) … … 218 214 CALL iom_rstput( 0, 0, inum4, 'gdepw', gdepw, ktype = jp_r4 ) 219 215 ELSE ! ! 2D bottom depth 220 DO jj = 1,jpj ; DO ji = 1,jpi221 ik = NINT( zprt(ji,jj) ) ! take care that mbathy is not what you think it is here !222 IF ( ik /= 0 ) THEN ; hdept(ji,jj) = gdept(ji,jj,ik) ; hdepw(ji,jj) = gdepw(ji,jj,ik+1)223 ELSE ; hdept(ji,jj) = 0. ; hdepw(ji,jj) = 0.224 END IF225 END DO ; END DO216 DO jj = 1,jpj 217 DO ji = 1,jpi 218 hdept(ji,jj) = gdept(ji,jj,mbkt(ji,jj) ) * tmask(ji,jj,1) 219 hdepw(ji,jj) = gdepw(ji,jj,mbkt(ji,jj)+1) * tmask(ji,jj,1) 220 END DO 221 END DO 226 222 CALL iom_rstput( 0, 0, inum4, 'hdept' , hdept, ktype = jp_r4 ) 227 223 CALL iom_rstput( 0, 0, inum4, 'hdepw' , hdepw, ktype = jp_r4 ) 228 224 ENDIF 229 225 ! 230 226 CALL iom_rstput( 0, 0, inum4, 'gdept_0', gdept_0 ) ! ! reference z-coord. 231 227 CALL iom_rstput( 0, 0, inum4, 'gdepw_0', gdepw_0 ) … … 233 229 CALL iom_rstput( 0, 0, inum4, 'e3w_0' , e3w_0 ) 234 230 ENDIF 235 236 231 237 232 IF( ln_zco ) THEN … … 256 251 CALL iom_close( inum4 ) 257 252 END SELECT 258 253 ! 259 254 END SUBROUTINE dom_wri 260 255 … … 268 263 !! ** Method : 1) aplly lbc_lnk on an array with different values for each element 269 264 !! 2) check which elements have been changed 270 !!271 265 !!---------------------------------------------------------------------- 272 266 CHARACTER(len=1) , INTENT(in ) :: cdgrd ! 273 267 REAL(wp), DIMENSION(jpi,jpj) :: puniq ! 274 268 ! 275 269 REAL(wp), DIMENSION(jpi,jpj ) :: ztstref ! array with different values for each element 276 270 REAL(wp) :: zshift ! shift value link to the process number … … 278 272 INTEGER :: ji ! dummy loop indices 279 273 !!---------------------------------------------------------------------- 280 274 ! 281 275 ! build an array with different values for each element 282 276 ! in mpp: make sure that these values are different even between process 283 277 ! -> apply a shift value according to the process number 284 278 zshift = jpi * jpj * ( narea - 1 ) 285 ztstref(:,:) = RESHAPE( (/ (zshift + REAL(ji, 286 279 ztstref(:,:) = RESHAPE( (/ (zshift + REAL(ji,wp), ji = 1, jpi*jpj) /), (/ jpi, jpj /) ) 280 ! 287 281 puniq(:,:) = ztstref(:,:) ! default definition 288 282 CALL lbc_lnk( puniq, cdgrd, 1. ) ! apply boundary conditions 289 283 lldbl(:,:,1) = puniq(:,:) == ztstref(:,:) ! check which values have been changed 290 284 ! 291 285 puniq(:,:) = 1. ! default definition 292 286 ! fill only the inner part of the cpu with llbl converted into real 293 puniq(nldi:nlei,nldj:nlej) = REAL( COUNT( lldbl(nldi:nlei,nldj:nlej,:), dim = 3 ), wp)294 287 puniq(nldi:nlei,nldj:nlej) = REAL( COUNT( lldbl(nldi:nlei,nldj:nlej,:), dim = 3 ) , wp ) 288 ! 295 289 END FUNCTION dom_uniq 296 297 290 298 291 !!====================================================================== -
branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/DTA/dtasal.F90
r2436 r2450 190 190 DO jj = 1, jpj ! interpolation of salinity at the last ocean level (i.e. the partial step) 191 191 DO ji = 1, jpi 192 ik = mb athy(ji,jj) - 1193 IF( ik > 2) THEN192 ik = mbkt(ji,jj) 193 IF( ik > 1 ) THEN 194 194 zl = ( gdept_0(ik) - fsdept_0(ji,jj,ik) ) / ( gdept_0(ik) - gdept_0(ik-1) ) 195 195 s_dta(ji,jj,ik) = (1.-zl) * s_dta(ji,jj,ik) + zl * s_dta(ji,jj,ik-1) -
branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/DTA/dtatem.F90
r2436 r2450 203 203 DO jj = 1, jpj ! interpolation of temperature at the last level 204 204 DO ji = 1, jpi 205 ik = mb athy(ji,jj) - 1206 IF( ik > 2) THEN205 ik = mbkt(ji,jj) 206 IF( ik > 1 ) THEN 207 207 zl = ( gdept_0(ik) - fsdept_0(ji,jj,ik) ) / ( gdept_0(ik) - gdept_0(ik-1) ) 208 208 t_dta(ji,jj,ik) = (1.-zl) * t_dta(ji,jj,ik) + zl * t_dta(ji,jj,ik-1) -
branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90
r2287 r2450 11 11 !! 8.5 ! 2002-08 (A. Bozec) hpg_zps: Original code 12 12 !! NEMO 1.0 ! 2005-10 (A. Beckmann, B.W. An) various s-coordinate options 13 !! 13 !! ! Original code for hpg_ctl, hpg_hel hpg_wdj, hpg_djc, hpg_rot 14 14 !! - ! 2005-11 (G. Madec) style & small optimisation 15 15 !! 3.3 ! 2010-10 (C. Ethe, G. Madec) reorganisation of initialisation phase … … 51 51 LOGICAL , PUBLIC :: ln_hpg_djc = .FALSE. !: s-coordinate (Density Jacobian with Cubic polynomial) 52 52 LOGICAL , PUBLIC :: ln_hpg_rot = .FALSE. !: s-coordinate (ROTated axes scheme) 53 REAL(wp), PUBLIC :: rn_gamma = 0. e0!: weighting coefficient53 REAL(wp), PUBLIC :: rn_gamma = 0._wp !: weighting coefficient 54 54 LOGICAL , PUBLIC :: ln_dynhpg_imp = .FALSE. !: semi-implicite hpg flag 55 55 INTEGER , PUBLIC :: nn_dynhpg_rst = 0 !: add dynhpg implicit variables in restart ot not … … 63 63 !! NEMO/OPA 3.3 , NEMO Consortium (2010) 64 64 !! $Id$ 65 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)65 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 66 66 !!---------------------------------------------------------------------- 67 68 67 CONTAINS 69 68 … … 82 81 REAL(wp), DIMENSION(jpi,jpj,jpk) :: ztrdu, ztrdv ! 3D temporary workspace 83 82 !!---------------------------------------------------------------------- 84 83 ! 85 84 IF( l_trddyn ) THEN ! Temporary saving of ua and va trends (l_trddyn) 86 85 ztrdu(:,:,:) = ua(:,:,:) 87 86 ztrdv(:,:,:) = va(:,:,:) 88 87 ENDIF 89 88 ! 90 89 SELECT CASE ( nhpg ) ! Hydrastatic pressure gradient computation 91 90 CASE ( 0 ) ; CALL hpg_zco ( kt ) ! z-coordinate … … 97 96 CASE ( 6 ) ; CALL hpg_rot ( kt ) ! s-coordinate (ROTated axes scheme) 98 97 END SELECT 99 98 ! 100 99 IF( l_trddyn ) THEN ! save the hydrostatic pressure gradient trends for momentum trend diagnostics 101 100 ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) … … 207 206 208 207 ! Local constant initialization 209 zcoef0 = - grav * 0.5 208 zcoef0 = - grav * 0.5_wp 210 209 211 210 ! Surface value … … 270 269 271 270 ! Local constant initialization 272 zcoef0 = - grav * 0.5 273 274 ! Surface value 271 zcoef0 = - grav * 0.5_wp 272 273 ! Surface value (also valid in partial step case) 275 274 DO jj = 2, jpjm1 276 275 DO ji = fs_2, fs_jpim1 ! vector opt. … … 305 304 END DO 306 305 307 ! partial steps correction at the last level ( new gradient with intgrd.F)306 ! partial steps correction at the last level (use gru & grv computed in zpshde.F90) 308 307 # if defined key_vectopt_loop 309 308 jj = 1 … … 313 312 DO ji = 2, jpim1 314 313 # endif 315 iku = MIN ( mbathy(ji,jj), mbathy(ji+1,jj) ) - 1316 ikv = MIN ( mbathy(ji,jj), mbathy(ji,jj+1) ) - 1314 iku = mbku(ji,jj) 315 ikv = mbkv(ji,jj) 317 316 zcoef2 = zcoef0 * MIN( fse3w(ji,jj,iku), fse3w(ji+1,jj ,iku) ) 318 317 zcoef3 = zcoef0 * MIN( fse3w(ji,jj,ikv), fse3w(ji ,jj+1,ikv) ) 319 ! on i-direction 320 IF ( iku > 2 ) THEN 321 ! subtract old value 322 ua(ji,jj,iku) = ua(ji,jj,iku) - zhpi(ji,jj,iku) 323 ! compute the new one 324 zhpi (ji,jj,iku) = zhpi(ji,jj,iku-1) & 325 + zcoef2 * ( rhd(ji+1,jj,iku-1) - rhd(ji,jj,iku-1) + gru(ji,jj) ) / e1u(ji,jj) 326 ! add the new one to the general momentum trend 327 ua(ji,jj,iku) = ua(ji,jj,iku) + zhpi(ji,jj,iku) 318 IF( iku > 1 ) THEN ! on i-direction (level 2 or more) 319 ua (ji,jj,iku) = ua(ji,jj,iku) - zhpi(ji,jj,iku) ! subtract old value 320 zhpi(ji,jj,iku) = zhpi(ji,jj,iku-1) & ! compute the new one 321 & + zcoef2 * ( rhd(ji+1,jj,iku-1) - rhd(ji,jj,iku-1) + gru(ji,jj) ) / e1u(ji,jj) 322 ua (ji,jj,iku) = ua(ji,jj,iku) + zhpi(ji,jj,iku) ! add the new one to the general momentum trend 328 323 ENDIF 329 ! on j-direction 330 IF ( ikv > 2 ) THEN 331 ! subtract old value 332 va(ji,jj,ikv) = va(ji,jj,ikv) - zhpj(ji,jj,ikv) 333 ! compute the new one 334 zhpj (ji,jj,ikv) = zhpj(ji,jj,ikv-1) & 335 + zcoef3 * ( rhd(ji,jj+1,ikv-1) - rhd(ji,jj,ikv-1) + grv(ji,jj) ) / e2v(ji,jj) 336 ! add the new one to the general momentum trend 337 va(ji,jj,ikv) = va(ji,jj,ikv) + zhpj(ji,jj,ikv) 324 IF( ikv > 1 ) THEN ! on j-direction (level 2 or more) 325 va (ji,jj,ikv) = va(ji,jj,ikv) - zhpj(ji,jj,ikv) ! subtract old value 326 zhpj(ji,jj,ikv) = zhpj(ji,jj,ikv-1) & ! compute the new one 327 & + zcoef3 * ( rhd(ji,jj+1,ikv-1) - rhd(ji,jj,ikv-1) + grv(ji,jj) ) / e2v(ji,jj) 328 va (ji,jj,ikv) = va(ji,jj,ikv) + zhpj(ji,jj,ikv) ! add the new one to the general momentum trend 338 329 ENDIF 339 330 # if ! defined key_vectopt_loop … … 379 370 380 371 ! Local constant initialization 381 zcoef0 = - grav * 0.5 372 zcoef0 = - grav * 0.5_wp 382 373 ! To use density and not density anomaly 383 IF ( lk_vvl ) THEN ; znad = 1. 384 ELSE ; znad = 0. e0! Fixed volume374 IF ( lk_vvl ) THEN ; znad = 1._wp ! Variable volume 375 ELSE ; znad = 0._wp ! Fixed volume 385 376 ENDIF 386 377 … … 394 385 & - fse3w(ji ,jj ,1) * ( znad + rhd(ji ,jj ,1) ) ) 395 386 ! s-coordinate pressure gradient correction 396 zuap = -zcoef0 * ( rhd (ji+1,jj,1) + rhd (ji,jj,1) + 2 *znad ) &387 zuap = -zcoef0 * ( rhd (ji+1,jj,1) + rhd (ji,jj,1) + 2._wp * znad ) & 397 388 & * ( fsde3w(ji+1,jj,1) - fsde3w(ji,jj,1) ) / e1u(ji,jj) 398 zvap = -zcoef0 * ( rhd (ji,jj+1,1) + rhd (ji,jj,1) + 2 *znad ) &389 zvap = -zcoef0 * ( rhd (ji,jj+1,1) + rhd (ji,jj,1) + 2._wp * znad ) & 399 390 & * ( fsde3w(ji,jj+1,1) - fsde3w(ji,jj,1) ) / e2v(ji,jj) 400 391 ! add to the general momentum trend … … 416 407 & - fse3w(ji,jj ,jk) * ( rhd(ji,jj, jk) + rhd(ji,jj ,jk-1) + 2*znad ) ) 417 408 ! s-coordinate pressure gradient correction 418 zuap = -zcoef0 * ( rhd (ji+1,jj ,jk) + rhd (ji,jj,jk) + 2 *znad ) &409 zuap = -zcoef0 * ( rhd (ji+1,jj ,jk) + rhd (ji,jj,jk) + 2._wp * znad ) & 419 410 & * ( fsde3w(ji+1,jj ,jk) - fsde3w(ji,jj,jk) ) / e1u(ji,jj) 420 zvap = -zcoef0 * ( rhd (ji ,jj+1,jk) + rhd (ji,jj,jk) + 2 *znad ) &411 zvap = -zcoef0 * ( rhd (ji ,jj+1,jk) + rhd (ji,jj,jk) + 2._wp * znad ) & 421 412 & * ( fsde3w(ji ,jj+1,jk) - fsde3w(ji,jj,jk) ) / e2v(ji,jj) 422 413 ! add to the general momentum trend … … 465 456 466 457 ! Local constant initialization 467 zcoef0 = - grav * 0.5 458 zcoef0 = - grav * 0.5_wp 468 459 469 460 ! Surface value … … 497 488 & -fse3t(ji ,jj,jk-1) * rhd(ji ,jj,jk-1) ) 498 489 zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) & 499 & + zcoef0 / e2v(ji,jj) * ( fse3t(ji,jj+1,jk ) * rhd(ji,jj+1,jk) &490 & + zcoef0 / e2v(ji,jj) * ( fse3t(ji,jj+1,jk ) * rhd(ji,jj+1,jk) & 500 491 & -fse3t(ji,jj ,jk ) * rhd(ji,jj, jk) ) & 501 & + zcoef0 / e2v(ji,jj) * ( fse3t(ji,jj+1,jk-1) * rhd(ji,jj+1,jk-1) &492 & + zcoef0 / e2v(ji,jj) * ( fse3t(ji,jj+1,jk-1) * rhd(ji,jj+1,jk-1) & 502 493 & -fse3t(ji,jj ,jk-1) * rhd(ji,jj, jk-1) ) 503 494 ! s-coordinate pressure gradient correction … … 543 534 544 535 ! Local constant initialization 545 zcoef0 = - grav * 0.5 546 zalph = 0.5 - rn_gamma ! weighting coefficients (alpha=0.5-rn_gamma547 zbeta = 0.5 + rn_gamma ! (beta =1-alpha=0.5+rn_gamma536 zcoef0 = - grav * 0.5_wp 537 zalph = 0.5_wp - rn_gamma ! weighting coefficients (alpha=0.5-rn_gamma 538 zbeta = 0.5_wp + rn_gamma ! (beta =1-alpha=0.5+rn_gamma 548 539 549 540 ! Surface value (no ponderation) … … 628 619 629 620 ! Local constant initialization 630 zcoef0 = - grav * 0.5 631 z1_10 = 1. 0 / 10.0632 z1_12 = 1. 0 / 12.0621 zcoef0 = - grav * 0.5_wp 622 z1_10 = 1._wp / 10._wp 623 z1_12 = 1._wp / 12._wp 633 624 634 625 !---------------------------------------------------------------------------------------- … … 662 653 DO jj = 2, jpjm1 663 654 DO ji = fs_2, fs_jpim1 ! vector opt. 664 cffw = 2. 0* drhoz(ji ,jj ,jk) * drhoz(ji,jj,jk-1)665 666 cffu = 2. 0* drhox(ji+1,jj ,jk) * drhox(ji,jj,jk )667 cffx = 2. 0* dzx (ji+1,jj ,jk) * dzx (ji,jj,jk )655 cffw = 2._wp * drhoz(ji ,jj ,jk) * drhoz(ji,jj,jk-1) 656 657 cffu = 2._wp * drhox(ji+1,jj ,jk) * drhox(ji,jj,jk ) 658 cffx = 2._wp * dzx (ji+1,jj ,jk) * dzx (ji,jj,jk ) 668 659 669 cffv = 2. 0* drhoy(ji ,jj+1,jk) * drhoy(ji,jj,jk )670 cffy = 2. 0* dzy (ji ,jj+1,jk) * dzy (ji,jj,jk )660 cffv = 2._wp * drhoy(ji ,jj+1,jk) * drhoy(ji,jj,jk ) 661 cffy = 2._wp * dzy (ji ,jj+1,jk) * dzy (ji,jj,jk ) 671 662 672 663 IF( cffw > zep) THEN 673 drhow(ji,jj,jk) = 2. 0* drhoz(ji,jj,jk) * drhoz(ji,jj,jk-1) &674 & / ( drhoz(ji,jj,jk) + drhoz(ji,jj,jk-1) )664 drhow(ji,jj,jk) = 2._wp * drhoz(ji,jj,jk) * drhoz(ji,jj,jk-1) & 665 & / ( drhoz(ji,jj,jk) + drhoz(ji,jj,jk-1) ) 675 666 ELSE 676 drhow(ji,jj,jk) = 0. e0667 drhow(ji,jj,jk) = 0._wp 677 668 ENDIF 678 669 679 dzw(ji,jj,jk) = 2. 0* dzz(ji,jj,jk) * dzz(ji,jj,jk-1) &680 & / ( dzz(ji,jj,jk) + dzz(ji,jj,jk-1) )670 dzw(ji,jj,jk) = 2._wp * dzz(ji,jj,jk) * dzz(ji,jj,jk-1) & 671 & / ( dzz(ji,jj,jk) + dzz(ji,jj,jk-1) ) 681 672 682 673 IF( cffu > zep ) THEN 683 drhou(ji,jj,jk) = 2. 0* drhox(ji+1,jj,jk) * drhox(ji,jj,jk) &684 & / ( drhox(ji+1,jj,jk) + drhox(ji,jj,jk) )674 drhou(ji,jj,jk) = 2._wp * drhox(ji+1,jj,jk) * drhox(ji,jj,jk) & 675 & / ( drhox(ji+1,jj,jk) + drhox(ji,jj,jk) ) 685 676 ELSE 686 drhou(ji,jj,jk ) = 0. e0677 drhou(ji,jj,jk ) = 0._wp 687 678 ENDIF 688 679 689 680 IF( cffx > zep ) THEN 690 dzu(ji,jj,jk) = 2. 0*dzx(ji+1,jj,jk)*dzx(ji,jj,jk) &691 & /(dzx(ji+1,jj,jk)+dzx(ji,jj,jk))681 dzu(ji,jj,jk) = 2._wp * dzx(ji+1,jj,jk) * dzx(ji,jj,jk) & 682 & / ( dzx(ji+1,jj,jk) + dzx(ji,jj,jk) ) 692 683 ELSE 693 dzu(ji,jj,jk) = 0. e0684 dzu(ji,jj,jk) = 0._wp 694 685 ENDIF 695 686 696 687 IF( cffv > zep ) THEN 697 drhov(ji,jj,jk) = 2. 0* drhoy(ji,jj+1,jk) * drhoy(ji,jj,jk) &698 & / ( drhoy(ji,jj+1,jk) + drhoy(ji,jj,jk) )688 drhov(ji,jj,jk) = 2._wp * drhoy(ji,jj+1,jk) * drhoy(ji,jj,jk) & 689 & / ( drhoy(ji,jj+1,jk) + drhoy(ji,jj,jk) ) 699 690 ELSE 700 drhov(ji,jj,jk) = 0. e0691 drhov(ji,jj,jk) = 0._wp 701 692 ENDIF 702 693 703 694 IF( cffy > zep ) THEN 704 dzv(ji,jj,jk) = 2. 0* dzy(ji,jj+1,jk) * dzy(ji,jj,jk) &705 & / ( dzy(ji,jj+1,jk) + dzy(ji,jj,jk) )695 dzv(ji,jj,jk) = 2._wp * dzy(ji,jj+1,jk) * dzy(ji,jj,jk) & 696 & / ( dzy(ji,jj+1,jk) + dzy(ji,jj,jk) ) 706 697 ELSE 707 dzv(ji,jj,jk) = 0. e0698 dzv(ji,jj,jk) = 0._wp 708 699 ENDIF 709 700 … … 715 706 ! apply boundary conditions at top and bottom using 5.36-5.37 716 707 !---------------------------------------------------------------------------------- 717 drhow(:,:, 1 ) = 1.5 * ( drhoz(:,:, 2 ) - drhoz(:,:, 1 ) ) - 0.5* drhow(:,:, 2 )718 drhou(:,:, 1 ) = 1.5 * ( drhox(:,:, 2 ) - drhox(:,:, 1 ) ) - 0.5* drhou(:,:, 2 )719 drhov(:,:, 1 ) = 1.5 * ( drhoy(:,:, 2 ) - drhoy(:,:, 1 ) ) - 0.5* drhov(:,:, 2 )720 721 drhow(:,:,jpk) = 1.5 * ( drhoz(:,:,jpk) - drhoz(:,:,jpkm1) ) - 0.5* drhow(:,:,jpkm1)722 drhou(:,:,jpk) = 1.5 * ( drhox(:,:,jpk) - drhox(:,:,jpkm1) ) - 0.5* drhou(:,:,jpkm1)723 drhov(:,:,jpk) = 1.5 * ( drhoy(:,:,jpk) - drhoy(:,:,jpkm1) ) - 0.5* drhov(:,:,jpkm1)708 drhow(:,:, 1 ) = 1.5_wp * ( drhoz(:,:, 2 ) - drhoz(:,:, 1 ) ) - 0.5_wp * drhow(:,:, 2 ) 709 drhou(:,:, 1 ) = 1.5_wp * ( drhox(:,:, 2 ) - drhox(:,:, 1 ) ) - 0.5_wp * drhou(:,:, 2 ) 710 drhov(:,:, 1 ) = 1.5_wp * ( drhoy(:,:, 2 ) - drhoy(:,:, 1 ) ) - 0.5_wp * drhov(:,:, 2 ) 711 712 drhow(:,:,jpk) = 1.5_wp * ( drhoz(:,:,jpk) - drhoz(:,:,jpkm1) ) - 0.5_wp * drhow(:,:,jpkm1) 713 drhou(:,:,jpk) = 1.5_wp * ( drhox(:,:,jpk) - drhox(:,:,jpkm1) ) - 0.5_wp * drhou(:,:,jpkm1) 714 drhov(:,:,jpk) = 1.5_wp * ( drhoy(:,:,jpk) - drhoy(:,:,jpkm1) ) - 0.5_wp * drhov(:,:,jpkm1) 724 715 725 716 … … 733 724 DO jj = 2, jpjm1 734 725 DO ji = fs_2, fs_jpim1 ! vector opt. 735 rho_k(ji,jj,1) = -grav * ( fse3w(ji,jj,1) - fsde3w(ji,jj,1) ) &736 & * ( rhd(ji,jj,1) &737 & + 0.5 * ( rhd(ji,jj,2) - rhd(ji,jj,1) ) &738 & * ( fse3w (ji,jj,1) - fsde3w(ji,jj,1) ) &739 & / ( fsde3w(ji,jj,2) - fsde3w(ji,jj,1) ) )726 rho_k(ji,jj,1) = -grav * ( fse3w(ji,jj,1) - fsde3w(ji,jj,1) ) & 727 & * ( rhd(ji,jj,1) & 728 & + 0.5_wp * ( rhd(ji,jj,2) - rhd(ji,jj,1) ) & 729 & * ( fse3w (ji,jj,1) - fsde3w(ji,jj,1) ) & 730 & / ( fsde3w(ji,jj,2) - fsde3w(ji,jj,1) ) ) 740 731 END DO 741 732 END DO … … 849 840 ! Local constant initialization 850 841 ! ------------------------------- 851 zcoef0 = - grav * 0.5 852 zforg = 0.95 e0853 zfrot = 1. e0- zforg842 zcoef0 = - grav * 0.5_wp 843 zforg = 0.95_wp 844 zfrot = 1._wp - zforg 854 845 855 846 ! inverse of the distance between 2 diagonal T-points (defined at F-point) (here zcoef0/distance) -
branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_imp.F90
r2287 r2450 69 69 REAL(wp) :: zzwi, zzws, zrhs ! temporary scalars 70 70 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zwi ! 3D workspace 71 #if defined key_zdfgls72 INTEGER :: ikbu, ikbv, ikbum1, ikbvm173 REAL(wp) :: zcbcu, zcbcv74 #endif75 71 !!---------------------------------------------------------------------- 76 72 … … 168 164 DO jj = 2, jpjm1 169 165 DO ji = fs_2, fs_jpim1 ! vector opt. 170 ikbu = MIN( mbathy(ji+1,jj ), mbathy(ji,jj) ) 171 ikbum1 = MAX( ikbu-1, 1 ) 172 wbotu(ji,jj) = bfrua(ji,jj) * ua(ji,jj,ikbum1) * umask(ji,jj,ikbum1) 166 wbotu(ji,jj) = bfrua(ji,jj) * ua(ji,jj,mbku(ji,jj)) * umask(ji,jj,mbku(ji,jj)) 173 167 END DO 174 168 END DO … … 269 263 DO jj = 2, jpjm1 270 264 DO ji = fs_2, fs_jpim1 ! vector opt. 271 ikbv = MIN( mbathy(ji,jj+1), mbathy(ji,jj) ) 272 ikbvm1 = MAX( ikbv-1, 1 ) 273 wbotv(ji,jj) = bfrva(ji,jj) * va(ji,jj,ikbvm1) * vmask(ji,jj,ikbvm1) 265 wbotv(ji,jj) = bfrva(ji,jj) * va(ji,jj,mbkv(ji,jj)) * vmask(ji,jj,mbku(ji,jj)) 274 266 END DO 275 267 END DO -
branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/LDF/ldfeiv.F90
r2371 r2450 4 4 !! Ocean physics: variable eddy induced velocity coefficients 5 5 !!====================================================================== 6 !! History : OPA ! 1999-03 (G. Madec, A. Jouzeau) Original code 7 !! NEMO 1.0 ! 2002-06 (G. Madec) Free form, F90 8 !!---------------------------------------------------------------------- 6 9 #if defined key_traldf_eiv && defined key_traldf_c2d 7 10 !!---------------------------------------------------------------------- … … 11 14 !! ldf_eiv : compute the eddy induced velocity coefficients 12 15 !!---------------------------------------------------------------------- 13 !! * Modules used14 16 USE oce ! ocean dynamics and tracers 15 17 USE dom_oce ! ocean space and time domain … … 22 24 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 23 25 USE prtctl ! Print control 24 USE iom 26 USE iom ! I/O library 25 27 26 28 IMPLICIT NONE 27 29 PRIVATE 28 30 29 !! * Routine accessibility 30 PUBLIC ldf_eiv ! routine called by step.F90 31 !!---------------------------------------------------------------------- 32 !! NEMO/OPA 3.3 , NEMO Consortium (2010) 33 !! $Id$ 34 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 35 !!---------------------------------------------------------------------- 31 PUBLIC ldf_eiv ! routine called by step.F90 32 36 33 !! * Substitutions 37 34 # include "domzgr_substitute.h90" 38 35 # include "vectopt_loop_substitute.h90" 39 36 !!---------------------------------------------------------------------- 40 37 !! NEMO/OPA 3.3 , NEMO Consortium (2010) 38 !! $Id$ 39 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 40 !!---------------------------------------------------------------------- 41 41 CONTAINS 42 42 … … 46 46 !! 47 47 !! ** Purpose : Compute the eddy induced velocity coefficient from the 48 !! growth rate of baroclinic instability.48 !! growth rate of baroclinic instability. 49 49 !! 50 50 !! ** Method : 51 51 !! 52 !! ** Action : - uslp(), : i- and j-slopes of neutral surfaces 53 !! - vslp() at u- and v-points, resp. 54 !! - wslpi(), : i- and j-slopes of neutral surfaces 55 !! - wslpj() at w-points. 56 !! 57 !! History : 58 !! 8.1 ! 99-03 (G. Madec, A. Jouzeau) Original code 59 !! 8.5 ! 02-06 (G. Madec) Free form, F90 52 !! ** Action : - uslp , vslp : i- and j-slopes of neutral surfaces at u- & v-points 53 !! - wslpi, wslpj : i- and j-slopes of neutral surfaces at w-points. 60 54 !!---------------------------------------------------------------------- 61 !! * Arguments 62 INTEGER, INTENT( in ) :: kt ! ocean time-step inedx 63 64 !! * Local declarations 65 INTEGER :: ji, jj, jk ! dummy loop indices 66 REAL(wp) :: & 67 zfw, ze3w, zn2, zf20, & ! temporary scalars 68 zaht, zaht_min 69 REAL(wp), DIMENSION(jpi,jpj) :: & 70 zn, zah, zhw, zross ! workspace 55 INTEGER, INTENT(in) :: kt ! ocean time-step inedx 56 !! 57 INTEGER :: ji, jj, jk ! dummy loop indices 58 REAL(wp) :: zfw, ze3w, zn2, zf20, zaht, zaht_min ! temporary scalars 59 REAL(wp), DIMENSION(jpi,jpj) :: zn, zah, zhw, zross ! 2D workspace 71 60 !!---------------------------------------------------------------------- 72 61 … … 79 68 ! 0. Local initialization 80 69 ! ----------------------- 81 zn (:,:) = 0. e082 zhw (:,:) = 5. e083 zah (:,:) = 0. e084 zross(:,:) = 0. e070 zn (:,:) = 0._wp 71 zhw (:,:) = 5._wp 72 zah (:,:) = 0._wp 73 zross(:,:) = 0._wp 85 74 86 75 87 76 ! 1. Compute lateral diffusive coefficient 88 77 ! ---------------------------------------- 89 IF (ln_traldf_grif) THEN78 IF( ln_traldf_grif ) THEN 90 79 DO jk = 1, jpk 91 80 # if defined key_vectopt_loop 92 81 !CDIR NOVERRCHK 93 82 DO ji = 1, jpij ! vector opt. 94 83 ! Take the max of N^2 and zero then take the vertical sum 95 84 ! of the square root of the resulting N^2 ( required to compute 96 85 ! internal Rossby radius Ro = .5 * sum_jpk(N) / f 97 zn2 = MAX( rn2b(ji,1,jk), 0. e0)86 zn2 = MAX( rn2b(ji,1,jk), 0._wp ) 98 87 zn(ji,1) = zn(ji,1) + SQRT( zn2 ) * fse3w(ji,1,jk) 99 88 ! Compute elements required for the inverse time scale of baroclinic … … 106 95 # else 107 96 DO jj = 2, jpjm1 108 97 !CDIR NOVERRCHK 109 98 DO ji = 2, jpim1 110 99 ! Take the max of N^2 and zero then take the vertical sum 111 100 ! of the square root of the resulting N^2 ( required to compute 112 101 ! internal Rossby radius Ro = .5 * sum_jpk(N) / f 113 zn2 = MAX( rn2b(ji,jj,jk), 0. e0)102 zn2 = MAX( rn2b(ji,jj,jk), 0._wp ) 114 103 zn(ji,jj) = zn(ji,jj) + SQRT( zn2 ) * fse3w(ji,jj,jk) 115 104 ! Compute elements required for the inverse time scale of baroclinic … … 126 115 DO jk = 1, jpk 127 116 # if defined key_vectopt_loop 128 117 !CDIR NOVERRCHK 129 118 DO ji = 1, jpij ! vector opt. 130 119 ! Take the max of N^2 and zero then take the vertical sum 131 120 ! of the square root of the resulting N^2 ( required to compute 132 121 ! internal Rossby radius Ro = .5 * sum_jpk(N) / f 133 zn2 = MAX( rn2b(ji,1,jk), 0. e0)122 zn2 = MAX( rn2b(ji,1,jk), 0._wp ) 134 123 zn(ji,1) = zn(ji,1) + SQRT( zn2 ) * fse3w(ji,1,jk) 135 124 ! Compute elements required for the inverse time scale of baroclinic … … 137 126 ! T^-1 = sqrt(m_jpk(N^2*(r1^2+r2^2)*e3w)) 138 127 ze3w = fse3w(ji,1,jk) * tmask(ji,1,jk) 139 zah(ji,1) = zah(ji,1) + zn2 & 140 * ( wslpi(ji,1,jk) * wslpi(ji,1,jk) & 141 + wslpj(ji,1,jk) * wslpj(ji,1,jk) ) & 142 * ze3w 128 zah(ji,1) = zah(ji,1) + zn2 * ( wslpi(ji,1,jk) * wslpi(ji,1,jk) & 129 & + wslpj(ji,1,jk) * wslpj(ji,1,jk) ) * ze3w 143 130 zhw(ji,1) = zhw(ji,1) + ze3w 144 131 END DO 145 132 # else 146 133 DO jj = 2, jpjm1 147 134 !CDIR NOVERRCHK 148 135 DO ji = 2, jpim1 149 136 ! Take the max of N^2 and zero then take the vertical sum 150 137 ! of the square root of the resulting N^2 ( required to compute 151 138 ! internal Rossby radius Ro = .5 * sum_jpk(N) / f 152 zn2 = MAX( rn2b(ji,jj,jk), 0. e0)139 zn2 = MAX( rn2b(ji,jj,jk), 0._wp ) 153 140 zn(ji,jj) = zn(ji,jj) + SQRT( zn2 ) * fse3w(ji,jj,jk) 154 141 ! Compute elements required for the inverse time scale of baroclinic … … 156 143 ! T^-1 = sqrt(m_jpk(N^2*(r1^2+r2^2)*e3w)) 157 144 ze3w = fse3w(ji,jj,jk) * tmask(ji,jj,jk) 158 zah(ji,jj) = zah(ji,jj) + zn2 & 159 * ( wslpi(ji,jj,jk) * wslpi(ji,jj,jk) & 160 + wslpj(ji,jj,jk) * wslpj(ji,jj,jk) ) & 161 * ze3w 145 zah(ji,jj) = zah(ji,jj) + zn2 * ( wslpi(ji,jj,jk) * wslpi(ji,jj,jk) & 146 & + wslpj(ji,jj,jk) * wslpj(ji,jj,jk) ) * ze3w 162 147 zhw(ji,jj) = zhw(ji,jj) + ze3w 163 148 END DO … … 178 163 END DO 179 164 180 IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN ! ORCA R 02165 IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN ! ORCA R2 181 166 DO jj = 2, jpjm1 182 167 !CDIR NOVERRCHK 183 168 DO ji = fs_2, fs_jpim1 ! vector opt. 184 ! Take the minimum between aeiw and aeiv0 for depth levels 185 ! lower than 20 (21 in w- point) 186 IF( mbathy(ji,jj) <= 21. ) aeiw(ji,jj) = MIN( aeiw(ji,jj), 1000. ) 169 ! Take the minimum between aeiw and 1000 m2/s over shelves (depth shallower than 650 m) 170 IF( mbkt(ji,jj) <= 20 ) aeiw(ji,jj) = MIN( aeiw(ji,jj), 1000. ) 187 171 END DO 188 172 END DO … … 190 174 191 175 ! Decrease the coefficient in the tropics (20N-20S) 192 zf20 = 2. * omega * sin( rad * 20.)176 zf20 = 2._wp * omega * sin( rad * 20._wp ) 193 177 DO jj = 2, jpjm1 194 178 DO ji = fs_2, fs_jpim1 ! vector opt. … … 205 189 END DO 206 190 ENDIF 207 208 ! lateral boundary condition on aeiw 209 CALL lbc_lnk( aeiw, 'W', 1. ) 191 CALL lbc_lnk( aeiw, 'W', 1. ) ! lateral boundary condition on aeiw 192 210 193 211 194 ! Average the diffusive coefficient at u- v- points 212 195 DO jj = 2, jpjm1 213 196 DO ji = fs_2, fs_jpim1 ! vector opt. 214 aeiu(ji,jj) = .5* ( aeiw(ji,jj) + aeiw(ji+1,jj ) )215 aeiv(ji,jj) = .5* ( aeiw(ji,jj) + aeiw(ji ,jj+1) )197 aeiu(ji,jj) = 0.5_wp * ( aeiw(ji,jj) + aeiw(ji+1,jj ) ) 198 aeiv(ji,jj) = 0.5_wp * ( aeiw(ji,jj) + aeiw(ji ,jj+1) ) 216 199 END DO 217 200 END DO 218 219 ! lateral boundary condition on aeiu, aeiv 220 CALL lbc_lnk( aeiu, 'U', 1. ) 221 CALL lbc_lnk( aeiv, 'V', 1. ) 222 223 IF(ln_ctl) THEN 201 CALL lbc_lnk( aeiu, 'U', 1. ) ; CALL lbc_lnk( aeiv, 'V', 1. ) ! lateral boundary condition 202 203 204 IF(ln_ctl) THEN 224 205 CALL prt_ctl(tab2d_1=aeiu, clinfo1=' eiv - u: ', ovlap=1) 225 206 CALL prt_ctl(tab2d_1=aeiv, clinfo1=' eiv - v: ', ovlap=1) … … 228 209 ! ORCA R05: add a space variation on aht (=aeiv except at the equator and river mouth) 229 210 IF( cp_cfg == "orca" .AND. jp_cfg == 05 ) THEN 230 zf20 = 2. * omega * SIN( rad * 20.)231 zaht_min = 100. 211 zf20 = 2._wp * omega * SIN( rad * 20._wp ) 212 zaht_min = 100._wp ! minimum value for aht 232 213 DO jj = 1, jpj 233 214 DO ji = 1, jpi 234 zaht = ( 1. - MIN( 1., ABS( ff(ji,jj) / zf20 ) ) ) * ( aht0 - zaht_min ) &215 zaht = ( 1._wp - MIN( 1._wp , ABS( ff(ji,jj) / zf20 ) ) ) * ( aht0 - zaht_min ) & 235 216 & + aht0 * rnfmsk(ji,jj) ! enhanced near river mouths 236 217 ahtu(ji,jj) = MAX( MAX( zaht_min, aeiu(ji,jj) ) + zaht, aht0 ) … … 246 227 ENDIF 247 228 248 IF( aeiv0 == 0. e0) THEN249 aeiu(:,:) = 0. e0250 aeiv(:,:) = 0. e0251 aeiw(:,:) = 0. e0229 IF( aeiv0 == 0._wp ) THEN 230 aeiu(:,:) = 0._wp 231 aeiv(:,:) = 0._wp 232 aeiw(:,:) = 0._wp 252 233 ENDIF 253 234 254 235 CALL iom_put( "aht2d" , ahtw ) ! lateral eddy diffusivity 255 236 CALL iom_put( "aht2d_eiv", aeiw ) ! EIV lateral eddy diffusivity 256 237 ! 257 238 END SUBROUTINE ldf_eiv 258 239 -
branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/LDF/ldfslp.F90
r2435 r2450 7 7 !! 8.0 ! 1997-06 (G. Madec) optimization, lbc 8 8 !! 8.1 ! 1999-10 (A. Jouzeau) NEW profile in the mixed layer 9 !! NEMO 0.5! 2002-10 (G. Madec) Free form, F9010 !! 1.0! 2005-10 (A. Beckmann) correction for s-coordinates9 !! NEMO 1.0 ! 2002-10 (G. Madec) Free form, F90 10 !! - ! 2005-10 (A. Beckmann) correction for s-coordinates 11 11 !! 3.3 ! 2010-10 (G. Nurser, C. Harris, G. Madec) add Griffies operator 12 12 !! - ! 2010-11 (F. Dupond, G. Madec) bug correction in slopes just below the ML … … 116 116 zm1_2g = -0.5_wp / grav 117 117 ! 118 zww(:,:,:) = 0. e0119 zwz(:,:,:) = 0. e0118 zww(:,:,:) = 0._wp 119 zwz(:,:,:) = 0._wp 120 120 ! 121 121 DO jk = 1, jpk !== i- & j-gradient of density ==! … … 135 135 DO ji = 1, jpim1 136 136 # endif 137 iku = MIN ( mbathy(ji,jj), mbathy(ji+1,jj) ) - 1 ! last ocean level 138 ikv = MIN ( mbathy(ji,jj), mbathy(ji,jj+1) ) - 1 139 zgru(ji,jj,iku) = gru(ji,jj) 140 zgrv(ji,jj,ikv) = grv(ji,jj) 137 zgru(ji,jj,mbku(ji,jj)) = gru(ji,jj) 138 zgrv(ji,jj,mbkv(ji,jj)) = grv(ji,jj) 141 139 END DO 142 140 END DO … … 329 327 ! ! Gibraltar Strait 330 328 ij0 = 50 ; ij1 = 53 331 ii0 = 69 ; ii1 = 71 ; uslp ( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0. e0329 ii0 = 69 ; ii1 = 71 ; uslp ( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0._wp 332 330 ij0 = 51 ; ij1 = 53 333 ii0 = 68 ; ii1 = 71 ; vslp ( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0. e0334 ii0 = 69 ; ii1 = 71 ; wslpi( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0. e0335 ii0 = 69 ; ii1 = 71 ; wslpj( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0. e0331 ii0 = 68 ; ii1 = 71 ; vslp ( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0._wp 332 ii0 = 69 ; ii1 = 71 ; wslpi( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0._wp 333 ii0 = 69 ; ii1 = 71 ; wslpj( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0._wp 336 334 ! 337 335 ! ! Mediterrannean Sea 338 336 ij0 = 49 ; ij1 = 56 339 ii0 = 71 ; ii1 = 90 ; uslp ( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0. e0337 ii0 = 71 ; ii1 = 90 ; uslp ( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0._wp 340 338 ij0 = 50 ; ij1 = 56 341 ii0 = 70 ; ii1 = 90 ; vslp ( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0. e0342 ii0 = 71 ; ii1 = 90 ; wslpi( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0. e0343 ii0 = 71 ; ii1 = 90 ; wslpj( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0. e0339 ii0 = 70 ; ii1 = 90 ; vslp ( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0._wp 340 ii0 = 71 ; ii1 = 90 ; wslpi( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0._wp 341 ii0 = 71 ; ii1 = 90 ; wslpj( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0._wp 344 342 ENDIF 345 343 … … 403 401 DO jj = 1, jpjm1 404 402 DO ji = 1, fs_jpim1 ! vector opt. 405 zdit(ji,jj,jk) = ( tb(ji+1,jj,jk) - tb(ji,jj,jk) ) * umask(ji,jj,jk) 403 zdit(ji,jj,jk) = ( tb(ji+1,jj,jk) - tb(ji,jj,jk) ) * umask(ji,jj,jk) ! i-gradient of T and S at jj 406 404 zdis(ji,jj,jk) = ( sb(ji+1,jj,jk) - sb(ji,jj,jk) ) * umask(ji,jj,jk) 407 zdjt(ji,jj,jk) = ( tb(ji,jj+1,jk) - tb(ji,jj,jk) ) * vmask(ji,jj,jk) 405 zdjt(ji,jj,jk) = ( tb(ji,jj+1,jk) - tb(ji,jj,jk) ) * vmask(ji,jj,jk) ! j-gradient of T and S at jj 408 406 zdjs(ji,jj,jk) = ( sb(ji,jj+1,jk) - sb(ji,jj,jk) ) * vmask(ji,jj,jk) 409 407 END DO … … 418 416 DO ji = 1, jpim1 419 417 # endif 420 iku = MAX( MIN( mbathy(ji,jj), mbathy(ji+1,jj ) ) - 1, 2 ) ! last ocean level 421 ikv = MAX( MIN( mbathy(ji,jj), mbathy(ji,jj+1 ) ) - 1, 2 ) 422 zdit (ji,jj,iku) = gtsu(ji,jj,jp_tem) ! i-gradient of T and S 423 zdis (ji,jj,iku) = gtsu(ji,jj,jp_sal) 424 zdjt (ji,jj,ikv) = gtsv(ji,jj,jp_tem) ! j-gradient of T and S 425 zdjs (ji,jj,ikv) = gtsv(ji,jj,jp_sal) 418 zdit(ji,jj,mbku(ji,jj)) = gtsu(ji,jj,jp_tem) ! i-gradient of T and S 419 zdis(ji,jj,mbku(ji,jj)) = gtsu(ji,jj,jp_sal) 420 zdjt(ji,jj,mbkv(ji,jj)) = gtsv(ji,jj,jp_tem) ! j-gradient of T and S 421 zdjs(ji,jj,mbkv(ji,jj)) = gtsv(ji,jj,jp_sal) 426 422 END DO 427 423 END DO … … 453 449 DO jj = 1, jpj ! NB: not masked due to the minimum value set 454 450 DO ji = 1, jpi ! vector opt. 455 451 zdzrho_raw = ( zalpha(ji,jj,jk) * zdkt(ji,jj,jk+kp) + zbeta(ji,jj,jk) * zdks(ji,jj,jk+kp) ) & 456 452 & / fse3w(ji,jj,jk+kp) 457 453 zdzrho(ji ,jj ,jk, kp) = - MIN( - repsln, zdzrho_raw ) ! force zdzrho >= repsln … … 462 458 ! 463 459 DO jj = 1, jpj !== Reciprocal depth of the w-point below ML base ==! 464 DO ji = 1, jpi ! i.e. 1 / (hmld+e3t(nmln)) where hmld=depw(nmln)465 jk = MIN( nmln(ji,jj), mb athy(ji,jj) - 1) + 1 ! MIN in case ML depth is the ocean depth460 DO ji = 1, jpi 461 jk = MIN( nmln(ji,jj), mbkt(ji,jj) ) + 1 ! MIN in case ML depth is the ocean depth 466 462 z1_mlbw(ji,jj) = 1._wp / fsdepw(ji,jj,jk) 467 463 END DO … … 471 467 ! 472 468 wslp2 (:,:,:) = 0._wp ! wslp2 will be cumulated 3D field set to zero 473 triadi_g(:,:,1,:,:) = 0._wp ; triadi_g(:,:,jpk,:,:) = 0._wp 469 triadi_g(:,:,1,:,:) = 0._wp ; triadi_g(:,:,jpk,:,:) = 0._wp ! set surface and bottom slope to zero 474 470 triadj_g(:,:,1,:,:) = 0._wp ; triadj_g(:,:,jpk,:,:) = 0._wp 475 471 !!gm _iso set to zero missing … … 486 482 DO ji = 1, fs_jpim1 487 483 ip = jl ; jp = jl 488 jk = MIN( nmln(ji+ip,jj) , mbathy(ji+ip,jj) - 1 ) + 1! ML level+1 (MIN in case ML depth is the ocean depth)484 jk = MIN( nmln(ji+ip,jj) , mbkt(ji+ip,jj) ) + 1 ! ML level+1 (MIN in case ML depth is the ocean depth) 489 485 zti_g_raw = ( zdxrho(ji+ip,jj,jk-kp,1-ip) / zdzrho(ji+ip,jj,jk-kp,kp) & 490 486 & + ( fsdept(ji+1,jj,jk-kp) - fsdept(ji,jj,jk-kp) ) / e1u(ji,jj) ) * umask(ji,jj,jk) 491 jk = MIN( nmln(ji,jj+jp) , mbathy(ji,jj+jp) - 1) + 1487 jk = MIN( nmln(ji,jj+jp) , mbkt(ji,jj+jp) ) + 1 492 488 ztj_g_raw = ( zdyrho(ji,jj+jp,jk-kp,1-jp) / zdzrho(ji,jj+jp,jk-kp,kp) & 493 489 & + ( fsdept(ji,jj+1,jk-kp) - fsdept(ji,jj,jk-kp) ) / e2v(ji,jj) ) * vmask(ji,jj,jk) … … 612 608 zm1_2g = -0.5_wp / grav 613 609 ! 614 uslpml (1,:) = 0. e0 ; uslpml (jpi,:) = 0.e0615 vslpml (1,:) = 0. e0 ; vslpml (jpi,:) = 0.e0616 wslpiml(1,:) = 0. e0 ; wslpiml(jpi,:) = 0.e0617 wslpjml(1,:) = 0. e0 ; wslpjml(jpi,:) = 0.e0610 uslpml (1,:) = 0._wp ; uslpml (jpi,:) = 0._wp 611 vslpml (1,:) = 0._wp ; vslpml (jpi,:) = 0._wp 612 wslpiml(1,:) = 0._wp ; wslpiml(jpi,:) = 0._wp 613 wslpjml(1,:) = 0._wp ; wslpjml(jpi,:) = 0._wp 618 614 ! 619 615 ! !== surface mixed layer mask ! … … 627 623 # endif 628 624 ik = nmln(ji,jj) - 1 629 IF( jk <= ik ) THEN ; omlmask(ji,jj,jk) = 1. e0630 ELSE ; omlmask(ji,jj,jk) = 0. e0625 IF( jk <= ik ) THEN ; omlmask(ji,jj,jk) = 1._wp 626 ELSE ; omlmask(ji,jj,jk) = 0._wp 631 627 ENDIF 632 628 END DO -
branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/SBC/sbcrnf.F90
r2384 r2450 19 19 USE phycst ! physical constants 20 20 USE sbc_oce ! surface boundary condition variables 21 USE fldread ! ???21 USE fldread ! read input field at current time step 22 22 USE in_out_manager ! I/O manager 23 23 USE iom ! I/O module 24 24 USE restart ! restart 25 USE closea 25 USE closea ! closed seas 26 26 27 27 IMPLICIT NONE 28 28 PRIVATE 29 29 30 PUBLIC sbc_rnf! routine call in sbcmod module31 PUBLIC sbc_rnf_div! routine called in sshwzv module30 PUBLIC sbc_rnf ! routine call in sbcmod module 31 PUBLIC sbc_rnf_div ! routine called in sshwzv module 32 32 33 33 ! !!* namsbc_rnf namelist * … … 43 43 TYPE(FLD_N) :: sn_dep_rnf !: information about the depth which river inflow affects 44 44 LOGICAL , PUBLIC :: ln_rnf_mouth = .false. !: specific treatment in mouths vicinity 45 REAL(wp) , PUBLIC :: rn_hrnf = 0.e0 !: runoffs, depth over which enhanced vertical mixing is used 46 REAL(wp) , PUBLIC :: rn_avt_rnf = 0.e0 !: runoffs, value of the additional vertical mixing coef. [m2/s] 47 REAL(wp) , PUBLIC :: rn_rfact = 1.e0 !: multiplicative factor for runoff 48 49 INTEGER , PUBLIC :: nkrnf = 0 !: number of levels over which Kz is increased at river mouths 50 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: rnfmsk !: river mouth mask (hori.) 51 REAL(wp), PUBLIC, DIMENSION(jpk) :: rnfmsk_z !: river mouth mask (vert.) 52 53 TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_rnf !: structure: river runoff (file information, fields read) 54 55 TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_s_rnf !: structure: river runoff salinity (file information, fields read) 56 TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_t_rnf !: structure: river runoff temperature (file information, fields read) 45 REAL(wp) , PUBLIC :: rn_hrnf = 0._wp !: runoffs, depth over which enhanced vertical mixing is used 46 REAL(wp) , PUBLIC :: rn_avt_rnf = 0._wp !: runoffs, value of the additional vertical mixing coef. [m2/s] 47 REAL(wp) , PUBLIC :: rn_rfact = 1._wp !: multiplicative factor for runoff 48 49 INTEGER , PUBLIC :: nkrnf = 0 !: nb of levels over which Kz is increased at river mouths 50 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: rnfmsk !: river mouth mask (hori.) 51 REAL(wp), PUBLIC, DIMENSION(jpk) :: rnfmsk_z !: river mouth mask (vert.) 52 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: h_rnf !: depth of runoff in m 53 INTEGER, PUBLIC, DIMENSION(jpi,jpj) :: nk_rnf !: depth of runoff in model levels 54 REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpts) :: rnf_tsc_b, rnf_tsc !: before and now T & S contents of runoffs [K.m/s & PSU.m/s] 55 56 REAL(wp) :: r1_rau0 ! = 1 / rau0 57 58 TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_rnf ! structure: river runoff (file information, fields read) 59 TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_s_rnf ! structure: river runoff salinity (file information, fields read) 60 TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_t_rnf ! structure: river runoff temperature (file information, fields read) 57 61 58 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: h_rnf !: depth of runoff in m59 INTEGER, PUBLIC, DIMENSION(jpi,jpj) :: nk_rnf !: depth of runoff in model levels60 61 REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpts) :: rnf_tsc_b, rnf_tsc !: before and now T & S contents of runoffs [K.m/s & PSU.m/s]62 63 62 !! * Substitutions 64 63 # include "domzgr_substitute.h90" … … 85 84 !! 86 85 INTEGER :: ji, jj ! dummy loop indices 87 REAL(wp) :: z1_rau0 ! local scalar88 86 !!---------------------------------------------------------------------- 89 87 ! … … 105 103 IF( ln_rnf_tem ) CALL fld_read ( kt, nn_fsbc, sf_t_rnf ) ! idem for runoffs temperature if required 106 104 IF( ln_rnf_sal ) CALL fld_read ( kt, nn_fsbc, sf_s_rnf ) ! idem for runoffs salinity if required 107 105 ! 108 106 ! Runoff reduction only associated to the ORCA2_LIM configuration 109 107 ! when reading the NetCDF file runoff_1m_nomask.nc 110 108 IF( cp_cfg == 'orca' .AND. jp_cfg == 2 ) THEN 111 DO jj = 1, jpj 112 DO ji = 1, jpi 113 IF( gphit(ji,jj) > 40 .AND. gphit(ji,jj) < 65 ) sf_rnf(1)%fnow(ji,jj,1) = 0.85 * sf_rnf(1)%fnow(ji,jj,1) 114 END DO 115 END DO 116 ENDIF 117 118 ! C a u t i o n : runoff is negative and in kg/m2/s 109 WHERE( 40._wp < gphit(:,:) .AND. gphit(:,:) < 65._wp ) 110 sf_rnf(1)%fnow(:,:,1) = 0.85 * sf_rnf(1)%fnow(:,:,1) 111 END WHERE 112 ENDIF 113 ! 119 114 IF( MOD( kt - 1, nn_fsbc ) == 0 ) THEN 120 rnf(:,:) 115 rnf(:,:) = rn_rfact * ( sf_rnf(1)%fnow(:,:,1) ) 121 116 ! 122 z1_rau0 = 1.e0/ rau0123 ! 124 IF( ln_rnf_tem ) THEN! use runoffs temperature data125 rnf_tsc(:,:,jp_tem) = ( sf_t_rnf(1)%fnow(:,:,1) ) * rnf(:,:) * z1_rau0126 WHERE( sf_t_rnf(1)%fnow(:,:,1) == -999 ) 127 rnf_tsc(:,:,jp_tem) = sst_m(:,:) * rnf(:,:) * z1_rau0128 END WHERE129 ELSE 130 rnf_tsc(:,:,jp_tem) = sst_m(:,:) * rnf(:,:) * z1_rau0117 r1_rau0 = 1._wp / rau0 118 ! ! set temperature & salinity content of runoffs 119 IF( ln_rnf_tem ) THEN ! use runoffs temperature data 120 rnf_tsc(:,:,jp_tem) = ( sf_t_rnf(1)%fnow(:,:,1) ) * rnf(:,:) * r1_rau0 121 WHERE( sf_t_rnf(1)%fnow(:,:,1) == -999 ) ! if missing data value use SST as runoffs temperature 122 rnf_tsc(:,:,jp_tem) = sst_m(:,:) * rnf(:,:) * r1_rau0 123 END WHERE 124 ELSE ! use SST as runoffs temperature 125 rnf_tsc(:,:,jp_tem) = sst_m(:,:) * rnf(:,:) * r1_rau0 131 126 ENDIF 132 ! 133 IF( ln_rnf_sal ) rnf_tsc(:,:,jp_sal) = ( sf_s_rnf(1)%fnow(:,:,1) ) * rnf(:,:) * z1_rau0134 ! 127 ! ! use runoffs salinity data 128 IF( ln_rnf_sal ) rnf_tsc(:,:,jp_sal) = ( sf_s_rnf(1)%fnow(:,:,1) ) * rnf(:,:) * r1_rau0 129 ! ! else use S=0 for runoffs (done one for all in the init) 135 130 ! 136 IF( ln_rnf_tem .OR. ln_rnf_sal ) THEN 137 WHERE( rnf(:,:) < 0. e0 )! example baltic model when flow is out of domain138 rnf_tsc(:,:,jp_tem) = sst_m(:,:) * rnf(:,:) * z1_rau0139 rnf_tsc(:,:,jp_sal) = sss_m(:,:) * rnf(:,:) * z1_rau0140 END WHERE131 IF( ln_rnf_tem .OR. ln_rnf_sal ) THEN ! runoffs as outflow: use ocean SST and SSS 132 WHERE( rnf(:,:) < 0._wp ) ! example baltic model when flow is out of domain 133 rnf_tsc(:,:,jp_tem) = sst_m(:,:) * rnf(:,:) * r1_rau0 134 rnf_tsc(:,:,jp_sal) = sss_m(:,:) * rnf(:,:) * r1_rau0 135 END WHERE 141 136 ENDIF 142 137 ! … … 190 185 !! 191 186 INTEGER :: ji, jj, jk ! dummy loop indices 192 REAL(wp) :: z1_rau0 ! local scalar187 REAL(wp) :: r1_rau0 ! local scalar 193 188 REAL(wp) :: zfact ! local scalar 194 189 !!---------------------------------------------------------------------- 195 190 ! 196 zfact = 0.5 e0197 ! 198 z1_rau0 = 1.e0/ rau0191 zfact = 0.5_wp 192 ! 193 r1_rau0 = 1._wp / rau0 199 194 IF( ln_rnf_depth ) THEN !== runoff distributed over several levels ==! 200 195 IF( lk_vvl ) THEN ! variable volume case 201 196 DO jj = 1, jpj ! update the depth over which runoffs are distributed 202 197 DO ji = 1, jpi 203 h_rnf(ji,jj) = 0. e0198 h_rnf(ji,jj) = 0._wp 204 199 DO jk = 1, nk_rnf(ji,jj) ! recalculates h_rnf to be the depth in metres 205 200 h_rnf(ji,jj) = h_rnf(ji,jj) + fse3t(ji,jj,jk) ! to the bottom of the relevant grid box … … 207 202 ! ! apply the runoff input flow 208 203 DO jk = 1, nk_rnf(ji,jj) 209 phdivn(ji,jj,jk) = phdivn(ji,jj,jk) - ( rnf(ji,jj) + rnf_b(ji,jj) ) * zfact * z1_rau0 / h_rnf(ji,jj)204 phdivn(ji,jj,jk) = phdivn(ji,jj,jk) - ( rnf(ji,jj) + rnf_b(ji,jj) ) * zfact * r1_rau0 / h_rnf(ji,jj) 210 205 END DO 211 206 END DO … … 215 210 DO ji = 1, jpi 216 211 DO jk = 1, nk_rnf(ji,jj) 217 phdivn(ji,jj,jk) = phdivn(ji,jj,jk) - ( rnf(ji,jj) + rnf_b(ji,jj) ) * zfact * z1_rau0 / h_rnf(ji,jj)212 phdivn(ji,jj,jk) = phdivn(ji,jj,jk) - ( rnf(ji,jj) + rnf_b(ji,jj) ) * zfact * r1_rau0 / h_rnf(ji,jj) 218 213 END DO 219 214 END DO … … 224 219 h_rnf(:,:) = fse3t(:,:,1) ! recalculate h_rnf to be depth of top box 225 220 ENDIF 226 phdivn(:,:,1) = phdivn(:,:,1) - ( rnf(:,:) + rnf_b(:,:) ) * zfact * z1_rau0 / fse3t(:,:,1)221 phdivn(:,:,1) = phdivn(:,:,1) - ( rnf(:,:) + rnf_b(:,:) ) * zfact * r1_rau0 / fse3t(:,:,1) 227 222 ENDIF 228 223 ! … … 303 298 CALL fld_fill( sf_rnf, (/ sn_rnf /), cn_dir, 'sbc_rnf_init', 'read runoffs data', 'namsbc_rnf' ) 304 299 ! 305 IF( ln_rnf_tem ) THEN ! Create (if required) sf_t_rnf structure300 IF( ln_rnf_tem ) THEN ! Create (if required) sf_t_rnf structure 306 301 IF(lwp) WRITE(numout,*) 307 302 IF(lwp) WRITE(numout,*) ' runoffs temperatures read in a file' … … 326 321 CALL fld_fill (sf_s_rnf, (/ sn_s_rnf /), cn_dir, 'sbc_rnf_init', 'read runoff salinity data', 'namsbc_rnf' ) 327 322 ENDIF 328 329 330 IF ( ln_rnf_depth ) THEN ! depth of runoffs set from a file 323 ! 324 IF( ln_rnf_depth ) THEN ! depth of runoffs set from a file 331 325 IF(lwp) WRITE(numout,*) 332 326 IF(lwp) WRITE(numout,*) ' runoffs depth read in a file' 333 327 rn_dep_file = TRIM( cn_dir )//TRIM( sn_dep_rnf%clname ) 334 328 CALL iom_open ( rn_dep_file, inum ) ! open file 335 CALL iom_get ( inum, jpdom_data, sn_dep_rnf%clvar, h_rnf ) 336 CALL iom_close( inum ) ! close file337 338 nk_rnf(:,:) = 0 ! set the number of level over which river runoffs are applied329 CALL iom_get ( inum, jpdom_data, sn_dep_rnf%clvar, h_rnf ) ! read the river mouth array 330 CALL iom_close( inum ) ! close file 331 ! 332 nk_rnf(:,:) = 0 ! set the number of level over which river runoffs are applied 339 333 DO jj = 1, jpj 340 DO ji = 1, jpi341 IF ( h_rnf(ji,jj) > 0.e0) THEN342 jk = 2343 DO WHILE ( jk /= ( mbathy(ji,jj) - 1 ) .AND. fsdept(ji,jj,jk) < h_rnf(ji,jj) ) ; jk = jk + 1 ; ENDDO344 nk_rnf(ji,jj) = jk345 ELSE IF( h_rnf(ji,jj) == -1 ) THEN ; nk_rnf(ji,jj) = 1346 ELSE IF ( h_rnf(ji,jj) == -999 ) THEN ; nk_rnf(ji,jj) = mbathy(ji,jj) - 1347 ELSE IF ( h_rnf(ji,jj) /= 0) THEN348 CALL ctl_stop( 'runoff depth not positive, and not -999 or -1, rnf value in file fort.999' )349 WRITE(999,*) 'ji, jj, rnf(ji,jj) :', ji, jj, rnf(ji,jj)350 ENDIF351 ENDDO352 END DO353 DO jj = 1, jpj ! set the associated depth354 DO ji = 1, jpi355 h_rnf(ji,jj) = 0.e0356 DO jk = 1, nk_rnf(ji,jj)357 h_rnf(ji,jj) = h_rnf(ji,jj)+fse3t(ji,jj,jk)358 ENDDO359 ENDDO360 END DO334 DO ji = 1, jpi 335 IF( h_rnf(ji,jj) > 0._wp ) THEN 336 jk = 2 337 DO WHILE ( jk /= mbkt(ji,jj) .AND. fsdept(ji,jj,jk) < h_rnf(ji,jj) ) ; jk = jk + 1 ; END DO 338 nk_rnf(ji,jj) = jk 339 ELSEIF( h_rnf(ji,jj) == -1 ) THEN ; nk_rnf(ji,jj) = 1 340 ELSEIF( h_rnf(ji,jj) == -999 ) THEN ; nk_rnf(ji,jj) = mbkt(ji,jj) 341 ELSEIF( h_rnf(ji,jj) /= 0 ) THEN 342 CALL ctl_stop( 'runoff depth not positive, and not -999 or -1, rnf value in file fort.999' ) 343 WRITE(999,*) 'ji, jj, rnf(ji,jj) :', ji, jj, rnf(ji,jj) 344 ENDIF 345 END DO 346 END DO 347 DO jj = 1, jpj ! set the associated depth 348 DO ji = 1, jpi 349 h_rnf(ji,jj) = 0._wp 350 DO jk = 1, nk_rnf(ji,jj) 351 h_rnf(ji,jj) = h_rnf(ji,jj) + fse3t(ji,jj,jk) 352 END DO 353 END DO 354 END DO 361 355 ELSE ! runoffs applied at the surface 362 356 nk_rnf(:,:) = 1 363 h_rnf (:,:)= fse3t(:,:,1)357 h_rnf (:,:) = fse3t(:,:,1) 364 358 ENDIF 365 ! 366 ENDIF 367 368 rnf_tsc(:,:,:) = 0.e0 ! runoffs temperature & salinty contents initilisation 359 ! 360 ENDIF 361 ! 362 rnf_tsc(:,:,:) = 0._wp ! runoffs temperature & salinty contents initilisation 363 ! 369 364 ! ! ======================== 370 365 ! ! River mouth vicinity … … 380 375 ! 381 376 nkrnf = 0 ! Number of level over which Kz increase 382 IF( rn_hrnf > 0. e0) THEN377 IF( rn_hrnf > 0._wp ) THEN 383 378 nkrnf = 2 384 379 DO WHILE( nkrnf /= jpkm1 .AND. gdepw_0(nkrnf+1) < rn_hrnf ) ; nkrnf = nkrnf + 1 ; END DO … … 398 393 IF(lwp) WRITE(numout,*) 399 394 IF(lwp) WRITE(numout,*) ' No specific treatment at river mouths' 400 rnfmsk (:,:) = 0. e0401 rnfmsk_z(:) = 0. e0395 rnfmsk (:,:) = 0._wp 396 rnfmsk_z(:) = 0._wp 402 397 nkrnf = 0 403 398 ENDIF … … 447 442 IF( nclosea == 1 ) CALL clo_rnf( rnfmsk ) ! closed sea inflow set as ruver mouth 448 443 449 rnfmsk_z(:) = 0. e0! vertical structure444 rnfmsk_z(:) = 0._wp ! vertical structure 450 445 rnfmsk_z(1) = 1.0 451 446 rnfmsk_z(2) = 1.0 ! ********** -
branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_bilap.F90
r2399 r2450 84 84 !! 85 85 INTEGER :: ji, jj, jk, jn ! dummy loop indices 86 INTEGER :: iku, ikv ! local integers87 86 REAL(wp) :: zbtr, ztra ! local scalars 88 87 REAL(wp), DIMENSION(jpi,jpj) :: zeeu, zeev, zlt ! 2D workspace … … 116 115 END DO 117 116 END DO 118 IF( ln_zps ) THEN ! set gradient at partial step level 117 IF( ln_zps ) THEN ! set gradient at partial step level (last ocean level) 119 118 DO jj = 1, jpjm1 120 119 DO ji = 1, jpim1 121 ! last level 122 iku = MIN ( mbathy(ji,jj), mbathy(ji+1,jj ) ) - 1 123 ikv = MIN ( mbathy(ji,jj), mbathy(ji ,jj+1) ) - 1 124 IF( iku == jk ) ztu(ji,jj,jk) = zeeu(ji,jj) * pgu(ji,jj,jn) 125 IF( ikv == jk ) ztv(ji,jj,jk) = zeev(ji,jj) * pgv(ji,jj,jn) 120 IF( mbku(ji,jj) == jk ) ztu(ji,jj,jk) = zeeu(ji,jj) * pgu(ji,jj,jn) 121 IF( mbkv(ji,jj) == jk ) ztv(ji,jj,jk) = zeev(ji,jj) * pgv(ji,jj,jn) 126 122 END DO 127 123 END DO -
branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_iso.F90
r2399 r2450 101 101 REAL(wp) , INTENT(in ) :: pahtb0 ! background diffusion coef 102 102 !! 103 INTEGER :: ji, jj, jk,jn ! dummy loop indices 104 INTEGER :: iku, ikv ! temporary integer 103 INTEGER :: ji, jj, jk, jn ! dummy loop indices 105 104 REAL(wp) :: zmsku, zabe1, zcof1, zcoef3 ! local scalars 106 105 REAL(wp) :: zmskv, zabe2, zcof2, zcoef4 ! - - … … 141 140 END DO 142 141 END DO 143 IF( ln_zps ) THEN ! partial steps correction at the last level142 IF( ln_zps ) THEN ! partial steps correction at the last ocean level 144 143 DO jj = 1, jpjm1 145 144 DO ji = 1, fs_jpim1 ! vector opt. 146 ! last level 147 iku = MIN( mbathy(ji,jj), mbathy(ji+1,jj ) ) - 1 148 ikv = MIN( mbathy(ji,jj), mbathy(ji ,jj+1) ) - 1 149 zdit(ji,jj,iku) = pgu(ji,jj,jn) 150 zdjt(ji,jj,ikv) = pgv(ji,jj,jn) 145 zdit(ji,jj,mbku(ji,jj)) = pgu(ji,jj,jn) 146 zdjt(ji,jj,mbkv(ji,jj)) = pgv(ji,jj,jn) 151 147 END DO 152 148 END DO … … 215 211 #if defined key_diaar5 216 212 IF( cdtype == 'TRA' .AND. jn == jp_tem ) THEN 217 z ztmp = 0.5 * rau0 * rcp218 z 2d(:,:) = 0.e0213 z2d(:,:) = 0._wp 214 zztmp = rau0 * rcp 219 215 DO jk = 1, jpkm1 220 216 DO jj = 2, jpjm1 221 217 DO ji = fs_2, fs_jpim1 ! vector opt. 222 z2d(ji,jj) = z2d(ji,jj) + zztmp * zftu(ji,jj,jk) & 223 & * ( ptn(ji,jj,jk,jn) + ptn(ji+1,jj,jk,jn) ) * e1u(ji,jj) * fse3u(ji,jj,jk) 218 z2d(ji,jj) = z2d(ji,jj) + zftu(ji,jj,jk) 224 219 END DO 225 220 END DO 226 221 END DO 222 z2d(:,:) = zztmp * z2d(:,:) 227 223 CALL lbc_lnk( z2d, 'U', -1. ) 228 224 CALL iom_put( "udiff_heattr", z2d ) ! heat transport in i-direction 229 z2d(:,:) = 0.e0230 225 DO jk = 1, jpkm1 231 226 DO jj = 2, jpjm1 232 227 DO ji = fs_2, fs_jpim1 ! vector opt. 233 z2d(ji,jj) = z2d(ji,jj) + zztmp * zftv(ji,jj,jk) & 234 & * ( ptn(ji,jj,jk,jn) + ptn(ji,jj+1,jk,jn) ) * e2v(ji,jj) * fse3v(ji,jj,jk) 228 z2d(ji,jj) = z2d(ji,jj) + zftv(ji,jj,jk) 235 229 END DO 236 230 END DO 237 231 END DO 232 z2d(:,:) = zztmp * z2d(:,:) 238 233 CALL lbc_lnk( z2d, 'V', -1. ) 239 234 CALL iom_put( "vdiff_heattr", z2d ) ! heat transport in i-direction -
branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_iso_grif.F90
r2399 r2450 53 53 SUBROUTINE tra_ldf_iso_grif( kt, cdtype, pgu, pgv, & 54 54 & ptb, pta, kjpt, pahtb0 ) 55 !!---------------------------------------------------------------------- 56 !! *** ROUTINE tra_ldf_iso_grif *** 57 !! 58 !! ** Purpose : Compute the before horizontal tracer (t & s) diffusive 59 !! trend for a laplacian tensor (ezxcept the dz[ dz[.] ] term) and 60 !! add it to the general trend of tracer equation. 61 !! 62 !! ** Method : The horizontal component of the lateral diffusive trends 63 !! is provided by a 2nd order operator rotated along neural or geopo- 64 !! tential surfaces to which an eddy induced advection can be added 65 !! It is computed using before fields (forward in time) and isopyc- 66 !! nal or geopotential slopes computed in routine ldfslp. 67 !! 68 !! 1st part : masked horizontal derivative of T ( di[ t ] ) 69 !! ======== with partial cell update if ln_zps=T. 70 !! 71 !! 2nd part : horizontal fluxes of the lateral mixing operator 72 !! ======== 73 !! zftu = (aht+ahtb0) e2u*e3u/e1u di[ tb ] 74 !! - aht e2u*uslp dk[ mi(mk(tb)) ] 75 !! zftv = (aht+ahtb0) e1v*e3v/e2v dj[ tb ] 76 !! - aht e2u*vslp dk[ mj(mk(tb)) ] 77 !! take the horizontal divergence of the fluxes: 78 !! difft = 1/(e1t*e2t*e3t) { di-1[ zftu ] + dj-1[ zftv ] } 79 !! Add this trend to the general trend (ta,sa): 80 !! ta = ta + difft 81 !! 82 !! 3rd part: vertical trends of the lateral mixing operator 83 !! ======== (excluding the vertical flux proportional to dk[t] ) 84 !! vertical fluxes associated with the rotated lateral mixing: 85 !! zftw =-aht { e2t*wslpi di[ mi(mk(tb)) ] 86 !! + e1t*wslpj dj[ mj(mk(tb)) ] } 87 !! take the horizontal divergence of the fluxes: 88 !! difft = 1/(e1t*e2t*e3t) dk[ zftw ] 89 !! Add this trend to the general trend (ta,sa): 90 !! pta = pta + difft 91 !! 92 !! ** Action : Update pta arrays with the before rotated diffusion 93 !!---------------------------------------------------------------------- 94 USE oce, zftu => ua ! use ua as workspace 95 USE oce, zftv => va ! use va as workspace 96 !! 97 INTEGER , INTENT(in ) :: kt ! ocean time-step index 98 CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator) 99 INTEGER , INTENT(in ) :: kjpt ! number of tracers 100 REAL(wp), DIMENSION(jpi,jpj ,kjpt), INTENT(in ) :: pgu, pgv ! tracer gradient at pstep levels 101 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: ptb ! before and now tracer fields 102 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pta ! tracer trend 103 REAL(wp) , INTENT(in ) :: pahtb0 ! background diffusion coef 104 !! 105 INTEGER :: ji, jj, jk,jn ! dummy loop indices 106 INTEGER :: ip,jp,kp ! dummy loop indices 107 INTEGER :: iku, ikv ! temporary integer 108 INTEGER :: ierr ! temporary integer 109 REAL(wp) :: zmsku, zabe1, zcof1, zcoef3 ! local scalars 110 REAL(wp) :: zmskv, zabe2, zcof2, zcoef4 ! - - 111 REAL(wp) :: zcoef0, zbtr ! - - 112 REAL(wp), DIMENSION(jpi,jpj,0:1) :: zdkt ! 2D+1 workspace 113 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdit, zdjt, ztfw ! 3D workspace 55 !!---------------------------------------------------------------------- 56 !! *** ROUTINE tra_ldf_iso_grif *** 57 !! 58 !! ** Purpose : Compute the before horizontal tracer (t & s) diffusive 59 !! trend for a laplacian tensor (ezxcept the dz[ dz[.] ] term) and 60 !! add it to the general trend of tracer equation. 61 !! 62 !! ** Method : The horizontal component of the lateral diffusive trends 63 !! is provided by a 2nd order operator rotated along neural or geopo- 64 !! tential surfaces to which an eddy induced advection can be added 65 !! It is computed using before fields (forward in time) and isopyc- 66 !! nal or geopotential slopes computed in routine ldfslp. 67 !! 68 !! 1st part : masked horizontal derivative of T ( di[ t ] ) 69 !! ======== with partial cell update if ln_zps=T. 70 !! 71 !! 2nd part : horizontal fluxes of the lateral mixing operator 72 !! ======== 73 !! zftu = (aht+ahtb0) e2u*e3u/e1u di[ tb ] 74 !! - aht e2u*uslp dk[ mi(mk(tb)) ] 75 !! zftv = (aht+ahtb0) e1v*e3v/e2v dj[ tb ] 76 !! - aht e2u*vslp dk[ mj(mk(tb)) ] 77 !! take the horizontal divergence of the fluxes: 78 !! difft = 1/(e1t*e2t*e3t) { di-1[ zftu ] + dj-1[ zftv ] } 79 !! Add this trend to the general trend (ta,sa): 80 !! ta = ta + difft 81 !! 82 !! 3rd part: vertical trends of the lateral mixing operator 83 !! ======== (excluding the vertical flux proportional to dk[t] ) 84 !! vertical fluxes associated with the rotated lateral mixing: 85 !! zftw =-aht { e2t*wslpi di[ mi(mk(tb)) ] 86 !! + e1t*wslpj dj[ mj(mk(tb)) ] } 87 !! take the horizontal divergence of the fluxes: 88 !! difft = 1/(e1t*e2t*e3t) dk[ zftw ] 89 !! Add this trend to the general trend (ta,sa): 90 !! pta = pta + difft 91 !! 92 !! ** Action : Update pta arrays with the before rotated diffusion 93 !!---------------------------------------------------------------------- 94 USE oce, zftu => ua ! use ua as workspace 95 USE oce, zftv => va ! use va as workspace 96 !! 97 INTEGER , INTENT(in ) :: kt ! ocean time-step index 98 CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator) 99 INTEGER , INTENT(in ) :: kjpt ! number of tracers 100 REAL(wp), DIMENSION(jpi,jpj ,kjpt), INTENT(in ) :: pgu, pgv ! tracer gradient at pstep levels 101 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: ptb ! before and now tracer fields 102 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pta ! tracer trend 103 REAL(wp) , INTENT(in ) :: pahtb0 ! background diffusion coef 104 !! 105 INTEGER :: ji, jj, jk,jn ! dummy loop indices 106 INTEGER :: ip,jp,kp ! dummy loop indices 107 INTEGER :: ierr ! temporary integer 108 REAL(wp) :: zmsku, zabe1, zcof1, zcoef3 ! local scalars 109 REAL(wp) :: zmskv, zabe2, zcof2, zcoef4 ! - - 110 REAL(wp) :: zcoef0, zbtr ! - - 111 REAL(wp), DIMENSION(jpi,jpj,0:1) :: zdkt ! 2D+1 workspace 112 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdit, zdjt, ztfw ! 3D workspace 114 113 ! 115 REAL(wp) :: zslope_skew,zslope_iso,zslope2, zbu, zbv116 REAL(wp) :: ze1ur,zdxt,ze2vr,ze3wr,zdyt,zdzt117 REAL(wp) :: ah,zah_slp,zaei_slp114 REAL(wp) :: zslope_skew,zslope_iso,zslope2, zbu, zbv 115 REAL(wp) :: ze1ur,zdxt,ze2vr,ze3wr,zdyt,zdzt 116 REAL(wp) :: ah,zah_slp,zaei_slp 118 117 #if defined key_diaar5 119 REAL(wp), DIMENSION(jpi,jpj) :: z2d ! 2D workspace120 REAL(wp) :: zztmp ! local scalar118 REAL(wp), DIMENSION(jpi,jpj) :: z2d ! 2D workspace 119 REAL(wp) :: zztmp ! local scalar 121 120 #endif 122 121 !!---------------------------------------------------------------------- 123 122 124 IF( kt == nit000 ) THEN 125 IF(lwp) WRITE(numout,*) 126 IF(lwp) WRITE(numout,*) 'tra_ldf_iso_grif : rotated laplacian diffusion operator on ', cdtype 127 IF(lwp) WRITE(numout,*) ' WARNING: STILL UNDER TEST, NOT RECOMMENDED. USE AT YOUR OWN PERIL' 128 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 129 ALLOCATE( ah_wslp2(jpi,jpj,jpk) , STAT=ierr ) 130 IF( ierr > 0 ) THEN 131 CALL ctl_stop( 'tra_ldf_iso_grif : unable to allocate Griffies operator ah_wslp2 ' ) ; RETURN 132 ENDIF 133 IF( ln_traldf_gdia ) THEN 134 ALLOCATE( psix_eiv(jpi,jpj,jpk) , psiy_eiv(jpi,jpj,jpk) , STAT=ierr ) 135 IF( ierr > 0 ) THEN 136 CALL ctl_stop( 'tra_ldf_iso_grif : unable to allocate Griffies operator diagnostics ' ) ; RETURN 137 ENDIF 138 ENDIF 139 ENDIF 140 141 ! 123 IF( kt == nit000 ) THEN 124 IF(lwp) WRITE(numout,*) 125 IF(lwp) WRITE(numout,*) 'tra_ldf_iso_grif : rotated laplacian diffusion operator on ', cdtype 126 IF(lwp) WRITE(numout,*) ' WARNING: STILL UNDER TEST, NOT RECOMMENDED. USE AT YOUR OWN PERIL' 127 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 128 ALLOCATE( ah_wslp2(jpi,jpj,jpk) , STAT=ierr ) 129 IF( ierr > 0 ) THEN 130 CALL ctl_stop( 'tra_ldf_iso_grif : unable to allocate Griffies operator ah_wslp2 ' ) ; RETURN 131 ENDIF 132 IF( ln_traldf_gdia ) THEN 133 ALLOCATE( psix_eiv(jpi,jpj,jpk) , psiy_eiv(jpi,jpj,jpk) , STAT=ierr ) 134 IF( ierr > 0 ) THEN 135 CALL ctl_stop( 'tra_ldf_iso_grif : unable to allocate Griffies operator diagnostics ' ) ; RETURN 136 ENDIF 137 ENDIF 138 ENDIF 139 142 140 !!---------------------------------------------------------------------- 143 141 !! 0 - calculate ah_wslp2, psix_eiv, psiy_eiv … … 146 144 !!gm Future development: consider using Ah defined at T-points and attached to the 4 t-point triads 147 145 148 ah_wslp2(:,:,:) = 0.0 149 IF( ln_traldf_gdia ) THEN 150 psix_eiv(:,:,:) = 0.0 151 psiy_eiv(:,:,:) = 0.0 152 ENDIF 153 154 DO ip=0,1 155 DO kp=0,1 156 DO jk=1,jpkm1 157 DO jj=1,jpjm1 158 DO ji=1,fs_jpim1 159 ze3wr=1.0_wp/fse3w(ji+ip,jj,jk+kp) 160 zbu = 0.25_wp*e1u(ji,jj)*e2u(ji,jj)*fse3u(ji,jj,jk) 161 ah = fsahtu(ji,jj,jk) ! fsaht(ji+ip,jj,jk) 162 zslope_skew = triadi_g(ji+ip,jj,jk,1-ip,kp) 163 zslope2=(zslope_skew & 164 & - umask(ji,jj,jk+kp)*(fsdept(ji+1,jj,jk) - fsdept(ji ,jj ,jk))*ze1ur & 146 ah_wslp2(:,:,:) = 0.0 147 IF( ln_traldf_gdia ) THEN 148 psix_eiv(:,:,:) = 0.0 149 psiy_eiv(:,:,:) = 0.0 150 ENDIF 151 152 DO ip=0,1 153 DO kp=0,1 154 DO jk=1,jpkm1 155 DO jj=1,jpjm1 156 DO ji=1,fs_jpim1 157 ze3wr=1.0_wp/fse3w(ji+ip,jj,jk+kp) 158 zbu = 0.25_wp*e1u(ji,jj)*e2u(ji,jj)*fse3u(ji,jj,jk) 159 ah = fsahtu(ji,jj,jk) ! fsaht(ji+ip,jj,jk) 160 zslope_skew = triadi_g(ji+ip,jj,jk,1-ip,kp) 161 zslope2 = ( zslope_skew & 162 & - umask(ji,jj,jk+kp)*(fsdept(ji+1,jj,jk) - fsdept(ji ,jj ,jk))*ze1ur )**2 163 ah_wslp2(ji+ip,jj,jk+kp)=ah_wslp2(ji+ip,jj,jk+kp) + & 164 & ah*((zbu*ze3wr)/(e1t(ji+ip,jj)*e2t(ji+ip,jj)))*zslope2 165 IF( ln_traldf_gdia ) THEN 166 zaei_slp = fsaeiw(ji+ip,jj,jk)*zslope_skew!fsaeit(ji+ip,jj,jk)*zslope_skew 167 psix_eiv(ji,jj,jk+kp) = psix_eiv(ji,jj,jk+kp) + 0.25_wp*zaei_slp 168 ENDIF 169 END DO 170 END DO 171 END DO 172 END DO 173 END DO 174 ! 175 DO jp=0,1 176 DO kp=0,1 177 DO jk=1,jpkm1 178 DO jj=1,jpjm1 179 DO ji=1,fs_jpim1 180 ze3wr=1.0_wp/fse3w(ji,jj+jp,jk+kp) 181 zbv = 0.25_wp*e1v(ji,jj)*e2v(ji,jj)*fse3v(ji,jj,jk) 182 ah = fsahtu(ji,jj,jk)!fsaht(ji,jj+jp,jk) 183 zslope_skew = triadj_g(ji,jj+jp,jk,1-jp,kp) 184 zslope2=(zslope_skew - vmask(ji,jj,jk+kp)*(fsdept(ji,jj+1,jk) - fsdept(ji ,jj ,jk))*ze2vr & 165 185 & )**2 166 ah_wslp2(ji+ip,jj,jk+kp)=ah_wslp2(ji+ip,jj,jk+kp) + & 167 & ah*((zbu*ze3wr)/(e1t(ji+ip,jj)*e2t(ji+ip,jj)))*zslope2 168 IF( ln_traldf_gdia ) THEN 169 zaei_slp = fsaeiw(ji+ip,jj,jk)*zslope_skew!fsaeit(ji+ip,jj,jk)*zslope_skew 170 psix_eiv(ji,jj,jk+kp) = psix_eiv(ji,jj,jk+kp) + 0.25_wp*zaei_slp 171 ENDIF 172 END DO 173 END DO 174 END DO 175 END DO 176 END DO 177 178 DO jp=0,1 179 DO kp=0,1 180 DO jk=1,jpkm1 181 DO jj=1,jpjm1 182 DO ji=1,fs_jpim1 183 ze3wr=1.0_wp/fse3w(ji,jj+jp,jk+kp) 184 zbv = 0.25_wp*e1v(ji,jj)*e2v(ji,jj)*fse3v(ji,jj,jk) 185 ah = fsahtu(ji,jj,jk)!fsaht(ji,jj+jp,jk) 186 zslope_skew = triadj_g(ji,jj+jp,jk,1-jp,kp) 187 zslope2=(zslope_skew - vmask(ji,jj,jk+kp)*(fsdept(ji,jj+1,jk) - fsdept(ji ,jj ,jk))*ze2vr & 188 & )**2 189 ah_wslp2(ji,jj+jp,jk+kp)=ah_wslp2(ji,jj+jp,jk+kp) + & 186 ah_wslp2(ji,jj+jp,jk+kp)=ah_wslp2(ji,jj+jp,jk+kp) + & 190 187 & ah*((zbv*ze3wr)/(e1t(ji,jj+jp)*e2t(ji,jj+jp)))*zslope2 191 IF( ln_traldf_gdia ) THEN192 zaei_slp = fsaeiw(ji,jj+jp,jk) * zslope_skew !fsaeit(ji,jj+jp,jk)*zslope_skew193 psiy_eiv(ji,jj,jk+kp) = psiy_eiv(ji,jj,jk+kp) + 0.25_wp*zaei_slp194 ENDIF195 END DO196 END DO197 END DO198 END DO188 IF( ln_traldf_gdia ) THEN 189 zaei_slp = fsaeiw(ji,jj+jp,jk) * zslope_skew !fsaeit(ji,jj+jp,jk)*zslope_skew 190 psiy_eiv(ji,jj,jk+kp) = psiy_eiv(ji,jj,jk+kp) + 0.25_wp*zaei_slp 191 ENDIF 192 END DO 193 END DO 194 END DO 195 END DO 199 196 END DO 200 201 197 ! 202 198 ! ! =========== … … 224 220 DO ji = 1, jpim1 225 221 # endif 226 iku = MIN( mbathy(ji,jj), mbathy(ji+1,jj ) ) - 1 ! last level 227 ikv = MIN( mbathy(ji,jj), mbathy(ji ,jj+1) ) - 1 228 zdit(ji,jj,iku) = pgu(ji,jj,jn) 229 zdjt(ji,jj,ikv) = pgv(ji,jj,jn) 222 zdit(ji,jj,mbku(ji,jj)) = pgu(ji,jj,jn) 223 zdjt(ji,jj,mbkv(ji,jj)) = pgv(ji,jj,jn) 230 224 END DO 231 225 END DO … … 294 288 END DO 295 289 296 ! !== divergence and add to the general trend ==!297 DO jj = 2 , jpjm1298 DO ji = fs_2, fs_jpim1 ! vector opt.299 zbtr = 1._wp / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )300 pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + zbtr * ( zftu(ji-1,jj,jk) - zftu(ji,jj,jk) &301 &+ zftv(ji,jj-1,jk) - zftv(ji,jj,jk) )302 END DO303 END DO304 !305 END DO306 !307 DO jk = 1, jpkm1 !== Divergence of vertical fluxes added to the general tracer trend308 DO jj = 2, jpjm1309 DO ji = fs_2, fs_jpim1 ! vector opt.310 pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ( ztfw(ji,jj,jk+1) - ztfw(ji,jj,jk) ) &311 &/ ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )312 END DO313 END DO314 END DO315 !316 ! ! "Poleward" diffusive heat or salt transports (T-S case only)317 IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nn_fptr ) == 0 ) ) THEN318 IF( jn == jp_tem) htr_ldf(:) = ptr_vj( zftv(:,:,:) ) ! 3.3 names319 IF( jn == jp_sal) str_ldf(:) = ptr_vj( zftv(:,:,:) )320 ENDIF290 ! !== divergence and add to the general trend ==! 291 DO jj = 2 , jpjm1 292 DO ji = fs_2, fs_jpim1 ! vector opt. 293 zbtr = 1._wp / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 294 pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + zbtr * ( zftu(ji-1,jj,jk) - zftu(ji,jj,jk) & 295 & + zftv(ji,jj-1,jk) - zftv(ji,jj,jk) ) 296 END DO 297 END DO 298 ! 299 END DO 300 ! 301 DO jk = 1, jpkm1 !== Divergence of vertical fluxes added to the general tracer trend 302 DO jj = 2, jpjm1 303 DO ji = fs_2, fs_jpim1 ! vector opt. 304 pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ( ztfw(ji,jj,jk+1) - ztfw(ji,jj,jk) ) & 305 & / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 306 END DO 307 END DO 308 END DO 309 ! 310 ! ! "Poleward" diffusive heat or salt transports (T-S case only) 311 IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nn_fptr ) == 0 ) ) THEN 312 IF( jn == jp_tem) htr_ldf(:) = ptr_vj( zftv(:,:,:) ) ! 3.3 names 313 IF( jn == jp_sal) str_ldf(:) = ptr_vj( zftv(:,:,:) ) 314 ENDIF 321 315 322 316 #if defined key_diaar5 323 IF( cdtype == 'TRA' .AND. jn == jp_tem ) THEN 324 zztmp = 0.5 * rau0 * rcp 325 z2d(:,:) = 0.e0 326 DO jk = 1, jpkm1 327 DO jj = 2, jpjm1 328 DO ji = fs_2, fs_jpim1 ! vector opt. 329 z2d(ji,jj) = z2d(ji,jj) + zztmp * zftu(ji,jj,jk) & 330 & * ( ptn(ji,jj,jk,jn) + ptn(ji+1,jj,jk,jn) ) * e1u(ji,jj) * fse3u(ji,jj,jk) 331 END DO 332 END DO 333 END DO 334 CALL lbc_lnk( z2d, 'U', -1. ) 335 CALL iom_put( "udiff_heattr", z2d ) ! heat transport in i-direction 336 z2d(:,:) = 0.e0 337 DO jk = 1, jpkm1 338 DO jj = 2, jpjm1 339 DO ji = fs_2, fs_jpim1 ! vector opt. 340 z2d(ji,jj) = z2d(ji,jj) + zztmp * zftv(ji,jj,jk) & 341 & * ( ptn(ji,jj,jk,jn) + ptn(ji,jj+1,jk,jn) ) * e2v(ji,jj) * fse3v(ji,jj,jk) 342 END DO 343 END DO 344 END DO 345 CALL lbc_lnk( z2d, 'V', -1. ) 346 CALL iom_put( "vdiff_heattr", z2d ) ! heat transport in i-direction 347 END IF 317 IF( cdtype == 'TRA' .AND. jn == jp_tem ) THEN 318 z2d(:,:) = 0._wp 319 zztmp = rau0 * rcp 320 DO jk = 1, jpkm1 321 DO jj = 2, jpjm1 322 DO ji = fs_2, fs_jpim1 ! vector opt. 323 z2d(ji,jj) = z2d(ji,jj) + zftu(ji,jj,jk) 324 END DO 325 END DO 326 END DO 327 z2d(:,:) = zztmp * z2d(:,:) 328 CALL lbc_lnk( z2d, 'U', -1. ) 329 CALL iom_put( "udiff_heattr", z2d ) ! heat transport in i-direction 330 DO jk = 1, jpkm1 331 DO jj = 2, jpjm1 332 DO ji = fs_2, fs_jpim1 ! vector opt. 333 z2d(ji,jj) = z2d(ji,jj) + zftv(ji,jj,jk) 334 END DO 335 END DO 336 END DO 337 z2d(:,:) = zztmp * z2d(:,:) 338 CALL lbc_lnk( z2d, 'V', -1. ) 339 CALL iom_put( "vdiff_heattr", z2d ) ! heat transport in i-direction 340 END IF 348 341 #endif 349 !350 END DO351 !342 ! 343 END DO 344 ! 352 345 END SUBROUTINE tra_ldf_iso_grif 353 346 -
branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_lap.F90
r2399 r2450 104 104 DO ji = 1, fs_jpim1 ! vector opt. 105 105 ! last level 106 iku = MIN ( mbathy(ji,jj), mbathy(ji+1,jj ) ) - 1107 ikv = MIN ( mbathy(ji,jj), mbathy(ji ,jj+1) ) - 1106 iku = mbku(ji,jj) 107 ikv = mbkv(ji,jj) 108 108 IF( iku == jk ) THEN 109 109 zabe1 = fsahtu(ji,jj,iku) * umask(ji,jj,iku) * e1ur(ji,jj) * fse3u(ji,jj,iku) -
branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/TRA/tranpc.F90
r2287 r2450 33 33 !! NEMO/OPA 3.3 , NEMO Consortium (2010) 34 34 !! $Id$ 35 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 36 !!---------------------------------------------------------------------- 37 35 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 36 !!---------------------------------------------------------------------- 38 37 CONTAINS 39 38 … … 131 130 IF( zwy(ji,1) /= 0.e0 ) THEN 132 131 ! 133 ikbot = mb athy(ji,jj) ! ikbot: ocean bottomlevel132 ikbot = mbkt(ji,jj) ! ikbot: ocean bottom T-level 134 133 ! 135 134 DO jiter = 1, jpk ! vertical iteration … … 140 139 220 CONTINUE 141 140 ik = ik + 1 142 IF( ik >= ikbot -1) GO TO 200141 IF( ik >= ikbot ) GO TO 200 143 142 zwx(ji,ik) = zrhop(ji,jj,ik) - zrhop(ji,jj,ik+1) 144 143 IF( zwx(ji,ik) <= 0.e0 ) GO TO 220 -
branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/TRA/zpshde.F90
r2287 r2450 1 1 MODULE zpshde 2 !!====================================================================== ========2 !!====================================================================== 3 3 !! *** MODULE zpshde *** 4 !! z-coordinate - partial step : Horizontal Derivative5 !!====================================================================== ========4 !! z-coordinate + partial step : Horizontal Derivative at ocean bottom level 5 !!====================================================================== 6 6 !! History : OPA ! 2002-04 (A. Bozec) Original code 7 !! 8.5! 2002-08 (G. Madec E. Durand) Optimization and Free form8 !! NEMO 1.0! 2004-03 (C. Ethe) adapted for passive tracers7 !! NEMO 1.0 ! 2002-08 (G. Madec E. Durand) Optimization and Free form 8 !! - ! 2004-03 (C. Ethe) adapted for passive tracers 9 9 !! 3.3 ! 2010-05 (C. Ethe, G. Madec) merge TRC-TRA 10 !!====================================================================== ========10 !!====================================================================== 11 11 12 12 !!---------------------------------------------------------------------- … … 14 14 !! ocean level (Z-coord. with Partial Steps) 15 15 !!---------------------------------------------------------------------- 16 USE dom_oce ! ocean space domainvariables17 USE oce ! ocean dynamics and tracersvariables16 USE oce ! ocean: dynamics and tracers variables 17 USE dom_oce ! domain: ocean variables 18 18 USE phycst ! physical constants 19 USE eosbn2 ! ocean equation of state 19 20 USE in_out_manager ! I/O manager 20 USE eosbn2 ! ocean equation of state21 21 USE lbclnk ! lateral boundary conditions (or mpp link) 22 22 … … 24 24 PRIVATE 25 25 26 PUBLIC zps_hde ! routine called by step.F90 27 PUBLIC zps_hde_init ! routine called by opa.F90 28 29 INTEGER, DIMENSION(jpi,jpj) :: mbatu, mbatv ! bottom ocean level index at U- and V-points 26 PUBLIC zps_hde ! routine called by step.F90 30 27 31 28 !! * Substitutions … … 35 32 !! NEMO/OPA 3.3 , NEMO Consortium (2010) 36 33 !! $Id$ 37 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)34 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 38 35 !!---------------------------------------------------------------------- 39 36 CONTAINS … … 44 41 !! *** ROUTINE zps_hde *** 45 42 !! 46 !! ** Purpose : Compute the horizontal derivative of T, S and r d43 !! ** Purpose : Compute the horizontal derivative of T, S and rho 47 44 !! at u- and v-points with a linear interpolation for z-coordinate 48 45 !! with partial steps. … … 78 75 !! formulation of the equation of state (eos). 79 76 !! Gradient formulation for rho : 80 !! di(rho) = rd~ - rd(i,j,k) or rd (i+1,j,k) - rd~ 81 !! 82 !! ** Action : - pgtu, pgtv: horizontal gradient of tracer at U/V-points 83 !! - pgru, pgrv: horizontal gradient of rd if present at U/V-points 84 !! and rd at V-points 77 !! di(rho) = rd~ - rd(i,j,k) or rd(i+1,j,k) - rd~ 78 !! 79 !! ** Action : - pgtu, pgtv: horizontal gradient of tracer at u- & v-points 80 !! - pgru, pgrv: horizontal gradient of rho (if present) at u- & v-points 85 81 !!---------------------------------------------------------------------- 86 82 INTEGER , INTENT(in ) :: kt ! ocean time-step index … … 92 88 !! 93 89 INTEGER :: ji, jj, jn ! Dummy loop indices 94 INTEGER :: iku, ikv ! partial step levelat u- and v-points90 INTEGER :: iku, ikv, ikum1, ikvm1 ! partial step level (ocean bottom level) at u- and v-points 95 91 REAL(wp), DIMENSION(jpi,jpj,kjpt) :: zti, ztj ! interpolated value of tracer 96 92 REAL(wp), DIMENSION(jpi,jpj) :: zri, zrj ! interpolated value of rd … … 99 95 !!---------------------------------------------------------------------- 100 96 101 102 ! Interpolation of tracers at the last ocean level 103 DO jn = 1, kjpt 97 DO jn = 1, kjpt !== Interpolation of tracers at the last ocean level ==! 104 98 ! 105 99 # if defined key_vectopt_loop … … 110 104 DO ji = 1, jpim1 111 105 # endif 112 ! last level 113 iku = mbatu(ji,jj) 114 ikv = mbatv(ji,jj) 115 ze3wu = fse3w(ji+1,jj ,iku) - fse3w(ji,jj,iku) 116 ze3wv = fse3w(ji ,jj+1,ikv) - fse3w(ji,jj,ikv) 117 106 iku = mbku(ji,jj) ; ikum1 = MAX( iku , 1 ) ! last and before last ocean level at u- & v-points 107 ikv = mbkv(ji,jj) ; ikvm1 = MAX( ikv , 1 ) ! if level first is a p-step, ik.m1=1 108 ze3wu = fse3w(ji+1,jj ,iku) - fse3w(ji,jj,iku) 109 ze3wv = fse3w(ji ,jj+1,ikv) - fse3w(ji,jj,ikv) 110 ! 118 111 ! i- direction 119 IF( ze3wu >= 0. ) THEN ! case 1112 IF( ze3wu >= 0._wp ) THEN ! case 1 120 113 zmaxu = ze3wu / fse3w(ji+1,jj,iku) 121 114 ! interpolated values of tracers 122 zti(ji,jj,jn) = pta(ji+1,jj,iku,jn) + zmaxu * ( pta(ji+1,jj,iku -1,jn) - pta(ji+1,jj,iku,jn) )115 zti(ji,jj,jn) = pta(ji+1,jj,iku,jn) + zmaxu * ( pta(ji+1,jj,ikum1,jn) - pta(ji+1,jj,iku,jn) ) 123 116 ! gradient of tracers 124 117 pgtu(ji,jj,jn) = umask(ji,jj,1) * ( zti(ji,jj,jn) - pta(ji,jj,iku,jn) ) 125 ELSE ! case 2118 ELSE ! case 2 126 119 zmaxu = -ze3wu / fse3w(ji,jj,iku) 127 120 ! interpolated values of tracers 128 zti(ji,jj,jn) = pta(ji,jj,iku,jn) + zmaxu * ( pta(ji,jj,iku -1,jn) - pta(ji,jj,iku,jn) )121 zti(ji,jj,jn) = pta(ji,jj,iku,jn) + zmaxu * ( pta(ji,jj,ikum1,jn) - pta(ji,jj,iku,jn) ) 129 122 ! gradient of tracers 130 123 pgtu(ji,jj,jn) = umask(ji,jj,1) * ( pta(ji+1,jj,iku,jn) - zti(ji,jj,jn) ) 131 124 ENDIF 132 125 ! 133 126 ! j- direction 134 IF( ze3wv >= 0. ) THEN ! case 1127 IF( ze3wv >= 0._wp ) THEN ! case 1 135 128 zmaxv = ze3wv / fse3w(ji,jj+1,ikv) 136 129 ! interpolated values of tracers 137 ztj(ji,jj,jn) = pta(ji,jj+1,ikv,jn) + zmaxv * ( pta(ji,jj+1,ikv -1,jn) - pta(ji,jj+1,ikv,jn) )130 ztj(ji,jj,jn) = pta(ji,jj+1,ikv,jn) + zmaxv * ( pta(ji,jj+1,ikvm1,jn) - pta(ji,jj+1,ikv,jn) ) 138 131 ! gradient of tracers 139 132 pgtv(ji,jj,jn) = vmask(ji,jj,1) * ( ztj(ji,jj,jn) - pta(ji,jj,ikv,jn) ) 140 ELSE ! case 2133 ELSE ! case 2 141 134 zmaxv = -ze3wv / fse3w(ji,jj,ikv) 142 135 ! interpolated values of tracers 143 ztj(ji,jj,jn) = pta(ji,jj,ikv,jn) + zmaxv * ( pta(ji,jj,ikv -1,jn) - pta(ji,jj,ikv,jn) )136 ztj(ji,jj,jn) = pta(ji,jj,ikv,jn) + zmaxv * ( pta(ji,jj,ikvm1,jn) - pta(ji,jj,ikv,jn) ) 144 137 ! gradient of tracers 145 138 pgtv(ji,jj,jn) = vmask(ji,jj,1) * ( pta(ji,jj+1,ikv,jn) - ztj(ji,jj,jn) ) … … 162 155 DO ji = 1, jpim1 163 156 # endif 164 iku = mb atu(ji,jj)165 ikv = mb atv(ji,jj)157 iku = mbku(ji,jj) 158 ikv = mbkv(ji,jj) 166 159 ze3wu = fse3w(ji+1,jj ,iku) - fse3w(ji,jj,iku) 167 160 ze3wv = fse3w(ji ,jj+1,ikv) - fse3w(ji,jj,ikv) 168 IF( ze3wu >= 0. ) THEN 169 zhi(ji,jj) = fsdept(ji ,jj,iku) 170 ELSE 171 zhi(ji,jj) = fsdept(ji+1,jj,iku) 172 ENDIF 173 IF( ze3wv >= 0. ) THEN 174 zhj(ji,jj) = fsdept(ji,jj ,ikv) 175 ELSE 176 zhj(ji,jj) = fsdept(ji,jj+1,ikv) 161 IF( ze3wu >= 0._wp ) THEN ; zhi(ji,jj) = fsdept(ji ,jj,iku) ! i-direction: case 1 162 ELSE ; zhi(ji,jj) = fsdept(ji+1,jj,iku) ! - - case 2 163 ENDIF 164 IF( ze3wv >= 0._wp ) THEN ; zhj(ji,jj) = fsdept(ji,jj ,ikv) ! j-direction: case 1 165 ELSE ; zhj(ji,jj) = fsdept(ji,jj+1,ikv) ! - - case 2 177 166 ENDIF 178 167 # if ! defined key_vectopt_loop … … 193 182 DO ji = 1, jpim1 194 183 # endif 195 iku = mb atu(ji,jj)196 ikv = mb atv(ji,jj)184 iku = mbku(ji,jj) 185 ikv = mbkv(ji,jj) 197 186 ze3wu = fse3w(ji+1,jj ,iku) - fse3w(ji,jj,iku) 198 187 ze3wv = fse3w(ji ,jj+1,ikv) - fse3w(ji,jj,ikv) 199 IF( ze3wu >= 0. ) THEN ! i-direction: case 1 200 pgru(ji,jj) = umask(ji,jj,1) * ( zri(ji,jj) - prd(ji,jj,iku) ) 201 ELSE ! i-direction: case 2 202 pgru(ji,jj) = umask(ji,jj,1) * ( prd(ji+1,jj,iku) - zri(ji,jj) ) 203 ENDIF 204 IF( ze3wv >= 0. ) THEN ! j-direction: case 1 205 pgrv(ji,jj) = vmask(ji,jj,1) * ( zrj(ji,jj) - prd(ji,jj,ikv) ) 206 ELSE ! j-direction: case 2 207 pgrv(ji,jj) = vmask(ji,jj,1) * ( prd(ji,jj+1,ikv) - zrj(ji,jj) ) 188 IF( ze3wu >= 0._wp ) THEN ; pgru(ji,jj) = umask(ji,jj,1) * ( zri(ji ,jj) - prd(ji,jj,iku) ) ! i: 1 189 ELSE ; pgru(ji,jj) = umask(ji,jj,1) * ( prd(ji+1,jj,iku) - zri(ji,jj) ) ! i: 2 190 ENDIF 191 IF( ze3wv >= 0._wp ) THEN ; pgrv(ji,jj) = vmask(ji,jj,1) * ( zrj(ji,jj ) - prd(ji,jj,ikv) ) ! j: 1 192 ELSE ; pgrv(ji,jj) = vmask(ji,jj,1) * ( prd(ji,jj+1,ikv) - zrj(ji,jj) ) ! j: 2 208 193 ENDIF 209 194 # if ! defined key_vectopt_loop … … 217 202 END SUBROUTINE zps_hde 218 203 219 220 SUBROUTINE zps_hde_init221 !!----------------------------------------------------------------------222 !! *** ROUTINE zps_hde_init ***223 !!224 !! ** Purpose : Computation of bottom ocean level index at U- and V-points225 !!226 !!----------------------------------------------------------------------227 INTEGER :: ji, jj ! Dummy loop indices228 REAL(wp), DIMENSION(jpi,jpj) :: zti, ztj ! 2D workspace229 !!----------------------------------------------------------------------230 !231 mbatu(:,:) = 0232 mbatv(:,:) = 0233 DO jj = 1, jpjm1234 DO ji = 1, fs_jpim1 ! vector opt.235 mbatu(ji,jj) = MAX( MIN( mbathy(ji,jj), mbathy(ji+1,jj ) ) - 1, 2 )236 mbatv(ji,jj) = MAX( MIN( mbathy(ji,jj), mbathy(ji ,jj+1) ) - 1, 2 )237 END DO238 END DO239 zti(:,:) = FLOAT( mbatu(:,:) )240 ztj(:,:) = FLOAT( mbatv(:,:) )241 ! lateral boundary conditions: T-point, sign unchanged242 CALL lbc_lnk( zti , 'U', 1. ) ; CALL lbc_lnk( ztj , 'V', 1. )243 mbatu(:,:) = MAX( INT( zti(:,:) ), 2 )244 mbatv(:,:) = MAX( INT( ztj(:,:) ), 2 )245 !246 END SUBROUTINE zps_hde_init247 204 !!====================================================================== 248 205 END MODULE zpshde -
branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/TRD/trdvor.F90
r2364 r2450 10 10 !!---------------------------------------------------------------------- 11 11 !! 'key_trdvor' : momentum trend diagnostics 12 !!----------------------------------------------------------------------13 12 !!---------------------------------------------------------------------- 14 13 !! trd_vor : momentum trends averaged over the depth … … 39 38 PUBLIC trd_vor_init ! routine called by opa.F90 40 39 41 INTEGER :: & 42 nh_t, nmoydpvor , & 43 nidvor, nhoridvor, & 44 ndexvor1(jpi*jpj), & 45 ndimvor1, icount, & 46 idebug ! (0/1) set it to 1 in case of problem to have more print 47 48 REAL(wp), DIMENSION(jpi,jpj) :: & 49 vor_avr , & ! average 50 vor_avrb , & ! before vorticity (kt-1) 51 vor_avrbb , & ! vorticity at begining of the nwrite-1 timestep averaging period 52 vor_avrbn , & ! after vorticity at time step after the 53 rotot , & ! begining of the NWRITE-1 timesteps 54 vor_avrtot , & 55 vor_avrres 56 57 REAL(wp), DIMENSION(jpi,jpj,jpltot_vor):: vortrd !: curl of trends 40 INTEGER :: nh_t, nmoydpvor, nidvor, nhoridvor, ndexvor1(jpi*jpj), ndimvor1, icount ! needs for IOIPSL output 41 INTEGER :: ndebug ! (0/1) set it to 1 in case of problem to have more print 42 43 REAL(wp), DIMENSION(jpi,jpj) :: vor_avr ! average 44 REAL(wp), DIMENSION(jpi,jpj) :: vor_avrb ! before vorticity (kt-1) 45 REAL(wp), DIMENSION(jpi,jpj) :: vor_avrbb ! vorticity at begining of the nwrite-1 timestep averaging period 46 REAL(wp), DIMENSION(jpi,jpj) :: vor_avrbn ! after vorticity at time step after the 47 REAL(wp), DIMENSION(jpi,jpj) :: rotot ! begining of the NWRITE-1 timesteps 48 REAL(wp), DIMENSION(jpi,jpj) :: vor_avrtot ! 49 REAL(wp), DIMENSION(jpi,jpj) :: vor_avrres ! 50 51 REAL(wp), DIMENSION(jpi,jpj,jpltot_vor) :: vortrd ! curl of trends 58 52 59 53 CHARACTER(len=12) :: cvort … … 66 60 !! NEMO/OPA 3.3 , NEMO Consortium (2010) 67 61 !! $Id$ 68 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 69 !!---------------------------------------------------------------------- 70 62 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 63 !!---------------------------------------------------------------------- 71 64 CONTAINS 72 65 … … 76 69 !! 77 70 !! ** Purpose : computation of vertically integrated vorticity budgets 78 !! from ocean surface down to control surface (NetCDF output) 79 !! 80 !! ** Method/usage : 81 !! integration done over nwrite-1 time steps 82 !! 83 !! 84 !! ** Action : 85 !! /comvor/ : 86 !! vor_avr average 87 !! vor_avrb vorticity at kt-1 88 !! vor_avrbb vorticity at begining of the NWRITE-1 89 !! time steps averaging period 90 !! vor_avrbn vorticity at time step after the 91 !! begining of the NWRITE-1 time 92 !! steps averaging period 93 !! 94 !! trends : 95 !! 71 !! from ocean surface down to control surface (NetCDF output) 72 !! 73 !! ** Method/usage : integration done over nwrite-1 time steps 74 !! 75 !! ** Action : trends : 96 76 !! vortrd (,, 1) = Pressure Gradient Trend 97 77 !! vortrd (,, 2) = KE Gradient Trend … … 104 84 !! vortrd (,, 9) = Beta V 105 85 !! vortrd (,,10) = forcing term 106 !! vortrd (,,11) = bottom friction term86 !! vortrd (,,11) = bottom friction term 107 87 !! rotot(,) : total cumulative trends over nwrite-1 time steps 108 88 !! vor_avrtot(,) : first membre of vrticity equation … … 110 90 !! 111 91 !! trends output in netCDF format using ioipsl 112 !! 113 !!---------------------------------------------------------------------- 114 INTEGER, INTENT( in ) :: ktrd ! ocean trend index 115 REAL(wp), DIMENSION(jpi,jpj), INTENT( inout ) :: & 116 putrdvor, & ! u vorticity trend 117 pvtrdvor ! v vorticity trend 118 !! 119 INTEGER :: ji, jj 120 INTEGER :: ikbu, ikbum1, ikbv, ikbvm1 92 !!---------------------------------------------------------------------- 93 INTEGER , INTENT(in ) :: ktrd ! ocean trend index 94 REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: putrdvor ! u vorticity trend 95 REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pvtrdvor ! v vorticity trend 96 !! 97 INTEGER :: ji, jj ! dummy loop indices 98 INTEGER :: ikbu, ikbv ! local integers 121 99 REAL(wp), DIMENSION(jpi,jpj) :: zudpvor, zvdpvor ! total cmulative trends 122 100 !!---------------------------------------------------------------------- 123 101 124 102 ! Initialization 125 zudpvor(:,:) = 0. e0126 zvdpvor(:,:) = 0. e0127 128 CALL lbc_lnk( putrdvor, 'U' , -1. ) 103 zudpvor(:,:) = 0._wp 104 zvdpvor(:,:) = 0._wp 105 ! 106 CALL lbc_lnk( putrdvor, 'U' , -1. ) ! lateral boundary condition on input momentum trends 129 107 CALL lbc_lnk( pvtrdvor, 'V' , -1. ) 130 108 … … 134 112 135 113 SELECT CASE (ktrd) 136 114 ! 137 115 CASE (jpvor_bfr) ! bottom friction 138 139 116 DO jj = 2, jpjm1 140 117 DO ji = fs_2, fs_jpim1 141 ikbu = min( mbathy(ji+1,jj), mbathy(ji,jj) ) 142 ikbum1 = max( ikbu-1, 1 ) 143 ikbv = min( mbathy(ji,jj+1), mbathy(ji,jj) ) 144 ikbvm1 = max( ikbv-1, 1 ) 145 146 zudpvor(ji,jj) = putrdvor(ji,jj) * fse3u(ji,jj,ikbum1) * e1u(ji,jj) * umask(ji,jj,ikbum1) 147 zvdpvor(ji,jj) = pvtrdvor(ji,jj) * fse3v(ji,jj,ikbvm1) * e2v(ji,jj) * vmask(ji,jj,ikbvm1) 118 ikbu = mbkv(ji,jj) 119 ikbv = mbkv(ji,jj) 120 zudpvor(ji,jj) = putrdvor(ji,jj) * fse3u(ji,jj,ikbu) * e1u(ji,jj) * umask(ji,jj,ikbu) 121 zvdpvor(ji,jj) = pvtrdvor(ji,jj) * fse3v(ji,jj,ikbv) * e2v(ji,jj) * vmask(ji,jj,ikbv) 148 122 END DO 149 123 END DO 150 124 ! 151 125 CASE (jpvor_swf) ! wind stress 152 153 126 zudpvor(:,:) = putrdvor(:,:) * fse3u(:,:,1) * e1u(:,:) * umask(:,:,1) 154 127 zvdpvor(:,:) = pvtrdvor(:,:) * fse3v(:,:,1) * e2v(:,:) * vmask(:,:,1) 155 128 ! 156 129 END SELECT 157 130 … … 163 136 DO ji=1,jpim1 164 137 DO jj=1,jpjm1 165 vortrd(ji,jj,ktrd) = ( zvdpvor(ji+1,jj) - zvdpvor(ji,jj) & 166 & - ( zudpvor(ji,jj+1) - zudpvor(ji,jj) ) ) & 167 & / ( e1f(ji,jj) * e2f(ji,jj) ) 138 vortrd(ji,jj,ktrd) = ( zvdpvor(ji+1,jj) - zvdpvor(ji,jj) & 139 & - ( zudpvor(ji,jj+1) - zudpvor(ji,jj) ) ) / ( e1f(ji,jj) * e2f(ji,jj) ) 168 140 END DO 169 141 END DO 170 171 ! Surface mask 172 vortrd(:,:,ktrd) = vortrd(:,:,ktrd) * fmask(:,:,1) 173 174 IF( idebug /= 0 ) THEN 142 vortrd(:,:,ktrd) = vortrd(:,:,ktrd) * fmask(:,:,1) ! Surface mask 143 144 IF( ndebug /= 0 ) THEN 175 145 IF(lwp) WRITE(numout,*) ' debuging trd_vor_zint: I done' 176 146 CALL FLUSH(numout) … … 185 155 !! 186 156 !! ** Purpose : computation of vertically integrated vorticity budgets 187 !! from ocean surface down to control surface (NetCDF output) 188 !! 189 !! ** Method/usage : 190 !! integration done over nwrite-1 time steps 191 !! 192 !! 193 !! ** Action : 194 !! /comvor/ : 195 !! vor_avr average 196 !! vor_avrb vorticity at kt-1 197 !! vor_avrbb vorticity at begining of the NWRITE-1 198 !! time steps averaging period 199 !! vor_avrbn vorticity at time step after the 200 !! begining of the NWRITE-1 time 201 !! steps averaging period 202 !! 203 !! trends : 204 !! 157 !! from ocean surface down to control surface (NetCDF output) 158 !! 159 !! ** Method/usage : integration done over nwrite-1 time steps 160 !! 161 !! ** Action : trends : 205 162 !! vortrd (,,1) = Pressure Gradient Trend 206 163 !! vortrd (,,2) = KE Gradient Trend … … 219 176 !! 220 177 !! trends output in netCDF format using ioipsl 221 !! 222 !!---------------------------------------------------------------------- 223 INTEGER, INTENT( in ) :: ktrd ! ocean trend index 224 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( inout ) :: & 225 putrdvor, & ! u vorticity trend 226 pvtrdvor ! v vorticity trend 178 !!---------------------------------------------------------------------- 179 INTEGER , INTENT(in ) :: ktrd ! ocean trend index 180 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: putrdvor ! u vorticity trend 181 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pvtrdvor ! v vorticity trend 227 182 !! 228 183 INTEGER :: ji, jj, jk 229 REAL(wp), DIMENSION(jpi,jpj) :: & 230 zubet, & ! u Beta.V case 231 zvbet, & ! v Beta.V case 232 zudpvor, & ! total cmulative trends 233 zvdpvor ! " " " 184 REAL(wp), DIMENSION(jpi,jpj) :: zubet , zvbet ! Beta.V 185 REAL(wp), DIMENSION(jpi,jpj) :: zudpvor, zvdpvor ! total cmulative trends 234 186 !!---------------------------------------------------------------------- 235 187 236 188 ! Initialization 237 zubet(:,:) = 0.e0 238 zvbet(:,:) = 0.e0 239 zudpvor(:,:) = 0.e0 240 zvdpvor(:,:) = 0.e0 189 zubet (:,:) = 0._wp 190 zvbet (:,:) = 0._wp 191 zudpvor(:,:) = 0._wp 192 zvdpvor(:,:) = 0._wp 193 ! 194 CALL lbc_lnk( putrdvor, 'U' , -1. ) ! lateral boundary condition on input momentum trends 195 CALL lbc_lnk( pvtrdvor, 'V' , -1. ) 241 196 242 197 ! ===================================== 243 198 ! I vertical integration of 3D trends 244 199 ! ===================================== 245 246 CALL lbc_lnk( putrdvor, 'U' , -1. )247 CALL lbc_lnk( pvtrdvor, 'V' , -1. )248 249 200 ! putrdvor and pvtrdvor terms 250 201 DO jk = 1,jpk … … 267 218 DO ji=1,jpim1 268 219 DO jj=1,jpjm1 269 vortrd(ji,jj,ktrd) = ( zvdpvor(ji+1,jj) - zvdpvor(ji,jj) - & 270 & ( zudpvor(ji,jj+1) - zudpvor(ji,jj) ) ) & 271 & / ( e1f(ji,jj) * e2f(ji,jj) ) 220 vortrd(ji,jj,ktrd) = ( zvdpvor(ji+1,jj) - zvdpvor(ji,jj) & 221 & - ( zudpvor(ji,jj+1) - zudpvor(ji,jj) ) ) / ( e1f(ji,jj) * e2f(ji,jj) ) 272 222 END DO 273 223 END DO … … 281 231 DO ji=1,jpim1 282 232 DO jj=1,jpjm1 283 vortrd(ji,jj,jpvor_bev) = ( zvbet(ji+1,jj) - zvbet(ji,jj) - & 284 & ( zubet(ji,jj+1) - zubet(ji,jj) ) ) & 285 & / ( e1f(ji,jj) * e2f(ji,jj) ) 233 vortrd(ji,jj,jpvor_bev) = ( zvbet(ji+1,jj) - zvbet(ji,jj) & 234 & - ( zubet(ji,jj+1) - zubet(ji,jj) ) ) / ( e1f(ji,jj) * e2f(ji,jj) ) 286 235 END DO 287 236 END DO … … 294 243 ENDIF 295 244 296 IF( idebug /= 0 ) THEN245 IF( ndebug /= 0 ) THEN 297 246 IF(lwp) WRITE(numout,*) ' debuging trd_vor_zint: I done' 298 247 CALL FLUSH(numout) … … 309 258 !! and make outputs (NetCDF or DIMG format) 310 259 !!---------------------------------------------------------------------- 311 INTEGER, INTENT( in ) :: kt ! ocean time-step index 312 !! 313 INTEGER :: ji, jj, jk, jl, it, itmod 314 REAL(wp) :: zmean 315 REAL(wp), DIMENSION(jpi,jpj) :: zun, zvn 260 INTEGER, INTENT(in) :: kt ! ocean time-step index 261 !! 262 INTEGER :: ji, jj, jk, jl ! dummy loop indices 263 INTEGER :: it, itmod ! local integers 264 REAL(wp) :: zmean ! local scalars 265 REAL(wp), DIMENSION(jpi,jpj) :: zun, zvn ! 2D workspace 316 266 !!---------------------------------------------------------------------- 317 267 … … 324 274 ! --------------------------------------------------- 325 275 326 IF( kt > nit000 ) THEN 327 vor_avrb(:,:) = vor_avr(:,:) 328 ENDIF 329 330 IF( idebug /= 0 ) THEN 276 IF( kt > nit000 ) vor_avrb(:,:) = vor_avr(:,:) 277 278 IF( ndebug /= 0 ) THEN 331 279 WRITE(numout,*) ' debuging trd_vor: I.1 done ' 332 280 CALL FLUSH(numout) … … 336 284 ! ---------------------------------- 337 285 338 vor_avr (:,:) = 0.339 zun (:,:)=0340 zvn (:,:)=0341 vor_avrtot(:,:) =0342 vor_avrres(:,:) =0286 vor_avr (:,:) = 0._wp 287 zun (:,:) = 0._wp 288 zvn (:,:) = 0._wp 289 vor_avrtot(:,:) = 0._wp 290 vor_avrres(:,:) = 0._wp 343 291 344 292 ! Vertically averaged velocity 345 293 DO jk = 1, jpk - 1 346 zun(:,:) =zun(:,:) + e1u(:,:)*un(:,:,jk)*fse3u(:,:,jk)347 zvn(:,:) =zvn(:,:) + e2v(:,:)*vn(:,:,jk)*fse3v(:,:,jk)294 zun(:,:) = zun(:,:) + e1u(:,:) * un(:,:,jk) * fse3u(:,:,jk) 295 zvn(:,:) = zvn(:,:) + e2v(:,:) * vn(:,:,jk) * fse3v(:,:,jk) 348 296 END DO 349 297 350 zun(:,:) =zun(:,:)*hur(:,:)351 zvn(:,:) =zvn(:,:)*hvr(:,:)298 zun(:,:) = zun(:,:) * hur(:,:) 299 zvn(:,:) = zvn(:,:) * hvr(:,:) 352 300 353 301 ! Curl 354 302 DO ji=1,jpim1 355 303 DO jj=1,jpjm1 356 vor_avr(ji,jj) = ((zvn(ji+1,jj)-zvn(ji,jj))- & 357 (zun(ji,jj+1)-zun(ji,jj))) & 358 /( e1f(ji,jj) * e2f(ji,jj) ) 359 vor_avr(ji,jj) = vor_avr(ji,jj)*fmask(ji,jj,1) 304 vor_avr(ji,jj) = ( ( zvn(ji+1,jj) - zvn(ji,jj) ) & 305 & - ( zun(ji,jj+1) - zun(ji,jj) ) ) / ( e1f(ji,jj) * e2f(ji,jj) ) * fmask(ji,jj,1) 360 306 END DO 361 307 END DO 362 308 363 IF( idebug /= 0) THEN309 IF( ndebug /= 0 ) THEN 364 310 WRITE(numout,*) ' debuging trd_vor: I.2 done' 365 311 CALL FLUSH(numout) … … 377 323 ENDIF 378 324 379 IF( idebug /= 0 ) THEN325 IF( ndebug /= 0 ) THEN 380 326 WRITE(numout,*) ' debuging trd_vor: I1.1 done' 381 327 CALL FLUSH(numout) … … 395 341 ENDIF 396 342 397 IF( idebug /= 0 ) THEN343 IF( ndebug /= 0 ) THEN 398 344 WRITE(numout,*) ' debuging trd_vor: II.2 done' 399 345 CALL FLUSH(numout) … … 405 351 406 352 ! define time axis 407 it = kt353 it = kt 408 354 itmod = kt - nit000 + 1 409 355 … … 412 358 ! III.1 compute total trend 413 359 ! ------------------------ 414 zmean = float(nmoydpvor) 415 416 vor_avrtot(:,:) = ( vor_avr(:,:) - vor_avrbn(:,:) + vor_avrb(:,:) - & 417 vor_avrbb(:,:) ) / (zmean * 2. * rdt) 418 419 IF( idebug /= 0 ) THEN 360 zmean = 1._wp / ( REAL( nmoydpvor, wp ) * 2._wp * rdt ) 361 vor_avrtot(:,:) = ( vor_avr(:,:) - vor_avrbn(:,:) + vor_avrb(:,:) - vor_avrbb(:,:) ) * zmean 362 363 IF( ndebug /= 0 ) THEN 420 364 WRITE(numout,*) ' zmean = ',zmean 421 365 WRITE(numout,*) ' debuging trd_vor: III.1 done' … … 425 369 ! III.2 compute residual 426 370 ! --------------------- 371 zmean = 1._wp / REAL( nmoydpvor, wp ) 427 372 vor_avrres(:,:) = vor_avrtot(:,:) - rotot(:,:) / zmean 428 373 … … 431 376 CALL lbc_lnk( vor_avrres, 'F', 1. ) 432 377 433 IF( idebug /= 0 ) THEN378 IF( ndebug /= 0 ) THEN 434 379 WRITE(numout,*) ' debuging trd_vor: III.2 done' 435 380 CALL FLUSH(numout) … … 439 384 ! ------------------------------ 440 385 vor_avrbb(:,:) = vor_avrb(:,:) 441 vor_avrbn(:,:) = vor_avr (:,:)442 443 IF( idebug /= 0 ) THEN386 vor_avrbn(:,:) = vor_avr (:,:) 387 388 IF( ndebug /= 0 ) THEN 444 389 WRITE(numout,*) ' debuging trd_vor: III.3 done' 445 390 CALL FLUSH(numout) 446 391 ENDIF 447 448 nmoydpvor =0449 392 ! 393 nmoydpvor = 0 394 ! 450 395 ENDIF 451 396 … … 475 420 CALL histwrite( nidvor,"sovorgap",it,vor_avrres ,ndimvor1,ndexvor1) ! gap between 1st and 2 nd mbre 476 421 ! 477 IF( idebug /= 0 ) THEN422 IF( ndebug /= 0 ) THEN 478 423 WRITE(numout,*) ' debuging trd_vor: III.4 done' 479 424 CALL FLUSH(numout) … … 508 453 509 454 ! Open specifier 510 idebug = 0 ! set it to 1 in case of problem to have more Print455 ndebug = 0 ! set it to 1 in case of problem to have more Print 511 456 512 457 IF(lwp) THEN … … 528 473 vor_avrres(:,:)=0 529 474 530 IF( idebug /= 0 ) THEN475 IF( ndebug /= 0 ) THEN 531 476 WRITE(numout,*) ' debuging trd_vor_init: I. done' 532 477 CALL FLUSH(numout) … … 600 545 CALL histend( nidvor, snc4set ) 601 546 602 IF( idebug /= 0 ) THEN547 IF( ndebug /= 0 ) THEN 603 548 WRITE(numout,*) ' debuging trd_vor_init: II. done' 604 549 CALL FLUSH(numout) … … 621 566 REAL, DIMENSION(:,:), INTENT( inout ) :: putrdvor, pvtrdvor 622 567 INTEGER, INTENT( in ) :: ktrd ! ocean trend index 623 WRITE(*,*) 'trd_vor_zint_2d: You should not have seen this print! error?', putrdvor(1,1) 624 WRITE(*,*) ' " " : You should not have seen this print! error?', pvtrdvor(1,1) 625 WRITE(*,*) ' " " : You should not have seen this print! error?', ktrd 568 WRITE(*,*) 'trd_vor_zint_2d: You should not have seen this print! error?', putrdvor(1,1), pvtrdvor(1,1), ktrd 626 569 END SUBROUTINE trd_vor_zint_2d 627 570 SUBROUTINE trd_vor_zint_3d( putrdvor, pvtrdvor, ktrd ) 628 571 REAL, DIMENSION(:,:,:), INTENT( inout ) :: putrdvor, pvtrdvor 629 572 INTEGER, INTENT( in ) :: ktrd ! ocean trend index 630 WRITE(*,*) 'trd_vor_zint_3d: You should not have seen this print! error?', putrdvor(1,1,1) 631 WRITE(*,*) ' " " : You should not have seen this print! error?', pvtrdvor(1,1,1) 632 WRITE(*,*) ' " " : You should not have seen this print! error?', ktrd 573 WRITE(*,*) 'trd_vor_zint_3d: You should not have seen this print! error?', putrdvor(1,1,1), pvtrdvor(1,1,1), ktrd 633 574 END SUBROUTINE trd_vor_zint_3d 634 575 SUBROUTINE trd_vor_init ! Empty routine -
branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfbfr.F90
r2346 r2450 40 40 REAL(wp) :: rn_bfri2 = 1.0e-3_wp ! bottom drag coefficient (non linear case) 41 41 REAL(wp) :: rn_bfeb2 = 2.5e-3_wp ! background bottom turbulent kinetic energy [m2/s2] 42 REAL(wp) :: rn_bfrien = 30 _wp! local factor to enhance coefficient bfri42 REAL(wp) :: rn_bfrien = 30._wp ! local factor to enhance coefficient bfri 43 43 LOGICAL :: ln_bfr2d = .false. ! logical switch for 2D enhancement 44 44 … … 74 74 INTEGER, INTENT( in ) :: kt ! ocean time-step index 75 75 !! 76 INTEGER :: ji, jj ! dummy loop indices 77 INTEGER :: ikbu, ikbum1 ! temporary integers 78 INTEGER :: ikbv, ikbvm1 ! - - 76 INTEGER :: ji, jj ! dummy loop indices 77 INTEGER :: ikbu, ikbv ! local integers 79 78 REAL(wp) :: zvu, zuv, zecu, zecv ! temporary scalars 80 79 !!---------------------------------------------------------------------- … … 96 95 DO ji = 2, jpim1 97 96 # endif 98 ikbu = MIN( mbathy(ji+1,jj ), mbathy(ji,jj) ) 99 ikbv = MIN( mbathy(ji ,jj+1), mbathy(ji,jj) ) 100 ikbum1 = MAX( ikbu-1, 1 ) 101 ikbvm1 = MAX( ikbv-1, 1 ) 97 ikbu = mbku(ji,jj) ! ocean bottom level at u- and v-points 98 ikbv = mbkv(ji,jj) ! (deepest ocean u- and v-points) 102 99 ! 103 zvu = 0.25 * ( vn(ji,jj ,ikbu m1) + vn(ji+1,jj ,ikbum1) &104 & + vn(ji,jj-1,ikbu m1) + vn(ji+1,jj-1,ikbum1) )105 zuv = 0.25 * ( un(ji,jj ,ikbv m1) + un(ji-1,jj ,ikbvm1) &106 & + un(ji,jj+1,ikbv m1) + un(ji-1,jj+1,ikbvm1) )100 zvu = 0.25 * ( vn(ji,jj ,ikbu) + vn(ji+1,jj ,ikbu) & 101 & + vn(ji,jj-1,ikbu) + vn(ji+1,jj-1,ikbu) ) 102 zuv = 0.25 * ( un(ji,jj ,ikbv) + un(ji-1,jj ,ikbv) & 103 & + un(ji,jj+1,ikbv) + un(ji-1,jj+1,ikbv) ) 107 104 ! 108 zecu = SQRT( un(ji,jj,ikbu m1) * un(ji,jj,ikbum1) + zvu*zvu + rn_bfeb2 )109 zecv = SQRT( vn(ji,jj,ikbv m1) * vn(ji,jj,ikbvm1) + zuv*zuv + rn_bfeb2 )105 zecu = SQRT( un(ji,jj,ikbu) * un(ji,jj,ikbu) + zvu*zvu + rn_bfeb2 ) 106 zecv = SQRT( vn(ji,jj,ikbv) * vn(ji,jj,ikbv) + zuv*zuv + rn_bfeb2 ) 110 107 ! 111 bfrua(ji,jj) = - 0.5 * ( bfrcoef2d(ji,jj) + bfrcoef2d(ji+1,jj ) ) * zecu112 bfrva(ji,jj) = - 0.5 * ( bfrcoef2d(ji,jj) + bfrcoef2d(ji ,jj+1) ) * zecv108 bfrua(ji,jj) = - 0.5_wp * ( bfrcoef2d(ji,jj) + bfrcoef2d(ji+1,jj ) ) * zecu 109 bfrva(ji,jj) = - 0.5_wp * ( bfrcoef2d(ji,jj) + bfrcoef2d(ji ,jj+1) ) * zecv 113 110 END DO 114 111 END DO … … 126 123 DO jj = 2, jpjm1 127 124 DO ji = fs_2, fs_jpim1 ! vector opt. 128 ! 129 ikbu = MIN( mbathy(ji+1,jj ), mbathy(ji,jj) ) 130 ikbum1 = MAX( ikbu-1, 1 ) 131 wbotu(ji,jj) = bfrua(ji,jj) * un(ji,jj,ikbum1) * umask(ji,jj,ikbum1) 132 ! 133 ikbv = MIN( mbathy(ji,jj+1), mbathy(ji,jj) ) 134 ikbvm1 = MAX( ikbv-1, 1 ) 135 wbotv(ji,jj) = bfrva(ji,jj) * vn(ji,jj,ikbvm1) * vmask(ji,jj,ikbvm1) 136 ! 125 wbotu(ji,jj) = bfrua(ji,jj) * un(ji,jj,mbku(ji,jj)) * umask(ji,jj,mbku(ji,jj)) 126 wbotv(ji,jj) = bfrva(ji,jj) * vn(ji,jj,mbkv(ji,jj)) * vmask(ji,jj,mbkv(ji,jj)) 137 127 END DO 138 128 END DO 139 CALL lbc_lnk( wbotu(:,:), 'U', -1. ) 140 CALL lbc_lnk( wbotv(:,:), 'V', -1. ) 129 CALL lbc_lnk( wbotu(:,:), 'U', -1. ) ; CALL lbc_lnk( wbotv(:,:), 'V', -1. ) ! lateral boundary condition 141 130 #endif 142 131 … … 242 231 DO ji = 2, jpim1 243 232 # endif 244 ikbu = MIN( mbathy(ji+1,jj ), mbathy(ji,jj) )245 ikbv = MIN( mbathy(ji ,jj+1), mbathy(ji,jj))233 ikbu = mbku(ji,jj) ! deepest ocean level at u- and v-points 234 ikbv = mbkv(ji,jj) 246 235 zfru = 0.5 * fse3u(ji,jj,ikbu) / rdt 247 236 zfrv = 0.5 * fse3v(ji,jj,ikbv) / rdt 248 237 IF( ABS( bfrcoef2d(ji,jj) ) > zfru ) THEN 249 238 IF( ln_ctl ) THEN 250 WRITE(numout,*) 'BFR ', narea,nimpp+ji,njmpp+jj,ikbu251 WRITE(numout,*) 'BFR ', ABS( bfrcoef2d(ji,jj) ), zfru239 WRITE(numout,*) 'BFR ', narea, nimpp+ji, njmpp+jj, ikbu 240 WRITE(numout,*) 'BFR ', ABS( bfrcoef2d(ji,jj) ), zfru 252 241 ENDIF 253 242 ictu = ictu + 1 … … 270 259 CALL mpp_max( zmaxbfr ) 271 260 ENDIF 272 IF( lwp .AND. ictu + ictv .GT.0 ) THEN261 IF( lwp .AND. ictu + ictv > 0 ) THEN 273 262 WRITE(numout,*) ' Bottom friction stability check failed at ', ictu, ' U-points ' 274 263 WRITE(numout,*) ' Bottom friction stability check failed at ', ictv, ' V-points ' -
branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfgls.F90
r2409 r2450 196 196 DO jj = 2, jpjm1 197 197 DO ji = fs_2, fs_jpim1 ! vector opt. 198 zup = mxln(ji,jj,jk) * fsdepw(ji,jj,mb athy(ji,jj))199 zdown = vkarmn * fsdepw(ji,jj,jk) * ( -fsdepw(ji,jj,jk) + fsdepw(ji,jj,mb athy(ji,jj)) )198 zup = mxln(ji,jj,jk) * fsdepw(ji,jj,mbkt(ji,jj)+1) 199 zdown = vkarmn * fsdepw(ji,jj,jk) * ( -fsdepw(ji,jj,jk) + fsdepw(ji,jj,mbkt(ji,jj)+1) ) 200 200 zcoef = ( zup / MAX( zdown, rsmall ) ) 201 201 zwall (ji,jj,jk) = ( 1._wp + re2 * zcoef*zcoef ) * tmask(ji,jj,jk) … … 360 360 !CDIR NOVERRCHK 361 361 DO ji = fs_2, fs_jpim1 ! vector opt. 362 ibot = mb athy(ji,jj)363 ibotm1 = ibot-1362 ibot = mbkt(ji,jj) + 1 ! k bottom level of w-point 363 ibotm1 = mbkt(ji,jj) ! k-1 bottom level of w-point but >=1 364 364 ! 365 365 ! Bottom level Dirichlet condition: … … 383 383 !CDIR NOVERRCHK 384 384 DO ji = fs_2, fs_jpim1 ! vector opt. 385 ibot = mb athy(ji,jj)386 ibotm1 = ibot-1385 ibot = mbkt(ji,jj) + 1 ! k bottom level of w-point 386 ibotm1 = mbkt(ji,jj) ! k-1 bottom level of w-point but >=1 387 387 ! 388 388 ! Bottom level Dirichlet condition: … … 623 623 !CDIR NOVERRCHK 624 624 DO ji = fs_2, fs_jpim1 ! vector opt. 625 ibot = mbathy(ji,jj)626 ibotm1 = ibot-1625 ibot = mbkt(ji,jj) + 1 ! k bottom level of w-point 626 ibotm1 = mbkt(ji,jj) ! k-1 bottom level of w-point but >=1 627 627 zdep(ji,jj) = vkarmn * rn_hbro 628 628 psi (ji,jj,ibot) = rc0**rpp * en(ji,jj,ibot)**rmm * zdep(ji,jj)**rnn … … 646 646 !CDIR NOVERRCHK 647 647 DO ji = fs_2, fs_jpim1 ! vector opt. 648 ibot = mb athy(ji,jj)649 ibotm1 = ibot-1648 ibot = mbkt(ji,jj) + 1 ! k bottom level of w-point 649 ibotm1 = mbkt(ji,jj) ! k-1 bottom level of w-point but >=1 650 650 ! 651 651 ! Bottom level Dirichlet condition: … … 825 825 DO jj = 2, jpjm1 826 826 DO ji = fs_2, fs_jpim1 ! vector opt. 827 ibot = mbathy(ji,jj) 828 avmv(ji,jj,ibot) = zcoef 827 avmv(ji,jj,mbkt(ji,jj)+1) = zcoef 829 828 END DO 830 829 END DO … … 1216 1215 DO jj = 2, jpjm1 1217 1216 DO ji = fs_2, fs_jpim1 ! vector opt. 1218 ikbu = MIN( mbathy(ji+1,jj ), mbathy(ji,jj) )1219 ikbum1 = MAX( ikbu-1, 1 )1220 ikbv = MIN( mbathy(ji,jj+1), mbathy(ji,jj) )1221 ikbvm1 = MAX( ikbv-1, 1)1217 ikbu = mbku(ji,jj) + 1 ! k bottom level of uw-point 1218 ikbum1 = mbku(ji,jj) ! k-1 bottom level of u -point, but >=1 1219 ikbv = mbkv(ji,jj) + 1 1220 ikbvm1 = mbkv(ji,jj) 1222 1221 cbx = avmu(ji,jj,ikbu) / fse3uw(ji,jj,ikbu) 1223 1222 cby = avmv(ji,jj,ikbv) / fse3vw(ji,jj,ikbv) 1224 wbotu(ji,jj) = -cbx * un(ji,jj,ikbum1) *umask(ji,jj,1)1225 wbotv(ji,jj) = -cby * vn(ji,jj,ikbvm1) *vmask(ji,jj,1)1223 wbotu(ji,jj) = -cbx * un(ji,jj,ikbum1) * umask(ji,jj,1) 1224 wbotv(ji,jj) = -cby * vn(ji,jj,ikbvm1) * vmask(ji,jj,1) 1226 1225 END DO 1227 1226 END DO … … 1235 1234 ! Initialize bottom stresses 1236 1235 DO jj = 2, jpjm1 1237 DO ji = fs_2, fs_jpim1 ! vector opt.1238 ikbu = MIN( mbathy(ji+1,jj ), mbathy(ji,jj) )1239 ikbum1 = MAX( ikbu-1, 1 )1240 ikbv = MIN( mbathy(ji,jj+1), mbathy(ji,jj) )1241 ikbvm1 = MAX( ikbv-1, 1)1242 cbx = avmu(ji,jj,ikbu) / fse3uw(ji,jj,ikbu)1243 cby = avmv(ji,jj,ikbv) / fse3vw(ji,jj,ikbv)1244 wbotu(ji,jj) = -cbx * un(ji,jj,ikbum1)*umask(ji,jj,1)1245 wbotv(ji,jj) = -cby * vn(ji,jj,ikbvm1)*vmask(ji,jj,1)1246 END DO1236 DO ji = fs_2, fs_jpim1 ! vector opt. 1237 ikbu = mbku(ji,jj) + 1 ! k bottom level of uw-point 1238 ikbum1 = mbku(ji,jj) ! k-1 bottom level of u -point, but >=1 1239 ikbv = mbkv(ji,jj) + 1 1240 ikbvm1 = mbkv(ji,jj) 1241 cbx = avmu(ji,jj,ikbu) / fse3uw(ji,jj,ikbu) 1242 cby = avmv(ji,jj,ikbv) / fse3vw(ji,jj,ikbv) 1243 wbotu(ji,jj) = -cbx * un(ji,jj,ikbum1) * umask(ji,jj,1) 1244 wbotv(ji,jj) = -cby * vn(ji,jj,ikbvm1) * vmask(ji,jj,1) 1245 END DO 1247 1246 END DO 1248 1247 ENDIF -
branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfkpp.F90
r2409 r2450 439 439 zria(ji ) = 0. 440 440 ! Maximum boundary layer depth 441 ikbot = mb athy(ji,jj) - 1! ikbot is the last T point in the water441 ikbot = mbkt(ji,jj) ! ikbot is the last T point in the water 442 442 zhmax(ji) = fsdept(ji,jj,ikbot) - 0.001 443 443 ! Compute Monin obukhov length scale at the surface and Ekman depth: -
branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfmxl.F90
r2287 r2450 69 69 70 70 ! w-level of the mixing and mixed layers 71 nmln(:,:) = mb athy(:,:) ! Initialization to the number of w ocean point mbathy72 imld(:,:) = mb athy(:,:)71 nmln(:,:) = mbkt(:,:) + 1 ! Initialization to the number of w ocean point 72 imld(:,:) = mbkt(:,:) + 1 73 73 DO jk = jpkm1, nlb10, -1 ! from the bottom to nlb10 74 74 DO jj = 1, jpj -
branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90
r2386 r2450 225 225 !CDIR NOVERRCHK 226 226 !! DO ji = fs_2, fs_jpim1 ! vector opt. 227 !! ikbu = MIN( mbathy(ji+1,jj), mbathy(ji,jj) ) 228 !! ikbv = MIN( mbathy(ji,jj+1), mbathy(ji,jj) ) 229 !! ikbum1 = MAX( ikbu-1, 1 ) 230 !! ikbvm1 = MAX( ikbv-1, 1 ) 231 !! ikbu = MIN( mbathy(ji,jj), mbathy(ji-1,jj) ) 232 !! ikbv = MIN( mbathy(ji,jj), mbathy(ji,jj-1) ) 233 !! ikbumm1 = MAX( ikbu-1, 1 ) 234 !! ikbvmm1 = MAX( ikbv-1, 1 ) 235 !! ikbt = MAX( mbathy(ji,jj), 1 ) 236 !! ztx2 = bfrua(ji-1,jj) * ub(ji-1,jj,ikbumm1) + & 237 !! bfrua(ji ,jj) * ub(ji ,jj ,ikbum1 ) 238 !! zty2 = bfrva(ji,jj ) * vb(ji ,jj ,ikbvm1) + & 239 !! bfrva(ji,jj-1) * vb(ji ,jj-1,ikbvmm1 ) 227 !! ztx2 = bfrua(ji-1,jj) * ub(ji-1,jj,mbku(ji-1,jj)) + & 228 !! bfrua(ji ,jj) * ub(ji ,jj,mbku(ji ,jj) ) 229 !! zty2 = bfrva(ji,jj ) * vb(ji,jj ,mbkv(ji,jj )) + & 230 !! bfrva(ji,jj-1) * vb(ji,jj-1,mbkv(ji,jj-1) ) 240 231 !! zebot = 0.001875_wp * SQRT( ztx2 * ztx2 + zty2 * zty2 ) ! where 0.001875 = (rn_ebb0/rau0) * 0.5 = 3.75*0.5/1000. 241 !! en (ji,jj, ikbt) = MAX( zebot, rn_emin ) * tmask(ji,jj,1)232 !! en (ji,jj,mbkt(ji,jj)+1) = MAX( zebot, rn_emin ) * tmask(ji,jj,1) 242 233 !! END DO 243 234 !! END DO … … 255 246 ! !* finite Langmuir Circulation depth 256 247 zcof = 0.5 * 0.016 * 0.016 / ( zrhoa * zcdrag ) 257 imlc(:,:) = mb athy(:,:) ! Initialization to the number of w ocean point mbathy248 imlc(:,:) = mbkt(:,:) + 1 ! Initialization to the number of w ocean point (=2 over land) 258 249 DO jk = jpkm1, 2, -1 259 250 DO jj = 1, jpj ! Last w-level at which zpelc>=0.5*us*us … … 509 500 DO ji = fs_2, fs_jpim1 ! vector opt. 510 501 zemxl = MIN( fsdepw(ji,jj,jk), zmxlm(ji,jj,jk), & 511 & fsdepw(ji,jj,mb athy(ji,jj)) - fsdepw(ji,jj,jk) )502 & fsdepw(ji,jj,mbkt(ji,jj)+1) - fsdepw(ji,jj,jk) ) 512 503 zmxlm(ji,jj,jk) = zemxl 513 504 zmxld(ji,jj,jk) = zemxl -
branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftmx.F90
r2287 r2450 371 371 DO jj = 1, jpj ! part independent of the level 372 372 DO ji = 1, jpi 373 zhdep(ji,jj) = fsdepw(ji,jj,mb athy(ji,jj)) ! depth of the ocean373 zhdep(ji,jj) = fsdepw(ji,jj,mbkt(ji,jj)+1) ! depth of the ocean 374 374 zfact(ji,jj) = rau0 * rn_htmx * ( 1. - EXP( -zhdep(ji,jj) / rn_htmx ) ) 375 375 IF( zfact(ji,jj) /= 0 ) zfact(ji,jj) = en_tmx(ji,jj) / zfact(ji,jj)
Note: See TracChangeset
for help on using the changeset viewer.