Changeset 1975 for branches/DEV_r1837_MLF
- Timestamp:
- 2010-06-28T19:22:14+02:00 (14 years ago)
- Location:
- branches/DEV_r1837_MLF/NEMO/OPA_SRC
- Files:
-
- 9 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/DEV_r1837_MLF/NEMO/OPA_SRC/DYN/dynnxt.F90
r1740 r1975 244 244 zve3b = vb(ji,jj,jk) * ze3v_b 245 245 ! 246 zuf = ( atfp 247 & / ( atfp 248 zvf = ( atfp 249 & / ( atfp 246 zuf = ( atfp * ( zue3b + zue3a ) + atfp1 * zue3n ) & 247 & / ( atfp * ( ze3u_b + ze3u_a ) + atfp1 * ze3u_n ) * umask(ji,jj,jk) 248 zvf = ( atfp * ( zve3b + zve3a ) + atfp1 * zve3n ) & 249 & / ( atfp * ( ze3v_b + ze3v_a ) + atfp1 * ze3v_n ) * vmask(ji,jj,jk) 250 250 ! 251 251 ub(ji,jj,jk) = zuf ! ub <-- filtered velocity -
branches/DEV_r1837_MLF/NEMO/OPA_SRC/DYN/sshwzv.F90
r1870 r1975 71 71 INTEGER :: ji, jj, jk ! dummy loop indices 72 72 REAL(wp) :: zcoefu, zcoefv, zcoeff ! temporary scalars 73 REAL(wp) :: z2dt, z raur! temporary scalars73 REAL(wp) :: z2dt, z1_2dt, z1_rau0 ! temporary scalars 74 74 REAL(wp), DIMENSION(jpi,jpj) :: zhdiv ! 2D workspace 75 75 REAL(wp), DIMENSION(jpi,jpj) :: z2d ! 2D workspace … … 94 94 sshv_b(ji,jj) = zcoefv * ( e1t(ji,jj ) * e2t(ji,jj ) * sshb(ji,jj ) & 95 95 & + e1t(ji,jj+1) * e2t(ji,jj+1) * sshb(ji,jj+1) ) 96 sshf_b(ji,jj) = zcoeff * ( sshb(ji ,jj) + sshb(ji ,jj+1) &97 & + sshb(ji+1,jj) + sshb(ji+1,jj+1) )98 96 sshu_n(ji,jj) = zcoefu * ( e1t(ji ,jj) * e2t(ji ,jj) * sshn(ji ,jj) & 99 97 & + e1t(ji+1,jj) * e2t(ji+1,jj) * sshn(ji+1,jj) ) 100 98 sshv_n(ji,jj) = zcoefv * ( e1t(ji,jj ) * e2t(ji,jj ) * sshn(ji,jj ) & 101 99 & + e1t(ji,jj+1) * e2t(ji,jj+1) * sshn(ji,jj+1) ) 102 sshf_n(ji,jj) = zcoeff * ( sshn(ji ,jj) + sshn(ji ,jj+1) &103 & + sshn(ji+1,jj) + sshn(ji+1,jj+1) )104 100 END DO 105 101 END DO 106 102 CALL lbc_lnk( sshu_b, 'U', 1. ) ; CALL lbc_lnk( sshu_n, 'U', 1. ) 107 103 CALL lbc_lnk( sshv_b, 'V', 1. ) ; CALL lbc_lnk( sshv_n, 'V', 1. ) 104 DO jj = 1, jpjm1 105 DO ji = 1, jpim1 ! NO Vector Opt. 106 sshf_b(ji,jj) = 0.5 * umask(ji,jj,1) * umask(ji,jj+1,1) & 107 & / ( e1f(ji,jj ) * e2f(ji,jj ) ) & 108 & * ( e1u(ji,jj ) * e2u(ji,jj ) * sshu_b(ji,jj ) & 109 & + e1u(ji,jj+1) * e2u(ji,jj+1) * sshu_b(ji,jj+1) ) 110 END DO 111 END DO 112 DO jj = 1, jpjm1 113 DO ji = 1, jpim1 ! NO Vector Opt. 114 sshf_n(ji,jj) = 0.5 * umask(ji,jj,1) * umask(ji,jj+1,1) & 115 & / ( e1f(ji,jj ) * e2f(ji,jj ) ) & 116 & * ( e1u(ji,jj ) * e2u(ji,jj ) * sshu_n(ji,jj ) & 117 & + e1u(ji,jj+1) * e2u(ji,jj+1) * sshu_n(ji,jj+1) ) 118 END DO 119 END DO 108 120 CALL lbc_lnk( sshf_b, 'F', 1. ) ; CALL lbc_lnk( sshf_n, 'F', 1. ) 109 121 ENDIF … … 145 157 ! ! After Sea Surface Height ! 146 158 ! !------------------------------! 147 z raur= 0.5 / rau0159 z1_rau0 = 0.5 / rau0 148 160 ! 149 161 zhdiv(:,:) = 0.e0 … … 152 164 END DO 153 165 ! ! Sea surface elevation time stepping 154 ssha(:,:) = ( sshb(:,:) - z2dt * ( zraur * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) ) ) * tmask(:,:,1) 166 ! In forward Euler time stepping case, the same formulation as in the leap-frog case can be used 167 ! because emp_b field is initialized with the vlaues of emp field. Hence, 0.5 * ( emp + emp_b ) = emp 168 ssha(:,:) = ( sshb(:,:) - z2dt * ( z1_rau0 * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) ) ) * tmask(:,:,1) 155 169 156 170 #if defined key_obc … … 170 184 & * ( e1t(ji,jj ) * e2t(ji,jj ) * ssha(ji,jj ) & 171 185 & + e1t(ji,jj+1) * e2t(ji,jj+1) * ssha(ji,jj+1) ) 172 sshf_a(ji,jj) = 0.25 * umask(ji,jj,1) * umask (ji,jj+1,1) &173 & * ( ssha(ji ,jj) + ssha(ji ,jj+1) &174 & + ssha(ji+1,jj) + ssha(ji+1,jj+1) )175 186 END DO 176 187 END DO 177 CALL lbc_lnk( sshu_a, 'U', 1. ) ! Boundaries conditions 188 ! Boundaries conditions 189 CALL lbc_lnk( sshu_a, 'U', 1. ) 178 190 CALL lbc_lnk( sshv_a, 'V', 1. ) 191 DO jj = 1, jpjm1 192 DO ji = 1, jpim1 ! NO Vector Opt. 193 sshf_a(ji,jj) = 0.5 * umask(ji,jj,1) * umask(ji,jj+1,1) & 194 & / ( e1f(ji,jj ) * e2f(ji,jj ) ) & 195 & * ( e1u(ji,jj ) * e2u(ji,jj ) * sshu_a(ji,jj ) & 196 & + e1u(ji,jj+1) * e2u(ji,jj+1) * sshu_a(ji,jj+1) ) 197 END DO 198 END DO 199 ! Boundaries conditions 179 200 CALL lbc_lnk( sshf_a, 'F', 1. ) 180 201 ENDIF … … 183 204 ! ! Now Vertical Velocity ! 184 205 ! !------------------------------! 206 z1_2dt = 1.e0 / z2dt 185 207 DO jk = jpkm1, 1, -1 ! integrate from the bottom the hor. divergence 186 wn(:,:,jk) = wn(:,:,jk+1) - fse3t_n(:,:,jk) * hdivn(:,:,jk) & 187 & - ( fse3t_a(:,:,jk) - fse3t_b(:,:,jk) ) * tmask(:,:,jk) / z2dt 208 ! - ML - need 3 lines here because replacement of fse3t by its expression yields too long lines otherwise 209 wn(:,:,jk) = wn(:,:,jk+1) - fse3t_n(:,:,jk) * hdivn(:,:,jk) & 210 & - ( fse3t_a(:,:,jk) - fse3t_b(:,:,jk) ) & 211 & * tmask(:,:,jk) * z1_2dt 188 212 END DO 189 213 … … 248 272 sshf_n(:,:) = sshf_a(:,:) 249 273 ELSE ! Leap-Frog time-stepping: Asselin filter + swap 274 zec = atfp * rdt / rau0 250 275 DO jj = 1, jpj 251 276 DO ji = 1, jpi ! before <-- now filtered 252 sshb (ji,jj) = sshn (ji,jj) + atfp * ( sshb (ji,jj) - 2 * sshn (ji,jj) + ssha (ji,jj) ) 253 sshu_b(ji,jj) = sshu_n(ji,jj) + atfp * ( sshu_b(ji,jj) - 2 * sshu_n(ji,jj) + sshu_a(ji,jj) ) 254 sshv_b(ji,jj) = sshv_n(ji,jj) + atfp * ( sshv_b(ji,jj) - 2 * sshv_n(ji,jj) + sshv_a(ji,jj) ) 255 sshf_b(ji,jj) = sshf_n(ji,jj) + atfp * ( sshf_b(ji,jj) - 2 * sshf_n(ji,jj) + sshf_a(ji,jj) ) 277 sshb (ji,jj) = sshn (ji,jj) + atfp * ( sshb (ji,jj) - 2 * sshn (ji,jj) + ssha (ji,jj) ) & 278 & - zec * ( emp_b(ji,jj) - emp(ji,jj) ) * tmask(ji,jj,1) 256 279 sshn (ji,jj) = ssha (ji,jj) ! now <-- after 257 280 sshu_n(ji,jj) = sshu_a(ji,jj) … … 260 283 END DO 261 284 END DO 285 DO jj = 1, jpjm1 286 DO ji = 1, jpim1 ! NO Vector Opt. 287 sshu_b(ji,jj) = 0.5 * umask(ji,jj,1) / ( e1u(ji ,jj) * e2u(ji ,jj) ) & 288 & * ( e1t(ji ,jj) * e2t(ji ,jj) * sshb(ji ,jj) & 289 & + e1t(ji+1,jj) * e2t(ji+1,jj) * sshb(ji+1,jj) ) 290 sshv_b(ji,jj) = 0.5 * vmask(ji,jj,1) / ( e1v(ji,jj ) * e2v(ji,jj ) ) & 291 & * ( e1t(ji,jj ) * e2t(ji,jj ) * sshb(ji,jj ) & 292 & + e1t(ji,jj+1) * e2t(ji,jj+1) * sshb(ji,jj+1) ) 293 END DO 294 END DO 295 ! Boundaries conditions 296 CALL lbc_lnk( sshu_b, 'U', 1. ) 297 CALL lbc_lnk( sshv_b, 'V', 1. ) 298 DO jj = 1, jpjm1 299 DO ji = 1, jpim1 ! NO Vector Opt. 300 sshf_b(ji,jj) = 0.5 * umask(ji,jj,1) * umask(ji,jj+1,1) & 301 & / ( e1f(ji,jj ) * e2f(ji,jj ) ) & 302 & * ( e1u(ji,jj ) * e2u(ji,jj ) * sshu_b(ji,jj ) & 303 & + e1u(ji,jj+1) * e2u(ji,jj+1) * sshu_b(ji,jj+1) ) 304 END DO 305 END DO 306 ! Boundaries conditions 307 CALL lbc_lnk( sshf_b, 'F', 1. ) 262 308 ENDIF 263 309 ! !--------------------------! … … 268 314 ! 269 315 ELSE ! Leap-Frog time-stepping: Asselin filter + swap 270 zec = atfp * rdt / rau0271 316 DO jj = 1, jpj 272 317 DO ji = 1, jpi ! before <-- now filtered 273 sshb(ji,jj) = sshn(ji,jj) + atfp * ( sshb(ji,jj) - 2 * sshn(ji,jj) + ssha(ji,jj) ) & 274 & - zec * ( emp_b(ji,jj) - emp(ji,jj) ) * tmask(ji,jj,1) 318 sshb(ji,jj) = sshn(ji,jj) + atfp * ( sshb(ji,jj) - 2 * sshn(ji,jj) + ssha(ji,jj) ) 275 319 sshn(ji,jj) = ssha(ji,jj) ! now <-- after 276 320 END DO … … 279 323 ! 280 324 ENDIF 281 282 ! time filter with global conservation correction and array swap283 sshbb(:,:) = sshb(:,:)284 sshb (:,:) = sshn(:,:) + atfp * ( sshb(:,:) - 2 * sshn(:,:) + zssha(:,:) ) &285 & - zfact *286 sshn (:,:) = zssha(:,:)287 empb (:,:) = emp(:,:)288 289 290 325 ! 291 326 IF(ln_ctl) CALL prt_ctl( tab2d_1=sshb, clinfo1=' sshb - : ', mask1=tmask, ovlap=1 ) -
branches/DEV_r1837_MLF/NEMO/OPA_SRC/IOM/restart.F90
r1613 r1975 217 217 hdivb(:,:,:) = hdivn(:,:,:) 218 218 sshb (:,:) = sshn (:,:) 219 ! - ML - sshbnc 219 220 ENDIF 220 221 ! -
branches/DEV_r1837_MLF/NEMO/OPA_SRC/SBC/sbc_oce.F90
r1870 r1975 52 52 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: emps , emps_b !: freshwater budget: concentration/dillution [Kg/m2/s] 53 53 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: emp_tot !: total E-P-R over ocean and ice [Kg/m2/s] 54 ! - ML - beginning 55 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: sbc_trd_hc_n !: sbc heat content trend now [K/m/s] 56 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: sbc_trd_hc_b !: " " " " before " 57 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: sbc_trd_sc_n !: sbc salt content trend now [psu/m/s] 58 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: sbc_trd_sc_b !: " " " " before " 59 REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) :: qsr_trd_hc_n !: heat content trend due to qsr flux now [K/m/s] 60 REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) :: qsr_trd_hc_b !: " " " " " " " before " 61 ! - ML - end 54 62 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: tprecip !: total precipitation [Kg/m2/s] 55 63 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: sprecip !: solid precipitation [Kg/m2/s] -
branches/DEV_r1837_MLF/NEMO/OPA_SRC/SBC/sbcmod.F90
r1870 r1975 193 193 ! ! ---------------------------------------- ! 194 194 utau_b(:,:) = utau(:,:) ! Swap the ocean forcing fields 195 utau_b(:,:) = utau(:,:) ! (except at nitOOOwhere before fields195 vtau_b(:,:) = vtau(:,:) ! (except at nit000 where before fields 196 196 qns_b (:,:) = qns (:,:) ! are set the end of the routine) 197 197 qsr_b (:,:) = qsr (:,:) 198 198 emp_b (:,:) = emp (:,:) 199 199 emps_b(:,:) = emps(:,:) 200 ! - ML - 201 sbc_trd_hc_b(:,:) = sbc_trd_hc_n(:,:) 202 qsr_trd_hc_b(:,:,:) = qsr_trd_hc_n(:,:,:) 203 IF ( .NOT. lk_vvl ) sbc_trd_sc_b(:,:) = sbc_trd_sc_n(:,:) 204 200 205 ENDIF 201 206 … … 256 261 IF(lwp) WRITE(numout,*) ' nit000-1 surface forcing fields red in the restart file' 257 262 CALL iom_get( numror, jpdom_autoglo, 'utau_b', utau_b ) ! before i-stress (U-point) 258 CALL iom_get( numror, jpdom_autoglo, ' utau_b', utau_b ) ! before j-stress (V-point)263 CALL iom_get( numror, jpdom_autoglo, 'vtau_b', vtau_b ) ! before j-stress (V-point) 259 264 CALL iom_get( numror, jpdom_autoglo, 'qns_b' , qns_b ) ! before non solar heat flux (T-point) 260 265 CALL iom_get( numror, jpdom_autoglo, 'qsr_b' , qsr_b ) ! before solar heat flux (T-point) 261 266 CALL iom_get( numror, jpdom_autoglo, 'emp_b' , emp_b ) ! before freshwater flux (T-point) 262 267 CALL iom_get( numror, jpdom_autoglo, 'emps_b', emp_b ) ! before C/D freshwater flux (T-point) 268 ! - ML - 269 CALL iom_get( numror, jpdom_autoglo, 'sbc_trd_hc_b', sbc_trd_hc_b ) ! before heat content sbc trend 270 CALL iom_get( numror, jpdom_autoglo, 'qsr_trd_hc_b', qsr_trd_hc_b ) ! before heat content trend due to Qsr flux 271 IF ( .NOT. lk_vvl ) THEN 272 CALL iom_get( numror, jpdom_autoglo, 'sbc_trd_sc_b', sbc_trd_sc_b ) ! before salt content sbc trend 273 ENDIF 263 274 ! 264 275 ELSE !* no restart: set from nit000 values 265 276 IF(lwp) WRITE(numout,*) ' nit000-1 surface forcing fields set to nit000' 266 277 utau_b(:,:) = utau(:,:) 267 utau_b(:,:) = utau(:,:)278 vtau_b(:,:) = vtau(:,:) 268 279 qns_b (:,:) = qns (:,:) 269 280 qsr_b (:,:) = qsr (:,:) … … 280 291 & 'at it= ', kt,' date= ', ndastp 281 292 IF(lwp) WRITE(numout,*) '~~~~' 282 CALL iom_rstput( kt, nitrst, numrow, 'utau_b' , utau ) ! 283 CALL iom_rstput( kt, nitrst, numrow, 'utau_b' , vtau ) 284 CALL iom_rstput( kt, nitrst, numrow, 'qns_b' , qns ) 285 CALL iom_rstput( kt, nitrst, numrow, 'qsr_b' , qsr ) 286 CALL iom_rstput( kt, nitrst, numrow, 'emp_b' , emp ) 287 CALL iom_rstput( kt, nitrst, numrow, 'emps_b' , emp ) 293 CALL iom_rstput( kt, nitrst, numrow, 'utau_b' , utau_b ) ! 294 CALL iom_rstput( kt, nitrst, numrow, 'vtau_b' , vtau_b ) 295 CALL iom_rstput( kt, nitrst, numrow, 'qns_b' , qns_b ) 296 CALL iom_rstput( kt, nitrst, numrow, 'qsr_b' , qsr_b ) 297 CALL iom_rstput( kt, nitrst, numrow, 'emp_b' , emp_b ) 298 CALL iom_rstput( kt, nitrst, numrow, 'emps_b' , emps_b ) 299 ! - ML - 300 CALL iom_rstput( kt, nitrst, numrow, 'sbc_trd_hc_b', sbc_trd_hc_b ) 301 CALL iom_rstput( kt, nitrst, numrow, 'qsr_trd_hc_b', qsr_trd_hc_b ) 302 IF ( .NOT. lk_vvl ) THEN 303 CALL iom_rstput( kt, nitrst, numrow, 'sbc_trd_sc_b', sbc_trd_sc_b ) 304 ENDIF 288 305 ! 289 306 ENDIF -
branches/DEV_r1837_MLF/NEMO/OPA_SRC/TRA/tranxt.F90
r1870 r1975 25 25 USE oce ! ocean dynamics and tracers variables 26 26 USE dom_oce ! ocean space and time domain variables 27 USE sbc_oce ! surface boundary condition: ocean 27 28 USE zdf_oce ! ??? 28 29 USE domvvl ! variable volume … … 37 38 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 38 39 USE prtctl ! Print control 40 USE traqsr ! penetrative solar radiation (needed for nksr) 39 41 USE agrif_opa_update 40 42 USE agrif_opa_interp … … 45 47 PUBLIC tra_nxt ! routine called by step.F90 46 48 47 REAL(wp) :: rbcp , r1_2bcp! Brown & Campana parameters for semi-implicit hpg49 REAL(wp) :: rbcp ! Brown & Campana parameters for semi-implicit hpg 48 50 REAL(wp), DIMENSION(jpk) :: r2dt_t ! vertical profile time step, =2*rdttra (leapfrog) or =rdttra (Euler) 49 51 … … 98 100 IF(lwp) WRITE(numout,*) '~~~~~~~' 99 101 ! 100 rbcp = 0.25 * (1 + atfp) * (1 + atfp * atfp) ! Brown & Campana parameter for semi-implicit hpg 101 r1_2bcp = 1.e0 - 2.e0 * rbcp 102 rbcp = 0.25 * (1. + atfp) * (1. + atfp) * ( 1. - atfp) ! Brown & Campana parameter for semi-implicit hpg 102 103 ENDIF 103 104 ! ! set time step size (Euler/Leapfrog) … … 164 165 !! - swap tracer fields to prepare the next time_step. 165 166 !! This can be summurized for temperature as: 166 !! ztm = (rbcp*ta+(2-rbcp)*tn+rbcp*tb)ln_dynhpg_imp = T167 !! ztm = tn + rbcp * [ta -2 tn + tb ] ln_dynhpg_imp = T 167 168 !! ztm = 0 otherwise 168 !! with rbcp= (1+atfp)*(1+atfp*atfp)/4169 !! with rbcp=1/4 * (1-atfp^4) / (1-atfp) 169 170 !! tb = tn + atfp*[ tb - 2 tn + ta ] 170 171 !! tn = ta … … 177 178 !! 178 179 INTEGER :: ji, jj, jk ! dummy loop indices 179 REAL(wp) :: zt m, ztf, zac! temporary scalars180 REAL(wp) :: z sm, zsf! - -180 REAL(wp) :: zt_m, zs_m ! temporary scalars 181 REAL(wp) :: ztn, zsn ! - - 181 182 !!---------------------------------------------------------------------- 182 183 … … 187 188 ENDIF 188 189 ! 189 ! ! ----------------------- ! 190 IF( ln_dynhpg_imp ) THEN ! semi-implicite hpg case ! 191 ! ! ----------------------- ! 192 ! 193 IF( neuler == 0 .AND. kt == nit000 ) THEN ! Euler time-stepping at first time-step 194 DO jk = 1, jpkm1 ! (only swap) 195 tn(:,:,jk) = ta(:,:,jk) ! tb <-- tn 196 sn(:,:,jk) = sa(:,:,jk) 197 END DO 198 ELSE ! general case (Leapfrog + Asselin filter) 199 DO jk = 1, jpkm1 200 DO jj = 1, jpj 201 DO ji = 1, jpi 202 ztm = rbcp * ( ta(ji,jj,jk) + tb(ji,jj,jk) ) + r1_2bcp * tn(ji,jj,jk) ! mean t 203 zsm = rbcp * ( sa(ji,jj,jk) + sb(ji,jj,jk) ) + r1_2bcp * sn(ji,jj,jk) 204 ztf = atfp * ( ta(ji,jj,jk) + tb(ji,jj,jk) - 2. * tn(ji,jj,jk) ) ! Asselin filter on t 205 zsf = atfp * ( sa(ji,jj,jk) + sb(ji,jj,jk) - 2. * sn(ji,jj,jk) ) 206 tb(ji,jj,jk) = tn(ji,jj,jk) + ztf ! tb <-- filtered tn 207 sb(ji,jj,jk) = sn(ji,jj,jk) + zsf 208 tn(ji,jj,jk) = ta(ji,jj,jk) ! tn <-- ta 209 sn(ji,jj,jk) = sa(ji,jj,jk) 210 ta(ji,jj,jk) = ztm ! ta <-- mean t 211 sa(ji,jj,jk) = zsm 212 END DO 190 IF( neuler == 0 .AND. kt == nit000 ) THEN ! Euler time-stepping at first time-step 191 DO jk = 1, jpkm1 ! (only swap) 192 tn(:,:,jk) = ta(:,:,jk) ! tn <-- ta 193 sn(:,:,jk) = sa(:,:,jk) ! sn <-- sa 194 END DO 195 ELSE ! General case (Leapfrog + Asselin filter) 196 DO jk = 1, jpkm1 197 DO jj = 1, jpj 198 DO ji = 1, jpi 199 IF( ln_dynhpg_imp ) THEN ! implicit hpg: keep tn, sn in memory 200 ztn = tn(ji,jj,jk) 201 zsn = sn(ji,jj,jk) 202 ENDIF 203 ! ! time laplacian on tracers 204 ! ! used for both Asselin and Brown & Campana filters 205 zt_m = ta(ji,jj,jk) - 2. * tn(ji,jj,jk) + tb(ji,jj,jk) 206 zs_m = sa(ji,jj,jk) - 2. * sn(ji,jj,jk) + sb(ji,jj,jk) 207 ! 208 ! ! swap of arrays 209 tb(ji,jj,jk) = tn(ji,jj,jk) + atfp * zt_m ! tb <-- tn filtered 210 sb(ji,jj,jk) = sn(ji,jj,jk) + atfp * zs_m ! sb <-- sn filtered 211 tn(ji,jj,jk) = ta(ji,jj,jk) ! tn <-- ta 212 sn(ji,jj,jk) = sa(ji,jj,jk) ! sn <-- sa 213 ! ! semi imlicit hpg computation (Brown & Campana) 214 IF( ln_dynhpg_imp ) THEN 215 ta(ji,jj,jk) = ztn + rbcp * zt_m ! ta <-- Brown & Campana average 216 sa(ji,jj,jk) = zsn + rbcp * zs_m ! sa <-- Brown & Campana average 217 ENDIF 213 218 END DO 214 219 END DO 215 ENDIF 216 ! ! ----------------------- ! 217 ELSE ! explicit hpg case ! 218 ! ! ----------------------- ! 219 ! 220 IF( neuler == 0 .AND. kt == nit000 ) THEN ! Euler time-stepping at first time-step 221 DO jk = 1, jpkm1 222 tn(:,:,jk) = ta(:,:,jk) ! tn <-- ta 223 sn(:,:,jk) = sa(:,:,jk) 224 END DO 225 ELSE ! general case (Leapfrog + Asselin filter) 226 DO jk = 1, jpkm1 227 DO jj = 1, jpj 228 DO ji = 1, jpi 229 ztf = atfp * ( ta(ji,jj,jk) - 2.* tn(ji,jj,jk) + tb(ji,jj,jk) ) ! Asselin filter on t 230 zsf = atfp * ( sa(ji,jj,jk) - 2.* sn(ji,jj,jk) + sb(ji,jj,jk) ) 231 tb(ji,jj,jk) = tn(ji,jj,jk) + ztf ! tb <-- filtered tn 232 sb(ji,jj,jk) = sn(ji,jj,jk) + zsf 233 tn(ji,jj,jk) = ta(ji,jj,jk) ! tn <-- ta 234 sn(ji,jj,jk) = sa(ji,jj,jk) 235 END DO 236 END DO 237 END DO 238 ENDIF 239 ENDIF 240 ! 241 !!gm from Matthieu : unchecked 242 IF( neuler /= 0 .OR. kt /= nit000 ) THEN ! remove the forcing from the Asselin filter 243 zac = atfp * rdttra(1) 244 tb(:,:,1) = tb(:,:,1) - zac * ( qns (:,:) - qns_b (:,:) ) ! non solar surface flux 245 sb(:,:,1) = sn(:,:,1) - zac * ( emps(:,:) - emps_b(:,:) ) ! surface salt flux 246 ! 247 IF( ln_traqsr ) ! solar penetrating flux 248 DO jk = 1, nksr 249 DO jj = 1, jpj 250 DO ji = 1, jpi 251 IF( ln_sco ) THEN 252 z_cor_qsr = etot3(ji,jj,jk) * ( qsr(ji,jj) - qsrb(ji,jj) ) 253 ELSEIF( ln_zps .OR. ln_zco ) THEN 254 z_cor_qsr = ( qsr(ji,jj) - qsrb(ji,jj) ) * & 255 & ( gdsr(jk) * tmask(ji,jj,jk) - gdsr(jk+1) * tmask(ji,jj,jk+1) ) 256 ENDIF 257 zt = zt - zfact1 * z_cor_qsr 258 END DO 259 END DO 220 END DO 260 221 ENDIF 261 222 ! … … 276 237 !! - swap tracer fields to prepare the next time_step. 277 238 !! This can be summurized for temperature as: 278 !! ztm = ( rbcp*e3t_a*ta + (2-rbtp)*e3t_n*tn + rbcp*e3t_b*tb) ln_dynhpg_imp = T279 !! /( rbcp*e3t_a + (2-rbcp)*e3t_n + rbcp*e3t_b)239 !! ztm = ( e3t_n*tn + rbcp*[ e3t_b*tb - 2 e3t_n*tn + e3t_a*ta ] ) ln_dynhpg_imp = T 240 !! /( e3t_n + rbcp*[ e3t_b - 2 e3t_n + e3t_a ] ) 280 241 !! ztm = 0 otherwise 281 242 !! tb = ( e3t_n*tn + atfp*[ e3t_b*tb - 2 e3t_n*tn + e3t_a*ta ] ) … … 287 248 !! - (ta,sa) time averaged (t,s) (ln_dynhpg_imp = T) 288 249 !!---------------------------------------------------------------------- 289 INTEGER, INTENT(in) :: kt ! ocean time-step index250 INTEGER, INTENT(in) :: kt ! ocean time-step index 290 251 !! 291 INTEGER :: ji, jj, jk ! dummy loop indices 292 REAL(wp) :: ztm , ztc_f , ztf , ztca, ztcn, ztcb ! temporary scalar 293 REAL(wp) :: zsm , zsc_f , zsf , zsca, zscn, zscb ! - - 294 REAL(wp) :: ze3mr, ze3fr ! - - 295 REAL(wp) :: ze3t_b, ze3t_n, ze3t_a, ze3t_f ! - - 296 !!---------------------------------------------------------------------- 297 252 INTEGER :: ji, jj, jk ! dummy loop indices 253 REAL :: ze3t_a, ze3t_n, ze3t_b ! temporary scalars 254 REAL :: ztc_a, ztc_n, ztc_b ! - - 255 REAL :: zsc_a, zsc_n, zsc_b ! - - 256 REAL :: ztc_f, zsc_f, ztc_m, zsc_m ! - - 257 REAL :: ze3t_f, ze3t_m ! - - 258 REAL :: zfact1, zfact2 ! - - 259 !!---------------------------------------------------------------------- 260 298 261 IF( kt == nit000 ) THEN 299 262 IF(lwp) WRITE(numout,*) … … 302 265 ENDIF 303 266 304 ! ! ----------------------- ! 305 IF( ln_dynhpg_imp ) THEN ! semi-implicite hpg case ! (mean tracer computation) 306 ! ! ----------------------- ! 307 ! 308 IF( neuler == 0 .AND. kt == nit000 ) THEN ! Euler time-stepping at first time-step 309 DO jk = 1, jpkm1 ! (only swap) 310 tn(:,:,jk) = ta(:,:,jk) ! tn <-- ta 311 sn(:,:,jk) = sa(:,:,jk) 312 END DO 313 ! ! general case (Leapfrog + Asselin filter) 314 ELSE ! apply filter on thickness weighted tracer and swap 315 DO jk = 1, jpkm1 316 DO jj = 1, jpj 317 DO ji = 1, jpi 318 ze3t_b = fse3t_b(ji,jj,jk) 319 ze3t_n = fse3t_n(ji,jj,jk) 320 ze3t_a = fse3t_a(ji,jj,jk) 321 ! ! tracer content at Before, now and after 322 ztcb = tb(ji,jj,jk) * ze3t_b ; zscb = sb(ji,jj,jk) * ze3t_b 323 ztcn = tn(ji,jj,jk) * ze3t_n ; zscn = sn(ji,jj,jk) * ze3t_n 324 ztca = ta(ji,jj,jk) * ze3t_a ; zsca = sa(ji,jj,jk) * ze3t_a 325 ! 326 ! ! Asselin filter on thickness and tracer content 327 ze3t_f = atfp * ( ze3t_a + ze3t_b - 2.* ze3t_n ) 328 ztc_f = atfp * ( ztca + ztcb - 2.* ztcn ) 329 zsc_f = atfp * ( zsca + zscb - 2.* zscn ) 330 ! 331 ! ! filtered tracer including the correction 332 ze3fr = 1.e0 / ( ze3t_n + ze3t_f ) 333 ztf = ze3fr * ( ztcn + ztc_f ) 334 zsf = ze3fr * ( zscn + zsc_f ) 335 ! ! mean thickness and tracer 336 ze3mr = 1.e0 / ( rbcp * ( ze3t_a + ze3t_b ) + r1_2bcp * ze3t_n ) 337 ztm = ze3mr * ( rbcp * ( ztca + ztcb ) + r1_2bcp * ztcn ) 338 zsm = ze3mr * ( rbcp * ( zsca + zscb ) + r1_2bcp * zscn ) 339 !!gm mean e3t have to be saved and used in dynhpg or it can be recomputed in dynhpg !! 340 !!gm e3t_m(ji,jj,jk) = 1.e0 / ze3mr 341 ! ! swap of arrays 342 tb(ji,jj,jk) = ztf ! tb <-- tn + filter 343 sb(ji,jj,jk) = zsf 344 tn(ji,jj,jk) = ta(ji,jj,jk) ! tn <-- ta 345 sn(ji,jj,jk) = sa(ji,jj,jk) 346 ta(ji,jj,jk) = ztm ! ta <-- mean t 347 sa(ji,jj,jk) = zsm 348 END DO 267 IF( neuler == 0 .AND. kt == nit000 ) THEN ! Euler time-stepping at first time-step 268 ! ! (only swap) 269 DO jk = 1, jpkm1 ! tn <-- ta 270 tn(:,:,jk) = ta(:,:,jk) ! sn <-- sa 271 sn(:,:,jk) = sa(:,:,jk) 272 END DO 273 ! ! General case (Leapfrog + Asselin filter) 274 ELSE ! apply filter on thickness weighted tracer and swap 275 DO jk = 1, jpkm1 276 zfact1 = atfp * r2dt_t(jk) 277 zfact2 = zfact1 / rau0 278 DO jj = 1, jpj 279 DO ji = 1, jpi 280 ! ! scale factors at Before, now and after 281 ze3t_b = fse3t_b(ji,jj,jk) 282 ze3t_n = fse3t_n(ji,jj,jk) 283 ze3t_a = fse3t_a(ji,jj,jk) 284 ! ! tracer content at Before, now and after 285 ztc_b = tb(ji,jj,jk) * ze3t_b ; zsc_b = sb(ji,jj,jk) * ze3t_b 286 ztc_n = tn(ji,jj,jk) * ze3t_n ; zsc_n = sn(ji,jj,jk) * ze3t_n 287 ztc_a = ta(ji,jj,jk) * ze3t_a ; zsc_a = sa(ji,jj,jk) * ze3t_a 288 ! 289 ! ! Time laplacian on thickness and tracer content 290 ! ! used for both Asselin and Brown & Campana filters 291 ze3t_m = ze3t_a - 2. * ze3t_n + ze3t_b 292 ztc_m = ztc_a - 2. * ztc_n + ztc_b 293 zsc_m = zsc_a - 2. * zsc_n + zsc_b 294 ! ! Asselin Filter + correction 295 ze3t_f = ze3t_n + atfp * ze3t_m 296 ztc_f = ztc_n + atfp * ztc_m 297 zsc_f = zsc_n + atfp * zsc_m 298 299 IF( jk == 1 ) THEN 300 ze3t_f = ze3t_f - zfact2 * ( emp_b (ji,jj) - emp (ji,jj) ) 301 ztc_f = ztc_f - zfact1 * ( sbc_trd_hc_n(ji,jj) - sbc_trd_hc_b(ji,jj) ) 302 ENDIF 303 IF( ln_traqsr .AND. ( jk .LE. nksr ) ) THEN 304 ztc_f = ztc_f - zfact1 * ( qsr_trd_hc_n(ji,jj,jk) - qsr_trd_hc_b(ji,jj,jk) ) 305 ENDIF 306 ! ! swap of arrays 307 ze3t_f = 1.e0 / ze3t_f 308 tb(ji,jj,jk) = ztc_f * ze3t_f ! tb <-- tn filtered 309 sb(ji,jj,jk) = zsc_f * ze3t_f ! sb <-- sn filtered 310 tn(ji,jj,jk) = ta(ji,jj,jk) ! tn <-- ta 311 sn(ji,jj,jk) = sa(ji,jj,jk) ! sn <-- sa 312 ! ! semi imlicit hpg computation (Brown & Campana) 313 IF( ln_dynhpg_imp ) THEN 314 ze3t_m = 1.e0 / ( ze3t_n + rbcp * ze3t_m ) 315 ta(ji,jj,jk) = ( ztc_n + rbcp * ztc_m ) * ze3t_m ! ta <-- Brown & Campana average 316 sa(ji,jj,jk) = ( zsc_n + rbcp * zsc_m ) * ze3t_m ! sa <-- Brown & Campana average 317 ENDIF 349 318 END DO 350 319 END DO 351 ENDIF 352 ! ! ----------------------- ! 353 ELSE ! explicit hpg case ! (NO mean tracer computation) 354 ! ! ----------------------- ! 355 ! 356 IF( neuler == 0 .AND. kt == nit000 ) THEN ! case of Euler time-stepping at first time-step 357 DO jk = 1, jpkm1 ! (only swap) 358 tn(:,:,jk) = ta(:,:,jk) ! tn <-- ta 359 sn(:,:,jk) = sa(:,:,jk) 360 END DO 361 ! ! general case (Leapfrog + Asselin filter) 362 ELSE ! apply filter on thickness weighted tracer and swap 363 DO jk = 1, jpkm1 364 DO jj = 1, jpj 365 DO ji = 1, jpi 366 ze3t_b = fse3t_b(ji,jj,jk) 367 ze3t_n = fse3t_n(ji,jj,jk) 368 ze3t_a = fse3t_a(ji,jj,jk) 369 ! ! tracer content at Before, now and after 370 ztcb = tb(ji,jj,jk) * ze3t_b ; zscb = sb(ji,jj,jk) * ze3t_b 371 ztcn = tn(ji,jj,jk) * ze3t_n ; zscn = sn(ji,jj,jk) * ze3t_n 372 ztca = ta(ji,jj,jk) * ze3t_a ; zsca = sa(ji,jj,jk) * ze3t_a 373 ! 374 ! ! Asselin filter on thickness and tracer content 375 ze3t_f = atfp * ( ze3t_a - 2.* ze3t_n + ze3t_b ) 376 ztc_f = atfp * ( ztca - 2.* ztcn + ztcb ) 377 zsc_f = atfp * ( zsca - 2.* zscn + zscb ) 378 ! 379 ! ! filtered tracer including the correction 380 ze3fr = 1.e0 / ( ze3t_n + ze3t_f ) 381 ztf = ( ztcn + ztc_f ) * ze3fr 382 zsf = ( zscn + zsc_f ) * ze3fr 383 ! ! swap of arrays 384 tb(ji,jj,jk) = ztf ! tb <-- tn filtered 385 sb(ji,jj,jk) = zsf 386 tn(ji,jj,jk) = ta(ji,jj,jk) ! tn <-- ta 387 sn(ji,jj,jk) = sa(ji,jj,jk) 388 END DO 389 END DO 390 END DO 391 ENDIF 392 ENDIF 393 !!gm from Matthieu : unchecked 394 ! 395 IF( neuler /= 0 .OR. kt /= nit000 ) THEN ! remove the forcing from the Asselin filter 396 IF( lk_vvl ) THEN ! Varying levels 397 DO jj = 1, jpj 398 DO ji = 1, jpi 399 ! - ML - modification for varaiance diagnostics 400 zssh1 = tmask(ji,jj,jk) / ( fse3t(ji,jj,jk) + atfp * ( zfse3ta(ji,jj,jk) - 2*fse3t(ji,jj,jk) & 401 & + fse3tb(ji,jj,jk) ) ) 402 zt = zssh1 * ( tn(ji,jj,jk) + atfp * ( ta(ji,jj,jk) - 2 * tn(ji,jj,jk) + tb(ji,jj,jk) ) ) 403 zs = zssh1 * ( sn(ji,jj,jk) + atfp * ( sa(ji,jj,jk) - 2 * sn(ji,jj,jk) + sb(ji,jj,jk) ) ) 404 ! filter correction for global conservation 405 !------------------------------------------ 406 zfact1 = atfp * rdttra(1) * zssh1 407 IF (jk == 1) THEN ! remove sbc trend from time filter 408 zt = zt - zfact1 * ( qn(ji,jj) - qb(ji,jj) ) 409 !!??? zs = zs - zfact1 * ( - emp( ji,jj) + empb(ji,jj) ) * zsrau * 35.5e0 410 ENDIF 411 ! remove qsr trend from time filter 412 IF (jk <= nksr) THEN 413 IF( ln_sco ) THEN 414 z_cor_qsr = etot3(ji,jj,jk) * ( qsr(ji,jj) - qsrb(ji,jj) ) 415 ELSEIF( ln_zps .OR. ln_zco ) THEN 416 z_cor_qsr = ( qsr(ji,jj) - qsrb(ji,jj) ) * & 417 & ( gdsr(jk) * tmask(ji,jj,jk) - gdsr(jk+1) * tmask(ji,jj,jk+1) ) 418 ENDIF 419 zt = zt - zfact1 * z_cor_qsr 420 IF (jk == nksr) qsrb(ji,jj) = qsr(ji,jj) 421 ENDIF 422 ! - ML - end of modification 423 zssh1 = tmask(ji,jj,1) / zfse3ta(ji,jj,jk) 424 tn(ji,jj,jk) = ta(ji,jj,jk) * zssh1 425 sn(ji,jj,jk) = sa(ji,jj,jk) * zssh1 426 END DO 427 END DO 428 ENDIF 429 !!gm end 320 END DO 321 ENDIF 430 322 ! 431 323 END SUBROUTINE tra_nxt_vvl -
branches/DEV_r1837_MLF/NEMO/OPA_SRC/TRA/traqsr.F90
r1870 r1975 47 47 TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_chl ! structure of input Chl (file informations, fields read) 48 48 49 INTEGER 49 INTEGER , PUBLIC :: nksr ! levels below which the light cannot penetrate (depth larger than 391 m) 50 50 REAL(wp), DIMENSION(3,61) :: rkrgb ! tabulated attenuation coefficients for RGB absorption 51 51 … … 97 97 REAL(wp) :: zchl, zcoef, zsi0r ! temporary scalars 98 98 REAL(wp) :: zc0, zc1, zc2, zc3 ! - - 99 REAL(wp) :: zqsr ! - -100 99 REAL(wp), DIMENSION(jpi,jpj) :: zekb, zekg, zekr ! 2D workspace 101 100 REAL(wp), DIMENSION(jpi,jpj,jpk) :: ze0, ze1 , ze2, ze3, zea ! 3D workspace … … 123 122 DO ji = fs_2, fs_jpim1 ! vector opt. 124 123 !!gm how to stecify the mean of time step here : TOP versus OPA time stepping strategy not obvious 125 ta(ji,jj,jk) = ta(ji,jj,jk) + ro0cpr * ( etot3(ji,jj,jk) - etot3(ji,jj,jk+1) ) / fse3t(ji,jj,jk)124 qsr_trd_hc_n(ji,jj,jk) = ro0cpr * ( etot3(ji,jj,jk) - etot3(ji,jj,jk+1) ) 126 125 END DO 127 126 END DO … … 155 154 zsi0r = 1.e0 / rn_si0 156 155 zcoef = ( 1. - rn_abs ) / 3.e0 ! equi-partition in R-G-B 157 !!gm bug !!! zqsr only a constant not an array 158 zqsr = 0.5 * ( qsr_b(:,:) + qsr(:,:) ) ! mean over 2 time steps 159 ze0(:,:,1) = rn_abs * zqsr 160 ze1(:,:,1) = zcoef * zqsr 161 ze2(:,:,1) = zcoef * zqsr 162 ze3(:,:,1) = zcoef * zqsr 163 zea(:,:,1) = zqsr 156 157 ze0(:,:,1) = rn_abs * qsr(:,:) 158 ze1(:,:,1) = zcoef * qsr(:,:) 159 ze2(:,:,1) = zcoef * qsr(:,:) 160 ze3(:,:,1) = zcoef * qsr(:,:) 161 zea(:,:,1) = qsr(:,:) 164 162 ! 165 163 DO jk = 2, nksr+1 ! deeper values … … 183 181 ! 184 182 DO jk = 1, nksr ! compute and add qsr trend to ta 185 ta(:,:,jk) = ta(:,:,jk) + ro0cpr * ( zea(:,:,jk) - zea(:,:,jk+1) ) / fse3t(:,:,jk)183 qsr_trd_hc_n(:,:,jk) = ro0cpr * ( zea(:,:,jk) - zea(:,:,jk+1) ) 186 184 END DO 187 185 zea(:,:,nksr+1:jpk) = 0.e0 ! below 400m set to zero … … 190 188 ELSE !* Constant Chlorophyll 191 189 DO jk = 1, nksr 192 ta(:,:,jk) = ta(:,:,jk) + etot3(:,:,jk) * 0.5 * ( qsr_b(:,:) + qsr(:,:))190 qsr_trd_hc_n(:,:,jk) = etot3(:,:,jk) * qsr(:,:) 193 191 END DO 194 192 ENDIF … … 201 199 DO jj = 2, jpjm1 202 200 DO ji = fs_2, fs_jpim1 ! vector opt. 203 ta(ji,jj,jk) = ta(ji,jj,jk) + etot3(ji,jj,jk) * 0.5 * ( qsr_b(:,:) + qsr(:,:))201 qsr_trd_hc_n(ji,jj,jk) = etot3(ji,jj,jk) * qsr(ji,jj) 204 202 END DO 205 203 END DO … … 208 206 ENDIF 209 207 ! 208 ENDIF 209 210 ! Add qsr trend to ta in all cases 211 IF( neuler == 0 .AND. kt == nit000 ) THEN 212 DO jk = 1, nksr 213 DO jj = 2, jpjm1 214 DO ji = fs_2, fs_jpim1 ! vector opt. 215 ta(ji,jj,jk) = ta(ji,jj,jk) + qsr_trd_hc_n(ji,jj,jk) / fse3t(ji,jj,jk) 216 END DO 217 END DO 218 END DO 219 ELSE 220 DO jk = 1, nksr 221 DO jj = 2, jpjm1 222 DO ji = fs_2, fs_jpim1 ! vector opt. 223 ta(ji,jj,jk) = ta(ji,jj,jk) + 0.5 * ( qsr_trd_hc_b(ji,jj,jk) + qsr_trd_hc_n(ji,jj,jk) ) / fse3t(ji,jj,jk) 224 END DO 225 END DO 226 END DO 210 227 ENDIF 211 228 … … 382 399 ! 383 400 DO jk = 1, nksr 384 etot3(:,:,jk) = ro0cpr * ( zea(:,:,jk) - zea(:,:,jk+1) ) / fse3t(:,:,jk)401 etot3(:,:,jk) = ro0cpr * ( zea(:,:,jk) - zea(:,:,jk+1) ) 385 402 END DO 386 403 etot3(:,:,nksr+1:jpk) = 0.e0 ! below 400m set to zero … … 404 421 zc0 = rn_abs * EXP( -fsdepw(ji,jj,jk )*zsi0r ) + (1.-rn_abs) * EXP( -fsdepw(ji,jj,jk )*zsi1r ) 405 422 zc1 = rn_abs * EXP( -fsdepw(ji,jj,jk+1)*zsi0r ) + (1.-rn_abs) * EXP( -fsdepw(ji,jj,jk+1)*zsi1r ) 406 etot3(ji,jj,jk) = ro0cpr * ( zc0 * tmask(ji,jj,jk) - zc1 * tmask(ji,jj,jk+1) ) / fse3t(ji,jj,jk)423 etot3(ji,jj,jk) = ro0cpr * ( zc0 * tmask(ji,jj,jk) - zc1 * tmask(ji,jj,jk+1) ) 407 424 END DO 408 425 END DO -
branches/DEV_r1837_MLF/NEMO/OPA_SRC/TRA/trasbc.F90
r1870 r1975 104 104 INTEGER, INTENT(in) :: kt ! ocean time-step index 105 105 !! 106 INTEGER :: ji, jj 107 REAL(wp) :: z ta, zsa, zsrau, zse3t! temporary scalars106 INTEGER :: ji, jj ! dummy loop indices 107 REAL(wp) :: zsrau, zse3t ! temporary scalars 108 108 !!---------------------------------------------------------------------- 109 109 … … 129 129 qsr(:,:) = 0.e0 ! qsr set to zero 130 130 IF( kt == nit000 ) THEN ! idem on before field at nit000 131 qns_b(:,:) = qns_b(:,:) + qsr_b(:,:)132 131 qsr_b(:,:) = 0.e0 133 132 ENDIF 134 133 ENDIF 135 136 ! Concentration dillution effect on (t,s)137 !! DO jj = 2, jpj138 !! DO ji = fs_2, fs_jpim1 ! vector opt.139 !!#if ! defined key_zco140 !! zse3t = 1. / fse3t(ji,jj,1)141 !!#endif142 !! IF( lk_vvl) THEN143 !! zta = ro0cpr * qns(ji,jj) * zse3t & ! temperature : heat flux144 !! & - emp(ji,jj) * zsrau * tn(ji,jj,1) * zse3t ! & cooling/heating effet of EMP flux145 !! zsa = 0.e0 ! No salinity concent./dilut. effect146 !! ELSE147 !! zta = ro0cpr * qns(ji,jj) * zse3t ! temperature : heat flux148 !! zsa = emps(ji,jj) * zsrau * sn(ji,jj,1) * zse3t ! salinity : concent./dilut. effect149 !! ENDIF150 !! ta(ji,jj,1) = ta(ji,jj,1) + zta ! add the trend to the general tracer trend151 !! sa(ji,jj,1) = sa(ji,jj,1) + zsa152 !! END DO153 !! END DO154 134 155 135 ! ! ---------------------- ! 156 136 IF( lk_vvl ) THEN ! Variable Volume case ! 157 137 ! ! ---------------------- ! 158 DO jj = 2, jpj 159 DO ji = fs_2, fs_jpim1 ! vector opt. 160 zse3t = 0.5 / fse3t(ji,jj,1) 161 ! ! temperature: heat flux + cooling/heating effet of EMP flux 162 ta(ji,jj,1) = ta(ji,jj,1) + ( ro0cpr * ( qns(ji,jj) + qns_b(ji,jj) ) & 163 & - zsrau * ( emp(ji,jj) + emp_b(ji,jj) ) * tn(ji,jj,1) ) * zse3t 164 ! ! salinity: salt flux 165 sa(ji,jj,1) = sa(ji,jj,1) + ( emps(ji,jj) + emps_b(ji,jj) ) * zse3t 166 167 !!gm BUG : in key_vvl emps must be modified to only include the salt flux due to with sea-ice freezing/melting 168 !!gm otherwise this flux will be missing ==> modification required in limsbc, limsbc_2 and CICE interface. 169 170 END DO 171 END DO 138 !!gm BUG : in key_vvl emps must be modified to only include the salt flux due to sea-ice freezing/melting 139 !!gm otherwise this flux will be missing ==> modification required in limsbc, limsbc_2 and CICE interface.s 140 IF ( neuler == 0 .AND. kt == nit000 ) THEN 141 DO jj = 2, jpj 142 DO ji = fs_2, fs_jpim1 ! vector opt. 143 ! temperature : heat flux + cooling/heating effet of EMP flux 144 sbc_trd_hc_n(ji,jj) = ro0cpr * qns(ji,jj) - zsrau * emp(ji,jj) * tn(ji,jj,1) 145 #if ! defined key_zco 146 zse3t = 1. / fse3t(ji,jj,1) 147 #endif 148 ta(ji,jj,1) = ta(ji,jj,1) + zse3t * sbc_trd_hc_n(ji,jj) 149 END DO 150 END DO 151 ELSE 152 DO jj = 2, jpj 153 DO ji = fs_2, fs_jpim1 ! vector opt. 154 ! temperature : heat flux + cooling/heating effet of EMP flux 155 sbc_trd_hc_n(ji,jj) = ro0cpr * qns(ji,jj) - zsrau * emp(ji,jj) * tn(ji,jj,1) 156 #if ! defined key_zco 157 zse3t = 1. / fse3t(ji,jj,1) 158 #endif 159 ta(ji,jj,1) = ta(ji,jj,1) + 0.5 * ( sbc_trd_hc_b(ji,jj) + sbc_trd_hc_n(ji,jj) ) * zse3t 160 END DO 161 END DO 162 ENDIF 172 163 ! ! ---------------------- ! 173 164 ELSE ! Constant Volume case ! 174 165 ! ! ---------------------- ! 175 DO jj = 2, jpj 176 DO ji = fs_2, fs_jpim1 ! vector opt. 177 zse3t = 0.5 / fse3t(ji,jj,1) 178 ! ! temperature: heat flux 179 ta(ji,jj,1) = ta(ji,jj,1) + ro0cpr * ( qns (ji,jj) + qns_b (ji,jj) ) * zse3t 180 ! ! salinity: salt flux + concent./dilut. effect (both in emps) 181 sa(ji,jj,1) = sa(ji,jj,1) + zsrau * ( emps(ji,jj) + emps_b(ji,jj) ) * sn(ji,jj,1) * zse3t 182 END DO 183 END DO 166 IF ( neuler == 0 .AND. kt == nit000 ) THEN 167 DO jj = 2, jpj 168 DO ji = fs_2, fs_jpim1 ! vector opt. 169 ! temperature : heat flux 170 sbc_trd_hc_n(ji,jj) = ro0cpr * qns(ji,jj) 171 ! salinity : salt flux + concent./dilut. effect (both in emps) 172 sbc_trd_sc_n(ji,jj) = zsrau * emps(ji,jj) * sn(ji,jj,1) 173 #if ! defined key_zco 174 zse3t = 1. / fse3t(ji,jj,1) 175 #endif 176 ta(ji,jj,1) = ta(ji,jj,1) + sbc_trd_hc_n(ji,jj) * zse3t 177 sa(ji,jj,1) = sa(ji,jj,1) + sbc_trd_sc_n(ji,jj) * zse3t 178 END DO 179 END DO 180 ELSE 181 DO jj = 2, jpj 182 DO ji = fs_2, fs_jpim1 ! vector opt. 183 ! temperature : heat flux 184 sbc_trd_hc_n(ji,jj) = ro0cpr * qns(ji,jj) 185 ! salinity : salt flux + concent./dilut. effect (both in emps) 186 sbc_trd_sc_n(ji,jj) = zsrau * emps(ji,jj) * sn(ji,jj,1) 187 #if ! defined key_zco 188 zse3t = 1. / fse3t(ji,jj,1) 189 #endif 190 ! temperature : heat flux 191 ta(ji,jj,1) = ta(ji,jj,1) + 0.5 * ( sbc_trd_hc_b(ji,jj) + sbc_trd_hc_n(ji,jj) ) * zse3t 192 sa(ji,jj,1) = sa(ji,jj,1) + 0.5 * ( sbc_trd_sc_b(ji,jj) + sbc_trd_sc_n(ji,jj) ) * zse3t 193 END DO 194 END DO 195 ENDIF 184 196 ! 185 197 ENDIF 186 198 187 IF( l_trdtra ) THEN ! save the sbc trends for diagnostic199 IF( l_trdtra ) THEN ! save the sbc trends for diagnostic 188 200 ztrdt(:,:,:) = ta(:,:,:) - ztrdt(:,:,:) 189 201 ztrds(:,:,:) = sa(:,:,:) - ztrds(:,:,:) 190 202 CALL trd_mod(ztrdt, ztrds, jptra_trd_nsr, 'TRA', kt) 191 203 ENDIF 192 ! ! control print204 ! ! control print 193 205 IF(ln_ctl) CALL prt_ctl( tab3d_1=ta, clinfo1=' sbc - Ta: ', mask1=tmask, & 194 206 & tab3d_2=sa, clinfo2= ' Sa: ', mask2=tmask, clinfo3='tra' ) -
branches/DEV_r1837_MLF/NEMO/OPA_SRC/step.F90
r1793 r1975 116 116 USE restart ! ocean restart (rst_wri routine) 117 117 USE prtctl ! Print control (prt_ctl routine) 118 ! - ML - 119 USE diatrb ! global tracer conservation (dia_trb routine) 118 120 119 121 #if defined key_agrif … … 317 319 318 320 CALL ssh_nxt( kstp ) ! sea surface height at next time step 321 CALL dia_trb( kstp ) ! - ML - global conservation diagnostics 319 322 320 323 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
Note: See TracChangeset
for help on using the changeset viewer.