- Timestamp:
- 2020-11-27T17:26:33+01:00 (4 years ago)
- Location:
- NEMO/branches/2020/tickets_icb_1900
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2020/tickets_icb_1900
- Property svn:externals
-
NEMO/branches/2020/tickets_icb_1900/src/OCE/ZDF/zdfosm.F90
r13237 r13899 300 300 zz0 = rn_abs ! surface equi-partition in 2-bands 301 301 zz1 = 1. - rn_abs 302 DO_2D _00_00302 DO_2D( 0, 0, 0, 0 ) 303 303 ! Surface downward irradiance (so always +ve) 304 304 zrad0(ji,jj) = qsr(ji,jj) * r1_rho0_rcp … … 310 310 END_2D 311 311 ! Turbulent surface fluxes and fluxes averaged over depth of the OSBL 312 DO_2D _00_00312 DO_2D( 0, 0, 0, 0 ) 313 313 zthermal = rab_n(ji,jj,1,jp_tem) 314 314 zbeta = rab_n(ji,jj,1,jp_sal) … … 337 337 ! Assume constant La#=0.3 338 338 CASE(0) 339 DO_2D _00_00339 DO_2D( 0, 0, 0, 0 ) 340 340 zus_x = zcos_wind(ji,jj) * zustar(ji,jj) / 0.3**2 341 341 zus_y = zsin_wind(ji,jj) * zustar(ji,jj) / 0.3**2 … … 345 345 ! Assume Pierson-Moskovitz wind-wave spectrum 346 346 CASE(1) 347 DO_2D _00_00347 DO_2D( 0, 0, 0, 0 ) 348 348 ! Use wind speed wndm included in sbc_oce module 349 349 zustke(ji,jj) = MAX ( 0.016 * wndm(ji,jj), 1.0e-8 ) … … 353 353 CASE(2) 354 354 zfac = 2.0_wp * rpi / 16.0_wp 355 DO_2D _00_00355 DO_2D( 0, 0, 0, 0 ) 356 356 ! The Langmur number from the ECMWF model appears to give La<0.3 for wind-driven seas. 357 357 ! The coefficient 0.8 gives La=0.3 in this situation. … … 366 366 ! Langmuir velocity scale (zwstrl), La # (zla) 367 367 ! mixed scale (zvstr), convective velocity scale (zwstrc) 368 DO_2D _00_00368 DO_2D( 0, 0, 0, 0 ) 369 369 ! Langmuir velocity scale (zwstrl), at T-point 370 370 zwstrl(ji,jj) = ( zustar(ji,jj) * zustar(ji,jj) * zustke(ji,jj) )**pthird … … 402 402 hbl(:,:) = MAX(hbl(:,:), gdepw(:,:,3,Kmm) ) 403 403 ibld(:,:) = 3 404 DO_3D _00_00(4, jpkm1 )404 DO_3D( 0, 0, 0, 0, 4, jpkm1 ) 405 405 IF ( hbl(ji,jj) >= gdepw(ji,jj,jk,Kmm) ) THEN 406 406 ibld(ji,jj) = MIN(mbkt(ji,jj), jk) … … 408 408 END_3D 409 409 410 DO_2D _00_00410 DO_2D( 0, 0, 0, 0 ) 411 411 zthermal = rab_n(ji,jj,1,jp_tem) !ideally use ibld not 1?? 412 412 zbeta = rab_n(ji,jj,1,jp_sal) … … 478 478 zdhdt(:,:) = MIN(zdhdt(:,:), (zhbl_t(:,:) - hbl(:,:))/rn_Dt + ww(ji,jj,ibld(ji,jj))) ! adjustment to represent limiting by ocean bottom 479 479 480 DO_3D _00_00(4, jpkm1 )480 DO_3D( 0, 0, 0, 0, 4, jpkm1 ) 481 481 IF ( zhbl_t(ji,jj) >= gdepw(ji,jj,jk,Kmm) ) THEN 482 482 ibld(ji,jj) = MIN(mbkt(ji,jj), jk) … … 487 487 ! Step through model levels taking account of buoyancy change to determine the effect on dhdt 488 488 ! 489 DO_2D _00_00489 DO_2D( 0, 0, 0, 0 ) 490 490 IF ( ibld(ji,jj) - imld(ji,jj) > 1 ) THEN 491 491 ! … … 552 552 ! Consider later combining this into the loop above and looking for columns 553 553 ! where the index for base of the boundary layer have changed 554 DO_2D _00_00554 DO_2D( 0, 0, 0, 0 ) 555 555 zthermal = rab_n(ji,jj,1,jp_tem) !ideally use ibld not 1?? 556 556 zbeta = rab_n(ji,jj,1,jp_sal) … … 635 635 ! Average over the depth of the mixed layer in the convective boundary layer 636 636 ! Also calculate entrainment fluxes for temperature and salinity 637 DO_2D _00_00637 DO_2D( 0, 0, 0, 0 ) 638 638 zthermal = rab_n(ji,jj,1,jp_tem) !ideally use ibld not 1?? 639 639 zbeta = rab_n(ji,jj,1,jp_sal) … … 705 705 ! 706 706 707 DO_2D _00_00707 DO_2D( 0, 0, 0, 0 ) 708 708 ztemp = zu_ml(ji,jj) 709 709 zu_ml(ji,jj) = zu_ml(ji,jj) * zcos_wind(ji,jj) + zv_ml(ji,jj) * zsin_wind(ji,jj) … … 723 723 zuw_bse = 0._wp 724 724 zvw_bse = 0._wp 725 DO_2D _00_00725 DO_2D( 0, 0, 0, 0 ) 726 726 727 727 IF ( lconv(ji,jj) ) THEN … … 740 740 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 741 741 742 DO_2D _00_00742 DO_2D( 0, 0, 0, 0 ) 743 743 ! 744 744 IF ( lconv (ji,jj) ) THEN … … 788 788 END_2D 789 789 ! 790 DO_2D _00_00790 DO_2D( 0, 0, 0, 0 ) 791 791 ! 792 792 IF ( lconv (ji,jj) ) THEN … … 832 832 ! zvisml_sc = zwstrl * zhbl * EXP ( -( zhol / 0.183_wp )**2 ) 833 833 ! ENDWHERE 834 DO_2D _00_00834 DO_2D( 0, 0, 0, 0 ) 835 835 IF ( lconv(ji,jj) ) THEN 836 836 zdifml_sc(ji,jj) = zhml(ji,jj) * ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird … … 846 846 END_2D 847 847 ! 848 DO_2D _00_00848 DO_2D( 0, 0, 0, 0 ) 849 849 IF ( lconv(ji,jj) ) THEN 850 850 DO jk = 2, imld(ji,jj) ! mixed layer diffusivity … … 896 896 897 897 898 DO_2D _00_00898 DO_2D( 0, 0, 0, 0 ) 899 899 IF ( lconv(ji,jj) ) THEN 900 900 DO jk = 2, imld(ji,jj) … … 929 929 ENDWHERE 930 930 931 DO_2D _00_00931 DO_2D( 0, 0, 0, 0 ) 932 932 IF ( lconv(ji,jj) ) THEN 933 933 DO jk = 2, imld(ji,jj) … … 961 961 ENDWHERE 962 962 963 DO_2D _00_00963 DO_2D( 0, 0, 0, 0 ) 964 964 IF (lconv(ji,jj) ) THEN 965 965 DO jk = 2, imld(ji,jj) … … 993 993 ENDWHERE 994 994 995 DO_2D _00_00995 DO_2D( 0, 0, 0, 0 ) 996 996 IF ( lconv(ji,jj) ) THEN 997 997 DO jk = 2 , imld(ji,jj) … … 1021 1021 ENDWHERE 1022 1022 1023 DO_2D _00_001023 DO_2D( 0, 0, 0, 0 ) 1024 1024 IF ( lconv(ji,jj) ) THEN 1025 1025 DO jk = 2, imld(ji,jj) … … 1058 1058 ENDWHERE 1059 1059 1060 DO_2D _00_001060 DO_2D( 0, 0, 0, 0 ) 1061 1061 IF ( lconv(ji,jj) ) THEN 1062 1062 DO jk = 2, imld(ji,jj) … … 1093 1093 ! Make surface forced velocity non-gradient terms go to zero at the base of the mixed layer. 1094 1094 1095 DO_2D _00_001095 DO_2D( 0, 0, 0, 0 ) 1096 1096 IF ( lconv(ji,jj) ) THEN 1097 1097 DO jk = 2, ibld(ji,jj) … … 1122 1122 ! Temporary fix to avoid instabilities when zdb_bl becomes very very small 1123 1123 zsc_uw_1 = 0._wp ! 50.0 * zla**(8.0/3.0) * zustar**2 * zhbl / ( zdb_bl + epsln ) 1124 DO_2D _00_001124 DO_2D( 0, 0, 0, 0 ) 1125 1125 DO jk= 2, ibld(ji,jj) 1126 1126 znd = gdepw(ji,jj,jk,Kmm) / zhbl(ji,jj) … … 1135 1135 ! Entrainment contribution. 1136 1136 1137 DO_2D _00_001137 DO_2D( 0, 0, 0, 0 ) 1138 1138 IF ( lconv(ji,jj) ) THEN 1139 1139 DO jk = 1, imld(ji,jj) - 1 … … 1170 1170 ! rotate non-gradient velocity terms back to model reference frame 1171 1171 1172 DO_2D _00_001172 DO_2D( 0, 0, 0, 0 ) 1173 1173 DO jk = 2, ibld(ji,jj) 1174 1174 ztemp = ghamu(ji,jj,jk) … … 1184 1184 ! KPP-style Ri# mixing 1185 1185 IF( ln_kpprimix) THEN 1186 DO_3D _10_10( 2, jpkm1)1186 DO_3D( 1, 0, 1, 0, 2, jpkm1 ) !* Shear production at uw- and vw-points (energy conserving form) 1187 1187 z3du(ji,jj,jk) = 0.5 * ( uu(ji,jj,jk-1,Kmm) - uu(ji ,jj,jk,Kmm) ) & 1188 1188 & * ( uu(ji,jj,jk-1,Kbb) - uu(ji ,jj,jk,Kbb) ) * wumask(ji,jj,jk) & … … 1193 1193 END_3D 1194 1194 ! 1195 DO_3D _00_00(2, jpkm1 )1195 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 1196 1196 ! ! shear prod. at w-point weightened by mask 1197 1197 zesh2 = ( z3du(ji-1,jj,jk) + z3du(ji,jj,jk) ) / MAX( 1._wp , umask(ji-1,jj,jk) + umask(ji,jj,jk) ) & … … 1204 1204 END_3D 1205 1205 1206 DO_2D _00_001206 DO_2D( 0, 0, 0, 0 ) 1207 1207 DO jk = ibld(ji,jj) + 1, jpkm1 1208 1208 zdiffut(ji,jj,jk) = zrimix(ji,jj,jk)*rn_difri … … 1215 1215 ! KPP-style set diffusivity large if unstable below BL 1216 1216 IF( ln_convmix) THEN 1217 DO_2D _00_001217 DO_2D( 0, 0, 0, 0 ) 1218 1218 DO jk = ibld(ji,jj) + 1, jpkm1 1219 1219 IF( MIN( rn2(ji,jj,jk), rn2b(ji,jj,jk) ) <= -1.e-12 ) zdiffut(ji,jj,jk) = rn_difconv … … 1227 1227 ! GN 25/8: need to change tmask --> wmask 1228 1228 1229 DO_3D _00_00(2, jpkm1 )1229 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 1230 1230 p_avt(ji,jj,jk) = MAX( zdiffut(ji,jj,jk), avtb(jk) ) * tmask(ji,jj,jk) 1231 1231 p_avm(ji,jj,jk) = MAX( zviscos(ji,jj,jk), avmb(jk) ) * tmask(ji,jj,jk) … … 1234 1234 CALL lbc_lnk_multi( 'zdfosm', p_avt, 'W', 1.0_wp , p_avm, 'W', 1.0_wp, & 1235 1235 & ghamu, 'W', 1.0_wp , ghamv, 'W', 1.0_wp ) 1236 DO_3D _00_00(2, jpkm1 )1236 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 1237 1237 ghamu(ji,jj,jk) = ( ghamu(ji,jj,jk) + ghamu(ji+1,jj,jk) ) & 1238 1238 & / MAX( 1., tmask(ji,jj,jk) + tmask (ji + 1,jj,jk) ) * umask(ji,jj,jk) … … 1348 1348 ENDIF 1349 1349 1350 1351 ! ! Check wave coupling settings ! 1352 ! ! Further work needed - see ticket #2447 ! 1353 IF( nn_osm_wave == 2 ) THEN 1354 IF (.NOT. ( ln_wave .AND. ln_sdw )) & 1355 & CALL ctl_stop( 'zdf_osm_init : ln_zdfosm and nn_osm_wave=2, ln_wave and ln_sdw must be true' ) 1356 END IF 1357 1350 1358 ! ! allocate zdfosm arrays 1351 1359 IF( zdf_osm_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'zdf_osm_init : unable to allocate arrays' ) … … 1387 1395 etmean(:,:,:) = 0.e0 1388 1396 1389 DO_3D _00_00(1, jpkm1 )1397 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 1390 1398 etmean(ji,jj,jk) = tmask(ji,jj,jk) & 1391 1399 & / MAX( 1., umask(ji-1,jj ,jk) + umask(ji,jj,jk) & … … 1401 1409 etmean(:,:,:) = 0.e0 1402 1410 1403 DO_3D _00_00(1, jpkm1 )1411 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 1404 1412 etmean(ji,jj,jk) = tmask(ji, jj,jk) & 1405 1413 & / MAX( 1., 2.* tmask(ji,jj,jk) & … … 1466 1474 id1 = iom_varid( numror, 'wn' , ldstop = .FALSE. ) 1467 1475 IF( id1 > 0 ) THEN ! 'wn' exists; read 1468 CALL iom_get( numror, jpdom_auto glo, 'wn', ww, ldxios = lrxios )1476 CALL iom_get( numror, jpdom_auto, 'wn', ww, ldxios = lrxios ) 1469 1477 WRITE(numout,*) ' ===>>>> : ww read from restart file' 1470 1478 ELSE … … 1475 1483 id2 = iom_varid( numror, 'hbli' , ldstop = .FALSE. ) 1476 1484 IF( id1 > 0 .AND. id2 > 0) THEN ! 'hbl' exists; read and return 1477 CALL iom_get( numror, jpdom_auto glo, 'hbl' , hbl , ldxios = lrxios )1478 CALL iom_get( numror, jpdom_auto glo, 'hbli', hbli, ldxios = lrxios )1485 CALL iom_get( numror, jpdom_auto, 'hbl' , hbl , ldxios = lrxios ) 1486 CALL iom_get( numror, jpdom_auto, 'hbli', hbli, ldxios = lrxios ) 1479 1487 WRITE(numout,*) ' ===>>>> : hbl & hbli read from restart file' 1480 1488 RETURN … … 1508 1516 ! 1509 1517 hbl(:,:) = 0._wp ! here hbl used as a dummy variable, integrating vertically N^2 1510 DO_3D _11_11( 1, jpkm1 )1518 DO_3D( 1, 1, 1, 1, 1, jpkm1 ) ! Mixed layer level: w-level 1511 1519 ikt = mbkt(ji,jj) 1512 1520 hbl(ji,jj) = hbl(ji,jj) + MAX( rn2(ji,jj,jk) , 0._wp ) * e3w(ji,jj,jk,Kmm) … … 1514 1522 END_3D 1515 1523 ! 1516 DO_2D _11_111524 DO_2D( 1, 1, 1, 1 ) 1517 1525 iiki = imld_rst(ji,jj) 1518 1526 hbl (ji,jj) = gdepw(ji,jj,iiki ,Kmm) * ssmask(ji,jj) ! Turbocline depth … … 1553 1561 1554 1562 ! add non-local temperature and salinity flux 1555 DO_3D _00_00(1, jpkm1 )1563 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 1556 1564 pts(ji,jj,jk,jp_tem,Krhs) = pts(ji,jj,jk,jp_tem,Krhs) & 1557 1565 & - ( ghamt(ji,jj,jk ) & … … 1621 1629 !code saving tracer trends removed, replace with trdmxl_oce 1622 1630 1623 DO_3D _00_00( 1, jpkm1 )1631 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) ! add non-local u and v fluxes 1624 1632 puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) & 1625 1633 & - ( ghamu(ji,jj,jk ) &
Note: See TracChangeset
for help on using the changeset viewer.