Changeset 12583 for NEMO/branches/2020
- Timestamp:
- 2020-03-21T15:40:52+01:00 (4 years ago)
- Location:
- NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src
- Files:
-
- 1 added
- 6 edited
-
ICE/iceistate.F90 (modified) (18 diffs)
-
OCE/DOM/domqe.F90 (modified) (3 diffs)
-
OCE/DOM/domzgr_substitute.h90 (added)
-
OCE/DYN/dynatfLF.F90 (modified) (1 diff)
-
OCE/SBC/sbcice_cice.F90 (modified) (36 diffs)
-
OCE/TRA/traatfLF.F90 (modified) (1 diff)
-
OCE/steplf.F90 (modified) (3 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/ICE/iceistate.F90
r12377 r12583 18 18 USE oce ! dynamics and tracers variables 19 19 USE dom_oce ! ocean domain 20 USE sbc_oce , ONLY : sst_m, sss_m, ln_ice_embd 20 USE sbc_oce , ONLY : sst_m, sss_m, ln_ice_embd 21 21 USE sbc_ice , ONLY : tn_ice, snwice_mass, snwice_mass_b 22 22 USE eosbn2 ! equation of state … … 60 60 INTEGER , PARAMETER :: jp_hpd = 9 ! index of pnd depth (m) 61 61 TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: si ! structure of input fields (file informations, fields read) 62 ! 62 ! 63 63 !! * Substitutions 64 64 # include "do_loop_substitute.h90" … … 77 77 !! 78 78 !! ** Method : This routine will put some ice where ocean 79 !! is at the freezing point, then fill in ice 80 !! state variables using prescribed initial 81 !! values in the namelist 79 !! is at the freezing point, then fill in ice 80 !! state variables using prescribed initial 81 !! values in the namelist 82 82 !! 83 83 !! ** Steps : 1) Set initial surface and basal temperatures … … 91 91 !! where there is no ice 92 92 !!-------------------------------------------------------------------- 93 INTEGER, INTENT(in) :: kt ! time step 93 INTEGER, INTENT(in) :: kt ! time step 94 94 INTEGER, INTENT(in) :: Kbb, Kmm, Kaa ! ocean time level indices 95 95 ! … … 117 117 ! basal temperature (considered at freezing point) [Kelvin] 118 118 CALL eos_fzp( sss_m(:,:), t_bo(:,:) ) 119 t_bo(:,:) = ( t_bo(:,:) + rt0 ) * tmask(:,:,1) 119 t_bo(:,:) = ( t_bo(:,:) + rt0 ) * tmask(:,:,1) 120 120 ! 121 121 ! surface temperature and conductivity … … 142 142 e_i (:,:,:,:) = 0._wp 143 143 e_s (:,:,:,:) = 0._wp 144 144 145 145 ! general fields 146 146 a_i (:,:,:) = 0._wp … … 215 215 IF( TRIM(si(jp_apd)%clrootname) == 'NOT USED' ) & 216 216 & si(jp_apd)%fnow(:,:,1) = ( rn_apd_ini_n * zswitch + rn_apd_ini_s * (1._wp - zswitch) ) * tmask(:,:,1) & ! rn_apd = pond fraction => rn_apnd * a_i = pond conc. 217 & * si(jp_ati)%fnow(:,:,1) 217 & * si(jp_ati)%fnow(:,:,1) 218 218 zapnd_ini(:,:) = si(jp_apd)%fnow(:,:,1) 219 219 ! … … 224 224 ! 225 225 ! change the switch for the following 226 WHERE( zat_i_ini(:,:) > 0._wp ) ; zswitch(:,:) = tmask(:,:,1) 226 WHERE( zat_i_ini(:,:) > 0._wp ) ; zswitch(:,:) = tmask(:,:,1) 227 227 ELSEWHERE ; zswitch(:,:) = 0._wp 228 228 END WHERE … … 231 231 ! !---------------! 232 232 ! no ice if (sst - Tfreez) >= thresold 233 WHERE( ( sst_m(:,:) - (t_bo(:,:) - rt0) ) * tmask(:,:,1) >= rn_thres_sst ) ; zswitch(:,:) = 0._wp 233 WHERE( ( sst_m(:,:) - (t_bo(:,:) - rt0) ) * tmask(:,:,1) >= rn_thres_sst ) ; zswitch(:,:) = 0._wp 234 234 ELSEWHERE ; zswitch(:,:) = tmask(:,:,1) 235 235 END WHERE … … 244 244 zt_su_ini(:,:) = rn_tsu_ini_n * zswitch(:,:) 245 245 ztm_s_ini(:,:) = rn_tms_ini_n * zswitch(:,:) 246 zapnd_ini(:,:) = rn_apd_ini_n * zswitch(:,:) * zat_i_ini(:,:) ! rn_apd = pond fraction => rn_apd * a_i = pond conc. 246 zapnd_ini(:,:) = rn_apd_ini_n * zswitch(:,:) * zat_i_ini(:,:) ! rn_apd = pond fraction => rn_apd * a_i = pond conc. 247 247 zhpnd_ini(:,:) = rn_hpd_ini_n * zswitch(:,:) 248 248 ELSEWHERE … … 265 265 zhpnd_ini(:,:) = 0._wp 266 266 ENDIF 267 267 268 268 !-------------! 269 269 ! fill fields ! … … 292 292 ALLOCATE( zhi_2d(npti,jpl), zhs_2d(npti,jpl), zai_2d (npti,jpl), & 293 293 & zti_2d(npti,jpl), zts_2d(npti,jpl), ztsu_2d(npti,jpl), zsi_2d(npti,jpl), zaip_2d(npti,jpl), zhip_2d(npti,jpl) ) 294 294 295 295 ! distribute 1-cat into jpl-cat: (jpi*jpj) -> (jpi*jpj,jpl) 296 296 CALL ice_var_itd( h_i_1d(1:npti) , h_s_1d(1:npti) , at_i_1d(1:npti), & … … 338 338 DO jl = 1, jpl 339 339 DO_3D_11_11( 1, nlay_i ) 340 t_i (ji,jj,jk,jl) = zti_3d(ji,jj,jl) 340 t_i (ji,jj,jk,jl) = zti_3d(ji,jj,jl) 341 341 ztmelts = - rTmlt * sz_i(ji,jj,jk,jl) + rt0 ! melting temperature in K 342 342 e_i(ji,jj,jk,jl) = zswitch(ji,jj) * v_i(ji,jj,jl) * r1_nlay_i * & … … 354 354 END WHERE 355 355 v_ip(:,:,:) = h_ip(:,:,:) * a_ip(:,:,:) 356 356 357 357 ! specific temperatures for coupled runs 358 358 tn_ice(:,:,:) = t_su(:,:,:) … … 374 374 ssh(:,:,Kbb) = ssh(:,:,Kbb) - snwice_mass(:,:) * r1_rau0 375 375 ! 376 IF( .NOT.ln_linssh ) THEN 377 ! 378 WHERE( ht_0(:,:) > 0 ) ; z2d(:,:) = 1._wp + ssh(:,:,Kmm)*tmask(:,:,1) / ht_0(:,:) 379 ELSEWHERE ; z2d(:,:) = 1._wp ; END WHERE 380 ! 381 DO jk = 1,jpkm1 ! adjust initial vertical scale factors 382 e3t(:,:,jk,Kmm) = e3t_0(:,:,jk) * z2d(:,:) 383 e3t(:,:,jk,Kbb) = e3t(:,:,jk,Kmm) 384 e3t(:,:,jk,Kaa) = e3t(:,:,jk,Kmm) 385 END DO 386 ! 387 ! Reconstruction of all vertical scale factors at now and before time-steps 388 ! ========================================================================= 389 ! Horizontal scale factor interpolations 390 ! -------------------------------------- 391 CALL dom_vvl_interpol( e3t(:,:,:,Kbb), e3u(:,:,:,Kbb), 'U' ) 392 CALL dom_vvl_interpol( e3t(:,:,:,Kbb), e3v(:,:,:,Kbb), 'V' ) 393 CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3u(:,:,:,Kmm), 'U' ) 394 CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3v(:,:,:,Kmm), 'V' ) 395 CALL dom_vvl_interpol( e3u(:,:,:,Kmm), e3f(:,:,:), 'F' ) 396 ! Vertical scale factor interpolations 397 ! ------------------------------------ 398 CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3w (:,:,:,Kmm), 'W' ) 399 CALL dom_vvl_interpol( e3u(:,:,:,Kmm), e3uw(:,:,:,Kmm), 'UW' ) 400 CALL dom_vvl_interpol( e3v(:,:,:,Kmm), e3vw(:,:,:,Kmm), 'VW' ) 401 CALL dom_vvl_interpol( e3u(:,:,:,Kbb), e3uw(:,:,:,Kbb), 'UW' ) 402 CALL dom_vvl_interpol( e3v(:,:,:,Kbb), e3vw(:,:,:,Kbb), 'VW' ) 403 ! t- and w- points depth 404 ! ---------------------- 405 !!gm not sure of that.... 406 gdept(:,:,1,Kmm) = 0.5_wp * e3w(:,:,1,Kmm) 407 gdepw(:,:,1,Kmm) = 0.0_wp 408 gde3w(:,:,1) = gdept(:,:,1,Kmm) - ssh(:,:,Kmm) 409 DO jk = 2, jpk 410 gdept(:,:,jk,Kmm) = gdept(:,:,jk-1,Kmm) + e3w(:,:,jk ,Kmm) 411 gdepw(:,:,jk,Kmm) = gdepw(:,:,jk-1,Kmm) + e3t(:,:,jk-1,Kmm) 412 gde3w(:,:,jk) = gdept(:,:,jk ,Kmm) - ssh (:,:,Kmm) 413 END DO 414 ENDIF 376 IF( .NOT.ln_linssh ) CALL dom_vvl_zgr( Kbb, Kmm, Kaa ) ! interpolation scale factor, depth and water column 377 ! !!st 378 ! IF( .NOT.ln_linssh ) THEN 379 ! ! 380 ! WHERE( ht_0(:,:) > 0 ) ; z2d(:,:) = 1._wp + ssh(:,:,Kmm)*tmask(:,:,1) / ht_0(:,:) 381 ! ELSEWHERE ; z2d(:,:) = 1._wp ; END WHERE 382 ! ! 383 ! DO jk = 1,jpkm1 ! adjust initial vertical scale factors 384 ! e3t(:,:,jk,Kmm) = e3t_0(:,:,jk) * z2d(:,:) 385 ! e3t(:,:,jk,Kbb) = e3t(:,:,jk,Kmm) 386 ! e3t(:,:,jk,Kaa) = e3t(:,:,jk,Kmm) 387 ! END DO 388 ! ! 389 ! ! Reconstruction of all vertical scale factors at now and before time-steps 390 ! ! ========================================================================= 391 ! ! Horizontal scale factor interpolations 392 ! ! -------------------------------------- 393 ! CALL dom_vvl_interpol( e3t(:,:,:,Kbb), e3u(:,:,:,Kbb), 'U' ) 394 ! CALL dom_vvl_interpol( e3t(:,:,:,Kbb), e3v(:,:,:,Kbb), 'V' ) 395 ! CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3u(:,:,:,Kmm), 'U' ) 396 ! CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3v(:,:,:,Kmm), 'V' ) 397 ! CALL dom_vvl_interpol( e3u(:,:,:,Kmm), e3f(:,:,:), 'F' ) 398 ! ! Vertical scale factor interpolations 399 ! ! ------------------------------------ 400 ! CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3w (:,:,:,Kmm), 'W' ) 401 ! CALL dom_vvl_interpol( e3u(:,:,:,Kmm), e3uw(:,:,:,Kmm), 'UW' ) 402 ! CALL dom_vvl_interpol( e3v(:,:,:,Kmm), e3vw(:,:,:,Kmm), 'VW' ) 403 ! CALL dom_vvl_interpol( e3u(:,:,:,Kbb), e3uw(:,:,:,Kbb), 'UW' ) 404 ! CALL dom_vvl_interpol( e3v(:,:,:,Kbb), e3vw(:,:,:,Kbb), 'VW' ) 405 ! ! t- and w- points depth 406 ! ! ---------------------- 407 ! !!gm not sure of that.... 408 ! gdept(:,:,1,Kmm) = 0.5_wp * e3w(:,:,1,Kmm) 409 ! gdepw(:,:,1,Kmm) = 0.0_wp 410 ! gde3w(:,:,1) = gdept(:,:,1,Kmm) - ssh(:,:,Kmm) 411 ! DO jk = 2, jpk 412 ! gdept(:,:,jk,Kmm) = gdept(:,:,jk-1,Kmm) + e3w(:,:,jk ,Kmm) 413 ! gdepw(:,:,jk,Kmm) = gdepw(:,:,jk-1,Kmm) + e3t(:,:,jk-1,Kmm) 414 ! gde3w(:,:,jk) = gdept(:,:,jk ,Kmm) - ssh (:,:,Kmm) 415 ! END DO 416 ! ENDIF 415 417 ENDIF 416 418 417 419 !------------------------------------ 418 420 ! 4) store fields at before time-step … … 429 431 v_ice_b(:,:) = v_ice(:,:) 430 432 ! total concentration is needed for Lupkes parameterizations 431 at_i_b (:,:) = at_i (:,:) 433 at_i_b (:,:) = at_i (:,:) 432 434 433 435 !!clem: output of initial state should be written here but it is impossible because … … 441 443 !!------------------------------------------------------------------- 442 444 !! *** ROUTINE ice_istate_init *** 443 !! 444 !! ** Purpose : Definition of initial state of the ice 445 !! 446 !! ** Method : Read the namini namelist and check the parameter 445 !! 446 !! ** Purpose : Definition of initial state of the ice 447 !! 448 !! ** Method : Read the namini namelist and check the parameter 447 449 !! values called at the first timestep (nit000) 448 450 !! … … 485 487 WRITE(numout,*) ' max ocean temp. above Tfreeze with initial ice rn_thres_sst = ', rn_thres_sst 486 488 IF( ln_iceini .AND. .NOT.ln_iceini_file ) THEN 487 WRITE(numout,*) ' initial snw thickness in the north-south rn_hts_ini = ', rn_hts_ini_n,rn_hts_ini_s 489 WRITE(numout,*) ' initial snw thickness in the north-south rn_hts_ini = ', rn_hts_ini_n,rn_hts_ini_s 488 490 WRITE(numout,*) ' initial ice thickness in the north-south rn_hti_ini = ', rn_hti_ini_n,rn_hti_ini_s 489 491 WRITE(numout,*) ' initial ice concentr in the north-south rn_ati_ini = ', rn_ati_ini_n,rn_ati_ini_s -
NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DOM/domqe.F90
r12581 r12583 13 13 14 14 !!---------------------------------------------------------------------- 15 !! dom_ vvl_init : define initial vertical scale factors, depths and column thickness16 !! dom_ vvl_sf_nxt : Compute next vertical scale factors17 !! dom_ vvl_sf_update : Swap vertical scale factors and update the vertical grid18 !! dom_ vvl_interpol : Interpolate vertical scale factors from one grid point to another19 !! dom_ vvl_r3c : Compute ssh/h_0 ratioat t-, u-, v-, and optionally f-points20 !! dom_ vvl_rst : read/write restart file21 !! dom_ vvl_ctl : Check the vvl options15 !! dom_qe_init : define initial vertical scale factors, depths and column thickness 16 !! dom_qe_sf_nxt : Compute next vertical scale factors 17 !! dom_qe_sf_update : Swap vertical scale factors and update the vertical grid 18 !! dom_qe_interpol : Interpolate vertical scale factors from one grid point to another 19 !! dom_qe_r3c : Compute ssh/h_0 ratioat t-, u-, v-, and optionally f-points 20 !! dom_qe_rst : read/write restart file 21 !! dom_qe_ctl : Check the vvl options 22 22 !!---------------------------------------------------------------------- 23 23 USE oce ! ocean dynamics and tracers … … 374 374 + MAX( 0._wp , gdepw_0(:,:,jk)-risfdep(:,:) ) * ( 1._wp + r3t(:,:,Kmm) ) 375 375 gde3w(:,:,jk) = gdept(:,:,jk,Kmm) - ssh(:,:,Kmm) 376 gdept(:,:,jk,Kbb) = MIN( risfdep(:,:) , gdept_0(:,:,jk) ) & 377 + MAX( 0._wp , gdept_0(:,:,jk)-risfdep(:,:) ) * ( 1._wp + r3t(:,:,Kbb) ) 378 gdepw(:,:,jk,Kbb) = MIN( risfdep(:,:) , gdepw_0(:,:,jk) ) & 379 + MAX( 0._wp , gdepw_0(:,:,jk)-risfdep(:,:) ) * ( 1._wp + r3t(:,:,Kbb) ) 376 380 END DO 377 381 ! … … 383 387 gdepw(:,:,jk,Kmm) = gdepw_0(:,:,jk) * ( 1._wp + r3t(:,:,Kmm) ) 384 388 gde3w(:,:,jk) = gdept (:,:,jk,Kmm) - ssh(:,:,Kmm) 389 gdept(:,:,jk,Kbb) = gdept_0(:,:,jk) * ( 1._wp + r3t(:,:,Kbb) ) 390 gdepw(:,:,jk,Kbb) = gdepw_0(:,:,jk) * ( 1._wp + r3t(:,:,Kbb) ) 385 391 END DO 386 392 ! -
NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DYN/dynatfLF.F90
r12581 r12583 116 116 ENDIF 117 117 118 IF ( ln_dynspg_ts ) THEN119 ! Ensure below that barotropic velocities match time splitting estimate120 ! Compute actual transport and replace it with ts estimate at "after" time step121 zue(:,:) = pe3u(:,:,1,Kaa) * puu(:,:,1,Kaa) * umask(:,:,1)122 zve(:,:) = pe3v(:,:,1,Kaa) * pvv(:,:,1,Kaa) * vmask(:,:,1)123 DO jk = 2, jpkm1124 zue(:,:) = zue(:,:) + pe3u(:,:,jk,Kaa) * puu(:,:,jk,Kaa) * umask(:,:,jk)125 zve(:,:) = zve(:,:) + pe3v(:,:,jk,Kaa) * pvv(:,:,jk,Kaa) * vmask(:,:,jk)126 END DO127 DO jk = 1, jpkm1128 puu(:,:,jk,Kaa) = ( puu(:,:,jk,Kaa) - zue(:,:) * r1_hu(:,:,Kaa) + uu_b(:,:,Kaa) ) * umask(:,:,jk)129 pvv(:,:,jk,Kaa) = ( pvv(:,:,jk,Kaa) - zve(:,:) * r1_hv(:,:,Kaa) + vv_b(:,:,Kaa) ) * vmask(:,:,jk)130 END DO131 !132 IF( .NOT.ln_bt_fw ) THEN133 ! Remove advective velocity from "now velocities"134 ! prior to asselin filtering135 ! In the forward case, this is done below after asselin filtering136 ! so that asselin contribution is removed at the same time137 DO jk = 1, jpkm1138 puu(:,:,jk,Kmm) = ( puu(:,:,jk,Kmm) - un_adv(:,:)*r1_hu(:,:,Kmm) + uu_b(:,:,Kmm) )*umask(:,:,jk)139 pvv(:,:,jk,Kmm) = ( pvv(:,:,jk,Kmm) - vn_adv(:,:)*r1_hv(:,:,Kmm) + vv_b(:,:,Kmm) )*vmask(:,:,jk)140 END DO141 ENDIF142 ENDIF143 144 ! Update after velocity on domain lateral boundaries145 ! --------------------------------------------------146 # if defined key_agrif147 CALL Agrif_dyn( kt ) !* AGRIF zoom boundaries148 # endif149 !150 CALL lbc_lnk_multi( 'dynatfLF', puu(:,:,:,Kaa), 'U', -1., pvv(:,:,:,Kaa), 'V', -1. ) !* local domain boundaries151 !152 ! !* BDY open boundaries153 IF( ln_bdy .AND. ln_dynspg_exp ) CALL bdy_dyn( kt, Kbb, puu, pvv, Kaa )154 IF( ln_bdy .AND. ln_dynspg_ts ) CALL bdy_dyn( kt, Kbb, puu, pvv, Kaa, dyn3d_only=.true. )118 ! IF ( ln_dynspg_ts ) THEN 119 ! ! Ensure below that barotropic velocities match time splitting estimate 120 ! ! Compute actual transport and replace it with ts estimate at "after" time step 121 ! zue(:,:) = pe3u(:,:,1,Kaa) * puu(:,:,1,Kaa) * umask(:,:,1) 122 ! zve(:,:) = pe3v(:,:,1,Kaa) * pvv(:,:,1,Kaa) * vmask(:,:,1) 123 ! DO jk = 2, jpkm1 124 ! zue(:,:) = zue(:,:) + pe3u(:,:,jk,Kaa) * puu(:,:,jk,Kaa) * umask(:,:,jk) 125 ! zve(:,:) = zve(:,:) + pe3v(:,:,jk,Kaa) * pvv(:,:,jk,Kaa) * vmask(:,:,jk) 126 ! END DO 127 ! DO jk = 1, jpkm1 128 ! puu(:,:,jk,Kaa) = ( puu(:,:,jk,Kaa) - zue(:,:) * r1_hu(:,:,Kaa) + uu_b(:,:,Kaa) ) * umask(:,:,jk) 129 ! pvv(:,:,jk,Kaa) = ( pvv(:,:,jk,Kaa) - zve(:,:) * r1_hv(:,:,Kaa) + vv_b(:,:,Kaa) ) * vmask(:,:,jk) 130 ! END DO 131 ! ! 132 ! IF( .NOT.ln_bt_fw ) THEN 133 ! ! Remove advective velocity from "now velocities" 134 ! ! prior to asselin filtering 135 ! ! In the forward case, this is done below after asselin filtering 136 ! ! so that asselin contribution is removed at the same time 137 ! DO jk = 1, jpkm1 138 ! puu(:,:,jk,Kmm) = ( puu(:,:,jk,Kmm) - un_adv(:,:)*r1_hu(:,:,Kmm) + uu_b(:,:,Kmm) )*umask(:,:,jk) 139 ! pvv(:,:,jk,Kmm) = ( pvv(:,:,jk,Kmm) - vn_adv(:,:)*r1_hv(:,:,Kmm) + vv_b(:,:,Kmm) )*vmask(:,:,jk) 140 ! END DO 141 ! ENDIF 142 ! ENDIF 143 ! 144 ! ! Update after velocity on domain lateral boundaries 145 ! ! -------------------------------------------------- 146 ! # if defined key_agrif 147 ! CALL Agrif_dyn( kt ) !* AGRIF zoom boundaries 148 ! # endif 149 ! ! 150 ! CALL lbc_lnk_multi( 'dynatfLF', puu(:,:,:,Kaa), 'U', -1., pvv(:,:,:,Kaa), 'V', -1. ) !* local domain boundaries 151 ! ! 152 ! ! !* BDY open boundaries 153 ! IF( ln_bdy .AND. ln_dynspg_exp ) CALL bdy_dyn( kt, Kbb, puu, pvv, Kaa ) 154 ! IF( ln_bdy .AND. ln_dynspg_ts ) CALL bdy_dyn( kt, Kbb, puu, pvv, Kaa, dyn3d_only=.true. ) 155 155 156 156 !!$ Do we need a call to bdy_vol here?? -
NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/SBC/sbcice_cice.F90
r12377 r12583 36 36 # if defined key_cice4 37 37 USE ice_flux, only: strax,stray,strocnx,strocny,frain,fsnow, & 38 strocnxT,strocnyT, & 38 strocnxT,strocnyT, & 39 39 sst,sss,uocn,vocn,ss_tltx,ss_tlty,fsalt_gbm, & 40 40 fresh_gbm,fhocn_gbm,fswthru_gbm,frzmlt, & … … 45 45 #else 46 46 USE ice_flux, only: strax,stray,strocnx,strocny,frain,fsnow, & 47 strocnxT,strocnyT, & 47 strocnxT,strocnyT, & 48 48 sst,sss,uocn,vocn,ss_tltx,ss_tlty,fsalt_ai, & 49 49 fresh_ai,fhocn_ai,fswthru_ai,frzmlt, & … … 70 70 INTEGER :: jj_off 71 71 72 INTEGER , PARAMETER :: jpfld = 13 ! maximum number of files to read 72 INTEGER , PARAMETER :: jpfld = 13 ! maximum number of files to read 73 73 INTEGER , PARAMETER :: jp_snow = 1 ! index of snow file 74 74 INTEGER , PARAMETER :: jp_rain = 2 ! index of rain file … … 109 109 !!--------------------------------------------------------------------- 110 110 !! *** ROUTINE sbc_ice_cice *** 111 !! 112 !! ** Purpose : update the ocean surface boundary condition via the 113 !! CICE Sea Ice Model time stepping 114 !! 115 !! ** Method : - Get any extra forcing fields for CICE 111 !! 112 !! ** Purpose : update the ocean surface boundary condition via the 113 !! CICE Sea Ice Model time stepping 114 !! 115 !! ** Method : - Get any extra forcing fields for CICE 116 116 !! - Prepare forcing fields 117 117 !! - CICE model time stepping 118 !! - call the routine that computes mass and 118 !! - call the routine that computes mass and 119 119 !! heat fluxes at the ice/ocean interface 120 120 !! … … 171 171 ! there is no restart file. 172 172 ! Values from a CICE restart file would overwrite this 173 IF( .NOT. ln_rstart ) THEN 174 CALL nemo2cice( ts(:,:,1,jp_tem,Kmm) , sst , 'T' , 1.) 175 ENDIF 173 IF( .NOT. ln_rstart ) THEN 174 CALL nemo2cice( ts(:,:,1,jp_tem,Kmm) , sst , 'T' , 1.) 175 ENDIF 176 176 #endif 177 177 … … 233 233 !!gm This should be put elsewhere.... (same remark for limsbc) 234 234 !!gm especially here it is assumed zstar coordinate, but it can be ztilde.... 235 IF( .NOT.ln_linssh ) THEN 236 ! 237 DO jk = 1,jpkm1 ! adjust initial vertical scale factors 238 e3t(:,:,jk,Kmm) = e3t_0(:,:,jk)*( 1._wp + ssh(:,:,Kmm)*tmask(:,:,1)/(ht_0(:,:) + 1.0 - tmask(:,:,1)) ) 239 e3t(:,:,jk,Kbb) = e3t_0(:,:,jk)*( 1._wp + ssh(:,:,Kbb)*tmask(:,:,1)/(ht_0(:,:) + 1.0 - tmask(:,:,1)) ) 240 ENDDO 241 e3t(:,:,:,Krhs) = e3t(:,:,:,Kbb) 242 ! Reconstruction of all vertical scale factors at now and before time-steps 243 ! ============================================================================= 244 ! Horizontal scale factor interpolations 245 ! -------------------------------------- 246 CALL dom_vvl_interpol( e3t(:,:,:,Kbb), e3u(:,:,:,Kbb), 'U' ) 247 CALL dom_vvl_interpol( e3t(:,:,:,Kbb), e3v(:,:,:,Kbb), 'V' ) 248 CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3u(:,:,:,Kmm), 'U' ) 249 CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3v(:,:,:,Kmm), 'V' ) 250 CALL dom_vvl_interpol( e3u(:,:,:,Kmm), e3f(:,:,:), 'F' ) 251 ! Vertical scale factor interpolations 252 ! ------------------------------------ 253 CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3w (:,:,:,Kmm), 'W' ) 254 CALL dom_vvl_interpol( e3u(:,:,:,Kmm), e3uw(:,:,:,Kmm), 'UW' ) 255 CALL dom_vvl_interpol( e3v(:,:,:,Kmm), e3vw(:,:,:,Kmm), 'VW' ) 256 CALL dom_vvl_interpol( e3u(:,:,:,Kbb), e3uw(:,:,:,Kbb), 'UW' ) 257 CALL dom_vvl_interpol( e3v(:,:,:,Kbb), e3vw(:,:,:,Kbb), 'VW' ) 258 ! t- and w- points depth 259 ! ---------------------- 260 gdept(:,:,1,Kmm) = 0.5_wp * e3w(:,:,1,Kmm) 261 gdepw(:,:,1,Kmm) = 0.0_wp 262 gde3w(:,:,1) = gdept(:,:,1,Kmm) - ssh(:,:,Kmm) 263 DO jk = 2, jpk 264 gdept(:,:,jk,Kmm) = gdept(:,:,jk-1,Kmm) + e3w(:,:,jk,Kmm) 265 gdepw(:,:,jk,Kmm) = gdepw(:,:,jk-1,Kmm) + e3t(:,:,jk-1,Kmm) 266 gde3w(:,:,jk) = gdept(:,:,jk ,Kmm) - sshn (:,:) 267 END DO 268 ENDIF 235 IF( .NOT.ln_linssh ) CALL dom_vvl_zgr( Kbb, Kmm, Kaa ) ! interpolation scale factor, depth and water column 236 ! IF( .NOT.ln_linssh ) THEN 237 ! ! 238 ! DO jk = 1,jpkm1 ! adjust initial vertical scale factors 239 ! e3t(:,:,jk,Kmm) = e3t_0(:,:,jk)*( 1._wp + ssh(:,:,Kmm)*tmask(:,:,1)/(ht_0(:,:) + 1.0 - tmask(:,:,1)) ) 240 ! e3t(:,:,jk,Kbb) = e3t_0(:,:,jk)*( 1._wp + ssh(:,:,Kbb)*tmask(:,:,1)/(ht_0(:,:) + 1.0 - tmask(:,:,1)) ) 241 ! ENDDO 242 ! e3t(:,:,:,Krhs) = e3t(:,:,:,Kbb) 243 ! ! Reconstruction of all vertical scale factors at now and before time-steps 244 ! ! ============================================================================= 245 ! ! Horizontal scale factor interpolations 246 ! ! -------------------------------------- 247 ! CALL dom_vvl_interpol( e3t(:,:,:,Kbb), e3u(:,:,:,Kbb), 'U' ) 248 ! CALL dom_vvl_interpol( e3t(:,:,:,Kbb), e3v(:,:,:,Kbb), 'V' ) 249 ! CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3u(:,:,:,Kmm), 'U' ) 250 ! CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3v(:,:,:,Kmm), 'V' ) 251 ! CALL dom_vvl_interpol( e3u(:,:,:,Kmm), e3f(:,:,:), 'F' ) 252 ! ! Vertical scale factor interpolations 253 ! ! ------------------------------------ 254 ! CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3w (:,:,:,Kmm), 'W' ) 255 ! CALL dom_vvl_interpol( e3u(:,:,:,Kmm), e3uw(:,:,:,Kmm), 'UW' ) 256 ! CALL dom_vvl_interpol( e3v(:,:,:,Kmm), e3vw(:,:,:,Kmm), 'VW' ) 257 ! CALL dom_vvl_interpol( e3u(:,:,:,Kbb), e3uw(:,:,:,Kbb), 'UW' ) 258 ! CALL dom_vvl_interpol( e3v(:,:,:,Kbb), e3vw(:,:,:,Kbb), 'VW' ) 259 ! ! t- and w- points depth 260 ! ! ---------------------- 261 ! gdept(:,:,1,Kmm) = 0.5_wp * e3w(:,:,1,Kmm) 262 ! gdepw(:,:,1,Kmm) = 0.0_wp 263 ! gde3w(:,:,1) = gdept(:,:,1,Kmm) - ssh(:,:,Kmm) 264 ! DO jk = 2, jpk 265 ! gdept(:,:,jk,Kmm) = gdept(:,:,jk-1,Kmm) + e3w(:,:,jk,Kmm) 266 ! gdepw(:,:,jk,Kmm) = gdepw(:,:,jk-1,Kmm) + e3t(:,:,jk-1,Kmm) 267 ! gde3w(:,:,jk) = gdept(:,:,jk ,Kmm) - sshn (:,:) 268 ! END DO 269 ! ENDIF 269 270 ENDIF 270 271 ENDIF … … 272 273 END SUBROUTINE cice_sbc_init 273 274 274 275 275 276 SUBROUTINE cice_sbc_in( kt, ksbc ) 276 277 !!--------------------------------------------------------------------- … … 281 282 INTEGER, INTENT(in ) :: ksbc ! surface forcing type 282 283 ! 283 INTEGER :: ji, jj, jl ! dummy loop indices 284 INTEGER :: ji, jj, jl ! dummy loop indices 284 285 REAL(wp), DIMENSION(jpi,jpj) :: ztmp, zpice 285 286 REAL(wp), DIMENSION(jpi,jpj,ncat) :: ztmpn … … 293 294 ztmp(:,:)=0.0 294 295 295 ! Aggregate ice concentration already set in cice_sbc_out (or cice_sbc_init on 296 ! Aggregate ice concentration already set in cice_sbc_out (or cice_sbc_init on 296 297 ! the first time-step) 297 298 298 ! forced and coupled case 299 ! forced and coupled case 299 300 300 301 IF( (ksbc == jp_flx).OR.(ksbc == jp_purecpl) ) THEN … … 356 357 ! Convert to GBM 357 358 IF(ksbc == jp_flx) THEN 358 ztmp(:,:) = (topmelt(:,:,jl)+botmelt(:,:,jl))*a_i(:,:,jl) 359 ztmp(:,:) = (topmelt(:,:,jl)+botmelt(:,:,jl))*a_i(:,:,jl) 359 360 ELSE 360 361 ztmp(:,:) = (topmelt(:,:,jl)+botmelt(:,:,jl)) … … 380 381 CALL nemo2cice(ztmp,Tair,'T', 1. ) ! Air temperature (K) 381 382 CALL nemo2cice(ztmp,potT,'T', 1. ) ! Potential temp (K) 382 ! Following line uses MAX(....) to avoid problems if tatm_ice has unset halo rows 383 ztmp(:,:) = 101000. / ( 287.04 * MAX(1.0,tatm_ice(:,:)) ) 383 ! Following line uses MAX(....) to avoid problems if tatm_ice has unset halo rows 384 ztmp(:,:) = 101000. / ( 287.04 * MAX(1.0,tatm_ice(:,:)) ) 384 385 ! Constant (101000.) atm pressure assumed 385 386 CALL nemo2cice(ztmp,rhoa,'T', 1. ) ! Air density (kg/m^3) … … 389 390 CALL nemo2cice(ztmp,zlvl,'T', 1. ) ! Atmos level height (m) 390 391 391 ! May want to check all values are physically realistic (as in CICE routine 392 ! May want to check all values are physically realistic (as in CICE routine 392 393 ! prepare_forcing)? 393 394 394 395 ! Divide shortwave into spectral bands (as in prepare_forcing) 395 396 ztmp(:,:)=qsr_ice(:,:,1)*frcvdr ! visible direct 396 CALL nemo2cice(ztmp,swvdr,'T', 1. ) 397 CALL nemo2cice(ztmp,swvdr,'T', 1. ) 397 398 ztmp(:,:)=qsr_ice(:,:,1)*frcvdf ! visible diffuse 398 CALL nemo2cice(ztmp,swvdf,'T', 1. ) 399 CALL nemo2cice(ztmp,swvdf,'T', 1. ) 399 400 ztmp(:,:)=qsr_ice(:,:,1)*frcidr ! near IR direct 400 401 CALL nemo2cice(ztmp,swidr,'T', 1. ) … … 406 407 ! Snowfall 407 408 ! Ensure fsnow is positive (as in CICE routine prepare_forcing) 408 IF( iom_use('snowpre') ) CALL iom_put('snowpre',MAX( (1.0-fr_i(:,:))*sprecip(:,:) ,0.0)) !!Joakim edit 409 ztmp(:,:)=MAX(fr_i(:,:)*sprecip(:,:),0.0) 410 CALL nemo2cice(ztmp,fsnow,'T', 1. ) 409 IF( iom_use('snowpre') ) CALL iom_put('snowpre',MAX( (1.0-fr_i(:,:))*sprecip(:,:) ,0.0)) !!Joakim edit 410 ztmp(:,:)=MAX(fr_i(:,:)*sprecip(:,:),0.0) 411 CALL nemo2cice(ztmp,fsnow,'T', 1. ) 411 412 412 413 ! Rainfall 413 414 IF( iom_use('precip') ) CALL iom_put('precip', (1.0-fr_i(:,:))*(tprecip(:,:)-sprecip(:,:)) ) !!Joakim edit 414 415 ztmp(:,:)=fr_i(:,:)*(tprecip(:,:)-sprecip(:,:)) 415 CALL nemo2cice(ztmp,frain,'T', 1. ) 416 CALL nemo2cice(ztmp,frain,'T', 1. ) 416 417 417 418 ! Freezing/melting potential … … 482 483 INTEGER, INTENT( in ) :: kt ! ocean time step 483 484 INTEGER, INTENT( in ) :: ksbc ! surface forcing type 484 485 485 486 INTEGER :: ji, jj, jl ! dummy loop indices 486 487 REAL(wp), DIMENSION(jpi,jpj) :: ztmp1, ztmp2 … … 490 491 IF(lwp) WRITE(numout,*)'cice_sbc_out' 491 492 ENDIF 492 493 ! x comp of ocean-ice stress 493 494 ! x comp of ocean-ice stress 494 495 CALL cice2nemo(strocnx,ztmp1,'F', -1. ) 495 496 ss_iou(:,:)=0.0 … … 500 501 CALL lbc_lnk( 'sbcice_cice', ss_iou , 'U', -1. ) 501 502 502 ! y comp of ocean-ice stress 503 ! y comp of ocean-ice stress 503 504 CALL cice2nemo(strocny,ztmp1,'F', -1. ) 504 505 ss_iov(:,:)=0.0 … … 513 514 ! Combine wind stress and ocean-ice stress 514 515 ! [Note that fr_iu hasn't yet been updated, so still from start of CICE timestep] 515 ! strocnx and strocny already weighted by ice fraction in CICE so not done here 516 ! strocnx and strocny already weighted by ice fraction in CICE so not done here 516 517 517 518 utau(:,:)=(1.0-fr_iu(:,:))*utau(:,:)-ss_iou(:,:) 518 vtau(:,:)=(1.0-fr_iv(:,:))*vtau(:,:)-ss_iov(:,:) 519 520 ! Also need ice/ocean stress on T points so that taum can be updated 521 ! This interpolation is already done in CICE so best to use those values 522 CALL cice2nemo(strocnxT,ztmp1,'T',-1.) 523 CALL cice2nemo(strocnyT,ztmp2,'T',-1.) 524 525 ! Update taum with modulus of ice-ocean stress 526 ! strocnxT and strocnyT are not weighted by ice fraction in CICE so must be done here 527 taum(:,:)=(1.0-fr_i(:,:))*taum(:,:)+fr_i(:,:)*SQRT(ztmp1*ztmp1 + ztmp2*ztmp2) 528 529 ! Freshwater fluxes 519 vtau(:,:)=(1.0-fr_iv(:,:))*vtau(:,:)-ss_iov(:,:) 520 521 ! Also need ice/ocean stress on T points so that taum can be updated 522 ! This interpolation is already done in CICE so best to use those values 523 CALL cice2nemo(strocnxT,ztmp1,'T',-1.) 524 CALL cice2nemo(strocnyT,ztmp2,'T',-1.) 525 526 ! Update taum with modulus of ice-ocean stress 527 ! strocnxT and strocnyT are not weighted by ice fraction in CICE so must be done here 528 taum(:,:)=(1.0-fr_i(:,:))*taum(:,:)+fr_i(:,:)*SQRT(ztmp1*ztmp1 + ztmp2*ztmp2) 529 530 ! Freshwater fluxes 530 531 531 532 IF(ksbc == jp_flx) THEN 532 533 ! Note that emp from the forcing files is evap*(1-aice)-(tprecip-aice*sprecip) 533 534 ! What we want here is evap*(1-aice)-tprecip*(1-aice) hence manipulation below 534 ! Not ideal since aice won't be the same as in the atmosphere. 535 ! Not ideal since aice won't be the same as in the atmosphere. 535 536 ! Better to use evap and tprecip? (but for now don't read in evap in this case) 536 537 emp(:,:) = emp(:,:)+fr_i(:,:)*(tprecip(:,:)-sprecip(:,:)) 537 538 ELSE IF(ksbc == jp_blk) THEN 538 emp(:,:) = (1.0-fr_i(:,:))*emp(:,:) 539 emp(:,:) = (1.0-fr_i(:,:))*emp(:,:) 539 540 ELSE IF(ksbc == jp_purecpl) THEN 540 ! emp_tot is set in sbc_cpl_ice_flx (called from cice_sbc_in above) 541 ! emp_tot is set in sbc_cpl_ice_flx (called from cice_sbc_in above) 541 542 ! This is currently as required with the coupling fields from the UM atmosphere 542 emp(:,:) = emp_tot(:,:)+tprecip(:,:)*fr_i(:,:) 543 emp(:,:) = emp_tot(:,:)+tprecip(:,:)*fr_i(:,:) 543 544 ENDIF 544 545 … … 560 561 emp(:,:)=emp(:,:)-ztmp1(:,:) 561 562 fmmflx(:,:) = ztmp1(:,:) !!Joakim edit 562 563 563 564 CALL lbc_lnk_multi( 'sbcice_cice', emp , 'T', 1., sfx , 'T', 1. ) 564 565 … … 634 635 !! *** ROUTINE cice_sbc_hadgam *** 635 636 !! ** Purpose: Prepare fields needed to pass to HadGAM3 atmosphere 636 !! 637 !! 637 638 !! 638 639 !!--------------------------------------------------------------------- … … 657 658 CALL cice2nemo(vvel,v_ice,'F', -1. ) 658 659 ! 659 ! Ice concentration (CO_1) = a_i calculated at end of cice_sbc_out 660 ! Ice concentration (CO_1) = a_i calculated at end of cice_sbc_out 660 661 ! 661 662 ! Snow and ice thicknesses (CO_2 and CO_3) … … 689 690 !!--------------------------------------------------------------------- 690 691 !! ** Method : READ monthly flux file in NetCDF files 691 !! 692 !! snowfall 693 !! rainfall 694 !! sublimation rate 692 !! 693 !! snowfall 694 !! rainfall 695 !! sublimation rate 695 696 !! topmelt (category) 696 697 !! botmelt (category) … … 709 710 TYPE(FLD_N) :: sn_snow, sn_rain, sn_sblm ! informations about the fields to be read 710 711 TYPE(FLD_N) :: sn_top1, sn_top2, sn_top3, sn_top4, sn_top5 711 TYPE(FLD_N) :: sn_bot1, sn_bot2, sn_bot3, sn_bot4, sn_bot5 712 TYPE(FLD_N) :: sn_bot1, sn_bot2, sn_bot3, sn_bot4, sn_bot5 712 713 !! 713 714 NAMELIST/namsbc_cice/ cn_dir, sn_snow, sn_rain, sn_sblm, & … … 727 728 ! ! file ! frequency ! variable ! time intep ! clim ! 'yearly' or ! weights ! rotation ! landmask 728 729 ! ! name ! (hours) ! name ! (T/F) ! (T/F) ! 'monthly' ! filename ! pairs ! file 729 sn_snow = FLD_N( 'snowfall_1m' , -1. , 'snowfall' , .true. , .true. , ' yearly' , '' , '' , '' ) 730 sn_rain = FLD_N( 'rainfall_1m' , -1. , 'rainfall' , .true. , .true. , ' yearly' , '' , '' , '' ) 730 sn_snow = FLD_N( 'snowfall_1m' , -1. , 'snowfall' , .true. , .true. , ' yearly' , '' , '' , '' ) 731 sn_rain = FLD_N( 'rainfall_1m' , -1. , 'rainfall' , .true. , .true. , ' yearly' , '' , '' , '' ) 731 732 sn_sblm = FLD_N( 'sublim_1m' , -1. , 'sublim' , .true. , .true. , ' yearly' , '' , '' , '' ) 732 733 sn_top1 = FLD_N( 'topmeltn1_1m' , -1. , 'topmeltn1' , .true. , .true. , ' yearly' , '' , '' , '' ) … … 754 755 slf_i(jp_bot2) = sn_bot2 ; slf_i(jp_bot3) = sn_bot3 ; slf_i(jp_bot4) = sn_bot4 755 756 slf_i(jp_bot5) = sn_bot5 756 757 757 758 ! set sf structure 758 759 ALLOCATE( sf(jpfld), STAT=ierror ) … … 792 793 ! control print (if less than 100 time-step asked) 793 794 IF( nitend-nit000 <= 100 .AND. lwp ) THEN 794 WRITE(numout,*) 795 WRITE(numout,*) 795 796 WRITE(numout,*) ' read forcing fluxes for CICE OK' 796 797 CALL FLUSH(numout) … … 802 803 !!--------------------------------------------------------------------- 803 804 !! *** ROUTINE nemo2cice *** 804 !! ** Purpose : Transfer field in NEMO array to field in CICE array. 805 !! ** Purpose : Transfer field in NEMO array to field in CICE array. 805 806 #if defined key_nemocice_decomp 806 !! 807 !! 807 808 !! NEMO and CICE PE sub domains are identical, hence 808 !! there is no need to gather or scatter data from 809 !! there is no need to gather or scatter data from 809 810 !! one PE configuration to another. 810 811 #else 811 !! Automatically gather/scatter between 812 !! Automatically gather/scatter between 812 813 !! different processors and blocks 813 814 !! ** Method : A. Ensure all haloes are filled in NEMO field (pn) 814 815 !! B. Gather pn into global array (png) 815 816 !! C. Map png into CICE global array (pcg) 816 !! D. Scatter pcg to CICE blocks (pc) + update haloes 817 !! D. Scatter pcg to CICE blocks (pc) + update haloes 817 818 #endif 818 819 !!--------------------------------------------------------------------- … … 858 859 IF( jpnij > 1) THEN 859 860 CALL mppsync 860 CALL mppgather (pn,0,png) 861 CALL mppgather (pn,0,png) 861 862 CALL mppsync 862 863 ELSE … … 869 870 ! (may be OK but not 100% sure) 870 871 871 IF(nproc==0) THEN 872 IF(nproc==0) THEN 872 873 ! pcg(:,:)=0.0 873 874 DO jn=1,jpnij … … 890 891 CASE ( 'T' ) 891 892 grid_loc=field_loc_center 892 CASE ( 'F' ) 893 CASE ( 'F' ) 893 894 grid_loc=field_loc_NEcorner 894 895 END SELECT … … 897 898 CASE ( -1 ) 898 899 field_type=field_type_vector 899 CASE ( 1 ) 900 CASE ( 1 ) 900 901 field_type=field_type_scalar 901 902 END SELECT … … 916 917 !! ** Purpose : Transfer field in CICE array to field in NEMO array. 917 918 #if defined key_nemocice_decomp 918 !! 919 !! 919 920 !! NEMO and CICE PE sub domains are identical, hence 920 !! there is no need to gather or scatter data from 921 !! there is no need to gather or scatter data from 921 922 !! one PE configuration to another. 922 #else 923 #else 923 924 !! Automatically deal with scatter/gather between 924 925 !! different processors and blocks … … 926 927 !! B. Map pcg into NEMO global array (png) 927 928 !! C. Scatter png into NEMO field (pn) for each processor 928 !! D. Ensure all haloes are filled in pn 929 !! D. Ensure all haloes are filled in pn 929 930 #endif 930 931 !!--------------------------------------------------------------------- … … 958 959 CASE ( 'T' ) 959 960 grid_loc=field_loc_center 960 CASE ( 'F' ) 961 CASE ( 'F' ) 961 962 grid_loc=field_loc_NEcorner 962 963 END SELECT … … 965 966 CASE ( -1 ) 966 967 field_type=field_type_vector 967 CASE ( 1 ) 968 CASE ( 1 ) 968 969 field_type=field_type_scalar 969 970 END SELECT … … 979 980 #else 980 981 981 ! A. Gather CICE blocks (pc) into global array (pcg) 982 ! A. Gather CICE blocks (pc) into global array (pcg) 982 983 983 984 CALL gather_global(pcg, pc, 0, distrb_info) … … 1005 1006 IF( jpnij > 1) THEN 1006 1007 CALL mppsync 1007 CALL mppscatter (png,0,pn) 1008 CALL mppscatter (png,0,pn) 1008 1009 CALL mppsync 1009 1010 ELSE -
NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/TRA/traatfLF.F90
r12581 r12583 103 103 ENDIF 104 104 105 ! Update after tracer on domain lateral boundaries106 !107 #if defined key_agrif108 CALL Agrif_tra ! AGRIF zoom boundaries109 #endif110 ! ! local domain boundaries (T-point, unchanged sign)111 CALL lbc_lnk_multi( 'traatfLF', pts(:,:,:,jp_tem,Kaa), 'T', 1., pts(:,:,:,jp_sal,Kaa), 'T', 1. )112 !113 IF( ln_bdy ) CALL bdy_tra( kt, Kbb, pts, Kaa ) ! BDY open boundaries105 ! ! Update after tracer on domain lateral boundaries 106 ! ! 107 ! #if defined key_agrif 108 ! CALL Agrif_tra ! AGRIF zoom boundaries 109 ! #endif 110 ! ! ! local domain boundaries (T-point, unchanged sign) 111 ! CALL lbc_lnk_multi( 'traatfLF', pts(:,:,:,jp_tem,Kaa), 'T', 1., pts(:,:,:,jp_sal,Kaa), 'T', 1. ) 112 ! ! 113 ! IF( ln_bdy ) CALL bdy_tra( kt, Kbb, pts, Kaa ) ! BDY open boundaries 114 114 115 115 ! set time step size (Euler/Leapfrog) -
NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/steplf.F90
r12581 r12583 43 43 USE traatfLF ! time filtering (tra_atf_lf routine) 44 44 USE dynatfLF ! time filtering (dyn_atf_lf routine) 45 USE dynspg_ts ! surface pressure gradient: split-explicit scheme (define un_adv) 46 USE bdydyn ! ocean open boundary conditions (define bdy_dyn) 45 47 46 48 IMPLICIT NONE … … 288 290 !! 289 291 !!jc2: dynnxt must be the latest call. e3t(:,:,:,Nbb) are indeed updated in that routine 292 CALL zdyn_ts ( Nnn, Naa, e3u, e3v, uu, vv ) ! barotrope ajustment 293 CALL finalize_sbc ( kstp, Nbb, Naa, uu, vv, ts ) ! boundary condifions 290 294 CALL ssh_atf ( kstp, Nbb, Nnn, Naa, ssh ) ! time filtering of "now" sea surface height 291 295 CALL dom_qe_r3c ( ssh(:,:,Nnn), r3t_f, r3u_f, r3v_f ) ! "now" ssh/h_0 ratio from filtrered ssh … … 352 356 ! 353 357 END SUBROUTINE stplf 358 359 360 SUBROUTINE zdyn_ts (Kmm, Kaa, pe3u, pe3v, puu, pvv) 361 !!---------------------------------------------------------------------- 362 !! *** ROUTINE zdyn_ts *** 363 !! 364 !! ** Purpose : Finalize after horizontal velocity. 365 !! 366 !! ** Method : * Ensure after velocities transport matches time splitting 367 !! estimate (ln_dynspg_ts=T) 368 !! 369 !! ** Action : puu(Kmm),pvv(Kmm),puu(Kaa),pvv(Kaa) now and after horizontal velocity 370 !!---------------------------------------------------------------------- 371 INTEGER , INTENT(in ) :: Kmm, Kaa ! before and after time level indices 372 REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv ! velocities 373 REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(in ) :: pe3u, pe3v ! scale factors 374 ! 375 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zue, zve 376 ! 377 INTEGER :: jk ! dummy loop indices 378 !!---------------------------------------------------------------------- 379 380 IF ( ln_dynspg_ts ) THEN 381 ALLOCATE( zue(jpi,jpj) , zve(jpi,jpj) ) 382 ! Ensure below that barotropic velocities match time splitting estimate 383 ! Compute actual transport and replace it with ts estimate at "after" time step 384 zue(:,:) = pe3u(:,:,1,Kaa) * puu(:,:,1,Kaa) * umask(:,:,1) 385 zve(:,:) = pe3v(:,:,1,Kaa) * pvv(:,:,1,Kaa) * vmask(:,:,1) 386 DO jk = 2, jpkm1 387 zue(:,:) = zue(:,:) + pe3u(:,:,jk,Kaa) * puu(:,:,jk,Kaa) * umask(:,:,jk) 388 zve(:,:) = zve(:,:) + pe3v(:,:,jk,Kaa) * pvv(:,:,jk,Kaa) * vmask(:,:,jk) 389 END DO 390 DO jk = 1, jpkm1 391 puu(:,:,jk,Kaa) = ( puu(:,:,jk,Kaa) - zue(:,:) * r1_hu(:,:,Kaa) + uu_b(:,:,Kaa) ) * umask(:,:,jk) 392 pvv(:,:,jk,Kaa) = ( pvv(:,:,jk,Kaa) - zve(:,:) * r1_hv(:,:,Kaa) + vv_b(:,:,Kaa) ) * vmask(:,:,jk) 393 END DO 394 ! 395 IF( .NOT.ln_bt_fw ) THEN 396 ! Remove advective velocity from "now velocities" 397 ! prior to asselin filtering 398 ! In the forward case, this is done below after asselin filtering 399 ! so that asselin contribution is removed at the same time 400 DO jk = 1, jpkm1 401 puu(:,:,jk,Kmm) = ( puu(:,:,jk,Kmm) - un_adv(:,:)*r1_hu(:,:,Kmm) + uu_b(:,:,Kmm) )*umask(:,:,jk) 402 pvv(:,:,jk,Kmm) = ( pvv(:,:,jk,Kmm) - vn_adv(:,:)*r1_hv(:,:,Kmm) + vv_b(:,:,Kmm) )*vmask(:,:,jk) 403 END DO 404 ENDIF 405 ! 406 DEALLOCATE( zue, zve ) 407 ! 408 ENDIF 409 ! 410 END SUBROUTINE zdyn_ts 411 412 413 SUBROUTINE finalize_sbc (kt, Kbb, Kaa, puu, pvv, pts) 414 !!---------------------------------------------------------------------- 415 !! *** ROUTINE finalize_sbc *** 416 !! 417 !! ** Purpose : Apply the boundary condition on the after velocity 418 !! 419 !! ** Method : * Apply lateral boundary conditions on after velocity 420 !! at the local domain boundaries through lbc_lnk call, 421 !! at the one-way open boundaries (ln_bdy=T), 422 !! at the AGRIF zoom boundaries (lk_agrif=T) 423 !! 424 !! ** Action : puu(Kaa),pvv(Kaa) after horizontal velocity and tracers 425 !!---------------------------------------------------------------------- 426 INTEGER , INTENT(in ) :: kt ! ocean time-step index 427 INTEGER , INTENT(in ) :: Kbb, Kaa ! before and after time level indices 428 REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv ! velocities to be time filtered 429 REAL(wp), DIMENSION(jpi,jpj,jpk,jpts,jpt), INTENT(inout) :: pts ! active tracers 430 ! 431 ! Update after tracer and velocity on domain lateral boundaries 432 ! 433 #if defined key_agrif 434 CALL Agrif_tra ! AGRIF zoom boundaries 435 CALL Agrif_dyn( kt ) !* AGRIF zoom boundaries 436 #endif 437 ! ! local domain boundaries (T-point, unchanged sign) 438 CALL lbc_lnk_multi( 'finalize_sbc', puu(:,:,:, Kaa), 'U', -1., pvv(:,:,: ,Kaa), 'V', -1. & 439 & , pts(:,:,:,jp_tem,Kaa), 'T', 1., pts(:,:,:,jp_sal,Kaa), 'T', 1. ) !* local domain boundaries 440 ! 441 ! !* BDY open boundaries 442 IF( ln_bdy ) THEN 443 CALL bdy_tra( kt, Kbb, pts, Kaa ) 444 IF( ln_dynspg_exp ) CALL bdy_dyn( kt, Kbb, puu, pvv, Kaa ) 445 IF( ln_dynspg_ts ) CALL bdy_dyn( kt, Kbb, puu, pvv, Kaa, dyn3d_only=.true. ) 446 ENDIF 447 ! 448 END SUBROUTINE finalize_sbc 449 354 450 ! 355 451 !!======================================================================
Note: See TracChangeset
for help on using the changeset viewer.
