Changeset 14015 for NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/src/OCE/ZDF
- Timestamp:
- 2020-12-02T16:26:43+01:00 (4 years ago)
- Location:
- NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source
- Files:
-
- 8 edited
- 1 copied
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source
- Property svn:externals
-
old new 8 8 9 9 # SETTE 10 ^/utils/CI/sette @13559sette10 ^/utils/CI/sette_wave@13990 sette
-
- Property svn:externals
-
NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/src/OCE/ZDF/zdf_oce.F90
r10425 r14015 40 40 LOGICAL , PUBLIC :: ln_zdfswm !: surface wave-induced mixing flag 41 41 LOGICAL , PUBLIC :: ln_zdfiwm !: internal wave-induced mixing flag 42 LOGICAL , PUBLIC :: ln_zdfmfc !: convection: eddy diffusivity Mass Flux Convection 42 43 ! ! coefficients 43 44 REAL(wp), PUBLIC :: rn_avm0 !: vertical eddy viscosity (m2/s) -
NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/src/OCE/ZDF/zdfgls.F90
r13886 r14015 1057 1057 CALL gls_rst( nit000, 'READ' ) ! (en, avt_k, avm_k, hmxl_n) 1058 1058 ! 1059 IF( lwxios ) THEN1060 CALL iom_set_rstw_var_active('en')1061 CALL iom_set_rstw_var_active('avt_k')1062 CALL iom_set_rstw_var_active('avm_k')1063 CALL iom_set_rstw_var_active('hmxl_n')1064 ENDIF1065 !1066 1059 END SUBROUTINE zdf_gls_init 1067 1060 … … 1097 1090 ! 1098 1091 IF( MIN( id1, id2, id3, id4 ) > 0 ) THEN ! all required arrays exist 1099 CALL iom_get( numror, jpdom_auto, 'en' , en , ldxios = lrxios)1100 CALL iom_get( numror, jpdom_auto, 'avt_k' , avt_k , ldxios = lrxios)1101 CALL iom_get( numror, jpdom_auto, 'avm_k' , avm_k , ldxios = lrxios)1102 CALL iom_get( numror, jpdom_auto, 'hmxl_n', hmxl_n , ldxios = lrxios)1092 CALL iom_get( numror, jpdom_auto, 'en' , en ) 1093 CALL iom_get( numror, jpdom_auto, 'avt_k' , avt_k ) 1094 CALL iom_get( numror, jpdom_auto, 'avm_k' , avm_k ) 1095 CALL iom_get( numror, jpdom_auto, 'hmxl_n', hmxl_n ) 1103 1096 ELSE 1104 1097 IF(lwp) WRITE(numout,*) … … 1119 1112 ! ! ------------------- 1120 1113 IF(lwp) WRITE(numout,*) '---- gls-rst ----' 1121 IF( lwxios ) CALL iom_swap( cwxios_context ) 1122 CALL iom_rstput( kt, nitrst, numrow, 'en' , en , ldxios = lwxios ) 1123 CALL iom_rstput( kt, nitrst, numrow, 'avt_k' , avt_k , ldxios = lwxios ) 1124 CALL iom_rstput( kt, nitrst, numrow, 'avm_k' , avm_k , ldxios = lwxios ) 1125 CALL iom_rstput( kt, nitrst, numrow, 'hmxl_n', hmxl_n, ldxios = lwxios ) 1126 IF( lwxios ) CALL iom_swap( cxios_context ) 1114 CALL iom_rstput( kt, nitrst, numrow, 'en' , en ) 1115 CALL iom_rstput( kt, nitrst, numrow, 'avt_k' , avt_k ) 1116 CALL iom_rstput( kt, nitrst, numrow, 'avm_k' , avm_k ) 1117 CALL iom_rstput( kt, nitrst, numrow, 'hmxl_n', hmxl_n ) 1127 1118 ! 1128 1119 ENDIF -
NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/src/OCE/ZDF/zdfosm.F90
r13497 r14015 1437 1437 ghamv(:,:,:) = 0. 1438 1438 ! 1439 IF( lwxios ) THEN1440 CALL iom_set_rstw_var_active('wn')1441 CALL iom_set_rstw_var_active('hbl')1442 CALL iom_set_rstw_var_active('hbli')1443 ENDIF1444 1439 END SUBROUTINE zdf_osm_init 1445 1440 … … 1474 1469 id1 = iom_varid( numror, 'wn' , ldstop = .FALSE. ) 1475 1470 IF( id1 > 0 ) THEN ! 'wn' exists; read 1476 CALL iom_get( numror, jpdom_auto, 'wn', ww , ldxios = lrxios)1471 CALL iom_get( numror, jpdom_auto, 'wn', ww ) 1477 1472 WRITE(numout,*) ' ===>>>> : ww read from restart file' 1478 1473 ELSE … … 1483 1478 id2 = iom_varid( numror, 'hbli' , ldstop = .FALSE. ) 1484 1479 IF( id1 > 0 .AND. id2 > 0) THEN ! 'hbl' exists; read and return 1485 CALL iom_get( numror, jpdom_auto, 'hbl' , hbl , ldxios = lrxios)1486 CALL iom_get( numror, jpdom_auto, 'hbli', hbli , ldxios = lrxios)1480 CALL iom_get( numror, jpdom_auto, 'hbl' , hbl ) 1481 CALL iom_get( numror, jpdom_auto, 'hbli', hbli ) 1487 1482 WRITE(numout,*) ' ===>>>> : hbl & hbli read from restart file' 1488 1483 RETURN … … 1496 1491 !!----------------------------------------------------------------------------- 1497 1492 IF( TRIM(cdrw) == 'WRITE') THEN !* Write hbli into the restart file, then return 1493 IF( ntile /= 0 .AND. ntile /= nijtile ) RETURN ! Do only on the last tile 1494 1498 1495 IF(lwp) WRITE(numout,*) '---- osm-rst ----' 1499 CALL iom_rstput( kt, nitrst, numrow, 'wn' , ww , ldxios = lwxios)1500 CALL iom_rstput( kt, nitrst, numrow, 'hbl' , hbl , ldxios = lwxios)1501 CALL iom_rstput( kt, nitrst, numrow, 'hbli' , hbli , ldxios = lwxios)1496 CALL iom_rstput( kt, nitrst, numrow, 'wn' , ww ) 1497 CALL iom_rstput( kt, nitrst, numrow, 'hbl' , hbl ) 1498 CALL iom_rstput( kt, nitrst, numrow, 'hbli' , hbli ) 1502 1499 RETURN 1503 1500 END IF … … 1550 1547 ! 1551 1548 IF( kt == nit000 ) THEN 1552 IF(lwp) WRITE(numout,*) 1553 IF(lwp) WRITE(numout,*) 'tra_osm : OSM non-local tracer fluxes' 1554 IF(lwp) WRITE(numout,*) '~~~~~~~ ' 1549 IF( ntile == 0 .OR. ntile == 1 ) THEN ! Do only on the first tile 1550 IF(lwp) WRITE(numout,*) 1551 IF(lwp) WRITE(numout,*) 'tra_osm : OSM non-local tracer fluxes' 1552 IF(lwp) WRITE(numout,*) '~~~~~~~ ' 1553 ENDIF 1555 1554 ENDIF 1556 1555 -
NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/src/OCE/ZDF/zdfphy.F90
r13886 r14015 21 21 USE zdfddm ! vertical physics: double diffusion mixing 22 22 USE zdfevd ! vertical physics: convection via enhanced vertical diffusion 23 USE zdfmfc ! vertical physics: Mass Flux Convection 23 24 USE zdfiwm ! vertical physics: internal wave-induced mixing 24 25 USE zdfswm ! vertical physics: surface wave-induced mixing … … 78 79 NAMELIST/namzdf/ ln_zdfcst, ln_zdfric, ln_zdftke, ln_zdfgls, & ! type of closure scheme 79 80 & ln_zdfosm, & ! type of closure scheme 81 & ln_zdfmfc, & ! convection : mass flux 80 82 & ln_zdfevd, nn_evdm, rn_evd , & ! convection : evd 81 83 & ln_zdfnpc, nn_npc , nn_npcp, & ! convection : npc … … 112 114 WRITE(numout,*) ' OSMOSIS-OBL closure (OSM) ln_zdfosm = ', ln_zdfosm 113 115 WRITE(numout,*) ' convection: ' 116 WRITE(numout,*) ' convection mass flux (mfc) ln_zdfmfc = ', ln_zdfmfc 114 117 WRITE(numout,*) ' enhanced vertical diffusion ln_zdfevd = ', ln_zdfevd 115 118 WRITE(numout,*) ' applied on momentum (=1/0) nn_evdm = ', nn_evdm … … 172 175 IF( ln_zdfnpc .AND. ln_zdfevd ) CALL ctl_stop( 'zdf_phy_init: chose between ln_zdfnpc and ln_zdfevd' ) 173 176 IF( ln_zdfosm .AND. ln_zdfevd ) CALL ctl_stop( 'zdf_phy_init: chose between ln_zdfosm and ln_zdfevd' ) 177 IF( ln_zdfmfc .AND. ln_zdfevd ) CALL ctl_stop( 'zdf_phy_init: chose between ln_zdfmfc and ln_zdfevd' ) 178 IF( ln_zdfmfc .AND. ln_zdfnpc ) CALL ctl_stop( 'zdf_phy_init: chose between ln_zdfmfc and ln_zdfnpc' ) 179 IF( ln_zdfmfc .AND. ln_zdfosm ) CALL ctl_stop( 'zdf_phy_init: chose between ln_zdfmfc and ln_zdfosm' ) 174 180 IF( lk_top .AND. ln_zdfnpc ) CALL ctl_stop( 'zdf_phy_init: npc scheme is not working with key_top' ) 175 181 IF( lk_top .AND. ln_zdfosm ) CALL ctl_stop( 'zdf_phy_init: osmosis scheme is not working with key_top' ) 182 IF( lk_top .AND. ln_zdfmfc ) CALL ctl_stop( 'zdf_phy_init: Mass Flux scheme is not working with key_top' ) 176 183 IF(lwp) THEN 177 184 WRITE(numout,*) 178 185 IF ( ln_zdfnpc ) THEN ; WRITE(numout,*) ' ==>>> convection: use non penetrative convective scheme' 179 186 ELSEIF( ln_zdfevd ) THEN ; WRITE(numout,*) ' ==>>> convection: use enhanced vertical diffusion scheme' 187 ELSEIF( ln_zdfmfc ) THEN ; WRITE(numout,*) ' ==>>> convection: use Mass Flux scheme' 180 188 ELSE ; WRITE(numout,*) ' ==>>> convection: no specific scheme used' 181 189 ENDIF … … 205 213 ELSE ; l_zdfsh2 = .TRUE. 206 214 ENDIF 207 215 ! !== Mass Flux Convectiive algorithm ==! 216 IF( ln_zdfmfc ) CALL zdf_mfc_init ! Convection computed with eddy diffusivity mass flux 217 ! 208 218 ! !== gravity wave-driven mixing ==! 209 219 IF( ln_zdfiwm ) CALL zdf_iwm_init ! internal wave-driven mixing -
NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/src/OCE/ZDF/zdfric.F90
r13497 r14015 103 103 CALL ric_rst( nit000, 'READ' ) !* read or initialize all required files 104 104 ! 105 IF( lwxios ) THEN106 CALL iom_set_rstw_var_active('avt_k')107 CALL iom_set_rstw_var_active('avm_k')108 ENDIF109 105 END SUBROUTINE zdf_ric_init 110 106 … … 214 210 ! 215 211 IF( MIN( id1, id2 ) > 0 ) THEN ! restart exists => read it 216 CALL iom_get( numror, jpdom_auto, 'avt_k', avt_k , ldxios = lrxios)217 CALL iom_get( numror, jpdom_auto, 'avm_k', avm_k , ldxios = lrxios)212 CALL iom_get( numror, jpdom_auto, 'avt_k', avt_k ) 213 CALL iom_get( numror, jpdom_auto, 'avm_k', avm_k ) 218 214 ENDIF 219 215 ENDIF … … 223 219 ! ! ------------------- 224 220 IF(lwp) WRITE(numout,*) '---- ric-rst ----' 225 IF( lwxios ) CALL iom_swap( cwxios_context ) 226 CALL iom_rstput( kt, nitrst, numrow, 'avt_k', avt_k, ldxios = lwxios ) 227 CALL iom_rstput( kt, nitrst, numrow, 'avm_k', avm_k, ldxios = lwxios) 228 IF( lwxios ) CALL iom_swap( cxios_context ) 221 CALL iom_rstput( kt, nitrst, numrow, 'avt_k', avt_k ) 222 CALL iom_rstput( kt, nitrst, numrow, 'avm_k', avm_k) 229 223 ! 230 224 ENDIF -
NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/src/OCE/ZDF/zdfsh2.F90
r13497 r14015 6 6 !! History : - ! 2014-10 (A. Barthelemy, G. Madec) original code 7 7 !! NEMO 4.0 ! 2017-04 (G. Madec) remove u-,v-pts avm 8 !! NEMO 4.2 ! 2020-12 (G. Madec, E. Clementi) add Stokes Drift Shear 9 ! ! for wave coupling 8 10 !!---------------------------------------------------------------------- 9 11 … … 13 15 USE oce 14 16 USE dom_oce ! domain: ocean 17 USE sbcwave ! Surface Waves (add Stokes shear) 18 USE sbc_oce , ONLY: ln_stshear !Stoked Drift shear contribution 15 19 ! 16 20 USE in_out_manager ! I/O manager … … 21 25 22 26 PUBLIC zdf_sh2 ! called by zdftke, zdfglf, and zdfric 23 27 24 28 !! * Substitutions 25 29 # include "do_loop_substitute.h90" … … 59 63 !!-------------------------------------------------------------------- 60 64 ! 61 DO jk = 2, jpkm1 62 DO_2D( 1, 0, 1, 0 ) !* 2 x shear production at uw- and vw-points (energy conserving form) 63 zsh2u(ji,jj) = ( p_avm(ji+1,jj,jk) + p_avm(ji,jj,jk) ) & 64 & * ( uu(ji,jj,jk-1,Kmm) - uu(ji,jj,jk,Kmm) ) & 65 & * ( uu(ji,jj,jk-1,Kbb) - uu(ji,jj,jk,Kbb) ) & 66 & / ( e3uw(ji,jj,jk ,Kmm) * e3uw(ji,jj,jk,Kbb) ) & 67 & * wumask(ji,jj,jk) 68 zsh2v(ji,jj) = ( p_avm(ji,jj+1,jk) + p_avm(ji,jj,jk) ) & 69 & * ( vv(ji,jj,jk-1,Kmm) - vv(ji,jj,jk,Kmm) ) & 70 & * ( vv(ji,jj,jk-1,Kbb) - vv(ji,jj,jk,Kbb) ) & 71 & / ( e3vw(ji,jj,jk ,Kmm) * e3vw(ji,jj,jk,Kbb) ) & 72 & * wvmask(ji,jj,jk) 73 END_2D 65 DO jk = 2, jpkm1 !* Shear production at uw- and vw-points (energy conserving form) 66 IF ( cpl_sdrftx .AND. ln_stshear ) THEN ! Surface Stokes Drift available ===>>> shear + stokes drift contibution 67 DO_2D( 1, 0, 1, 0 ) 68 zsh2u(ji,jj) = ( p_avm(ji+1,jj,jk) + p_avm(ji,jj,jk) ) & 69 & * ( uu (ji,jj,jk-1,Kmm) - uu (ji,jj,jk,Kmm) & 70 & + usd(ji,jj,jk-1) - usd(ji,jj,jk) ) & 71 & * ( uu (ji,jj,jk-1,Kbb) - uu (ji,jj,jk,Kbb) ) & 72 & / ( e3uw(ji,jj,jk,Kmm) * e3uw(ji,jj,jk,Kbb) ) * wumask(ji,jj,jk) 73 zsh2v(ji,jj) = ( p_avm(ji,jj+1,jk) + p_avm(ji,jj,jk) ) & 74 & * ( vv (ji,jj,jk-1,Kmm) - vv (ji,jj,jk,Kmm) & 75 & + vsd(ji,jj,jk-1) - vsd(ji,jj,jk) ) & 76 & * ( vv (ji,jj,jk-1,Kbb) - vv (ji,jj,jk,Kbb) ) & 77 &/ ( e3vw(ji,jj,jk,Kmm) * e3vw(ji,jj,jk,Kbb) ) * wvmask(ji,jj,jk) 78 END_2D 79 ELSE 80 DO_2D( 1, 0, 1, 0 ) !* 2 x shear production at uw- and vw-points (energy conserving form) 81 zsh2u(ji,jj) = ( p_avm(ji+1,jj,jk) + p_avm(ji,jj,jk) ) & 82 & * ( uu(ji,jj,jk-1,Kmm) - uu(ji,jj,jk,Kmm) ) & 83 & * ( uu(ji,jj,jk-1,Kbb) - uu(ji,jj,jk,Kbb) ) & 84 & / ( e3uw(ji,jj,jk ,Kmm) * e3uw(ji,jj,jk,Kbb) ) & 85 & * wumask(ji,jj,jk) 86 zsh2v(ji,jj) = ( p_avm(ji,jj+1,jk) + p_avm(ji,jj,jk) ) & 87 & * ( vv(ji,jj,jk-1,Kmm) - vv(ji,jj,jk,Kmm) ) & 88 & * ( vv(ji,jj,jk-1,Kbb) - vv(ji,jj,jk,Kbb) ) & 89 & / ( e3vw(ji,jj,jk ,Kmm) * e3vw(ji,jj,jk,Kbb) ) & 90 & * wvmask(ji,jj,jk) 91 END_2D 92 ENDIF 74 93 DO_2D( 0, 0, 0, 0 ) !* shear production at w-point ! coast mask: =2 at the coast ; =1 otherwise (NB: wmask useless as zsh2 are masked) 75 94 p_sh2(ji,jj,jk) = 0.25 * ( ( zsh2u(ji-1,jj) + zsh2u(ji,jj) ) * ( 2. - umask(ji-1,jj,jk) * umask(ji,jj,jk) ) & 76 95 & + ( zsh2v(ji,jj-1) + zsh2v(ji,jj) ) * ( 2. - vmask(ji,jj-1,jk) * vmask(ji,jj,jk) ) ) 77 96 END_2D 78 END DO 97 END DO 79 98 ! 80 99 END SUBROUTINE zdf_sh2 -
NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/src/OCE/ZDF/zdftke.F90
r13886 r14015 29 29 !! 4.0 ! 2017-04 (G. Madec) remove CPP ddm key & avm at t-point only 30 30 !! - ! 2017-05 (G. Madec) add top/bottom friction as boundary condition 31 !! 4.2 ! 2020-12 (G. Madec, E. Clementi) add wave coupling 32 ! ! following Couvelard et al., 2019 31 33 !!---------------------------------------------------------------------- 32 34 … … 58 60 USE prtctl ! Print control 59 61 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 62 USE sbcwave ! Surface boundary waves 60 63 61 64 IMPLICIT NONE … … 68 71 ! !!** Namelist namzdf_tke ** 69 72 LOGICAL :: ln_mxl0 ! mixing length scale surface value as function of wind stress or not 73 LOGICAL :: ln_mxhsw ! mixing length scale surface value as a fonction of wave height 70 74 INTEGER :: nn_mxlice ! type of scaling under sea-ice (=0/1/2/3) 71 75 REAL(wp) :: rn_mxlice ! ice thickness value when scaling under sea-ice … … 81 85 INTEGER :: nn_etau ! type of depth penetration of surface tke (=0/1/2/3) 82 86 INTEGER :: nn_htau ! type of tke profile of penetration (=0/1) 87 INTEGER :: nn_bc_surf! surface condition (0/1=Dir/Neum) ! Only applicable for wave coupling 88 INTEGER :: nn_bc_bot ! surface condition (0/1=Dir/Neum) ! Only applicable for wave coupling 83 89 REAL(wp) :: rn_efr ! fraction of TKE surface value which penetrates in the ocean 84 90 LOGICAL :: ln_lc ! Langmuir cells (LC) as a source term of TKE or not … … 209 215 REAL(wp) :: zus , zwlc , zind ! - - 210 216 REAL(wp) :: zzd_up, zzd_lw ! - - 217 REAL(wp) :: ztaui, ztauj, z1_norm 211 218 INTEGER , DIMENSION(jpi,jpj) :: imlc 212 REAL(wp), DIMENSION(jpi,jpj) :: zice_fra, zhlc, zus3 219 REAL(wp), DIMENSION(jpi,jpj) :: zice_fra, zhlc, zus3, zWlc2 213 220 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zpelc, zdiag, zd_up, zd_lw 214 221 !!-------------------------------------------------------------------- … … 219 226 zfact2 = 1.5_wp * rn_Dt * rn_ediss 220 227 zfact3 = 0.5_wp * rn_ediss 228 ! 229 zpelc(:,:,:) = 0._wp ! need to be initialised in case ln_lc is not used 221 230 ! 222 231 ! ice fraction considered for attenuation of langmuir & wave breaking … … 232 241 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 233 242 ! 234 DO_2D( 0, 0, 0, 0 ) ! en(1) = rn_ebb taum / rau0 (min value rn_emin0) 235 !! clem: this should be the right formulation but it makes the model unstable unless drags are calculated implicitly 236 !! one way around would be to increase zbbirau 237 !! en(ji,jj,1) = MAX( rn_emin0, ( ( 1._wp - fr_i(ji,jj) ) * zbbrau + & 238 !! & fr_i(ji,jj) * zbbirau ) * taum(ji,jj) ) * tmask(ji,jj,1) 243 DO_2D( 0, 0, 0, 0 ) 239 244 en(ji,jj,1) = MAX( rn_emin0, zbbrau * taum(ji,jj) ) * tmask(ji,jj,1) 245 zdiag(ji,jj,1) = 1._wp/en(ji,jj,1) 246 zd_lw(ji,jj,1) = 1._wp 247 zd_up(ji,jj,1) = 0._wp 240 248 END_2D 241 249 ! … … 274 282 ! 275 283 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 276 IF( ln_lc ) THEN ! Langmuir circulation source term added to tke !(Axell JGR 2002)284 IF( ln_lc ) THEN ! Langmuir circulation source term added to tke (Axell JGR 2002) 277 285 ! !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 278 286 ! 279 ! !* total energy produce by LC : cumulative sum over jk 287 ! !* Langmuir velocity scale 288 ! 289 IF ( cpl_sdrftx ) THEN ! Surface Stokes Drift available 290 ! ! Craik-Leibovich velocity scale Wlc = ( u* u_s )^1/2 with u* = (taum/rho0)^1/2 291 ! ! associated kinetic energy : 1/2 (Wlc)^2 = u* u_s 292 ! ! more precisely, it is the dot product that must be used : 293 ! ! 1/2 (W_lc)^2 = MAX( u* u_s + v* v_s , 0 ) only the positive part 294 !!gm ! PS: currently we don't have neither the 2 stress components at t-point !nor the angle between u* and u_s 295 !!gm ! so we will overestimate the LC velocity.... !!gm I will do the work if !LC have an effect ! 296 DO_2D( 0, 0, 0, 0 ) 297 !!XC zWlc2(ji,jj) = 0.5_wp * SQRT( taum(ji,jj) * r1_rho0 * ( ut0sd(ji,jj)**2 +vt0sd(ji,jj)**2 ) ) 298 zWlc2(ji,jj) = 0.5_wp * ( ut0sd(ji,jj)**2 +vt0sd(ji,jj)**2 ) 299 END_2D 300 ! 301 ! Projection of Stokes drift in the wind stress direction 302 ! 303 DO_2D( 0, 0, 0, 0 ) 304 ztaui = 0.5_wp * ( utau(ji,jj) + utau(ji-1,jj) ) 305 ztauj = 0.5_wp * ( vtau(ji,jj) + vtau(ji,jj-1) ) 306 z1_norm = 1._wp / MAX( SQRT(ztaui*ztaui+ztauj*ztauj), 1.e-12 ) * tmask(ji,jj,1) 307 zWlc2(ji,jj) = 0.5_wp * z1_norm * ( MAX( ut0sd(ji,jj)*ztaui + vt0sd(ji,jj)*ztauj, 0._wp ) )**2 308 END_2D 309 CALL lbc_lnk ( 'zdftke', zWlc2, 'T', 1. ) 310 ! 311 ELSE ! Surface Stokes drift deduced from surface stress 312 ! ! Wlc = u_s with u_s = 0.016*U_10m, the surface stokes drift (Axell 2002, Eq.44) 313 ! ! using |tau| = rho_air Cd |U_10m|^2 , it comes: 314 ! ! Wlc = 0.016 * [|tau|/(rho_air Cdrag) ]^1/2 and thus: 315 ! ! 1/2 Wlc^2 = 0.5 * 0.016 * 0.016 |tau| /( rho_air Cdrag ) 316 zcof = 0.5 * 0.016 * 0.016 / ( zrhoa * zcdrag ) ! to convert stress in 10m wind using a constant drag 317 DO_2D( 1, 1, 1, 1 ) 318 zWlc2(ji,jj) = zcof * taum(ji,jj) 319 END_2D 320 ! 321 ENDIF 322 ! 323 ! !* Depth of the LC circulation (Axell 2002, Eq.47) 324 ! !- LHS of Eq.47 280 325 zpelc(:,:,1) = MAX( rn2b(:,:,1), 0._wp ) * gdepw(:,:,1,Kmm) * e3w(:,:,1,Kmm) 281 326 DO jk = 2, jpk … … 283 328 & MAX( rn2b(:,:,jk), 0._wp ) * gdepw(:,:,jk,Kmm) * e3w(:,:,jk,Kmm) 284 329 END DO 285 ! !* finite Langmuir Circulation depth286 zcof = 0.5 * 0.016 * 0.016 / ( zrhoa * zcdrag )330 ! 331 ! !- compare LHS to RHS of Eq.47 287 332 imlc(:,:) = mbkt(:,:) + 1 ! Initialization to the number of w ocean point (=2 over land) 288 DO_3DS( 1, 1, 1, 1, jpkm1, 2, -1 ) ! Last w-level at which zpelc>=0.5*us*us 289 zus = zcof * taum(ji,jj) ! with us=0.016*wind(starting from jpk-1) 290 IF( zpelc(ji,jj,jk) > zus ) imlc(ji,jj) = jk 333 DO_3DS( 1, 1, 1, 1, jpkm1, 2, -1 ) 334 IF( zpelc(ji,jj,jk) > zWlc2(ji,jj) ) imlc(ji,jj) = jk 291 335 END_3D 292 336 ! ! finite LC depth … … 294 338 zhlc(ji,jj) = gdepw(ji,jj,imlc(ji,jj),Kmm) 295 339 END_2D 340 ! 296 341 zcof = 0.016 / SQRT( zrhoa * zcdrag ) 297 342 DO_2D( 0, 0, 0, 0 ) 298 zus = zcof * SQRT( taum(ji,jj) )! Stokes drift343 zus = SQRT( 2. * zWlc2(ji,jj) ) ! Stokes drift 299 344 zus3(ji,jj) = MAX( 0._wp, 1._wp - zice_fra(ji,jj) ) * zus * zus * zus * tmask(ji,jj,1) ! zus > 0. ok 300 345 END_2D … … 351 396 & ) * wmask(ji,jj,jk) 352 397 END_3D 398 ! 399 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 400 ! ! Surface boundary condition on tke if 401 ! ! coupling with waves 402 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 403 ! 404 IF ( cpl_phioc .and. ln_phioc ) THEN 405 SELECT CASE (nn_bc_surf) ! Boundary Condition using surface TKE flux from waves 406 407 CASE ( 0 ) ! Dirichlet BC 408 DO_2D( 0, 0, 0, 0 ) ! en(1) = rn_ebb taum / rho0 (min value rn_emin0) 409 IF ( phioc(ji,jj) < 0 ) phioc(ji,jj) = 0._wp 410 en(ji,jj,1) = MAX( rn_emin0, .5 * ( 15.8 * phioc(ji,jj) / rho0 )**(2./3.) ) * tmask(ji,jj,1) 411 zdiag(ji,jj,1) = 1._wp/en(ji,jj,1) ! choose to keep coherence with former estimation of 412 END_2D 413 414 CASE ( 1 ) ! Neumann BC 415 DO_2D( 0, 0, 0, 0 ) 416 IF ( phioc(ji,jj) < 0 ) phioc(ji,jj) = 0._wp 417 en(ji,jj,2) = en(ji,jj,2) + ( rn_Dt * phioc(ji,jj) / rho0 ) /e3w(ji,jj,2,Kmm) 418 en(ji,jj,1) = en(ji,jj,2) + (2 * e3t(ji,jj,1,Kmm) * phioc(ji,jj)/rho0) / ( p_avm(ji,jj,1) + p_avm(ji,jj,2) ) 419 zdiag(ji,jj,2) = zdiag(ji,jj,2) + zd_lw(ji,jj,2) 420 zdiag(ji,jj,1) = 1._wp 421 zd_lw(ji,jj,2) = 0._wp 422 END_2D 423 424 END SELECT 425 426 ENDIF 427 ! 353 428 ! !* Matrix inversion from level 2 (tke prescribed at level 1) 354 DO_3D( 0, 0, 0, 0, 3, jpkm1 ) ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1429 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 355 430 zdiag(ji,jj,jk) = zdiag(ji,jj,jk) - zd_lw(ji,jj,jk) * zd_up(ji,jj,jk-1) / zdiag(ji,jj,jk-1) 356 431 END_3D 357 DO_2D( 0, 0, 0, 0 ) ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1 358 zd_lw(ji,jj,2) = en(ji,jj,2) - zd_lw(ji,jj,2) * en(ji,jj,1) ! Surface boudary conditions on tke 359 END_2D 360 DO_3D( 0, 0, 0, 0, 3, jpkm1 ) 432 !XC : commented to allow for neumann boundary condition 433 ! DO_2D( 0, 0, 0, 0 ) 434 ! zd_lw(ji,jj,2) = en(ji,jj,2) - zd_lw(ji,jj,2) * en(ji,jj,1) ! Surface boudary conditions on tke 435 ! END_2D 436 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 361 437 zd_lw(ji,jj,jk) = en(ji,jj,jk) - zd_lw(ji,jj,jk) / zdiag(ji,jj,jk-1) *zd_lw(ji,jj,jk-1) 362 438 END_3D … … 460 536 zmxlm(:,:,:) = rmxl_min 461 537 zmxld(:,:,:) = rmxl_min 538 ! 539 IF(ln_sdw .AND. ln_mxhsw) THEN 540 zmxlm(:,:,1)= vkarmn * MAX ( 1.6 * hsw(:,:) , 0.02 ) ! surface mixing length = F(wave height) 541 ! from terray et al 1999 and mellor and blumberg 2004 it should be 0.85 and not 1.6 542 zcoef = vkarmn * ( (rn_ediff*rn_ediss)**0.25 ) / rn_ediff 543 zmxlm(:,:,1)= zcoef * MAX ( 1.6 * hsw(:,:) , 0.02 ) ! surface mixing length = F(wave height) 544 ELSE 462 545 ! 463 IF( ln_mxl0 ) THEN ! surface mixing length = F(stress) : l=vkarmn*2.e5*taum/(rho0*g)464 ! 465 zraug = vkarmn * 2.e5_wp / ( rho0 * grav )546 IF( ln_mxl0 ) THEN ! surface mixing length = F(stress) : l=vkarmn*2.e5*taum/(rho0*g) 547 ! 548 zraug = vkarmn * 2.e5_wp / ( rho0 * grav ) 466 549 #if ! defined key_si3 && ! defined key_cice 467 DO_2D( 0, 0, 0, 0 ) ! No sea-ice468 zmxlm(ji,jj,1) = zraug * taum(ji,jj) * tmask(ji,jj,1)469 END_2D550 DO_2D( 0, 0, 0, 0 ) ! No sea-ice 551 zmxlm(ji,jj,1) = zraug * taum(ji,jj) * tmask(ji,jj,1) 552 END_2D 470 553 #else 471 SELECT CASE( nn_mxlice ) ! Type of scaling under sea-ice 472 ! 473 CASE( 0 ) ! No scaling under sea-ice 554 SELECT CASE( nn_mxlice ) ! Type of scaling under sea-ice 555 ! 556 CASE( 0 ) ! No scaling under sea-ice 557 DO_2D( 0, 0, 0, 0 ) 558 zmxlm(ji,jj,1) = zraug * taum(ji,jj) * tmask(ji,jj,1) 559 END_2D 560 ! 561 CASE( 1 ) ! scaling with constant sea-ice thickness 562 DO_2D( 0, 0, 0, 0 ) 563 zmxlm(ji,jj,1) = ( ( 1._wp - fr_i(ji,jj) ) * zraug * taum(ji,jj) + & 564 & fr_i(ji,jj) * rn_mxlice ) * tmask(ji,jj,1) 565 END_2D 566 ! 567 CASE( 2 ) ! scaling with mean sea-ice thickness 568 DO_2D( 0, 0, 0, 0 ) 569 #if defined key_si3 570 zmxlm(ji,jj,1) = ( ( 1._wp - fr_i(ji,jj) ) * zraug * taum(ji,jj) + & 571 & fr_i(ji,jj) * hm_i(ji,jj) * 2._wp ) * tmask(ji,jj,1) 572 #elif defined key_cice 573 zmaxice = MAXVAL( h_i(ji,jj,:) ) 574 zmxlm(ji,jj,1) = ( ( 1._wp - fr_i(ji,jj) ) * zraug * taum(ji,jj) + & 575 & fr_i(ji,jj) * zmaxice ) * tmask(ji,jj,1) 576 #endif 577 END_2D 578 ! 579 CASE( 3 ) ! scaling with max sea-ice thickness 580 DO_2D( 0, 0, 0, 0 ) 581 zmaxice = MAXVAL( h_i(ji,jj,:) ) 582 zmxlm(ji,jj,1) = ( ( 1._wp - fr_i(ji,jj) ) * zraug * taum(ji,jj) + & 583 & fr_i(ji,jj) * zmaxice ) * tmask(ji,jj,1) 584 END_2D 585 ! 586 END SELECT 587 #endif 588 ! 474 589 DO_2D( 0, 0, 0, 0 ) 475 zmxlm(ji,jj,1) = zraug * taum(ji,jj) * tmask(ji,jj,1)590 zmxlm(ji,jj,1) = MAX( rn_mxl0, zmxlm(ji,jj,1) ) 476 591 END_2D 477 592 ! 478 CASE( 1 ) ! scaling with constant sea-ice thickness 479 DO_2D( 0, 0, 0, 0 ) 480 zmxlm(ji,jj,1) = ( ( 1._wp - fr_i(ji,jj) ) * zraug * taum(ji,jj) + & 481 & fr_i(ji,jj) * rn_mxlice ) * tmask(ji,jj,1) 482 END_2D 483 ! 484 CASE( 2 ) ! scaling with mean sea-ice thickness 485 DO_2D( 0, 0, 0, 0 ) 486 #if defined key_si3 487 zmxlm(ji,jj,1) = ( ( 1._wp - fr_i(ji,jj) ) * zraug * taum(ji,jj) + & 488 & fr_i(ji,jj) * hm_i(ji,jj) * 2._wp ) * tmask(ji,jj,1) 489 #elif defined key_cice 490 zmaxice = MAXVAL( h_i(ji,jj,:) ) 491 zmxlm(ji,jj,1) = ( ( 1._wp - fr_i(ji,jj) ) * zraug * taum(ji,jj) + & 492 & fr_i(ji,jj) * zmaxice ) * tmask(ji,jj,1) 493 #endif 494 END_2D 495 ! 496 CASE( 3 ) ! scaling with max sea-ice thickness 497 DO_2D( 0, 0, 0, 0 ) 498 zmaxice = MAXVAL( h_i(ji,jj,:) ) 499 zmxlm(ji,jj,1) = ( ( 1._wp - fr_i(ji,jj) ) * zraug * taum(ji,jj) + & 500 & fr_i(ji,jj) * zmaxice ) * tmask(ji,jj,1) 501 END_2D 502 ! 503 END SELECT 504 #endif 505 ! 506 DO_2D( 0, 0, 0, 0 ) 507 zmxlm(ji,jj,1) = MAX( rn_mxl0, zmxlm(ji,jj,1) ) 508 END_2D 509 ! 510 ELSE 511 zmxlm(:,:,1) = rn_mxl0 512 ENDIF 513 593 ELSE 594 zmxlm(:,:,1) = rn_mxl0 595 ENDIF 596 ENDIF 514 597 ! 515 598 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) … … 624 707 & rn_mxl0 , nn_mxlice, rn_mxlice, & 625 708 & nn_pdl , ln_lc , rn_lc , & 626 & nn_etau , nn_htau , rn_efr , nn_eice 709 & nn_etau , nn_htau , rn_efr , nn_eice , & 710 & nn_bc_surf, nn_bc_bot, ln_mxhsw 627 711 !!---------------------------------------------------------------------- 628 712 ! … … 666 750 WRITE(numout,*) ' Langmuir cells parametrization ln_lc = ', ln_lc 667 751 WRITE(numout,*) ' coef to compute vertical velocity of LC rn_lc = ', rn_lc 752 IF ( cpl_phioc .and. ln_phioc ) THEN 753 SELECT CASE( nn_bc_surf) ! Type of scaling under sea-ice 754 CASE( 0 ) ; WRITE(numout,*) ' nn_bc_surf=0 ==>>> DIRICHLET SBC using surface TKE flux from waves' 755 CASE( 1 ) ; WRITE(numout,*) ' nn_bc_surf=1 ==>>> NEUMANN SBC using surface TKE flux from waves' 756 END SELECT 757 ENDIF 668 758 WRITE(numout,*) ' test param. to add tke induced by wind nn_etau = ', nn_etau 669 759 WRITE(numout,*) ' type of tke penetration profile nn_htau = ', nn_htau … … 721 811 CALL tke_rst( nit000, 'READ' ) ! (en, avt_k, avm_k, dissl) 722 812 ! 723 IF( lwxios ) THEN724 CALL iom_set_rstw_var_active('en')725 CALL iom_set_rstw_var_active('avt_k')726 CALL iom_set_rstw_var_active('avm_k')727 CALL iom_set_rstw_var_active('dissl')728 ENDIF729 813 END SUBROUTINE zdf_tke_init 730 814 … … 758 842 ! 759 843 IF( MIN( id1, id2, id3, id4 ) > 0 ) THEN ! fields exist 760 CALL iom_get( numror, jpdom_auto, 'en' , en , ldxios = lrxios)761 CALL iom_get( numror, jpdom_auto, 'avt_k', avt_k , ldxios = lrxios)762 CALL iom_get( numror, jpdom_auto, 'avm_k', avm_k , ldxios = lrxios)763 CALL iom_get( numror, jpdom_auto, 'dissl', dissl , ldxios = lrxios)844 CALL iom_get( numror, jpdom_auto, 'en' , en ) 845 CALL iom_get( numror, jpdom_auto, 'avt_k', avt_k ) 846 CALL iom_get( numror, jpdom_auto, 'avm_k', avm_k ) 847 CALL iom_get( numror, jpdom_auto, 'dissl', dissl ) 764 848 ELSE ! start TKE from rest 765 849 IF(lwp) WRITE(numout,*) … … 780 864 ! ! ------------------- 781 865 IF(lwp) WRITE(numout,*) '---- tke_rst ----' 782 IF( lwxios ) CALL iom_swap( cwxios_context ) 783 CALL iom_rstput( kt, nitrst, numrow, 'en' , en , ldxios = lwxios ) 784 CALL iom_rstput( kt, nitrst, numrow, 'avt_k', avt_k, ldxios = lwxios ) 785 CALL iom_rstput( kt, nitrst, numrow, 'avm_k', avm_k, ldxios = lwxios ) 786 CALL iom_rstput( kt, nitrst, numrow, 'dissl', dissl, ldxios = lwxios ) 787 IF( lwxios ) CALL iom_swap( cxios_context ) 866 CALL iom_rstput( kt, nitrst, numrow, 'en' , en ) 867 CALL iom_rstput( kt, nitrst, numrow, 'avt_k', avt_k ) 868 CALL iom_rstput( kt, nitrst, numrow, 'avm_k', avm_k ) 869 CALL iom_rstput( kt, nitrst, numrow, 'dissl', dissl ) 788 870 ! 789 871 ENDIF
Note: See TracChangeset
for help on using the changeset viewer.