- Timestamp:
- 2016-11-28T17:04:10+01:00 (7 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2016/dev_INGV_UKMO_2016/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90
r5407 r7351 23 23 !! lim_sbc_tau : update i- and j-stresses, and its modulus at the ocean surface 24 24 !!---------------------------------------------------------------------- 25 USE par_oce ! ocean parameters 26 USE phycst ! physical constants 27 USE dom_oce ! ocean domain 28 USE ice ! LIM sea-ice variables 29 USE sbc_ice ! Surface boundary condition: sea-ice fields 30 USE sbc_oce ! Surface boundary condition: ocean fields 31 USE sbccpl 32 USE oce , ONLY : sshn, sshb, snwice_mass, snwice_mass_b, snwice_fmass 33 USE albedo ! albedo parameters 34 USE lbclnk ! ocean lateral boundary condition - MPP exchanges 35 USE lib_mpp ! MPP library 36 USE wrk_nemo ! work arrays 37 USE in_out_manager ! I/O manager 38 USE prtctl ! Print control 39 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 40 USE traqsr ! add penetration of solar flux in the calculation of heat budget 41 USE iom 42 USE domvvl ! Variable volume 43 USE limctl 44 USE limcons 25 USE par_oce ! ocean parameters 26 USE oce , ONLY : sshn, sshb, snwice_mass, snwice_mass_b, snwice_fmass 27 USE phycst ! physical constants 28 USE dom_oce ! ocean domain 29 USE ice ! LIM sea-ice variables 30 USE sbc_ice ! Surface boundary condition: sea-ice fields 31 USE sbc_oce ! Surface boundary condition: ocean fields 32 USE sbccpl ! Surface boundary condition: coupled interface 33 USE albedo ! albedo parameters 34 USE traqsr ! add penetration of solar flux in the calculation of heat budget 35 USE domvvl ! Variable volume 36 USE limctl ! 37 USE limcons ! 38 ! 39 USE in_out_manager ! I/O manager 40 USE iom ! xIO server 41 USE lbclnk ! ocean lateral boundary condition - MPP exchanges 42 USE lib_mpp ! MPP library 43 USE wrk_nemo ! work arrays 44 USE prtctl ! Print control 45 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 45 46 46 47 IMPLICIT NONE 47 48 PRIVATE 48 49 49 PUBLIC lim_sbc_init ! called by sbc _lim_init50 PUBLIC lim_sbc_init ! called by sbcice_lim 50 51 PUBLIC lim_sbc_flx ! called by sbc_ice_lim 51 52 PUBLIC lim_sbc_tau ! called by sbc_ice_lim … … 57 58 !! * Substitutions 58 59 # include "vectopt_loop_substitute.h90" 59 # include "domzgr_substitute.h90"60 60 !!---------------------------------------------------------------------- 61 61 !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) … … 101 101 !! The ref should be Rousset et al., 2015 102 102 !!--------------------------------------------------------------------- 103 INTEGER, INTENT(in) :: kt ! number of iteration 104 INTEGER :: ji, jj, jl, jk ! dummy loop indices 105 REAL(wp) :: zqmass ! Heat flux associated with mass exchange ice->ocean (W.m-2) 106 REAL(wp) :: zqsr ! New solar flux received by the ocean 107 ! 108 REAL(wp), POINTER, DIMENSION(:,:,:) :: zalb_cs, zalb_os ! 2D/3D workspace 103 INTEGER, INTENT(in) :: kt ! number of iteration 104 ! 105 INTEGER :: ji, jj, jl, jk ! dummy loop indices 106 REAL(wp) :: zqmass ! Heat flux associated with mass exchange ice->ocean (W.m-2) 107 REAL(wp) :: zqsr ! New solar flux received by the ocean 108 REAL(wp), POINTER, DIMENSION(:,:,:) :: zalb_cs, zalb_os ! 3D workspace 109 REAL(wp), POINTER, DIMENSION(:,:) :: zalb ! 2D workspace 109 110 !!--------------------------------------------------------------------- 110 111 ! 111 112 ! make calls for heat fluxes before it is modified 113 ! pfrld is the lead fraction at the previous time step (actually between TRP and THD) 112 114 IF( iom_use('qsr_oce') ) CALL iom_put( "qsr_oce" , qsr_oce(:,:) * pfrld(:,:) ) ! solar flux at ocean surface 113 115 IF( iom_use('qns_oce') ) CALL iom_put( "qns_oce" , qns_oce(:,:) * pfrld(:,:) + qemp_oce(:,:) ) ! non-solar flux at ocean surface … … 118 120 IF( iom_use('qt_ice' ) ) CALL iom_put( "qt_ice" , SUM( ( qns_ice(:,:,:) + qsr_ice(:,:,:) ) & 119 121 & * a_i_b(:,:,:), dim=3 ) + qemp_ice(:,:) ) 120 IF( iom_use('qemp_oce' ) ) CALL iom_put( "qemp_oce" , qemp_oce(:,:) ) 121 IF( iom_use('qemp_ice' ) ) CALL iom_put( "qemp_ice" , qemp_ice(:,:) ) 122 123 ! pfrld is the lead fraction at the previous time step (actually between TRP and THD) 122 IF( iom_use('qemp_oce') ) CALL iom_put( "qemp_oce" , qemp_oce(:,:) ) 123 IF( iom_use('qemp_ice') ) CALL iom_put( "qemp_ice" , qemp_ice(:,:) ) 124 IF( iom_use('emp_oce' ) ) CALL iom_put( "emp_oce" , emp_oce(:,:) ) ! emp over ocean (taking into account the snow blown away from the ice) 125 IF( iom_use('emp_ice' ) ) CALL iom_put( "emp_ice" , emp_ice(:,:) ) ! emp over ice (taking into account the snow blown away from the ice) 126 127 ! albedo output 128 CALL wrk_alloc( jpi,jpj, zalb ) 129 130 zalb(:,:) = 0._wp 131 WHERE ( SUM( a_i_b, dim=3 ) <= epsi06 ) ; zalb(:,:) = 0.066_wp 132 ELSEWHERE ; zalb(:,:) = SUM( alb_ice * a_i_b, dim=3 ) / SUM( a_i_b, dim=3 ) 133 END WHERE 134 IF( iom_use('alb_ice' ) ) CALL iom_put( "alb_ice" , zalb(:,:) ) ! ice albedo output 135 136 zalb(:,:) = SUM( alb_ice * a_i_b, dim=3 ) + 0.066_wp * ( 1._wp - SUM( a_i_b, dim=3 ) ) 137 IF( iom_use('albedo' ) ) CALL iom_put( "albedo" , zalb(:,:) ) ! ice albedo output 138 139 CALL wrk_dealloc( jpi,jpj, zalb ) 140 124 141 DO jj = 1, jpj 125 142 DO ji = 1, jpi … … 140 157 hfx_out(ji,jj) = hfx_out(ji,jj) + zqmass + zqsr 141 158 142 ! Add the residual from heat diffusion equation (W.m-2) 143 !------------------------------------------------------- 144 hfx_out(ji,jj) = hfx_out(ji,jj) + hfx_err_dif(ji,jj) 159 ! Add the residual from heat diffusion equation and sublimation (W.m-2) 160 !---------------------------------------------------------------------- 161 hfx_out(ji,jj) = hfx_out(ji,jj) + hfx_err_dif(ji,jj) + & 162 & ( hfx_sub(ji,jj) - SUM( qevap_ice(ji,jj,:) * a_i_b(ji,jj,:) ) ) 145 163 146 164 ! New qsr and qns used to compute the oceanic heat flux at the next time step 147 !--------------------------------------------------- 165 !---------------------------------------------------------------------------- 148 166 qsr(ji,jj) = zqsr 149 167 qns(ji,jj) = hfx_out(ji,jj) - zqsr … … 165 183 166 184 ! mass flux at the ocean/ice interface 167 fmmflx(ji,jj) = - ( wfx_ice(ji,jj) + wfx_snw(ji,jj) ) * r1_rdtice ! F/M mass flux save at least for biogeochemical model 168 emp(ji,jj) = emp_oce(ji,jj) - wfx_ice(ji,jj) - wfx_snw(ji,jj) ! mass flux + F/M mass flux (always ice/ocean mass exchange) 169 185 fmmflx(ji,jj) = - ( wfx_ice(ji,jj) + wfx_snw(ji,jj) + wfx_err_sub(ji,jj) ) ! F/M mass flux save at least for biogeochemical model 186 emp(ji,jj) = emp_oce(ji,jj) - wfx_ice(ji,jj) - wfx_snw(ji,jj) - wfx_err_sub(ji,jj) ! mass flux + F/M mass flux (always ice/ocean mass exchange) 170 187 END DO 171 188 END DO … … 175 192 !------------------------------------------! 176 193 sfx(:,:) = sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) + sfx_opw(:,:) & 177 & + sfx_res(:,:) + sfx_dyn(:,:) + sfx_bri(:,:) 194 & + sfx_res(:,:) + sfx_dyn(:,:) + sfx_bri(:,:) + sfx_sub(:,:) 178 195 179 196 !-------------------------------------------------------------! … … 198 215 ! Snow/ice albedo (only if sent to coupler, useless in forced mode) ! 199 216 !------------------------------------------------------------------------! 200 CALL wrk_alloc( jpi, jpj, jpl,zalb_cs, zalb_os )217 CALL wrk_alloc( jpi,jpj,jpl, zalb_cs, zalb_os ) 201 218 CALL albedo_ice( t_su, ht_i, ht_s, zalb_cs, zalb_os ) ! cloud-sky and overcast-sky ice albedos 202 219 alb_ice(:,:,:) = ( 1. - cldf_ice ) * zalb_cs(:,:,:) + cldf_ice * zalb_os(:,:,:) 203 CALL wrk_dealloc( jpi, jpj, jpl,zalb_cs, zalb_os )220 CALL wrk_dealloc( jpi,jpj,jpl, zalb_cs, zalb_os ) 204 221 205 222 ! conservation test 206 IF( ln_limdiahsb ) CALL lim_cons_final( 'limsbc' )223 IF( ln_limdiahsb ) CALL lim_cons_final( 'limsbc' ) 207 224 208 225 ! control prints 209 226 IF( ln_icectl ) CALL lim_prt( kt, iiceprt, jiceprt, 3, ' - Final state lim_sbc - ' ) 210 227 ! 211 228 IF(ln_ctl) THEN 212 229 CALL prt_ctl( tab2d_1=qsr , clinfo1=' lim_sbc: qsr : ', tab2d_2=qns , clinfo2=' qns : ' ) … … 215 232 CALL prt_ctl( tab3d_1=tn_ice, clinfo1=' lim_sbc: tn_ice : ', kdim=jpl ) 216 233 ENDIF 217 234 ! 218 235 END SUBROUTINE lim_sbc_flx 219 236 … … 246 263 INTEGER , INTENT(in) :: kt ! ocean time-step index 247 264 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pu_oce, pv_oce ! surface ocean currents 248 ! !265 ! 249 266 INTEGER :: ji, jj ! dummy loop indices 250 267 REAL(wp) :: zat_u, zutau_ice, zu_t, zmodt ! local scalar … … 303 320 !! ** input : Namelist namicedia 304 321 !!------------------------------------------------------------------- 305 INTEGER :: ji, jj, jk ! dummy loop indices 306 REAL(wp) :: zcoefu, zcoefv, zcoeff ! local scalar 322 INTEGER :: ji, jj, jk ! dummy loop indices 323 REAL(wp) :: zcoefu, zcoefv, zcoeff ! local scalar 324 !!------------------------------------------------------------------- 325 ! 307 326 IF(lwp) WRITE(numout,*) 308 327 IF(lwp) WRITE(numout,*) 'lim_sbc_init : LIM-3 sea-ice - surface boundary condition' … … 335 354 sshn(:,:) = sshn(:,:) - snwice_mass(:,:) * r1_rau0 336 355 sshb(:,:) = sshb(:,:) - snwice_mass(:,:) * r1_rau0 337 #if defined key_vvl 338 ! key_vvl necessary? clem: yes for compilation purpose 339 DO jk = 1,jpkm1 ! adjust initial vertical scale factors 340 fse3t_n(:,:,jk) = e3t_0(:,:,jk)*( 1._wp + sshn(:,:)*tmask(:,:,1)/(ht_0(:,:) + 1.0 - tmask(:,:,1)) ) 341 fse3t_b(:,:,jk) = e3t_0(:,:,jk)*( 1._wp + sshb(:,:)*tmask(:,:,1)/(ht_0(:,:) + 1.0 - tmask(:,:,1)) ) 342 ENDDO 343 fse3t_a(:,:,:) = fse3t_b(:,:,:) 344 ! Reconstruction of all vertical scale factors at now and before time 345 ! steps 346 ! ============================================================================= 347 ! Horizontal scale factor interpolations 348 ! -------------------------------------- 349 CALL dom_vvl_interpol( fse3t_b(:,:,:), fse3u_b(:,:,:), 'U' ) 350 CALL dom_vvl_interpol( fse3t_b(:,:,:), fse3v_b(:,:,:), 'V' ) 351 CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3u_n(:,:,:), 'U' ) 352 CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3v_n(:,:,:), 'V' ) 353 CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3f_n(:,:,:), 'F' ) 354 ! Vertical scale factor interpolations 355 ! ------------------------------------ 356 CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3w_n (:,:,:), 'W' ) 357 CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3uw_n(:,:,:), 'UW' ) 358 CALL dom_vvl_interpol( fse3v_n(:,:,:), fse3vw_n(:,:,:), 'VW' ) 359 CALL dom_vvl_interpol( fse3u_b(:,:,:), fse3uw_b(:,:,:), 'UW' ) 360 CALL dom_vvl_interpol( fse3v_b(:,:,:), fse3vw_b(:,:,:), 'VW' ) 361 ! t- and w- points depth 362 ! ---------------------- 363 fsdept_n(:,:,1) = 0.5_wp * fse3w_n(:,:,1) 364 fsdepw_n(:,:,1) = 0.0_wp 365 fsde3w_n(:,:,1) = fsdept_n(:,:,1) - sshn(:,:) 366 DO jk = 2, jpk 367 fsdept_n(:,:,jk) = fsdept_n(:,:,jk-1) + fse3w_n(:,:,jk) 368 fsdepw_n(:,:,jk) = fsdepw_n(:,:,jk-1) + fse3t_n(:,:,jk-1) 369 fsde3w_n(:,:,jk) = fsdept_n(:,:,jk ) - sshn (:,:) 370 END DO 371 #endif 356 357 !!gm I really don't like this staff here... Find a way to put that elsewhere or differently 358 !!gm 359 IF( .NOT.ln_linssh ) THEN 360 DO jk = 1,jpkm1 ! adjust initial vertical scale factors 361 e3t_n(:,:,jk) = e3t_0(:,:,jk)*( 1._wp + sshn(:,:)*tmask(:,:,1)/(ht_0(:,:) + 1.0 - tmask(:,:,1)) ) 362 e3t_b(:,:,jk) = e3t_0(:,:,jk)*( 1._wp + sshb(:,:)*tmask(:,:,1)/(ht_0(:,:) + 1.0 - tmask(:,:,1)) ) 363 END DO 364 e3t_a(:,:,:) = e3t_b(:,:,:) 365 ! Reconstruction of all vertical scale factors at now and before time-steps 366 ! ========================================================================= 367 ! Horizontal scale factor interpolations 368 ! -------------------------------------- 369 CALL dom_vvl_interpol( e3t_b(:,:,:), e3u_b(:,:,:), 'U' ) 370 CALL dom_vvl_interpol( e3t_b(:,:,:), e3v_b(:,:,:), 'V' ) 371 CALL dom_vvl_interpol( e3t_n(:,:,:), e3u_n(:,:,:), 'U' ) 372 CALL dom_vvl_interpol( e3t_n(:,:,:), e3v_n(:,:,:), 'V' ) 373 CALL dom_vvl_interpol( e3u_n(:,:,:), e3f_n(:,:,:), 'F' ) 374 ! Vertical scale factor interpolations 375 ! ------------------------------------ 376 CALL dom_vvl_interpol( e3t_n(:,:,:), e3w_n (:,:,:), 'W' ) 377 CALL dom_vvl_interpol( e3u_n(:,:,:), e3uw_n(:,:,:), 'UW' ) 378 CALL dom_vvl_interpol( e3v_n(:,:,:), e3vw_n(:,:,:), 'VW' ) 379 CALL dom_vvl_interpol( e3u_b(:,:,:), e3uw_b(:,:,:), 'UW' ) 380 CALL dom_vvl_interpol( e3v_b(:,:,:), e3vw_b(:,:,:), 'VW' ) 381 ! t- and w- points depth 382 ! ---------------------- 383 !!gm not sure of that.... 384 gdept_n(:,:,1) = 0.5_wp * e3w_n(:,:,1) 385 gdepw_n(:,:,1) = 0.0_wp 386 gde3w_n(:,:,1) = gdept_n(:,:,1) - sshn(:,:) 387 DO jk = 2, jpk 388 gdept_n(:,:,jk) = gdept_n(:,:,jk-1) + e3w_n(:,:,jk) 389 gdepw_n(:,:,jk) = gdepw_n(:,:,jk-1) + e3t_n(:,:,jk-1) 390 gde3w_n(:,:,jk) = gdept_n(:,:,jk ) - sshn (:,:) 391 END DO 392 ENDIF 372 393 ENDIF 373 394 ENDIF ! .NOT. ln_rstart 374 395 ! 375 376 396 END SUBROUTINE lim_sbc_init 377 397
Note: See TracChangeset
for help on using the changeset viewer.
