Changeset 8540
- Timestamp:
- 2017-09-19T10:01:20+02:00 (7 years ago)
- Location:
- branches/UKMO/dev_r5518_GO6_package_lim3/NEMOGCM
- Files:
-
- 32 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/UKMO/dev_r5518_GO6_package_lim3/NEMOGCM/CONFIG/SHARED/field_def.xml
r8308 r8540 373 373 <field id="vfxsnw" long_name="snw melt/growth" unit="m/day" /> 374 374 <field id="vfxsub" long_name="snw sublimation" unit="m/day" /> 375 <field id="vfxsub_err" long_name="excess of snw sublimation sent to ocean" unit="m/day" /> 375 376 <field id="vfxspr" long_name="snw precipitation on ice" unit="m/day" /> 376 377 <field id="vfxthin" long_name="daily thermo ice prod. for thin ice(<20cm) + open water" unit="m/day" /> … … 599 600 <!-- LIM3 scalar variables --> 600 601 601 <field_group id="SBC_scalar" grid_ref="scalar" >602 <field_group id="SBC_scalar" domain_ref="1point" > 602 603 <!-- available with ln_limdiaout --> 603 <field id="ibgvoltot" long_name="global mean ice volume" unit="km3" /> 604 <field id="sbgvoltot" long_name="global mean snow volume" unit="km3" /> 605 <field id="ibgarea" long_name="global mean ice area" unit="km2" /> 606 <field id="ibgsaline" long_name="global mean ice salinity" unit="0.001" /> 607 <field id="ibgtemper" long_name="global mean ice temperature" unit="degree_C" /> 608 <field id="ibgheatco" long_name="global mean ice heat content" unit="10^20J" /> 609 <field id="sbgheatco" long_name="global mean snow heat content" unit="10^20J" /> 610 <field id="ibgsaltco" long_name="global mean ice salt content" unit="0.001*km3" /> 611 612 <field id="ibgvfx" long_name="global mean volume flux (emp)" unit="m/day" /> 613 <field id="ibgvfxbog" long_name="global mean volume flux (bottom growth)" unit="m/day" /> 614 <field id="ibgvfxopw" long_name="global mean volume flux (open water growth)" unit="m/day" /> 615 <field id="ibgvfxsni" long_name="global mean volume flux (snow-ice growth)" unit="m/day" /> 616 <field id="ibgvfxdyn" long_name="global mean volume flux (dynamic growth)" unit="m/day" /> 617 <field id="ibgvfxbom" long_name="global mean volume flux (bottom melt)" unit="m/day" /> 618 <field id="ibgvfxsum" long_name="global mean volume flux (surface melt)" unit="m/day" /> 619 <field id="ibgvfxres" long_name="global mean volume flux (resultant)" unit="m/day" /> 620 <field id="ibgvfxspr" long_name="global mean volume flux (snow precip)" unit="m/day" /> 621 <field id="ibgvfxsnw" long_name="global mean volume flux (snow melt)" unit="m/day" /> 622 <field id="ibgvfxsub" long_name="global mean volume flux (snow sublimation)" unit="m/day" /> 623 624 <field id="ibgsfx" long_name="global mean salt flux (total)" unit="0.001*m/day" /> 625 <field id="ibgsfxbri" long_name="global mean salt flux (brines)" unit="0.001*m/day" /> 626 <field id="ibgsfxdyn" long_name="global mean salt flux (dynamic)" unit="0.001*m/day" /> 627 <field id="ibgsfxres" long_name="global mean salt flux (resultant)" unit="0.001*m/day" /> 628 <field id="ibgsfxbog" long_name="global mean salt flux (thermo)" unit="0.001*m/day" /> 629 <field id="ibgsfxopw" long_name="global mean salt flux (thermo)" unit="0.001*m/day" /> 630 <field id="ibgsfxsni" long_name="global mean salt flux (thermo)" unit="0.001*m/day" /> 631 <field id="ibgsfxbom" long_name="global mean salt flux (thermo)" unit="0.001*m/day" /> 632 <field id="ibgsfxsum" long_name="global mean salt flux (thermo)" unit="0.001*m/day" /> 633 <field id="ibgsfxsub" long_name="global mean salt flux (thermo)" unit="0.001*m/day" /> 634 635 <field id="ibghfxdhc" long_name="Heat content variation in snow and ice" unit="W" /> 636 <field id="ibghfxspr" long_name="Heat content of snow precip" unit="W" /> 637 638 <field id="ibghfxthd" long_name="heat fluxes from ice-ocean exchange during thermo" unit="W" /> 639 <field id="ibghfxsum" long_name="heat fluxes causing surface ice melt" unit="W" /> 640 <field id="ibghfxbom" long_name="heat fluxes causing bottom ice melt" unit="W" /> 641 <field id="ibghfxbog" long_name="heat fluxes causing bottom ice growth" unit="W" /> 642 <field id="ibghfxdif" long_name="heat fluxes causing ice temperature change" unit="W" /> 643 <field id="ibghfxopw" long_name="heat fluxes causing open water ice formation" unit="W" /> 644 <field id="ibghfxdyn" long_name="heat fluxes from ice-ocean exchange during dynamic" unit="W" /> 645 <field id="ibghfxres" long_name="heat fluxes from ice-ocean exchange during resultant" unit="W" /> 646 <field id="ibghfxsub" long_name="heat fluxes from sublimation" unit="W" /> 647 <field id="ibghfxsnw" long_name="heat fluxes from snow-ocean exchange" unit="W" /> 648 <field id="ibghfxout" long_name="non solar heat fluxes received by the ocean" unit="W" /> 649 <field id="ibghfxin" long_name="total heat fluxes at the ice surface" unit="W" /> 650 651 <field id="ibgfrcvol" long_name="global mean forcing volume (emp)" unit="km3" /> 652 <field id="ibgfrcsfx" long_name="global mean forcing salt (sfx)" unit="0.001*km3" /> 653 <field id="ibgvolgrm" long_name="global mean ice growth+melt volume" unit="km3" /> 604 <field id="ibgfrcvoltop" long_name="global mean ice/snow forcing at interface ice/snow-atm (volume equivalent ocean volume)" unit="km3" /> 605 <field id="ibgfrcvolbot" long_name="global mean ice/snow forcing at interface ice/snow-ocean (volume equivalent ocean volume)" unit="km3" /> 606 <field id="ibgfrctemtop" long_name="global mean heat on top of ice/snw/ocean-atm " unit="1e20J" /> 607 <field id="ibgfrctembot" long_name="global mean heat below ice (on top of ocean) " unit="1e20J" /> 608 <field id="ibgfrcsal" long_name="global mean ice/snow forcing (salt equivalent ocean volume)" unit="pss*km3" /> 609 <field id="ibgfrchfxtop" long_name="global mean heat flux on top of ice/snw/ocean-atm " unit="W/m2" /> 610 <field id="ibgfrchfxbot" long_name="global mean heat flux below ice (on top of ocean) " unit="W/m2" /> 611 612 <field id="ibgvolume" long_name="drift in ice/snow volume (equivalent ocean volume)" unit="km3" /> 613 <field id="ibgsaltco" long_name="drift in ice salt content (equivalent ocean volume)" unit="pss*km3" /> 614 <field id="ibgheatco" long_name="drift in ice/snow heat content" unit="1e20J" /> 615 <field id="ibgheatfx" long_name="drift in ice/snow heat flux" unit="W/m2" /> 616 617 <field id="ibgvol_tot" long_name="global mean ice volume" unit="km3" /> 618 <field id="sbgvol_tot" long_name="global mean snow volume" unit="km3" /> 619 <field id="ibgarea_tot" long_name="global mean ice area" unit="km2" /> 620 <field id="ibgsalt_tot" long_name="global mean ice salt content" unit="1e-3*km3" /> 621 <field id="ibgheat_tot" long_name="global mean ice heat content" unit="1e20J" /> 622 <field id="sbgheat_tot" long_name="global mean snow heat content" unit="1e20J" /> 654 623 </field_group> 655 624 -
branches/UKMO/dev_r5518_GO6_package_lim3/NEMOGCM/CONFIG/SHARED/namelist_ice_lim3_ref
r6498 r8540 69 69 rn_relast = 0.333 ! ratio of elastic timescale to ice time step: Telast = dt_ice * rn_relast 70 70 ! advised value: 1/3 (rn_nevp=120) or 1/9 (rn_nevp=300) 71 nn_ahi0 = 2 ! horizontal diffusivity computation72 ! 0: use rn_ahi0_ref73 ! 1: use rn_ahi0_ref x mean grid cell length / ( 2deg mean grid cell length )74 ! 2: use rn_ahi0_ref x grid cell length / ( 2deg mean grid cell length )75 rn_ahi0_ref = 350.0 ! horizontal sea ice diffusivity (m2/s)76 ! if nn_ahi0 > 0, rn_ahi0_ref is the reference value at a nominal 2 deg resolution77 71 / 78 72 !------------------------------------------------------------------------------ 79 73 &namicehdf ! Ice horizontal diffusion 80 74 !------------------------------------------------------------------------------ 75 nn_ahi0 = -1 ! horizontal diffusivity computation 76 ! -1: no diffusion (bypass limhdf) 77 ! 0: use rn_ahi0_ref 78 ! 1: use rn_ahi0_ref x mean grid cell length / ( 2deg mean grid cell length ) 79 ! 2: use rn_ahi0_ref x grid cell length / ( 2deg mean grid cell length ) 80 rn_ahi0_ref = 350.0 ! horizontal sea ice diffusivity (m2/s) 81 ! if nn_ahi0 > 0, rn_ahi0_ref is the reference value at a nominal 2 deg resolution 81 82 nn_convfrq = 5 ! convergence check frequency of the Crant-Nicholson scheme (perf. optimization) 82 83 / … … 98 99 ! 0: k = k0 + beta.S/T (Untersteiner, 1964) 99 100 ! 1: k = k0 + beta1.S/T - beta2.T (Pringle et al., 2007) 101 rn_cdsn = 0.31 ! thermal conductivity of the snow (0.31 W/m/K, Maykut and Untersteiner, 1971) 102 ! Obs: 0.1-0.5 (Lecomte et al, JAMES 2013) 100 103 nn_monocat = 0 ! virtual ITD mono-category parameterizations (1, jpl = 1 only) or not (0) 101 104 ! 2: simple piling instead of ridging --- temporary option -
branches/UKMO/dev_r5518_GO6_package_lim3/NEMOGCM/NEMO/LIM_SRC_3/dom_ice.F90
r6486 r8540 25 25 !!---------------------------------------------------------------------- 26 26 !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) 27 !! $Id $27 !! $Id: dom_ice.F90 5123 2015-03-04 16:06:03Z clem $ 28 28 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 29 29 !!---------------------------------------------------------------------- -
branches/UKMO/dev_r5518_GO6_package_lim3/NEMOGCM/NEMO/LIM_SRC_3/ice.F90
r7993 r8540 212 212 REAL(wp), PUBLIC :: rn_betas !: coef. for partitioning of snowfall between leads and sea ice 213 213 REAL(wp), PUBLIC :: rn_kappa_i !: coef. for the extinction of radiation Grenfell et al. (2006) [1/m] 214 REAL(wp), PUBLIC :: rn_cdsn !: thermal conductivity of the snow [W/m/K] 214 215 REAL(wp), PUBLIC :: nn_conv_dif !: maximal number of iterations for heat diffusion 215 216 REAL(wp), PUBLIC :: rn_terr_dif !: maximal tolerated error (C) for heat diffusion … … 243 244 ! 244 245 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: sist !: Average Sea-Ice Surface Temperature [Kelvin] 245 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: icethi !: total ice thickness (for all categories) (diag only)246 246 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: t_bo !: Sea-Ice bottom temperature [Kelvin] 247 247 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: frld !: Leads fraction = 1 - ice fraction … … 320 320 ! ! this is an extensive variable that has to be transported 321 321 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: o_i !: Sea-Ice Age (days) 322 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: ov_i !: Sea-Ice Age times volume per area (days.m)323 322 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: oa_i !: Sea-Ice Age times ice area (days) 323 324 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: bv_i !: brine volume 324 325 325 326 !! Variables summed over all categories, or associated to all the ice in a single grid cell 326 327 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: u_ice, v_ice !: components of the ice velocity (m/s) 327 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: tio_u, tio_v !: components of the ice-ocean stress (N/m2)328 328 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: vt_i , vt_s !: ice and snow total volume per unit area (m) 329 329 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: at_i !: ice total fractional area (ice concentration) 330 330 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ato_i !: =1-at_i ; total open water fractional area 331 331 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: et_i , et_s !: ice and snow total heat content 332 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ot_i !: mean age over all categories 333 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: tm_i !: mean ice temperature over all categories 334 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: bv_i !: brine volume averaged over all categories 335 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: smt_i !: mean sea ice salinity averaged over all categories [PSU] 332 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: tm_i !: mean ice temperature over all categories 333 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: bvm_i !: brine volume averaged over all categories 334 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: smt_i !: mean sea ice salinity averaged over all categories [PSU] 335 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: tm_su !: mean surface temperature over all categories 336 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: htm_i !: mean ice thickness over all categories 337 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: htm_s !: mean snow thickness over all categories 338 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: om_i !: mean ice age over all categories 336 339 337 340 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: t_s !: Snow temperatures [K] … … 405 408 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: diag_vice !: ice volume variation [m/s] 406 409 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: diag_vsnw !: snw volume variation [m/s] 410 407 411 ! 408 412 !!---------------------------------------------------------------------- 409 413 !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2010) 410 !! $Id $414 !! $Id: ice.F90 7814 2017-03-20 16:21:42Z clem $ 411 415 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 412 416 !!---------------------------------------------------------------------- … … 415 419 FUNCTION ice_alloc() 416 420 !!----------------------------------------------------------------- 417 !! *** Routine ice_alloc _2***421 !! *** Routine ice_alloc *** 418 422 !!----------------------------------------------------------------- 419 423 INTEGER :: ice_alloc … … 435 439 436 440 ii = ii + 1 437 ALLOCATE( sist (jpi,jpj) , icethi (jpi,jpj) ,t_bo (jpi,jpj) , &441 ALLOCATE( sist (jpi,jpj) , t_bo (jpi,jpj) , & 438 442 & frld (jpi,jpj) , pfrld (jpi,jpj) , phicif (jpi,jpj) , & 439 443 & wfx_snw(jpi,jpj) , wfx_ice(jpi,jpj) , wfx_sub(jpi,jpj) , & … … 456 460 & v_s (jpi,jpj,jpl) , ht_s (jpi,jpj,jpl) , t_su (jpi,jpj,jpl) , & 457 461 & sm_i (jpi,jpj,jpl) , smv_i(jpi,jpj,jpl) , o_i (jpi,jpj,jpl) , & 458 & o v_i (jpi,jpj,jpl) , oa_i (jpi,jpj,jpl), STAT=ierr(ii) )459 ii = ii + 1 460 ALLOCATE( u_ice(jpi,jpj) , v_ice(jpi,jpj) , tio_u(jpi,jpj) , tio_v(jpi,jpj) ,&462 & oa_i (jpi,jpj,jpl) , bv_i (jpi,jpj,jpl) , STAT=ierr(ii) ) 463 ii = ii + 1 464 ALLOCATE( u_ice(jpi,jpj) , v_ice(jpi,jpj) , & 461 465 & vt_i (jpi,jpj) , vt_s (jpi,jpj) , at_i (jpi,jpj) , ato_i(jpi,jpj) , & 462 & et_i (jpi,jpj) , et_s (jpi,jpj) , ot_i (jpi,jpj) , tm_i (jpi,jpj) , & 463 & bv_i (jpi,jpj) , smt_i(jpi,jpj) , STAT=ierr(ii) ) 466 & et_i (jpi,jpj) , et_s (jpi,jpj) , tm_i (jpi,jpj) , bvm_i(jpi,jpj) , & 467 & smt_i(jpi,jpj) , tm_su(jpi,jpj) , htm_i(jpi,jpj) , htm_s(jpi,jpj) , & 468 & om_i (jpi,jpj) , STAT=ierr(ii) ) 464 469 ii = ii + 1 465 470 ALLOCATE( t_s(jpi,jpj,nlay_s,jpl) , e_s(jpi,jpj,nlay_s,jpl) , STAT=ierr(ii) ) … … 501 506 502 507 ice_alloc = MAXVAL( ierr(:) ) 503 IF( ice_alloc /= 0 ) CALL ctl_warn('ice_alloc _2: failed to allocate arrays.')508 IF( ice_alloc /= 0 ) CALL ctl_warn('ice_alloc: failed to allocate arrays.') 504 509 ! 505 510 END FUNCTION ice_alloc -
branches/UKMO/dev_r5518_GO6_package_lim3/NEMOGCM/NEMO/LIM_SRC_3/limadv.F90
r6486 r8540 35 35 !!---------------------------------------------------------------------- 36 36 !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) 37 !! $Id $37 !! $Id: limadv.F90 5429 2015-06-16 09:57:07Z smasson $ 38 38 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 39 39 !!---------------------------------------------------------------------- -
branches/UKMO/dev_r5518_GO6_package_lim3/NEMOGCM/NEMO/LIM_SRC_3/limcons.F90
r6498 r8540 37 37 !!---------------------------------------------------------------------- 38 38 !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) 39 !! $Id $39 !! $Id: limcons.F90 6963 2016-09-30 12:40:04Z clem $ 40 40 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 41 41 !!---------------------------------------------------------------------- … … 288 288 #if ! defined key_bdy 289 289 ! heat flux 290 zhfx = glob_sum( ( hfx_in - hfx_out - diag_heat - diag_trp_ei - diag_trp_es - SUM( qevap_ice * a_i_b, dim=3 ) ) & 291 & * e12t * tmask(:,:,1) * zconv ) 290 zhfx = glob_sum( ( hfx_in - hfx_out - diag_heat - diag_trp_ei - diag_trp_es & 291 ! & - SUM( qevap_ice * a_i_b, dim=3 ) & !!clem: I think this line must be commented (but need check) 292 & ) * e12t * tmask(:,:,1) * zconv ) 292 293 ! salt flux 293 294 zsfx = glob_sum( ( sfx + diag_smvi ) * e12t * tmask(:,:,1) * zconv ) * rday -
branches/UKMO/dev_r5518_GO6_package_lim3/NEMOGCM/NEMO/LIM_SRC_3/limctl.F90
r6486 r8540 40 40 !!---------------------------------------------------------------------- 41 41 !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) 42 !! $Id $42 !! $Id: limctl.F90 5043 2015-01-28 16:44:18Z clem $ 43 43 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 44 44 !!---------------------------------------------------------------------- -
branches/UKMO/dev_r5518_GO6_package_lim3/NEMOGCM/NEMO/LIM_SRC_3/limdiahsb.F90
r6498 r8540 31 31 32 32 PUBLIC lim_diahsb ! routine called by ice_step.F90 33 34 real(wp) :: frc_sal, frc_vol ! global forcing trends 35 real(wp) :: bg_grme ! global ice growth+melt trends 36 33 PUBLIC lim_diahsb_init ! routine called in sbcice_lim.F90 34 35 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: vol_loc_ini, sal_loc_ini, tem_loc_ini ! initial volume, salt and heat contents 36 REAL(wp) :: frc_sal, frc_voltop, frc_volbot, frc_temtop, frc_tembot ! global forcing trends 37 37 38 !! * Substitutions 38 39 # include "vectopt_loop_substitute.h90" … … 40 41 !!---------------------------------------------------------------------- 41 42 !! NEMO/OPA 3.4 , NEMO Consortium (2012) 42 !! $Id $43 !! $Id: limdiahsb.F90 6963 2016-09-30 12:40:04Z clem $ 43 44 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 44 45 !!---------------------------------------------------------------------- … … 46 47 CONTAINS 47 48 48 SUBROUTINE lim_diahsb 49 SUBROUTINE lim_diahsb( kt ) 49 50 !!--------------------------------------------------------------------------- 50 51 !! *** ROUTINE lim_diahsb *** … … 53 54 !! 54 55 !!--------------------------------------------------------------------------- 56 INTEGER, INTENT(in) :: kt ! number of iteration 55 57 !! 56 real(wp) :: zbg_ivo, zbg_svo, zbg_are, zbg_sal ,zbg_tem ,zbg_ihc ,zbg_shc 57 real(wp) :: zbg_sfx, zbg_sfx_bri, zbg_sfx_bog, zbg_sfx_bom, zbg_sfx_sum, zbg_sfx_sni, & 58 & zbg_sfx_opw, zbg_sfx_res, zbg_sfx_dyn, zbg_sfx_sub 59 real(wp) :: zbg_vfx, zbg_vfx_bog, zbg_vfx_opw, zbg_vfx_sni, zbg_vfx_dyn 60 real(wp) :: zbg_vfx_bom, zbg_vfx_sum, zbg_vfx_res, zbg_vfx_spr, zbg_vfx_snw, zbg_vfx_sub 61 real(wp) :: zbg_hfx_dhc, zbg_hfx_spr 62 real(wp) :: zbg_hfx_res, zbg_hfx_sub, zbg_hfx_dyn, zbg_hfx_thd, zbg_hfx_snw, zbg_hfx_out, zbg_hfx_in 63 real(wp) :: zbg_hfx_sum, zbg_hfx_bom, zbg_hfx_bog, zbg_hfx_dif, zbg_hfx_opw 64 real(wp) :: z_frc_vol, z_frc_sal, z_bg_grme 65 real(wp) :: z1_area ! - - 66 REAL(wp) :: ztmp 58 real(wp) :: zbg_ivol, zbg_svol, zbg_area, zbg_isal, zbg_item ,zbg_stem 59 REAL(wp) :: z_frc_voltop, z_frc_volbot, z_frc_sal, z_frc_temtop, z_frc_tembot 60 REAL(wp) :: zdiff_vol, zdiff_sal, zdiff_tem 67 61 !!--------------------------------------------------------------------------- 68 62 IF( nn_timing == 1 ) CALL timing_start('lim_diahsb') 69 63 70 IF( numit == nstart ) CALL lim_diahsb_init 71 72 ! 1/area 73 z1_area = 1._wp / MAX( glob_sum( e12t(:,:) * tmask(:,:,1) ), epsi06 ) 74 75 rswitch = MAX( 0._wp , SIGN( 1._wp , glob_sum( e12t(:,:) * tmask(:,:,1) ) - epsi06 ) ) 76 ! ----------------------- ! 77 ! 1 - Content variations ! 78 ! ----------------------- ! 79 zbg_ivo = glob_sum( vt_i(:,:) * e12t(:,:) * tmask(:,:,1) ) ! volume ice 80 zbg_svo = glob_sum( vt_s(:,:) * e12t(:,:) * tmask(:,:,1) ) ! volume snow 81 zbg_are = glob_sum( at_i(:,:) * e12t(:,:) * tmask(:,:,1) ) ! area 82 zbg_sal = glob_sum( SUM( smv_i(:,:,:), dim=3 ) * e12t(:,:) * tmask(:,:,1) ) ! mean salt content 83 zbg_tem = glob_sum( ( tm_i(:,:) - rt0 ) * vt_i(:,:) * e12t(:,:) * tmask(:,:,1) ) ! mean temp content 84 85 !zbg_ihc = glob_sum( et_i(:,:) * e12t(:,:) * tmask(:,:,1) ) / MAX( zbg_ivo,epsi06 ) ! ice heat content 86 !zbg_shc = glob_sum( et_s(:,:) * e12t(:,:) * tmask(:,:,1) ) / MAX( zbg_svo,epsi06 ) ! snow heat content 87 88 ! Volume 89 ztmp = rswitch * z1_area * r1_rau0 * rday 90 zbg_vfx = ztmp * glob_sum( emp(:,:) * e12t(:,:) * tmask(:,:,1) ) 91 zbg_vfx_bog = ztmp * glob_sum( wfx_bog(:,:) * e12t(:,:) * tmask(:,:,1) ) 92 zbg_vfx_opw = ztmp * glob_sum( wfx_opw(:,:) * e12t(:,:) * tmask(:,:,1) ) 93 zbg_vfx_sni = ztmp * glob_sum( wfx_sni(:,:) * e12t(:,:) * tmask(:,:,1) ) 94 zbg_vfx_dyn = ztmp * glob_sum( wfx_dyn(:,:) * e12t(:,:) * tmask(:,:,1) ) 95 zbg_vfx_bom = ztmp * glob_sum( wfx_bom(:,:) * e12t(:,:) * tmask(:,:,1) ) 96 zbg_vfx_sum = ztmp * glob_sum( wfx_sum(:,:) * e12t(:,:) * tmask(:,:,1) ) 97 zbg_vfx_res = ztmp * glob_sum( wfx_res(:,:) * e12t(:,:) * tmask(:,:,1) ) 98 zbg_vfx_spr = ztmp * glob_sum( wfx_spr(:,:) * e12t(:,:) * tmask(:,:,1) ) 99 zbg_vfx_snw = ztmp * glob_sum( wfx_snw(:,:) * e12t(:,:) * tmask(:,:,1) ) 100 zbg_vfx_sub = ztmp * glob_sum( wfx_sub(:,:) * e12t(:,:) * tmask(:,:,1) ) 101 102 ! Salt 103 zbg_sfx = ztmp * glob_sum( sfx(:,:) * e12t(:,:) * tmask(:,:,1) ) 104 zbg_sfx_bri = ztmp * glob_sum( sfx_bri(:,:) * e12t(:,:) * tmask(:,:,1) ) 105 zbg_sfx_res = ztmp * glob_sum( sfx_res(:,:) * e12t(:,:) * tmask(:,:,1) ) 106 zbg_sfx_dyn = ztmp * glob_sum( sfx_dyn(:,:) * e12t(:,:) * tmask(:,:,1) ) 107 108 zbg_sfx_bog = ztmp * glob_sum( sfx_bog(:,:) * e12t(:,:) * tmask(:,:,1) ) 109 zbg_sfx_opw = ztmp * glob_sum( sfx_opw(:,:) * e12t(:,:) * tmask(:,:,1) ) 110 zbg_sfx_sni = ztmp * glob_sum( sfx_sni(:,:) * e12t(:,:) * tmask(:,:,1) ) 111 zbg_sfx_bom = ztmp * glob_sum( sfx_bom(:,:) * e12t(:,:) * tmask(:,:,1) ) 112 zbg_sfx_sum = ztmp * glob_sum( sfx_sum(:,:) * e12t(:,:) * tmask(:,:,1) ) 113 zbg_sfx_sub = ztmp * glob_sum( sfx_sub(:,:) * e12t(:,:) * tmask(:,:,1) ) 114 115 ! Heat budget 116 zbg_ihc = glob_sum( et_i(:,:) * e12t(:,:) * 1.e-20 ) ! ice heat content [1.e20 J] 117 zbg_shc = glob_sum( et_s(:,:) * e12t(:,:) * 1.e-20 ) ! snow heat content [1.e20 J] 118 zbg_hfx_dhc = glob_sum( diag_heat(:,:) * e12t(:,:) * tmask(:,:,1) ) ! [in W] 119 zbg_hfx_spr = glob_sum( hfx_spr(:,:) * e12t(:,:) * tmask(:,:,1) ) ! [in W] 120 121 zbg_hfx_thd = glob_sum( hfx_thd(:,:) * e12t(:,:) * tmask(:,:,1) ) ! [in W] 122 zbg_hfx_dyn = glob_sum( hfx_dyn(:,:) * e12t(:,:) * tmask(:,:,1) ) ! [in W] 123 zbg_hfx_res = glob_sum( hfx_res(:,:) * e12t(:,:) * tmask(:,:,1) ) ! [in W] 124 zbg_hfx_sub = glob_sum( hfx_sub(:,:) * e12t(:,:) * tmask(:,:,1) ) ! [in W] 125 zbg_hfx_snw = glob_sum( hfx_snw(:,:) * e12t(:,:) * tmask(:,:,1) ) ! [in W] 126 zbg_hfx_sum = glob_sum( hfx_sum(:,:) * e12t(:,:) * tmask(:,:,1) ) ! [in W] 127 zbg_hfx_bom = glob_sum( hfx_bom(:,:) * e12t(:,:) * tmask(:,:,1) ) ! [in W] 128 zbg_hfx_bog = glob_sum( hfx_bog(:,:) * e12t(:,:) * tmask(:,:,1) ) ! [in W] 129 zbg_hfx_dif = glob_sum( hfx_dif(:,:) * e12t(:,:) * tmask(:,:,1) ) ! [in W] 130 zbg_hfx_opw = glob_sum( hfx_opw(:,:) * e12t(:,:) * tmask(:,:,1) ) ! [in W] 131 zbg_hfx_out = glob_sum( hfx_out(:,:) * e12t(:,:) * tmask(:,:,1) ) ! [in W] 132 zbg_hfx_in = glob_sum( hfx_in(:,:) * e12t(:,:) * tmask(:,:,1) ) ! [in W] 133 134 ! --------------------------------------------- ! 135 ! 2 - Trends due to forcing and ice growth/melt ! 136 ! --------------------------------------------- ! 137 z_frc_vol = r1_rau0 * glob_sum( - emp(:,:) * e12t(:,:) * tmask(:,:,1) ) ! volume fluxes 138 z_frc_sal = r1_rau0 * glob_sum( sfx(:,:) * e12t(:,:) * tmask(:,:,1) ) ! salt fluxes 139 z_bg_grme = glob_sum( - ( wfx_bog(:,:) + wfx_opw(:,:) + wfx_sni(:,:) + wfx_dyn(:,:) + & 140 & wfx_bom(:,:) + wfx_sum(:,:) + wfx_res(:,:) + wfx_snw(:,:) + & 141 & wfx_sub(:,:) ) * e12t(:,:) * tmask(:,:,1) ) ! volume fluxes 142 ! 143 frc_vol = frc_vol + z_frc_vol * rdt_ice 144 frc_sal = frc_sal + z_frc_sal * rdt_ice 145 bg_grme = bg_grme + z_bg_grme * rdt_ice 64 ! ----------------------- ! 65 ! 1 - Contents ! 66 ! ----------------------- ! 67 zbg_ivol = glob_sum( vt_i(:,:) * e12t(:,:) * tmask(:,:,1) * 1.e-9 ) ! ice volume (km3) 68 zbg_svol = glob_sum( vt_s(:,:) * e12t(:,:) * tmask(:,:,1) * 1.e-9 ) ! snow volume (km3) 69 zbg_area = glob_sum( at_i(:,:) * e12t(:,:) * tmask(:,:,1) * 1.e-6 ) ! area (km2) 70 zbg_isal = glob_sum( SUM( smv_i(:,:,:), dim=3 ) * e12t(:,:) * tmask(:,:,1) * 1.e-9 ) ! salt content (pss*km3) 71 zbg_item = glob_sum( et_i * e12t(:,:) * tmask(:,:,1) * 1.e-20 ) ! heat content (1.e20 J) 72 zbg_stem = glob_sum( et_s * e12t(:,:) * tmask(:,:,1) * 1.e-20 ) ! heat content (1.e20 J) 146 73 147 ! difference 148 !frc_vol = zbg_ivo - frc_vol 149 !frc_sal = zbg_sal - frc_sal 150 151 ! ----------------------- ! 152 ! 3 - Diagnostics writing ! 153 ! ----------------------- ! 154 rswitch = MAX( 0._wp , SIGN( 1._wp , zbg_ivo - epsi06 ) ) 155 ! 156 IF( iom_use('ibgvoltot') ) & 157 CALL iom_put( 'ibgvoltot' , zbg_ivo * rhoic * r1_rau0 * 1.e-9 ) ! ice volume (km3 equivalent liquid) 158 IF( iom_use('sbgvoltot') ) & 159 CALL iom_put( 'sbgvoltot' , zbg_svo * rhosn * r1_rau0 * 1.e-9 ) ! snw volume (km3 equivalent liquid) 160 IF( iom_use('ibgarea') ) & 161 CALL iom_put( 'ibgarea' , zbg_are * 1.e-6 ) ! ice area (km2) 162 IF( iom_use('ibgsaline') ) & 163 CALL iom_put( 'ibgsaline' , rswitch * zbg_sal / MAX( zbg_ivo, epsi06 ) ) ! ice saline (psu) 164 IF( iom_use('ibgtemper') ) & 165 CALL iom_put( 'ibgtemper' , rswitch * zbg_tem / MAX( zbg_ivo, epsi06 ) ) ! ice temper (C) 166 CALL iom_put( 'ibgheatco' , zbg_ihc ) ! ice heat content (1.e20 J) 167 CALL iom_put( 'sbgheatco' , zbg_shc ) ! snw heat content (1.e20 J) 168 IF( iom_use('ibgsaltco') ) & 169 CALL iom_put( 'ibgsaltco' , zbg_sal * rhoic * r1_rau0 * 1.e-9 ) ! ice salt content (psu*km3 equivalent liquid) 170 171 CALL iom_put( 'ibgvfx' , zbg_vfx ) ! volume flux emp (m/day liquid) 172 CALL iom_put( 'ibgvfxbog' , zbg_vfx_bog ) ! volume flux bottom growth -(m/day equivalent liquid) 173 CALL iom_put( 'ibgvfxopw' , zbg_vfx_opw ) ! volume flux open water growth - 174 CALL iom_put( 'ibgvfxsni' , zbg_vfx_sni ) ! volume flux snow ice growth - 175 CALL iom_put( 'ibgvfxdyn' , zbg_vfx_dyn ) ! volume flux dynamic growth - 176 CALL iom_put( 'ibgvfxbom' , zbg_vfx_bom ) ! volume flux bottom melt - 177 CALL iom_put( 'ibgvfxsum' , zbg_vfx_sum ) ! volume flux surface melt - 178 CALL iom_put( 'ibgvfxres' , zbg_vfx_res ) ! volume flux resultant - 179 CALL iom_put( 'ibgvfxspr' , zbg_vfx_spr ) ! volume flux from snow precip - 180 CALL iom_put( 'ibgvfxsnw' , zbg_vfx_snw ) ! volume flux from snow melt - 181 CALL iom_put( 'ibgvfxsub' , zbg_vfx_sub ) ! volume flux from sublimation - 182 183 CALL iom_put( 'ibgsfx' , zbg_sfx ) ! salt flux -(psu*m/day equivalent liquid) 184 CALL iom_put( 'ibgsfxbri' , zbg_sfx_bri ) ! salt flux brines - 185 CALL iom_put( 'ibgsfxdyn' , zbg_sfx_dyn ) ! salt flux dynamic - 186 CALL iom_put( 'ibgsfxres' , zbg_sfx_res ) ! salt flux result - 187 CALL iom_put( 'ibgsfxbog' , zbg_sfx_bog ) ! salt flux bottom growth 188 CALL iom_put( 'ibgsfxopw' , zbg_sfx_opw ) ! salt flux open water growth - 189 CALL iom_put( 'ibgsfxsni' , zbg_sfx_sni ) ! salt flux snow ice growth - 190 CALL iom_put( 'ibgsfxbom' , zbg_sfx_bom ) ! salt flux bottom melt - 191 CALL iom_put( 'ibgsfxsum' , zbg_sfx_sum ) ! salt flux surface melt - 192 CALL iom_put( 'ibgsfxsub' , zbg_sfx_sub ) ! salt flux sublimation - 193 194 CALL iom_put( 'ibghfxdhc' , zbg_hfx_dhc ) ! Heat content variation in snow and ice [W] 195 CALL iom_put( 'ibghfxspr' , zbg_hfx_spr ) ! Heat content of snow precip [W] 196 197 CALL iom_put( 'ibghfxres' , zbg_hfx_res ) ! 198 CALL iom_put( 'ibghfxsub' , zbg_hfx_sub ) ! 199 CALL iom_put( 'ibghfxdyn' , zbg_hfx_dyn ) ! 200 CALL iom_put( 'ibghfxthd' , zbg_hfx_thd ) ! 201 CALL iom_put( 'ibghfxsnw' , zbg_hfx_snw ) ! 202 CALL iom_put( 'ibghfxsum' , zbg_hfx_sum ) ! 203 CALL iom_put( 'ibghfxbom' , zbg_hfx_bom ) ! 204 CALL iom_put( 'ibghfxbog' , zbg_hfx_bog ) ! 205 CALL iom_put( 'ibghfxdif' , zbg_hfx_dif ) ! 206 CALL iom_put( 'ibghfxopw' , zbg_hfx_opw ) ! 207 CALL iom_put( 'ibghfxout' , zbg_hfx_out ) ! 208 CALL iom_put( 'ibghfxin' , zbg_hfx_in ) ! 209 210 CALL iom_put( 'ibgfrcvol' , frc_vol * 1.e-9 ) ! vol - forcing (km3 equivalent liquid) 211 CALL iom_put( 'ibgfrcsfx' , frc_sal * 1.e-9 ) ! sal - forcing (psu*km3 equivalent liquid) 212 IF( iom_use('ibgvolgrm') ) & 213 CALL iom_put( 'ibgvolgrm' , bg_grme * r1_rau0 * 1.e-9 ) ! vol growth + melt (km3 equivalent liquid) 214 74 ! ---------------------------! 75 ! 2 - Trends due to forcing ! 76 ! ---------------------------! 77 z_frc_volbot = r1_rau0 * glob_sum( - ( wfx_ice(:,:) + wfx_snw(:,:) + wfx_err_sub(:,:) ) * e12t(:,:) * tmask(:,:,1) * 1.e-9 ) ! freshwater flux ice/snow-ocean 78 z_frc_voltop = r1_rau0 * glob_sum( - ( wfx_sub(:,:) + wfx_spr(:,:) ) * e12t(:,:) * tmask(:,:,1) * 1.e-9 ) ! freshwater flux ice/snow-atm 79 z_frc_sal = r1_rau0 * glob_sum( - sfx(:,:) * e12t(:,:) * tmask(:,:,1) * 1.e-9 ) ! salt fluxes ice/snow-ocean 80 z_frc_tembot = glob_sum( hfx_out(:,:) * e12t(:,:) * tmask(:,:,1) * 1.e-20 ) ! heat on top of ocean (and below ice) 81 z_frc_temtop = glob_sum( hfx_in (:,:) * e12t(:,:) * tmask(:,:,1) * 1.e-20 ) ! heat on top of ice-coean 82 ! 83 frc_voltop = frc_voltop + z_frc_voltop * rdt_ice ! km3 84 frc_volbot = frc_volbot + z_frc_volbot * rdt_ice ! km3 85 frc_sal = frc_sal + z_frc_sal * rdt_ice ! km3*pss 86 frc_temtop = frc_temtop + z_frc_temtop * rdt_ice ! 1.e20 J 87 frc_tembot = frc_tembot + z_frc_tembot * rdt_ice ! 1.e20 J 88 89 ! ----------------------- ! 90 ! 3 - Content variations ! 91 ! ----------------------- ! 92 zdiff_vol = r1_rau0 * glob_sum( ( rhoic * vt_i(:,:) + rhosn * vt_s(:,:) - vol_loc_ini(:,:) & ! freshwater trend (km3) 93 & ) * e12t(:,:) * tmask(:,:,1) * 1.e-9 ) 94 zdiff_sal = r1_rau0 * glob_sum( ( rhoic * SUM( smv_i(:,:,:), dim=3 ) - sal_loc_ini(:,:) & ! salt content trend (km3*pss) 95 & ) * e12t(:,:) * tmask(:,:,1) * 1.e-9 ) 96 zdiff_tem = glob_sum( ( et_i(:,:) + et_s(:,:) - tem_loc_ini(:,:) & ! heat content trend (1.e20 J) 97 ! & + SUM( qevap_ice * a_i_b, dim=3 ) & !! clem: I think this line should be commented (but needs a check) 98 & ) * e12t(:,:) * tmask(:,:,1) * 1.e-20 ) 99 100 ! ----------------------- ! 101 ! 4 - Drifts ! 102 ! ----------------------- ! 103 zdiff_vol = zdiff_vol - ( frc_voltop + frc_volbot ) 104 zdiff_sal = zdiff_sal - frc_sal 105 zdiff_tem = zdiff_tem - ( frc_tembot - frc_temtop ) 106 107 ! ----------------------- ! 108 ! 5 - Diagnostics writing ! 109 ! ----------------------- ! 110 ! 111 IF( iom_use('ibgvolume') ) CALL iom_put( 'ibgvolume' , zdiff_vol ) ! ice/snow volume drift (km3 equivalent ocean water) 112 IF( iom_use('ibgsaltco') ) CALL iom_put( 'ibgsaltco' , zdiff_sal ) ! ice salt content drift (psu*km3 equivalent ocean water) 113 IF( iom_use('ibgheatco') ) CALL iom_put( 'ibgheatco' , zdiff_tem ) ! ice/snow heat content drift (1.e20 J) 114 IF( iom_use('ibgheatfx') ) CALL iom_put( 'ibgheatfx' , zdiff_tem / & ! ice/snow heat flux drift (W/m2) 115 & glob_sum( e12t(:,:) * tmask(:,:,1) * 1.e-20 * kt*rdt ) ) 116 117 IF( iom_use('ibgfrcvoltop') ) CALL iom_put( 'ibgfrcvoltop' , frc_voltop ) ! vol forcing ice/snw-atm (km3 equivalent ocean water) 118 IF( iom_use('ibgfrcvolbot') ) CALL iom_put( 'ibgfrcvolbot' , frc_volbot ) ! vol forcing ice/snw-ocean (km3 equivalent ocean water) 119 IF( iom_use('ibgfrcsal') ) CALL iom_put( 'ibgfrcsal' , frc_sal ) ! sal - forcing (psu*km3 equivalent ocean water) 120 IF( iom_use('ibgfrctemtop') ) CALL iom_put( 'ibgfrctemtop' , frc_temtop ) ! heat on top of ice/snw/ocean (1.e20 J) 121 IF( iom_use('ibgfrctembot') ) CALL iom_put( 'ibgfrctembot' , frc_tembot ) ! heat on top of ocean(below ice) (1.e20 J) 122 IF( iom_use('ibgfrchfxtop') ) CALL iom_put( 'ibgfrchfxtop' , frc_temtop / & ! heat on top of ice/snw/ocean (W/m2) 123 & glob_sum( e12t(:,:) * tmask(:,:,1) * 1.e-20 * kt*rdt ) ) 124 IF( iom_use('ibgfrchfxbot') ) CALL iom_put( 'ibgfrchfxbot' , frc_tembot / & ! heat on top of ocean(below ice) (W/m2) 125 & glob_sum( e12t(:,:) * tmask(:,:,1) * 1.e-20 * kt*rdt ) ) 126 127 IF( iom_use('ibgvol_tot' ) ) CALL iom_put( 'ibgvol_tot' , zbg_ivol ) ! ice volume (km3) 128 IF( iom_use('sbgvol_tot' ) ) CALL iom_put( 'sbgvol_tot' , zbg_svol ) ! snow volume (km3) 129 IF( iom_use('ibgarea_tot') ) CALL iom_put( 'ibgarea_tot' , zbg_area ) ! ice area (km2) 130 IF( iom_use('ibgsalt_tot') ) CALL iom_put( 'ibgsalt_tot' , zbg_isal ) ! ice salinity content (pss*km3) 131 IF( iom_use('ibgheat_tot') ) CALL iom_put( 'ibgheat_tot' , zbg_item ) ! ice heat content (1.e20 J) 132 IF( iom_use('sbgheat_tot') ) CALL iom_put( 'sbgheat_tot' , zbg_stem ) ! snow heat content (1.e20 J) 215 133 ! 216 134 IF( lrst_ice ) CALL lim_diahsb_rst( numit, 'WRITE' ) 217 135 ! 218 136 IF( nn_timing == 1 ) CALL timing_stop('lim_diahsb') 219 !137 ! 220 138 END SUBROUTINE lim_diahsb 221 139 … … 233 151 !! - Compute coefficients for conversion 234 152 !!--------------------------------------------------------------------------- 235 INTEGER :: jk ! dummy loop indice236 153 INTEGER :: ierror ! local integer 237 154 !! … … 247 164 WRITE(numout,*) '~~~~~~~~~~~~' 248 165 ENDIF 249 ! 166 ! 167 ALLOCATE( vol_loc_ini(jpi,jpj), sal_loc_ini(jpi,jpj), tem_loc_ini(jpi,jpj), STAT=ierror ) 168 IF( ierror > 0 ) THEN 169 CALL ctl_stop( 'lim_diahsb: unable to allocate vol_loc_ini' ) 170 RETURN 171 ENDIF 172 250 173 CALL lim_diahsb_rst( nstart, 'READ' ) !* read or initialize all required files 251 174 ! … … 263 186 CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag 264 187 ! 265 INTEGER :: id1, id2, id3 ! local integers266 188 !!---------------------------------------------------------------------- 267 189 ! 268 190 IF( TRIM(cdrw) == 'READ' ) THEN ! Read/initialise 269 191 IF( ln_rstart ) THEN !* Read the restart file 270 !id1 = iom_varid( numrir, 'frc_vol' , ldstop = .TRUE. )271 192 ! 272 193 IF(lwp) WRITE(numout,*) '~~~~~~~' 273 IF(lwp) WRITE(numout,*) ' lim_diahsb_rst at it= ', kt,' date= ', ndastp 274 IF(lwp) WRITE(numout,*) '~~~~~~~' 275 CALL iom_get( numrir, 'frc_vol', frc_vol ) 276 CALL iom_get( numrir, 'frc_sal', frc_sal ) 277 CALL iom_get( numrir, 'bg_grme', bg_grme ) 194 IF(lwp) WRITE(numout,*) ' lim_diahsb_rst read at it= ', kt,' date= ', ndastp 195 IF(lwp) WRITE(numout,*) '~~~~~~~' 196 CALL iom_get( numrir, 'frc_voltop' , frc_voltop ) 197 CALL iom_get( numrir, 'frc_volbot' , frc_volbot ) 198 CALL iom_get( numrir, 'frc_temtop' , frc_temtop ) 199 CALL iom_get( numrir, 'frc_tembot' , frc_tembot ) 200 CALL iom_get( numrir, 'frc_sal' , frc_sal ) 201 CALL iom_get( numrir, jpdom_autoglo, 'vol_loc_ini', vol_loc_ini ) 202 CALL iom_get( numrir, jpdom_autoglo, 'tem_loc_ini', tem_loc_ini ) 203 CALL iom_get( numrir, jpdom_autoglo, 'sal_loc_ini', sal_loc_ini ) 278 204 ELSE 279 205 IF(lwp) WRITE(numout,*) '~~~~~~~' 280 206 IF(lwp) WRITE(numout,*) ' lim_diahsb at initial state ' 281 207 IF(lwp) WRITE(numout,*) '~~~~~~~' 282 frc_vol = 0._wp 283 frc_sal = 0._wp 284 bg_grme = 0._wp 208 ! set trends to 0 209 frc_voltop = 0._wp 210 frc_volbot = 0._wp 211 frc_temtop = 0._wp 212 frc_tembot = 0._wp 213 frc_sal = 0._wp 214 ! record initial ice volume, salt and temp 215 vol_loc_ini(:,:) = rhoic * vt_i(:,:) + rhosn * vt_s(:,:) ! ice/snow volume (kg/m2) 216 tem_loc_ini(:,:) = et_i(:,:) + et_s(:,:) ! ice/snow heat content (J) 217 sal_loc_ini(:,:) = rhoic * SUM( smv_i(:,:,:), dim=3 ) ! ice salt content (pss*kg/m2) 218 285 219 ENDIF 286 220 … … 288 222 ! ! ------------------- 289 223 IF(lwp) WRITE(numout,*) '~~~~~~~' 290 IF(lwp) WRITE(numout,*) ' lim_diahsb_rst at it= ', kt,' date= ', ndastp224 IF(lwp) WRITE(numout,*) ' lim_diahsb_rst write at it= ', kt,' date= ', ndastp 291 225 IF(lwp) WRITE(numout,*) '~~~~~~~' 292 CALL iom_rstput( kt, nitrst, numriw, 'frc_vol' , frc_vol ) 293 CALL iom_rstput( kt, nitrst, numriw, 'frc_sal' , frc_sal ) 294 CALL iom_rstput( kt, nitrst, numriw, 'bg_grme' , bg_grme ) 226 CALL iom_rstput( kt, nitrst, numriw, 'frc_voltop' , frc_voltop ) 227 CALL iom_rstput( kt, nitrst, numriw, 'frc_volbot' , frc_volbot ) 228 CALL iom_rstput( kt, nitrst, numriw, 'frc_temtop' , frc_temtop ) 229 CALL iom_rstput( kt, nitrst, numriw, 'frc_tembot' , frc_tembot ) 230 CALL iom_rstput( kt, nitrst, numriw, 'frc_sal' , frc_sal ) 231 CALL iom_rstput( kt, nitrst, numriw, 'vol_loc_ini', vol_loc_ini ) 232 CALL iom_rstput( kt, nitrst, numriw, 'tem_loc_ini', tem_loc_ini ) 233 CALL iom_rstput( kt, nitrst, numriw, 'sal_loc_ini', sal_loc_ini ) 295 234 ! 296 235 ENDIF -
branches/UKMO/dev_r5518_GO6_package_lim3/NEMOGCM/NEMO/LIM_SRC_3/limdyn.F90
r6486 r8540 41 41 !!---------------------------------------------------------------------- 42 42 !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) 43 !! $Id $43 !! $Id: limdyn.F90 7607 2017-01-25 15:37:31Z cetlod $ 44 44 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 45 45 !!---------------------------------------------------------------------- … … 244 244 INTEGER :: ios ! Local integer output status for namelist read 245 245 NAMELIST/namicedyn/ nn_icestr, ln_icestr_bvf, rn_pe_rdg, rn_pstar, rn_crhg, rn_cio, rn_creepl, rn_ecc, & 246 & nn_nevp, rn_relast, nn_ahi0, rn_ahi0_ref 247 INTEGER :: ji, jj 248 REAL(wp) :: za00, zd_max 246 & nn_nevp, rn_relast 249 247 !!------------------------------------------------------------------- 250 248 … … 272 270 WRITE(numout,*) ' number of iterations for subcycling nn_nevp = ', nn_nevp 273 271 WRITE(numout,*) ' ratio of elastic timescale over ice time step rn_relast = ', rn_relast 274 WRITE(numout,*) ' horizontal diffusivity calculation nn_ahi0 = ', nn_ahi0275 WRITE(numout,*) ' horizontal diffusivity coeff. (orca2 grid) rn_ahi0_ref = ', rn_ahi0_ref276 272 ENDIF 277 273 ! … … 279 275 rhoco = rau0 * rn_cio 280 276 ! 281 ! Diffusion coefficients282 SELECT CASE( nn_ahi0 )283 284 CASE( 0 )285 ahiu(:,:) = rn_ahi0_ref286 ahiv(:,:) = rn_ahi0_ref287 288 IF(lwp) WRITE(numout,*) ''289 IF(lwp) WRITE(numout,*) ' laplacian operator: ahim constant = rn_ahi0_ref'290 291 CASE( 1 )292 293 zd_max = MAX( MAXVAL( e1t(:,:) ), MAXVAL( e2t(:,:) ) )294 IF( lk_mpp ) CALL mpp_max( zd_max ) ! max over the global domain295 296 ahiu(:,:) = rn_ahi0_ref * zd_max * 1.e-05_wp ! 1.e05 = 100km = max grid space at 60° latitude in orca2297 ! (60° = min latitude for ice cover)298 ahiv(:,:) = rn_ahi0_ref * zd_max * 1.e-05_wp299 300 IF(lwp) WRITE(numout,*) ''301 IF(lwp) WRITE(numout,*) ' laplacian operator: ahim proportional to max of e1 e2 over the domain (', zd_max, ')'302 IF(lwp) WRITE(numout,*) ' value for ahim = ', rn_ahi0_ref * zd_max * 1.e-05_wp303 304 CASE( 2 )305 306 zd_max = MAX( MAXVAL( e1t(:,:) ), MAXVAL( e2t(:,:) ) )307 IF( lk_mpp ) CALL mpp_max( zd_max ) ! max over the global domain308 309 za00 = rn_ahi0_ref * 1.e-05_wp ! 1.e05 = 100km = max grid space at 60° latitude in orca2310 ! (60° = min latitude for ice cover)311 DO jj = 1, jpj312 DO ji = 1, jpi313 ahiu(ji,jj) = za00 * MAX( e1t(ji,jj), e2t(ji,jj) ) * umask(ji,jj,1)314 ahiv(ji,jj) = za00 * MAX( e1f(ji,jj), e2f(ji,jj) ) * vmask(ji,jj,1)315 END DO316 END DO317 !318 IF(lwp) WRITE(numout,*) ''319 IF(lwp) WRITE(numout,*) ' laplacian operator: ahim proportional to e1'320 IF(lwp) WRITE(numout,*) ' maximum grid-spacing = ', zd_max, ' maximum value for ahim = ', za00*zd_max321 322 END SELECT323 324 277 END SUBROUTINE lim_dyn_init 325 278 -
branches/UKMO/dev_r5518_GO6_package_lim3/NEMOGCM/NEMO/LIM_SRC_3/limhdf.F90
r7993 r8540 39 39 !!---------------------------------------------------------------------- 40 40 !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2010) 41 !! $Id $41 !! $Id: limhdf.F90 7621 2017-01-31 08:43:47Z cetlod $ 42 42 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 43 43 !!---------------------------------------------------------------------- … … 202 202 ! ----------------------- 203 203 204 IF(ln_ctl) THEN205 DO jk = 1 , isize206 zrlx(:,:,jk) = ptab(:,:,jk) - ztab0(:,:,jk)207 WRITE(charout,FMT="(' lim_hdf : zconv =',D23.16, ' iter =',I4,2X)") zconv, iter208 CALL prt_ctl( tab2d_1=zrlx(:,:,jk), clinfo1=charout )209 END DO210 ENDIF204 ! IF(ln_ctl) THEN 205 ! DO jk = 1 , isize 206 ! zrlx(:,:,jk) = ptab(:,:,jk) - ztab0(:,:,jk) 207 ! WRITE(charout,FMT="('lim_hdf : zconv =',D23.16, ' iter =',I4)") zconv, iter 208 ! CALL prt_ctl( tab2d_1=zrlx(:,:,jk), clinfo1=charout ) 209 ! END DO 210 ! ENDIF 211 211 ! 212 212 CALL wrk_dealloc( jpi, jpj, isize, zrlx, zdiv0, ztab0 ) … … 233 233 !!------------------------------------------------------------------- 234 234 INTEGER :: ios ! Local integer output status for namelist read 235 NAMELIST/namicehdf/ nn_convfrq 236 !!------------------------------------------------------------------- 237 ! 238 IF(lwp) THEN 239 WRITE(numout,*) 240 WRITE(numout,*) 'lim_hdf : Ice horizontal diffusion' 241 WRITE(numout,*) '~~~~~~~' 242 ENDIF 235 NAMELIST/namicehdf/ nn_ahi0, rn_ahi0_ref, nn_convfrq 236 INTEGER :: ji, jj 237 REAL(wp) :: za00, zd_max 238 !!------------------------------------------------------------------- 243 239 ! 244 240 REWIND( numnam_ice_ref ) ! Namelist namicehdf in reference namelist : Ice horizontal diffusion … … 253 249 IF(lwp) THEN ! control print 254 250 WRITE(numout,*) 255 WRITE(numout,*)' Namelist of ice parameters for ice horizontal diffusion computation ' 256 WRITE(numout,*)' convergence check frequency of the Crant-Nicholson scheme nn_convfrq = ', nn_convfrq 251 WRITE(numout,*) 'lim_hdf_init : Ice horizontal diffusion' 252 WRITE(numout,*) '~~~~~~~~~~~' 253 WRITE(numout,*) ' horizontal diffusivity calculation nn_ahi0 = ', nn_ahi0 254 WRITE(numout,*) ' horizontal diffusivity coeff. (orca2 grid) rn_ahi0_ref = ', rn_ahi0_ref 255 WRITE(numout,*) ' convergence check frequency of the Crant-Nicholson scheme nn_convfrq = ', nn_convfrq 257 256 ENDIF 257 ! 258 ! Diffusion coefficients 259 SELECT CASE( nn_ahi0 ) 260 261 CASE( -1 ) 262 ahiu(:,:) = 0._wp 263 ahiv(:,:) = 0._wp 264 265 IF(lwp) WRITE(numout,*) '' 266 IF(lwp) WRITE(numout,*) ' No sea-ice diffusion applied' 267 268 CASE( 0 ) 269 ahiu(:,:) = rn_ahi0_ref 270 ahiv(:,:) = rn_ahi0_ref 271 272 IF(lwp) WRITE(numout,*) '' 273 IF(lwp) WRITE(numout,*) ' laplacian operator: ahim constant = rn_ahi0_ref' 274 275 CASE( 1 ) 276 277 zd_max = MAX( MAXVAL( e1t(:,:) ), MAXVAL( e2t(:,:) ) ) 278 IF( lk_mpp ) CALL mpp_max( zd_max ) ! max over the global domain 279 280 ahiu(:,:) = rn_ahi0_ref * zd_max * 1.e-05_wp ! 1.e05 = 100km = max grid space at 60deg latitude in orca2 281 ! (60deg = min latitude for ice cover) 282 ahiv(:,:) = rn_ahi0_ref * zd_max * 1.e-05_wp 283 284 IF(lwp) WRITE(numout,*) '' 285 IF(lwp) WRITE(numout,*) ' laplacian operator: ahim proportional to max of e1 e2 over the domain (', zd_max, ')' 286 IF(lwp) WRITE(numout,*) ' value for ahim = ', rn_ahi0_ref * zd_max * 1.e-05_wp 287 288 CASE( 2 ) 289 290 zd_max = MAX( MAXVAL( e1t(:,:) ), MAXVAL( e2t(:,:) ) ) 291 IF( lk_mpp ) CALL mpp_max( zd_max ) ! max over the global domain 292 293 za00 = rn_ahi0_ref * 1.e-05_wp ! 1.e05 = 100km = max grid space at 60deg latitude in orca2 294 ! (60deg = min latitude for ice cover) 295 DO jj = 1, jpj 296 DO ji = 1, jpi 297 ahiu(ji,jj) = za00 * MAX( e1t(ji,jj), e2t(ji,jj) ) * umask(ji,jj,1) 298 ahiv(ji,jj) = za00 * MAX( e1f(ji,jj), e2f(ji,jj) ) * vmask(ji,jj,1) 299 END DO 300 END DO 301 ! 302 IF(lwp) WRITE(numout,*) '' 303 IF(lwp) WRITE(numout,*) ' laplacian operator: ahim proportional to e1' 304 IF(lwp) WRITE(numout,*) ' maximum grid-spacing = ', zd_max, ' maximum value for ahim = ', za00*zd_max 305 306 END SELECT 258 307 ! 259 308 END SUBROUTINE lim_hdf_init -
branches/UKMO/dev_r5518_GO6_package_lim3/NEMOGCM/NEMO/LIM_SRC_3/limistate.F90
r7993 r8540 51 51 !!---------------------------------------------------------------------- 52 52 !! LIM 3.0, UCL-LOCEAN-IPSL (2008) 53 !! $Id $53 !! $Id: limistate.F90 6696 2016-06-13 14:55:42Z clem $ 54 54 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 55 55 !!---------------------------------------------------------------------- -
branches/UKMO/dev_r5518_GO6_package_lim3/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90
r7993 r8540 61 61 !!---------------------------------------------------------------------- 62 62 !! NEMO/LIM3 3.3 , UCL - NEMO Consortium (2010) 63 !! $Id $63 !! $Id: limitd_me.F90 7814 2017-03-20 16:21:42Z clem $ 64 64 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 65 65 !!---------------------------------------------------------------------- … … 259 259 closing_net(ji,jj) = 0._wp 260 260 opning (ji,jj) = 0._wp 261 ato_i (ji,jj) = MAX( 0._wp, 1._wp - SUM( a_i(ji,jj,:) ) ) 261 262 ELSE 262 263 iterate_ridging = 1 … … 844 845 END DO 845 846 846 strength(:,:) = rn_pe_rdg * Cp * strength(:,:) / aksum(:,:) 847 strength(:,:) = rn_pe_rdg * Cp * strength(:,:) / aksum(:,:) * tmask(:,:,1) 847 848 ! where Cp = (g/2)*(rhow-rhoi)*(rhoi/rhow) and rn_pe_rdg accounts for frictional dissipation 848 849 ksmooth = 1 … … 853 854 ELSE ! kstrngth ne 1: Hibler (1979) form 854 855 ! 855 strength(:,:) = rn_pstar * vt_i(:,:) * EXP( - rn_crhg * ( 1._wp - at_i(:,:) ) ) 856 strength(:,:) = rn_pstar * vt_i(:,:) * EXP( - rn_crhg * ( 1._wp - at_i(:,:) ) ) * tmask(:,:,1) 856 857 ! 857 858 ksmooth = 1 … … 866 867 DO jj = 1, jpj 867 868 DO ji = 1, jpi 868 strength(ji,jj) = strength(ji,jj) * exp(-5.88*SQRT(MAX(bv _i(ji,jj),0.0)))869 strength(ji,jj) = strength(ji,jj) * exp(-5.88*SQRT(MAX(bvm_i(ji,jj),0.0))) 869 870 END DO 870 871 END DO -
branches/UKMO/dev_r5518_GO6_package_lim3/NEMOGCM/NEMO/LIM_SRC_3/limitd_th.F90
r6486 r8540 40 40 !!---------------------------------------------------------------------- 41 41 !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2010) 42 !! $Id $42 !! $Id: limitd_th.F90 5407 2015-06-11 19:13:22Z smasson $ 43 43 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 44 44 !!---------------------------------------------------------------------- -
branches/UKMO/dev_r5518_GO6_package_lim3/NEMOGCM/NEMO/LIM_SRC_3/limmsh.F90
r6486 r8540 27 27 !!---------------------------------------------------------------------- 28 28 !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) 29 !! $Id $29 !! $Id: limmsh.F90 5123 2015-03-04 16:06:03Z clem $ 30 30 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 31 31 !!---------------------------------------------------------------------- -
branches/UKMO/dev_r5518_GO6_package_lim3/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90
r6487 r8540 10 10 !! 3.4 ! 2011-01 (A. Porter) dynamical allocation 11 11 !! 3.5 ! 2012-08 (R. Benshila) AGRIF 12 !! 3.6 ! 2016-06 (C. Rousset) Rewriting (conserves energy) 12 13 !!---------------------------------------------------------------------- 13 14 #if defined key_lim3 || ( defined key_lim2 && ! defined key_lim2_vp ) … … 54 55 !!---------------------------------------------------------------------- 55 56 !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) 56 !! $Id $57 !! $Id: limrhg.F90 6964 2016-09-30 12:41:39Z clem $ 57 58 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 58 59 !!---------------------------------------------------------------------- … … 95 96 !! coriolis terms of the momentum equation 96 97 !! 3) Solve the momentum equation (iterative procedure) 97 !! 4) Prevent high velocities if the ice is thin 98 !! 5) Recompute invariants of the strain rate tensor 98 !! 4) Recompute invariants of the strain rate tensor 99 99 !! which are inputs of the ITD, store stress 100 100 !! for the next time step 101 !! 6) Control prints of residual (convergence)101 !! 5) Control prints of residual (convergence) 102 102 !! and charge ellipse. 103 103 !! The user should make sure that the parameters … … 106 106 !! e.g. in the Canadian Archipelago 107 107 !! 108 !! ** Notes : Boundary condition for ice is chosen no-slip 109 !! but can be adjusted with param rn_shlat 110 !! 108 111 !! References : Hunke and Dukowicz, JPO97 109 112 !! Bouillon et al., Ocean Modelling 2009 … … 115 118 INTEGER :: jter ! local integers 116 119 CHARACTER (len=50) :: charout 117 REAL(wp) :: zt11, zt12, zt21, zt22, ztagnx, ztagny, delta ! 118 REAL(wp) :: za, zstms ! local scalars 119 REAL(wp) :: zc1, zc2, zc3 ! ice mass 120 121 REAL(wp) :: dtevp , z1_dtevp ! time step for subcycling 122 REAL(wp) :: dtotel, z1_dtotel, ecc2, ecci ! square of yield ellipse eccenticity 123 REAL(wp) :: z0, zr, zcca, zccb ! temporary scalars 124 REAL(wp) :: zu_ice2, zv_ice1 ! 125 REAL(wp) :: zddc, zdtc ! delta on corners and on centre 126 REAL(wp) :: zdst ! shear at the center of the grid point 127 REAL(wp) :: zdsshx, zdsshy ! term for the gradient of ocean surface 128 REAL(wp) :: sigma1, sigma2 ! internal ice stress 129 130 REAL(wp) :: zresm ! Maximal error on ice velocity 131 REAL(wp) :: zintb, zintn ! dummy argument 132 133 REAL(wp), POINTER, DIMENSION(:,:) :: zpresh ! temporary array for ice strength 134 REAL(wp), POINTER, DIMENSION(:,:) :: zpreshc ! Ice strength on grid cell corners (zpreshc) 135 REAL(wp), POINTER, DIMENSION(:,:) :: zfrld1, zfrld2 ! lead fraction on U/V points 136 REAL(wp), POINTER, DIMENSION(:,:) :: zmass1, zmass2 ! ice/snow mass on U/V points 137 REAL(wp), POINTER, DIMENSION(:,:) :: zcorl1, zcorl2 ! coriolis parameter on U/V points 138 REAL(wp), POINTER, DIMENSION(:,:) :: za1ct , za2ct ! temporary arrays 139 REAL(wp), POINTER, DIMENSION(:,:) :: v_oce1 ! ocean u/v component on U points 140 REAL(wp), POINTER, DIMENSION(:,:) :: u_oce2 ! ocean u/v component on V points 141 REAL(wp), POINTER, DIMENSION(:,:) :: u_ice2, v_ice1 ! ice u/v component on V/U point 142 REAL(wp), POINTER, DIMENSION(:,:) :: zf1 , zf2 ! arrays for internal stresses 143 REAL(wp), POINTER, DIMENSION(:,:) :: zmask ! mask ocean grid points 120 121 REAL(wp) :: zdtevp, z1_dtevp ! time step for subcycling 122 REAL(wp) :: ecc2, z1_ecc2 ! square of yield ellipse eccenticity 123 REAL(wp) :: zbeta, zalph1, z1_alph1, zalph2, z1_alph2 ! alpha and beta from Bouillon 2009 and 2013 124 REAL(wp) :: zm1, zm2, zm3, zmassU, zmassV ! ice/snow mass 125 REAL(wp) :: zdelta, zp_delf, zds2, zdt, zdt2, zdiv, zdiv2 ! temporary scalars 126 REAL(wp) :: zTauO, zTauE, zCor ! temporary scalars 127 128 REAL(wp) :: zsig1, zsig2 ! internal ice stress 129 REAL(wp) :: zresm ! Maximal error on ice velocity 130 REAL(wp) :: zintb, zintn ! dummy argument 144 131 145 REAL(wp), POINTER, DIMENSION(:,:) :: zdt ! tension at centre of grid cells 146 REAL(wp), POINTER, DIMENSION(:,:) :: zds ! Shear on northeast corner of grid cells 147 REAL(wp), POINTER, DIMENSION(:,:) :: zs1 , zs2 ! Diagonal stress tensor components zs1 and zs2 148 REAL(wp), POINTER, DIMENSION(:,:) :: zs12 ! Non-diagonal stress tensor component zs12 149 REAL(wp), POINTER, DIMENSION(:,:) :: zu_ice, zv_ice, zresr ! Local error on velocity 150 REAL(wp), POINTER, DIMENSION(:,:) :: zpice ! array used for the calculation of ice surface slope: 151 ! ocean surface (ssh_m) if ice is not embedded 152 ! ice top surface if ice is embedded 153 154 REAL(wp), PARAMETER :: zepsi = 1.0e-20_wp ! tolerance parameter 155 REAL(wp), PARAMETER :: zvmin = 1.0e-03_wp ! ice volume below which ice velocity equals ocean velocity 132 REAL(wp), POINTER, DIMENSION(:,:) :: zpresh ! temporary array for ice strength 133 REAL(wp), POINTER, DIMENSION(:,:) :: z1_e1t0, z1_e2t0 ! scale factors 134 REAL(wp), POINTER, DIMENSION(:,:) :: zp_delt ! P/delta at T points 135 ! 136 REAL(wp), POINTER, DIMENSION(:,:) :: zaU , zaV ! ice fraction on U/V points 137 REAL(wp), POINTER, DIMENSION(:,:) :: zmU_t, zmV_t ! ice/snow mass/dt on U/V points 138 REAL(wp), POINTER, DIMENSION(:,:) :: zmf ! coriolis parameter at T points 139 REAL(wp), POINTER, DIMENSION(:,:) :: zTauU_ia , ztauV_ia ! ice-atm. stress at U-V points 140 REAL(wp), POINTER, DIMENSION(:,:) :: zspgU , zspgV ! surface pressure gradient at U/V points 141 REAL(wp), POINTER, DIMENSION(:,:) :: v_oceU, u_oceV, v_iceU, u_iceV ! ocean/ice u/v component on V/U points 142 REAL(wp), POINTER, DIMENSION(:,:) :: zfU , zfV ! internal stresses 143 144 REAL(wp), POINTER, DIMENSION(:,:) :: zds ! shear 145 REAL(wp), POINTER, DIMENSION(:,:) :: zs1, zs2, zs12 ! stress tensor components 146 REAL(wp), POINTER, DIMENSION(:,:) :: zu_ice, zv_ice, zresr ! check convergence 147 REAL(wp), POINTER, DIMENSION(:,:) :: zpice ! array used for the calculation of ice surface slope: 148 ! ocean surface (ssh_m) if ice is not embedded 149 ! ice top surface if ice is embedded 150 REAL(wp), POINTER, DIMENSION(:,:) :: zswitchU, zswitchV ! dummy arrays 151 REAL(wp), POINTER, DIMENSION(:,:) :: zmaskU, zmaskV ! mask for ice presence 152 REAL(wp), POINTER, DIMENSION(:,:) :: zfmask, zwf ! mask at F points for the ice 153 154 REAL(wp), PARAMETER :: zepsi = 1.0e-20_wp ! tolerance parameter 155 REAL(wp), PARAMETER :: zmmin = 1._wp ! ice mass (kg/m2) below which ice velocity equals ocean velocity 156 REAL(wp), PARAMETER :: zshlat = 2._wp ! boundary condition for sea-ice velocity (2=no slip ; 0=free slip) 156 157 !!------------------------------------------------------------------- 157 158 158 CALL wrk_alloc( jpi,jpj, zpresh, zfrld1, zmass1, zcorl1, za1ct , zpreshc, zfrld2, zmass2, zcorl2, za2ct ) 159 CALL wrk_alloc( jpi,jpj, u_oce2, u_ice2, v_oce1 , v_ice1 , zmask ) 160 CALL wrk_alloc( jpi,jpj, zf1 , zu_ice, zf2 , zv_ice , zdt , zds ) 161 CALL wrk_alloc( jpi,jpj, zs1 , zs2 , zs12 , zresr , zpice ) 159 CALL wrk_alloc( jpi,jpj, zpresh, z1_e1t0, z1_e2t0, zp_delt ) 160 CALL wrk_alloc( jpi,jpj, zaU, zaV, zmU_t, zmV_t, zmf, zTauU_ia, ztauV_ia ) 161 CALL wrk_alloc( jpi,jpj, zspgU, zspgV, v_oceU, u_oceV, v_iceU, u_iceV, zfU, zfV ) 162 CALL wrk_alloc( jpi,jpj, zds, zs1, zs2, zs12, zu_ice, zv_ice, zresr, zpice ) 163 CALL wrk_alloc( jpi,jpj, zswitchU, zswitchV, zmaskU, zmaskV, zfmask, zwf ) 162 164 163 165 #if defined key_lim2 && ! defined key_lim2_vp … … 176 178 ! 177 179 !------------------------------------------------------------------------------! 178 ! 1) Ice strength (zpresh) ! 179 !------------------------------------------------------------------------------! 180 ! 181 ! Put every vector to 0 182 delta_i(:,:) = 0._wp ; 183 zpresh (:,:) = 0._wp ; 184 zpreshc(:,:) = 0._wp 185 u_ice2 (:,:) = 0._wp ; v_ice1(:,:) = 0._wp 186 divu_i (:,:) = 0._wp ; zdt (:,:) = 0._wp ; zds(:,:) = 0._wp 187 shear_i(:,:) = 0._wp 188 180 ! 0) mask at F points for the ice (on the whole domain, not only k_j1,k_jpj) 181 !------------------------------------------------------------------------------! 182 ! ocean/land mask 183 DO jj = 1, jpjm1 184 DO ji = 1, jpim1 ! NO vector opt. 185 zfmask(ji,jj) = tmask(ji,jj,1) * tmask(ji+1,jj,1) * tmask(ji,jj+1,1) * tmask(ji+1,jj+1,1) 186 END DO 187 END DO 188 CALL lbc_lnk( zfmask, 'F', 1._wp ) 189 190 ! Lateral boundary conditions on velocity (modify zfmask) 191 zwf(:,:) = zfmask(:,:) 192 DO jj = 2, jpjm1 193 DO ji = fs_2, fs_jpim1 ! vector opt. 194 IF( zfmask(ji,jj) == 0._wp ) THEN 195 zfmask(ji,jj) = zshlat * MIN( 1._wp , MAX( zwf(ji+1,jj), zwf(ji,jj+1), zwf(ji-1,jj), zwf(ji,jj-1) ) ) 196 ENDIF 197 END DO 198 END DO 199 DO jj = 2, jpjm1 200 IF( zfmask(1,jj) == 0._wp ) THEN 201 zfmask(1 ,jj) = zshlat * MIN( 1._wp , MAX( zwf(2,jj), zwf(1,jj+1), zwf(1,jj-1) ) ) 202 ENDIF 203 IF( zfmask(jpi,jj) == 0._wp ) THEN 204 zfmask(jpi,jj) = zshlat * MIN( 1._wp , MAX( zwf(jpi,jj+1), zwf(jpim1,jj), zwf(jpi,jj-1) ) ) 205 ENDIF 206 END DO 207 DO ji = 2, jpim1 208 IF( zfmask(ji,1) == 0._wp ) THEN 209 zfmask(ji,1 ) = zshlat * MIN( 1._wp , MAX( zwf(ji+1,1), zwf(ji,2), zwf(ji-1,1) ) ) 210 ENDIF 211 IF( zfmask(ji,jpj) == 0._wp ) THEN 212 zfmask(ji,jpj) = zshlat * MIN( 1._wp , MAX( zwf(ji+1,jpj), zwf(ji-1,jpj), zwf(ji,jpjm1) ) ) 213 ENDIF 214 END DO 215 CALL lbc_lnk( zfmask, 'F', 1._wp ) 216 217 !------------------------------------------------------------------------------! 218 ! 1) define some variables and initialize arrays 219 !------------------------------------------------------------------------------! 220 ! ecc2: square of yield ellipse eccenticrity 221 ecc2 = rn_ecc * rn_ecc 222 z1_ecc2 = 1._wp / ecc2 223 224 ! Time step for subcycling 225 zdtevp = rdt_ice / REAL( nn_nevp ) 226 z1_dtevp = 1._wp / zdtevp 227 228 ! alpha parameters (Bouillon 2009) 189 229 #if defined key_lim3 190 CALL lim_itd_me_icestrength( nn_icestr ) ! LIM-3: Ice strength on T-points 191 #endif 192 193 DO jj = k_j1 , k_jpj ! Ice mass and temp variables 194 DO ji = 1 , jpi 230 zalph1 = ( 2._wp * rn_relast * rdt_ice ) * z1_dtevp 231 #else 232 zalph1 = ( 2._wp * telast ) * z1_dtevp 233 #endif 234 zalph2 = zalph1 * z1_ecc2 235 236 z1_alph1 = 1._wp / ( zalph1 + 1._wp ) 237 z1_alph2 = 1._wp / ( zalph2 + 1._wp ) 238 239 ! Initialise stress tensor 240 zs1 (:,:) = stress1_i (:,:) 241 zs2 (:,:) = stress2_i (:,:) 242 zs12(:,:) = stress12_i(:,:) 243 244 ! Ice strength 195 245 #if defined key_lim3 196 zpresh(ji,jj) = tmask(ji,jj,1) * strength(ji,jj) 197 #endif 198 #if defined key_lim2 199 zpresh(ji,jj) = tmask(ji,jj,1) * pstar * vt_i(ji,jj) * EXP( -c_rhg * (1. - at_i(ji,jj) ) ) 200 #endif 201 ! zmask = 1 where there is ice or on land 202 zmask(ji,jj) = 1._wp - ( 1._wp - MAX( 0._wp , SIGN ( 1._wp , vt_i(ji,jj) - zepsi ) ) ) * tmask(ji,jj,1) 246 CALL lim_itd_me_icestrength( nn_icestr ) 247 zpresh(:,:) = tmask(:,:,1) * strength(:,:) 248 #else 249 zpresh(:,:) = tmask(:,:,1) * pstar * vt_i(:,:) * EXP( -c_rhg * (1. - at_i(:,:) ) ) 250 #endif 251 252 ! scale factors 253 DO jj = k_j1+1, k_jpj-1 254 DO ji = fs_2, fs_jpim1 255 z1_e1t0(ji,jj) = 1._wp / ( e1t(ji+1,jj ) + e1t(ji,jj ) ) 256 z1_e2t0(ji,jj) = 1._wp / ( e2t(ji ,jj+1) + e2t(ji,jj ) ) 203 257 END DO 204 258 END DO 205 206 ! Ice strength on grid cell corners (zpreshc) 207 ! needed for calculation of shear stress 208 DO jj = k_j1+1, k_jpj-1 209 DO ji = 2, jpim1 !RB caution no fs_ (ji+1,jj+1) 210 zstms = tmask(ji+1,jj+1,1) * wght(ji+1,jj+1,2,2) + tmask(ji,jj+1,1) * wght(ji+1,jj+1,1,2) + & 211 & tmask(ji+1,jj,1) * wght(ji+1,jj+1,2,1) + tmask(ji,jj,1) * wght(ji+1,jj+1,1,1) 212 zpreshc(ji,jj) = ( zpresh(ji+1,jj+1) * wght(ji+1,jj+1,2,2) + zpresh(ji,jj+1) * wght(ji+1,jj+1,1,2) + & 213 & zpresh(ji+1,jj) * wght(ji+1,jj+1,2,1) + zpresh(ji,jj) * wght(ji+1,jj+1,1,1) & 214 & ) / MAX( zstms, zepsi ) 215 END DO 216 END DO 217 CALL lbc_lnk( zpreshc(:,:), 'F', 1. ) 259 218 260 ! 219 261 !------------------------------------------------------------------------------! 220 262 ! 2) Wind / ocean stress, mass terms, coriolis terms 221 263 !------------------------------------------------------------------------------! 222 !223 ! Wind stress, coriolis and mass terms on the sides of the squares224 ! zfrld1: lead fraction on U-points225 ! zfrld2: lead fraction on V-points226 ! zmass1: ice/snow mass on U-points227 ! zmass2: ice/snow mass on V-points228 ! zcorl1: Coriolis parameter on U-points229 ! zcorl2: Coriolis parameter on V-points230 ! (ztagnx,ztagny): wind stress on U/V points231 ! v_oce1: ocean v component on u points232 ! u_oce2: ocean u component on v points233 264 234 265 IF( nn_ice_embd == 2 ) THEN !== embedded sea ice: compute representative ice top surface ==! … … 242 273 zintb = REAL( nn_fsbc + 1 ) / REAL( nn_fsbc ) * 0.5_wp 243 274 ! 244 zpice(:,:) = ssh_m(:,:) + ( zintn * snwice_mass(:,:) + zintb * snwice_mass_b(:,:)) * r1_rau0275 zpice(:,:) = ssh_m(:,:) + ( zintn * snwice_mass(:,:) + zintb * snwice_mass_b(:,:) ) * r1_rau0 245 276 ! 246 277 ELSE !== non-embedded sea ice: use ocean surface for slope calculation ==! … … 251 282 DO ji = fs_2, fs_jpim1 252 283 253 zc1 = tmask(ji ,jj ,1) * ( rhosn * vt_s(ji ,jj ) + rhoic * vt_i(ji ,jj ) ) 254 zc2 = tmask(ji+1,jj ,1) * ( rhosn * vt_s(ji+1,jj ) + rhoic * vt_i(ji+1,jj ) ) 255 zc3 = tmask(ji ,jj+1,1) * ( rhosn * vt_s(ji ,jj+1) + rhoic * vt_i(ji ,jj+1) ) 256 257 zt11 = tmask(ji ,jj,1) * e1t(ji ,jj) 258 zt12 = tmask(ji+1,jj,1) * e1t(ji+1,jj) 259 zt21 = tmask(ji,jj ,1) * e2t(ji,jj ) 260 zt22 = tmask(ji,jj+1,1) * e2t(ji,jj+1) 261 262 ! Leads area. 263 zfrld1(ji,jj) = ( zt12 * ( 1.0 - at_i(ji,jj) ) + zt11 * ( 1.0 - at_i(ji+1,jj) ) ) / ( zt11 + zt12 + zepsi ) 264 zfrld2(ji,jj) = ( zt22 * ( 1.0 - at_i(ji,jj) ) + zt21 * ( 1.0 - at_i(ji,jj+1) ) ) / ( zt21 + zt22 + zepsi ) 265 266 ! Mass, coriolis coeff. and currents 267 zmass1(ji,jj) = ( zt12 * zc1 + zt11 * zc2 ) / ( zt11 + zt12 + zepsi ) 268 zmass2(ji,jj) = ( zt22 * zc1 + zt21 * zc3 ) / ( zt21 + zt22 + zepsi ) 269 zcorl1(ji,jj) = zmass1(ji,jj) * ( e1t(ji+1,jj) * fcor(ji,jj) + e1t(ji,jj) * fcor(ji+1,jj) ) & 270 & / ( e1t(ji,jj) + e1t(ji+1,jj) + zepsi ) 271 zcorl2(ji,jj) = zmass2(ji,jj) * ( e2t(ji,jj+1) * fcor(ji,jj) + e2t(ji,jj) * fcor(ji,jj+1) ) & 272 & / ( e2t(ji,jj+1) + e2t(ji,jj) + zepsi ) 273 ! 274 ! Ocean has no slip boundary condition 275 v_oce1(ji,jj) = 0.5 * ( ( v_oce(ji ,jj) + v_oce(ji ,jj-1) ) * e1t(ji,jj) & 276 & + ( v_oce(ji+1,jj) + v_oce(ji+1,jj-1) ) * e1t(ji+1,jj) ) & 277 & / ( e1t(ji+1,jj) + e1t(ji,jj) ) * umask(ji,jj,1) 278 279 u_oce2(ji,jj) = 0.5 * ( ( u_oce(ji,jj ) + u_oce(ji-1,jj ) ) * e2t(ji,jj) & 280 & + ( u_oce(ji,jj+1) + u_oce(ji-1,jj+1) ) * e2t(ji,jj+1) ) & 281 & / ( e2t(ji,jj+1) + e2t(ji,jj) ) * vmask(ji,jj,1) 282 283 ! Wind stress at U,V-point 284 ztagnx = ( 1. - zfrld1(ji,jj) ) * utau_ice(ji,jj) 285 ztagny = ( 1. - zfrld2(ji,jj) ) * vtau_ice(ji,jj) 286 287 ! Computation of the velocity field taking into account the ice internal interaction. 288 ! Terms that are independent of the velocity field. 289 290 ! SB On utilise maintenant le gradient de la pente de l'ocean 291 ! include it later 292 293 zdsshx = ( zpice(ji+1,jj) - zpice(ji,jj) ) * r1_e1u(ji,jj) 294 zdsshy = ( zpice(ji,jj+1) - zpice(ji,jj) ) * r1_e2v(ji,jj) 295 296 za1ct(ji,jj) = ztagnx - zmass1(ji,jj) * grav * zdsshx 297 za2ct(ji,jj) = ztagny - zmass2(ji,jj) * grav * zdsshy 284 ! ice fraction at U-V points 285 zaU(ji,jj) = 0.5_wp * ( at_i(ji,jj) * e12t(ji,jj) + at_i(ji+1,jj) * e12t(ji+1,jj) ) * r1_e12u(ji,jj) * umask(ji,jj,1) 286 zaV(ji,jj) = 0.5_wp * ( at_i(ji,jj) * e12t(ji,jj) + at_i(ji,jj+1) * e12t(ji,jj+1) ) * r1_e12v(ji,jj) * vmask(ji,jj,1) 287 288 ! Ice/snow mass at U-V points 289 zm1 = ( rhosn * vt_s(ji ,jj ) + rhoic * vt_i(ji ,jj ) ) 290 zm2 = ( rhosn * vt_s(ji+1,jj ) + rhoic * vt_i(ji+1,jj ) ) 291 zm3 = ( rhosn * vt_s(ji ,jj+1) + rhoic * vt_i(ji ,jj+1) ) 292 zmassU = 0.5_wp * ( zm1 * e12t(ji,jj) + zm2 * e12t(ji+1,jj) ) * r1_e12u(ji,jj) * umask(ji,jj,1) 293 zmassV = 0.5_wp * ( zm1 * e12t(ji,jj) + zm3 * e12t(ji,jj+1) ) * r1_e12v(ji,jj) * vmask(ji,jj,1) 294 295 ! Ocean currents at U-V points 296 v_oceU(ji,jj) = 0.5_wp * ( ( v_oce(ji ,jj) + v_oce(ji ,jj-1) ) * e1t(ji+1,jj) & 297 & + ( v_oce(ji+1,jj) + v_oce(ji+1,jj-1) ) * e1t(ji ,jj) ) * z1_e1t0(ji,jj) * umask(ji,jj,1) 298 299 u_oceV(ji,jj) = 0.5_wp * ( ( u_oce(ji,jj ) + u_oce(ji-1,jj ) ) * e2t(ji,jj+1) & 300 & + ( u_oce(ji,jj+1) + u_oce(ji-1,jj+1) ) * e2t(ji,jj ) ) * z1_e2t0(ji,jj) * vmask(ji,jj,1) 301 302 ! Coriolis at T points (m*f) 303 zmf(ji,jj) = zm1 * fcor(ji,jj) 304 305 ! m/dt 306 zmU_t(ji,jj) = zmassU * z1_dtevp 307 zmV_t(ji,jj) = zmassV * z1_dtevp 308 309 ! Drag ice-atm. 310 zTauU_ia(ji,jj) = zaU(ji,jj) * utau_ice(ji,jj) 311 zTauV_ia(ji,jj) = zaV(ji,jj) * vtau_ice(ji,jj) 312 313 ! Surface pressure gradient (- m*g*GRAD(ssh)) at U-V points 314 zspgU(ji,jj) = - zmassU * grav * ( zpice(ji+1,jj) - zpice(ji,jj) ) * r1_e1u(ji,jj) 315 zspgV(ji,jj) = - zmassV * grav * ( zpice(ji,jj+1) - zpice(ji,jj) ) * r1_e2v(ji,jj) 316 317 ! masks 318 zmaskU(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassU ) ) ! 0 if no ice 319 zmaskV(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassV ) ) ! 0 if no ice 320 321 ! switches 322 zswitchU(ji,jj) = MAX( 0._wp, SIGN( 1._wp, zmassU - zmmin ) ) ! 0 if ice mass < zmmin 323 zswitchV(ji,jj) = MAX( 0._wp, SIGN( 1._wp, zmassV - zmmin ) ) ! 0 if ice mass < zmmin 298 324 299 325 END DO 300 326 END DO 301 327 CALL lbc_lnk( zmf, 'T', 1. ) 302 328 ! 303 329 !------------------------------------------------------------------------------! … … 305 331 !------------------------------------------------------------------------------! 306 332 ! 307 ! Time step for subcycling308 dtevp = rdt_ice / nn_nevp309 #if defined key_lim3310 dtotel = dtevp / ( 2._wp * rn_relast * rdt_ice )311 #else312 dtotel = dtevp / ( 2._wp * telast )313 #endif314 z1_dtotel = 1._wp / ( 1._wp + dtotel )315 z1_dtevp = 1._wp / dtevp316 !-ecc2: square of yield ellipse eccenticrity (reminder: must become a namelist parameter)317 ecc2 = rn_ecc * rn_ecc318 ecci = 1. / ecc2319 320 !-Initialise stress tensor321 zs1 (:,:) = stress1_i (:,:)322 zs2 (:,:) = stress2_i (:,:)323 zs12(:,:) = stress12_i(:,:)324 325 333 ! !----------------------! 326 334 DO jter = 1 , nn_nevp ! loop over jter ! 327 335 ! !----------------------! 328 DO jj = k_j1, k_jpj-1 329 zu_ice(:,jj) = u_ice(:,jj) ! velocity at previous time step 330 zv_ice(:,jj) = v_ice(:,jj) 331 END DO 332 333 DO jj = k_j1+1, k_jpj-1 334 DO ji = fs_2, fs_jpim1 !RB bug no vect opt due to zmask 335 336 ! 337 !- Divergence, tension and shear (Section a. Appendix B of Hunke & Dukowicz, 2002) 338 !- divu_i(:,:), zdt(:,:): divergence and tension at centre of grid cells 339 !- zds(:,:): shear on northeast corner of grid cells 340 ! 341 !- IMPORTANT REMINDER: Dear Gurvan, note that, the way these terms are coded, 342 ! there are many repeated calculations. 343 ! Speed could be improved by regrouping terms. For 344 ! the moment, however, the stress is on clarity of coding to avoid 345 ! bugs (Martin, for Miguel). 346 ! 347 !- ALSO: arrays zdt, zds and delta could 348 ! be removed in the future to minimise memory demand. 349 ! 350 !- MORE NOTES: Note that we are calculating deformation rates and stresses on the corners of 351 ! grid cells, exactly as in the B grid case. For simplicity, the indexation on 352 ! the corners is the same as in the B grid. 353 ! 354 ! 355 divu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj) & 356 & + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1) & 357 & ) * r1_e12t(ji,jj) 358 359 zdt(ji,jj) = ( ( u_ice(ji,jj) * r1_e2u(ji,jj) - u_ice(ji-1,jj) * r1_e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj) & 360 & - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj) & 361 & ) * r1_e12t(ji,jj) 362 363 ! 336 IF(ln_ctl) THEN ! Convergence test 337 DO jj = k_j1, k_jpj-1 338 zu_ice(:,jj) = u_ice(:,jj) ! velocity at previous time step 339 zv_ice(:,jj) = v_ice(:,jj) 340 END DO 341 ENDIF 342 343 ! --- divergence, tension & shear (Appendix B of Hunke & Dukowicz, 2002) --- ! 344 DO jj = k_j1, k_jpj-1 ! loops start at 1 since there is no boundary condition (lbc_lnk) at i=1 and j=1 for F points 345 DO ji = 1, jpim1 346 347 ! shear at F points 364 348 zds(ji,jj) = ( ( u_ice(ji,jj+1) * r1_e1u(ji,jj+1) - u_ice(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj) & 365 349 & + ( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj) & 366 & ) * r1_e12f(ji,jj) * ( 2._wp - fmask(ji,jj,1) ) & 367 & * zmask(ji,jj) * zmask(ji,jj+1) * zmask(ji+1,jj) * zmask(ji+1,jj+1) 368 369 370 v_ice1(ji,jj) = 0.5_wp * ( ( v_ice(ji ,jj) + v_ice(ji ,jj-1) ) * e1t(ji+1,jj) & 371 & + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji ,jj) ) & 372 & / ( e1t(ji+1,jj) + e1t(ji,jj) ) * umask(ji,jj,1) 373 374 u_ice2(ji,jj) = 0.5_wp * ( ( u_ice(ji,jj ) + u_ice(ji-1,jj ) ) * e2t(ji,jj+1) & 375 & + ( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj ) ) & 376 & / ( e2t(ji,jj+1) + e2t(ji,jj) ) * vmask(ji,jj,1) 377 END DO 378 END DO 379 380 CALL lbc_lnk_multi( v_ice1, 'U', -1., u_ice2, 'V', -1. ) ! lateral boundary cond. 381 350 & ) * r1_e12f(ji,jj) * zfmask(ji,jj) 351 352 END DO 353 END DO 354 CALL lbc_lnk( zds, 'F', 1. ) 355 382 356 DO jj = k_j1+1, k_jpj-1 383 DO ji = fs_2, fs_jpim1 384 385 !- Calculate Delta at centre of grid cells 386 zdst = ( e2u(ji,jj) * v_ice1(ji,jj) - e2u(ji-1,jj ) * v_ice1(ji-1,jj ) & 387 & + e1v(ji,jj) * u_ice2(ji,jj) - e1v(ji ,jj-1) * u_ice2(ji ,jj-1) & 388 & ) * r1_e12t(ji,jj) 389 390 delta = SQRT( divu_i(ji,jj)**2 + ( zdt(ji,jj)**2 + zdst**2 ) * usecc2 ) 391 delta_i(ji,jj) = delta + rn_creepl 392 393 !- Calculate Delta on corners 394 zddc = ( ( v_ice1(ji,jj+1) * r1_e1u(ji,jj+1) - v_ice1(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj) & 395 & + ( u_ice2(ji+1,jj) * r1_e2v(ji+1,jj) - u_ice2(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj) & 396 & ) * r1_e12f(ji,jj) 397 398 zdtc = (- ( v_ice1(ji,jj+1) * r1_e1u(ji,jj+1) - v_ice1(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj) & 399 & + ( u_ice2(ji+1,jj) * r1_e2v(ji+1,jj) - u_ice2(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj) & 400 & ) * r1_e12f(ji,jj) 401 402 zddc = SQRT( zddc**2 + ( zdtc**2 + zds(ji,jj)**2 ) * usecc2 ) + rn_creepl 403 404 !-Calculate stress tensor components zs1 and zs2 at centre of grid cells (see section 3.5 of CICE user's guide). 405 zs1(ji,jj) = ( zs1 (ji,jj) + dtotel * ( divu_i(ji,jj) - delta ) / delta_i(ji,jj) * zpresh(ji,jj) & 406 & ) * z1_dtotel 407 zs2(ji,jj) = ( zs2 (ji,jj) + dtotel * ecci * zdt(ji,jj) / delta_i(ji,jj) * zpresh(ji,jj) & 408 & ) * z1_dtotel 409 !-Calculate stress tensor component zs12 at corners 410 zs12(ji,jj) = ( zs12(ji,jj) + dtotel * ecci * zds(ji,jj) / ( 2._wp * zddc ) * zpreshc(ji,jj) & 411 & ) * z1_dtotel 412 413 END DO 414 END DO 415 416 CALL lbc_lnk_multi( zs1 , 'T', 1., zs2, 'T', 1., zs12, 'F', 1. ) 357 DO ji = 2, jpim1 ! no vector loop 358 359 ! shear**2 at T points (doc eq. A16) 360 zds2 = ( zds(ji,jj ) * zds(ji,jj ) * e12f(ji,jj ) + zds(ji-1,jj ) * zds(ji-1,jj ) * e12f(ji-1,jj ) & 361 & + zds(ji,jj-1) * zds(ji,jj-1) * e12f(ji,jj-1) + zds(ji-1,jj-1) * zds(ji-1,jj-1) * e12f(ji-1,jj-1) & 362 & ) * 0.25_wp * r1_e12t(ji,jj) 363 364 ! divergence at T points 365 zdiv = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj) & 366 & + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1) & 367 & ) * r1_e12t(ji,jj) 368 zdiv2 = zdiv * zdiv 369 370 ! tension at T points 371 zdt = ( ( u_ice(ji,jj) * r1_e2u(ji,jj) - u_ice(ji-1,jj) * r1_e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj) & 372 & - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj) & 373 & ) * r1_e12t(ji,jj) 374 zdt2 = zdt * zdt 375 376 ! delta at T points 377 zdelta = SQRT( zdiv2 + ( zdt2 + zds2 ) * usecc2 ) 378 379 ! P/delta at T points 380 zp_delt(ji,jj) = zpresh(ji,jj) / ( zdelta + rn_creepl ) 381 382 ! stress at T points 383 zs1(ji,jj) = ( zs1(ji,jj) * zalph1 + zp_delt(ji,jj) * ( zdiv - zdelta ) ) * z1_alph1 384 zs2(ji,jj) = ( zs2(ji,jj) * zalph2 + zp_delt(ji,jj) * ( zdt * z1_ecc2 ) ) * z1_alph2 385 386 END DO 387 END DO 388 CALL lbc_lnk( zp_delt, 'T', 1. ) 389 390 DO jj = k_j1, k_jpj-1 391 DO ji = 1, jpim1 392 393 ! P/delta at F points 394 zp_delf = 0.25_wp * ( zp_delt(ji,jj) + zp_delt(ji+1,jj) + zp_delt(ji,jj+1) + zp_delt(ji+1,jj+1) ) 395 396 ! stress at F points 397 zs12(ji,jj)= ( zs12(ji,jj) * zalph2 + zp_delf * ( zds(ji,jj) * z1_ecc2 ) * 0.5_wp ) * z1_alph2 398 399 END DO 400 END DO 401 CALL lbc_lnk_multi( zs1, 'T', 1., zs2, 'T', 1., zs12, 'F', 1. ) 417 402 418 ! Ice internal stresses (Appendix C of Hunke and Dukowicz, 2002)403 ! --- Ice internal stresses (Appendix C of Hunke and Dukowicz, 2002) --- ! 419 404 DO jj = k_j1+1, k_jpj-1 420 DO ji = fs_2, fs_jpim1 421 !- contribution of zs1, zs2 and zs12 to zf1 422 zf1(ji,jj) = 0.5 * ( ( zs1(ji+1,jj) - zs1(ji,jj) ) * e2u(ji,jj) & 423 & + ( zs2(ji+1,jj) * e2t(ji+1,jj)**2 - zs2(ji,jj) * e2t(ji,jj)**2 ) * r1_e2u(ji,jj) & 424 & + 2.0 * ( zs12(ji,jj) * e1f(ji,jj)**2 - zs12(ji,jj-1) * e1f(ji,jj-1)**2 ) * r1_e1u(ji,jj) & 425 & ) * r1_e12u(ji,jj) 426 ! contribution of zs1, zs2 and zs12 to zf2 427 zf2(ji,jj) = 0.5 * ( ( zs1(ji,jj+1) - zs1(ji,jj) ) * e1v(ji,jj) & 428 & - ( zs2(ji,jj+1) * e1t(ji,jj+1)**2 - zs2(ji,jj) * e1t(ji,jj)**2 ) * r1_e1v(ji,jj) & 429 & + 2.0 * ( zs12(ji,jj) * e2f(ji,jj)**2 - zs12(ji-1,jj) * e2f(ji-1,jj)**2 ) * r1_e2v(ji,jj) & 430 & ) * r1_e12v(ji,jj) 405 DO ji = fs_2, fs_jpim1 406 407 ! U points 408 zfU(ji,jj) = 0.5_wp * ( ( zs1(ji+1,jj) - zs1(ji,jj) ) * e2u(ji,jj) & 409 & + ( zs2(ji+1,jj) * e2t(ji+1,jj) * e2t(ji+1,jj) - zs2(ji,jj) * e2t(ji,jj) * e2t(ji,jj) & 410 & ) * r1_e2u(ji,jj) & 411 & + ( zs12(ji,jj) * e1f(ji,jj) * e1f(ji,jj) - zs12(ji,jj-1) * e1f(ji,jj-1) * e1f(ji,jj-1) & 412 & ) * 2._wp * r1_e1u(ji,jj) & 413 & ) * r1_e12u(ji,jj) 414 415 ! V points 416 zfV(ji,jj) = 0.5_wp * ( ( zs1(ji,jj+1) - zs1(ji,jj) ) * e1v(ji,jj) & 417 & - ( zs2(ji,jj+1) * e1t(ji,jj+1) * e1t(ji,jj+1) - zs2(ji,jj) * e1t(ji,jj) * e1t(ji,jj) & 418 & ) * r1_e1v(ji,jj) & 419 & + ( zs12(ji,jj) * e2f(ji,jj) * e2f(ji,jj) - zs12(ji-1,jj) * e2f(ji-1,jj) * e2f(ji-1,jj) & 420 & ) * 2._wp * r1_e2v(ji,jj) & 421 & ) * r1_e12v(ji,jj) 422 423 ! u_ice at V point 424 u_iceV(ji,jj) = 0.5_wp * ( ( u_ice(ji,jj ) + u_ice(ji-1,jj ) ) * e2t(ji,jj+1) & 425 & + ( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj ) ) * z1_e2t0(ji,jj) * vmask(ji,jj,1) 426 427 ! v_ice at U point 428 v_iceU(ji,jj) = 0.5_wp * ( ( v_ice(ji ,jj) + v_ice(ji ,jj-1) ) * e1t(ji+1,jj) & 429 & + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji ,jj) ) * z1_e1t0(ji,jj) * umask(ji,jj,1) 430 431 431 END DO 432 432 END DO 433 433 ! 434 ! Computation of ice velocity 435 ! 436 ! Both the Coriolis term and the ice-ocean drag are solved semi-implicitly. 437 ! 438 IF (MOD(jter,2).eq.0) THEN 439 434 ! --- Computation of ice velocity --- ! 435 ! Bouillon et al. 2013 (eq 47-48) => unstable unless alpha, beta are chosen wisely and large nn_nevp 436 ! Bouillon et al. 2009 (eq 34-35) => stable 437 IF( MOD(jter,2) .EQ. 0 ) THEN ! even iterations 438 440 439 DO jj = k_j1+1, k_jpj-1 441 440 DO ji = fs_2, fs_jpim1 442 rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zmass1(ji,jj) ) ) ) * umask(ji,jj,1) 443 z0 = zmass1(ji,jj) * z1_dtevp 444 445 ! SB modif because ocean has no slip boundary condition 446 zv_ice1 = 0.5 * ( ( v_ice(ji ,jj) + v_ice(ji ,jj-1) ) * e1t(ji ,jj) & 447 & + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji+1,jj) ) & 448 & / ( e1t(ji+1,jj) + e1t(ji,jj) ) * umask(ji,jj,1) 449 za = rhoco * SQRT( ( u_ice(ji,jj) - u_oce(ji,jj) )**2 + & 450 & ( zv_ice1 - v_oce1(ji,jj) )**2 ) * ( 1.0 - zfrld1(ji,jj) ) 451 zr = z0 * u_ice(ji,jj) + zf1(ji,jj) + za1ct(ji,jj) + za * u_oce(ji,jj) 452 zcca = z0 + za 453 zccb = zcorl1(ji,jj) 454 u_ice(ji,jj) = ( zr + zccb * zv_ice1 ) / ( zcca + zepsi ) * rswitch 441 442 ! tau_io/(v_oce - v_ice) 443 zTauO = zaV(ji,jj) * rhoco * SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) ) & 444 & + ( u_iceV(ji,jj) - u_oceV(ji,jj) ) * ( u_iceV(ji,jj) - u_oceV(ji,jj) ) ) 445 446 ! Coriolis at V-points (energy conserving formulation) 447 zCor = - 0.25_wp * r1_e2v(ji,jj) * & 448 & ( zmf(ji,jj ) * ( e2u(ji,jj ) * u_ice(ji,jj ) + e2u(ji-1,jj ) * u_ice(ji-1,jj ) ) & 449 & + zmf(ji,jj+1) * ( e2u(ji,jj+1) * u_ice(ji,jj+1) + e2u(ji-1,jj+1) * u_ice(ji-1,jj+1) ) ) 450 451 ! Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io 452 zTauE = zfV(ji,jj) + zTauV_ia(ji,jj) + zCor + zspgV(ji,jj) + zTauO * ( v_oce(ji,jj) - v_ice(ji,jj) ) 453 454 ! ice velocity using implicit formulation (cf Madec doc & Bouillon 2009) 455 v_ice(ji,jj) = ( ( zmV_t(ji,jj) * v_ice(ji,jj) + zTauE + zTauO * v_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 456 & ) / MAX( zepsi, zmV_t(ji,jj) + zTauO ) * zswitchV(ji,jj) & ! m/dt + tau_io(only ice part) 457 & + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) ) & ! v_ice = v_oce if mass < zmmin 458 & ) * zmaskV(ji,jj) 455 459 END DO 456 460 END DO 457 458 CALL lbc_lnk( u_ice(:,:), 'U', -1. ) 461 CALL lbc_lnk( v_ice, 'V', -1. ) 462 463 #if defined key_agrif && defined key_lim2 464 CALL agrif_rhg_lim2( jter, nn_nevp, 'V' ) 465 #endif 466 #if defined key_bdy 467 CALL bdy_ice_lim_dyn( 'V' ) 468 #endif 469 470 DO jj = k_j1+1, k_jpj-1 471 DO ji = fs_2, fs_jpim1 472 473 ! tau_io/(u_oce - u_ice) 474 zTauO = zaU(ji,jj) * rhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) ) & 475 & + ( v_iceU(ji,jj) - v_oceU(ji,jj) ) * ( v_iceU(ji,jj) - v_oceU(ji,jj) ) ) 476 477 ! Coriolis at U-points (energy conserving formulation) 478 zCor = 0.25_wp * r1_e1u(ji,jj) * & 479 & ( zmf(ji ,jj) * ( e1v(ji ,jj) * v_ice(ji ,jj) + e1v(ji ,jj-1) * v_ice(ji ,jj-1) ) & 480 & + zmf(ji+1,jj) * ( e1v(ji+1,jj) * v_ice(ji+1,jj) + e1v(ji+1,jj-1) * v_ice(ji+1,jj-1) ) ) 481 482 ! Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io 483 zTauE = zfU(ji,jj) + zTauU_ia(ji,jj) + zCor + zspgU(ji,jj) + zTauO * ( u_oce(ji,jj) - u_ice(ji,jj) ) 484 485 ! ice velocity using implicit formulation (cf Madec doc & Bouillon 2009) 486 u_ice(ji,jj) = ( ( zmU_t(ji,jj) * u_ice(ji,jj) + zTauE + zTauO * u_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 487 & ) / MAX( zepsi, zmU_t(ji,jj) + zTauO ) * zswitchU(ji,jj) & ! m/dt + tau_io(only ice part) 488 & + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) ) & ! v_ice = v_oce if mass < zmmin 489 & ) * zmaskU(ji,jj) 490 END DO 491 END DO 492 CALL lbc_lnk( u_ice, 'U', -1. ) 493 459 494 #if defined key_agrif && defined key_lim2 460 495 CALL agrif_rhg_lim2( jter, nn_nevp, 'U' ) 461 496 #endif 462 497 #if defined key_bdy 463 CALL bdy_ice_lim_dyn( 'U' )498 CALL bdy_ice_lim_dyn( 'U' ) 464 499 #endif 500 501 ELSE ! odd iterations 465 502 466 503 DO jj = k_j1+1, k_jpj-1 467 504 DO ji = fs_2, fs_jpim1 468 469 rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zmass2(ji,jj) ) ) ) * vmask(ji,jj,1) 470 z0 = zmass2(ji,jj) * z1_dtevp 471 ! SB modif because ocean has no slip boundary condition 472 zu_ice2 = 0.5 * ( ( u_ice(ji,jj ) + u_ice(ji-1,jj ) ) * e2t(ji,jj) & 473 & + ( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj+1) ) & 474 & / ( e2t(ji,jj+1) + e2t(ji,jj) ) * vmask(ji,jj,1) 475 za = rhoco * SQRT( ( zu_ice2 - u_oce2(ji,jj) )**2 + & 476 & ( v_ice(ji,jj) - v_oce(ji,jj))**2 ) * ( 1.0 - zfrld2(ji,jj) ) 477 zr = z0 * v_ice(ji,jj) + zf2(ji,jj) + za2ct(ji,jj) + za * v_oce(ji,jj) 478 zcca = z0 + za 479 zccb = zcorl2(ji,jj) 480 v_ice(ji,jj) = ( zr - zccb * zu_ice2 ) / ( zcca + zepsi ) * rswitch 505 506 ! tau_io/(u_oce - u_ice) 507 zTauO = zaU(ji,jj) * rhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) ) & 508 & + ( v_iceU(ji,jj) - v_oceU(ji,jj) ) * ( v_iceU(ji,jj) - v_oceU(ji,jj) ) ) 509 510 ! Coriolis at U-points (energy conserving formulation) 511 zCor = 0.25_wp * r1_e1u(ji,jj) * & 512 & ( zmf(ji ,jj) * ( e1v(ji ,jj) * v_ice(ji ,jj) + e1v(ji ,jj-1) * v_ice(ji ,jj-1) ) & 513 & + zmf(ji+1,jj) * ( e1v(ji+1,jj) * v_ice(ji+1,jj) + e1v(ji+1,jj-1) * v_ice(ji+1,jj-1) ) ) 514 515 ! Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io 516 zTauE = zfU(ji,jj) + zTauU_ia(ji,jj) + zCor + zspgU(ji,jj) + zTauO * ( u_oce(ji,jj) - u_ice(ji,jj) ) 517 518 ! ice velocity using implicit formulation (cf Madec doc & Bouillon 2009) 519 u_ice(ji,jj) = ( ( zmU_t(ji,jj) * u_ice(ji,jj) + zTauE + zTauO * u_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 520 & ) / MAX( zepsi, zmU_t(ji,jj) + zTauO ) * zswitchU(ji,jj) & ! m/dt + tau_io(only ice part) 521 & + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) ) & ! v_ice = v_oce if mass < zmmin 522 & ) * zmaskU(ji,jj) 481 523 END DO 482 524 END DO 483 484 CALL lbc_lnk( v_ice(:,:), 'V', -1. ) 525 CALL lbc_lnk( u_ice, 'U', -1. ) 526 527 #if defined key_agrif && defined key_lim2 528 CALL agrif_rhg_lim2( jter, nn_nevp, 'U' ) 529 #endif 530 #if defined key_bdy 531 CALL bdy_ice_lim_dyn( 'U' ) 532 #endif 533 534 DO jj = k_j1+1, k_jpj-1 535 DO ji = fs_2, fs_jpim1 536 537 ! tau_io/(v_oce - v_ice) 538 zTauO = zaV(ji,jj) * rhoco * SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) ) & 539 & + ( u_iceV(ji,jj) - u_oceV(ji,jj) ) * ( u_iceV(ji,jj) - u_oceV(ji,jj) ) ) 540 541 ! Coriolis at V-points (energy conserving formulation) 542 zCor = - 0.25_wp * r1_e2v(ji,jj) * & 543 & ( zmf(ji,jj ) * ( e2u(ji,jj ) * u_ice(ji,jj ) + e2u(ji-1,jj ) * u_ice(ji-1,jj ) ) & 544 & + zmf(ji,jj+1) * ( e2u(ji,jj+1) * u_ice(ji,jj+1) + e2u(ji-1,jj+1) * u_ice(ji-1,jj+1) ) ) 545 546 ! Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io 547 zTauE = zfV(ji,jj) + zTauV_ia(ji,jj) + zCor + zspgV(ji,jj) + zTauO * ( v_oce(ji,jj) - v_ice(ji,jj) ) 548 549 ! ice velocity using implicit formulation (cf Madec doc & Bouillon 2009) 550 v_ice(ji,jj) = ( ( zmV_t(ji,jj) * v_ice(ji,jj) + zTauE + zTauO * v_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 551 & ) / MAX( zepsi, zmV_t(ji,jj) + zTauO ) * zswitchV(ji,jj) & ! m/dt + tau_io(only ice part) 552 & + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) ) & ! v_ice = v_oce if mass < zmmin 553 & ) * zmaskV(ji,jj) 554 END DO 555 END DO 556 CALL lbc_lnk( v_ice, 'V', -1. ) 557 485 558 #if defined key_agrif && defined key_lim2 486 559 CALL agrif_rhg_lim2( jter, nn_nevp, 'V' ) 487 560 #endif 488 561 #if defined key_bdy 489 CALL bdy_ice_lim_dyn( 'V' )562 CALL bdy_ice_lim_dyn( 'V' ) 490 563 #endif 491 564 492 ELSE 565 ENDIF 566 567 IF(ln_ctl) THEN ! Convergence test 493 568 DO jj = k_j1+1, k_jpj-1 494 DO ji = fs_2, fs_jpim1495 rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zmass2(ji,jj) ) ) ) * vmask(ji,jj,1)496 z0 = zmass2(ji,jj) * z1_dtevp497 ! SB modif because ocean has no slip boundary condition498 zu_ice2 = 0.5 * ( ( u_ice(ji,jj ) + u_ice(ji-1,jj ) ) * e2t(ji,jj) &499 & +( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj+1) ) &500 & / ( e2t(ji,jj+1) + e2t(ji,jj) ) * vmask(ji,jj,1)501 502 za = rhoco * SQRT( ( zu_ice2 - u_oce2(ji,jj) )**2 + &503 & ( v_ice(ji,jj) - v_oce(ji,jj) )**2 ) * ( 1.0 - zfrld2(ji,jj) )504 zr = z0 * v_ice(ji,jj) + zf2(ji,jj) + za2ct(ji,jj) + za * v_oce(ji,jj)505 zcca = z0 + za506 zccb = zcorl2(ji,jj)507 v_ice(ji,jj) = ( zr - zccb * zu_ice2 ) / ( zcca + zepsi ) * rswitch508 END DO509 END DO510 511 CALL lbc_lnk( v_ice(:,:), 'V', -1. )512 #if defined key_agrif && defined key_lim2513 CALL agrif_rhg_lim2( jter, nn_nevp, 'V' )514 #endif515 #if defined key_bdy516 CALL bdy_ice_lim_dyn( 'V' )517 #endif518 519 DO jj = k_j1+1, k_jpj-1520 DO ji = fs_2, fs_jpim1521 rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zmass1(ji,jj) ) ) ) * umask(ji,jj,1)522 z0 = zmass1(ji,jj) * z1_dtevp523 zv_ice1 = 0.5 * ( ( v_ice(ji ,jj) + v_ice(ji ,jj-1) ) * e1t(ji,jj) &524 & + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji+1,jj) ) &525 & / ( e1t(ji+1,jj) + e1t(ji,jj) ) * umask(ji,jj,1)526 527 za = rhoco * SQRT( ( u_ice(ji,jj) - u_oce(ji,jj) )**2 + &528 & ( zv_ice1 - v_oce1(ji,jj) )**2 ) * ( 1.0 - zfrld1(ji,jj) )529 zr = z0 * u_ice(ji,jj) + zf1(ji,jj) + za1ct(ji,jj) + za * u_oce(ji,jj)530 zcca = z0 + za531 zccb = zcorl1(ji,jj)532 u_ice(ji,jj) = ( zr + zccb * zv_ice1 ) / ( zcca + zepsi ) * rswitch533 END DO534 END DO535 536 CALL lbc_lnk( u_ice(:,:), 'U', -1. )537 #if defined key_agrif && defined key_lim2538 CALL agrif_rhg_lim2( jter, nn_nevp, 'U' )539 #endif540 #if defined key_bdy541 CALL bdy_ice_lim_dyn( 'U' )542 #endif543 544 ENDIF545 546 IF(ln_ctl) THEN547 !--- Convergence test.548 DO jj = k_j1+1 , k_jpj-1549 569 zresr(:,jj) = MAX( ABS( u_ice(:,jj) - zu_ice(:,jj) ), ABS( v_ice(:,jj) - zv_ice(:,jj) ) ) 550 570 END DO … … 552 572 IF( lk_mpp ) CALL mpp_max( zresm ) ! max over the global domain 553 573 ENDIF 554 574 ! 555 575 ! ! ==================== ! 556 576 END DO ! end loop over jter ! … … 558 578 ! 559 579 !------------------------------------------------------------------------------! 560 ! 4) Prevent ice velocities when the ice is thin 561 !------------------------------------------------------------------------------! 562 ! If the ice volume is below zvmin then ice velocity should equal the 563 ! ocean velocity. This prevents high velocity when ice is thin 564 DO jj = k_j1+1, k_jpj-1 565 DO ji = fs_2, fs_jpim1 566 IF ( vt_i(ji,jj) <= zvmin ) THEN 567 u_ice(ji,jj) = u_oce(ji,jj) 568 v_ice(ji,jj) = v_oce(ji,jj) 569 ENDIF 580 ! 4) Recompute delta, shear and div (inputs for mechanical redistribution) 581 !------------------------------------------------------------------------------! 582 DO jj = k_j1, k_jpj-1 583 DO ji = 1, jpim1 584 585 ! shear at F points 586 zds(ji,jj) = ( ( u_ice(ji,jj+1) * r1_e1u(ji,jj+1) - u_ice(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj) & 587 & + ( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj) & 588 & ) * r1_e12f(ji,jj) * zfmask(ji,jj) 589 590 END DO 591 END DO 592 CALL lbc_lnk( zds, 'F', 1. ) 593 594 DO jj = k_j1+1, k_jpj-1 595 DO ji = 2, jpim1 ! no vector loop 596 597 ! tension**2 at T points 598 zdt = ( ( u_ice(ji,jj) * r1_e2u(ji,jj) - u_ice(ji-1,jj) * r1_e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj) & 599 & - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj) & 600 & ) * r1_e12t(ji,jj) 601 zdt2 = zdt * zdt 602 603 ! shear**2 at T points (doc eq. A16) 604 zds2 = ( zds(ji,jj ) * zds(ji,jj ) * e12f(ji,jj ) + zds(ji-1,jj ) * zds(ji-1,jj ) * e12f(ji-1,jj ) & 605 & + zds(ji,jj-1) * zds(ji,jj-1) * e12f(ji,jj-1) + zds(ji-1,jj-1) * zds(ji-1,jj-1) * e12f(ji-1,jj-1) & 606 & ) * 0.25_wp * r1_e12t(ji,jj) 607 608 ! shear at T points 609 shear_i(ji,jj) = SQRT( zdt2 + zds2 ) 610 611 ! divergence at T points 612 divu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj) & 613 & + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1) & 614 & ) * r1_e12t(ji,jj) 615 616 ! delta at T points 617 zdelta = SQRT( divu_i(ji,jj) * divu_i(ji,jj) + ( zdt2 + zds2 ) * usecc2 ) 618 rswitch = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zdelta ) ) ! 0 if delta=0 619 delta_i(ji,jj) = zdelta + rn_creepl * rswitch 620 570 621 END DO 571 622 END DO 572 573 CALL lbc_lnk_multi( u_ice(:,:), 'U', -1., v_ice(:,:), 'V', -1. ) 574 575 #if defined key_agrif && defined key_lim2 576 CALL agrif_rhg_lim2( nn_nevp , nn_nevp, 'U' ) 577 CALL agrif_rhg_lim2( nn_nevp , nn_nevp, 'V' ) 578 #endif 579 #if defined key_bdy 580 CALL bdy_ice_lim_dyn( 'U' ) 581 CALL bdy_ice_lim_dyn( 'V' ) 582 #endif 583 584 DO jj = k_j1+1, k_jpj-1 585 DO ji = fs_2, fs_jpim1 586 IF ( vt_i(ji,jj) <= zvmin ) THEN 587 v_ice1(ji,jj) = 0.5_wp * ( ( v_ice(ji ,jj) + v_ice(ji, jj-1) ) * e1t(ji+1,jj) & 588 & + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji ,jj) ) & 589 & / ( e1t(ji+1,jj) + e1t(ji,jj) ) * umask(ji,jj,1) 590 591 u_ice2(ji,jj) = 0.5_wp * ( ( u_ice(ji,jj ) + u_ice(ji-1,jj ) ) * e2t(ji,jj+1) & 592 & + ( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj ) ) & 593 & / ( e2t(ji,jj+1) + e2t(ji,jj) ) * vmask(ji,jj,1) 594 ENDIF 595 END DO 596 END DO 597 598 CALL lbc_lnk_multi( u_ice2(:,:), 'V', -1., v_ice1(:,:), 'U', -1. ) 599 600 ! Recompute delta, shear and div, inputs for mechanical redistribution 601 DO jj = k_j1+1, k_jpj-1 602 DO ji = fs_2, jpim1 !RB bug no vect opt due to zmask 603 !- divu_i(:,:), zdt(:,:): divergence and tension at centre 604 !- zds(:,:): shear on northeast corner of grid cells 605 IF ( vt_i(ji,jj) <= zvmin ) THEN 606 607 divu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj ) * u_ice(ji-1,jj ) & 608 & + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji ,jj-1) * v_ice(ji ,jj-1) & 609 & ) * r1_e12t(ji,jj) 610 611 zdt(ji,jj) = ( ( u_ice(ji,jj) * r1_e2u(ji,jj) - u_ice(ji-1,jj) * r1_e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj) & 612 & -( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj) & 613 & ) * r1_e12t(ji,jj) 614 ! 615 ! SB modif because ocean has no slip boundary condition 616 zds(ji,jj) = ( ( u_ice(ji,jj+1) * r1_e1u(ji,jj+1) - u_ice(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj) & 617 & +( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj) & 618 & ) * r1_e12f(ji,jj) * ( 2.0 - fmask(ji,jj,1) ) & 619 & * zmask(ji,jj) * zmask(ji,jj+1) * zmask(ji+1,jj) * zmask(ji+1,jj+1) 620 621 zdst = ( e2u(ji,jj) * v_ice1(ji,jj) - e2u(ji-1,jj ) * v_ice1(ji-1,jj ) & 622 & + e1v(ji,jj) * u_ice2(ji,jj) - e1v(ji ,jj-1) * u_ice2(ji ,jj-1) ) * r1_e12t(ji,jj) 623 624 delta = SQRT( divu_i(ji,jj)**2 + ( zdt(ji,jj)**2 + zdst**2 ) * usecc2 ) 625 delta_i(ji,jj) = delta + rn_creepl 626 627 ENDIF 628 END DO 629 END DO 630 ! 631 !------------------------------------------------------------------------------! 632 ! 5) Store stress tensor and its invariants 633 !------------------------------------------------------------------------------! 634 ! * Invariants of the stress tensor are required for limitd_me 635 ! (accelerates convergence and improves stability) 636 DO jj = k_j1+1, k_jpj-1 637 DO ji = fs_2, fs_jpim1 638 zdst = ( e2u(ji,jj) * v_ice1(ji,jj) - e2u( ji-1, jj ) * v_ice1(ji-1,jj) & 639 & + e1v(ji,jj) * u_ice2(ji,jj) - e1v( ji , jj-1 ) * u_ice2(ji,jj-1) ) * r1_e12t(ji,jj) 640 shear_i(ji,jj) = SQRT( zdt(ji,jj) * zdt(ji,jj) + zdst * zdst ) 641 END DO 642 END DO 643 644 ! Lateral boundary condition 645 CALL lbc_lnk_multi( divu_i (:,:), 'T', 1., delta_i(:,:), 'T', 1., shear_i(:,:), 'T', 1. ) 646 647 ! * Store the stress tensor for the next time step 623 CALL lbc_lnk_multi( shear_i, 'T', 1., divu_i, 'T', 1., delta_i, 'T', 1. ) 624 625 ! --- Store the stress tensor for the next time step --- ! 648 626 stress1_i (:,:) = zs1 (:,:) 649 627 stress2_i (:,:) = zs2 (:,:) … … 652 630 ! 653 631 !------------------------------------------------------------------------------! 654 ! 6) Control prints of residual and charge ellipse632 ! 5) Control prints of residual and charge ellipse 655 633 !------------------------------------------------------------------------------! 656 634 ! … … 675 653 DO ji = 2, jpim1 676 654 IF (zpresh(ji,jj) > 1.0) THEN 677 sigma1 = ( zs1(ji,jj) + (zs2(ji,jj)**2 + 4*zs12(ji,jj)**2 )**0.5 ) / ( 2*zpresh(ji,jj) )678 sigma2 = ( zs1(ji,jj) - (zs2(ji,jj)**2 + 4*zs12(ji,jj)**2 )**0.5 ) / ( 2*zpresh(ji,jj) )655 zsig1 = ( zs1(ji,jj) + (zs2(ji,jj)**2 + 4*zs12(ji,jj)**2 )**0.5 ) / ( 2*zpresh(ji,jj) ) 656 zsig2 = ( zs1(ji,jj) - (zs2(ji,jj)**2 + 4*zs12(ji,jj)**2 )**0.5 ) / ( 2*zpresh(ji,jj) ) 679 657 WRITE(charout,FMT="('lim_rhg :', I4, I4, D23.16, D23.16, D23.16, D23.16, A10)") 680 658 CALL prt_ctl_info(charout) … … 687 665 ENDIF 688 666 ! 689 CALL wrk_dealloc( jpi,jpj, zpresh, zfrld1, zmass1, zcorl1, za1ct , zpreshc, zfrld2, zmass2, zcorl2, za2ct ) 690 CALL wrk_dealloc( jpi,jpj, u_oce2, u_ice2, v_oce1 , v_ice1 , zmask ) 691 CALL wrk_dealloc( jpi,jpj, zf1 , zu_ice, zf2 , zv_ice , zdt , zds ) 692 CALL wrk_dealloc( jpi,jpj, zs1 , zs2 , zs12 , zresr , zpice ) 667 CALL wrk_dealloc( jpi,jpj, zpresh, z1_e1t0, z1_e2t0, zp_delt ) 668 CALL wrk_dealloc( jpi,jpj, zaU, zaV, zmU_t, zmV_t, zmf, zTauU_ia, ztauV_ia ) 669 CALL wrk_dealloc( jpi,jpj, zspgU, zspgV, v_oceU, u_oceV, v_iceU, u_iceV, zfU, zfV ) 670 CALL wrk_dealloc( jpi,jpj, zds, zs1, zs2, zs12, zu_ice, zv_ice, zresr, zpice ) 671 CALL wrk_dealloc( jpi,jpj, zswitchU, zswitchV, zmaskU, zmaskV, zfmask, zwf ) 693 672 694 673 END SUBROUTINE lim_rhg -
branches/UKMO/dev_r5518_GO6_package_lim3/NEMOGCM/NEMO/LIM_SRC_3/limrst.F90
r6486 r8540 40 40 !!---------------------------------------------------------------------- 41 41 !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) 42 !! $Id $42 !! $Id: limrst.F90 7814 2017-03-20 16:21:42Z clem $ 43 43 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 44 44 !!---------------------------------------------------------------------- … … 108 108 INTEGER :: iter 109 109 CHARACTER(len=15) :: znam 110 CHARACTER(len= 1) :: zchar, zchar1110 CHARACTER(len=2) :: zchar, zchar1 111 111 REAL(wp), POINTER, DIMENSION(:,:) :: z2d 112 112 !!---------------------------------------------------------------------- … … 130 130 ! Prognostic variables 131 131 DO jl = 1, jpl 132 WRITE(zchar,'(I 1)') jl133 znam = 'v_i'//'_htc'// zchar132 WRITE(zchar,'(I2)') jl 133 znam = 'v_i'//'_htc'//TRIM(ADJUSTL(zchar)) 134 134 z2d(:,:) = v_i(:,:,jl) 135 135 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 136 znam = 'v_s'//'_htc'// zchar136 znam = 'v_s'//'_htc'//TRIM(ADJUSTL(zchar)) 137 137 z2d(:,:) = v_s(:,:,jl) 138 138 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 139 znam = 'smv_i'//'_htc'// zchar139 znam = 'smv_i'//'_htc'//TRIM(ADJUSTL(zchar)) 140 140 z2d(:,:) = smv_i(:,:,jl) 141 141 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 142 znam = 'oa_i'//'_htc'// zchar142 znam = 'oa_i'//'_htc'//TRIM(ADJUSTL(zchar)) 143 143 z2d(:,:) = oa_i(:,:,jl) 144 144 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 145 znam = 'a_i'//'_htc'// zchar145 znam = 'a_i'//'_htc'//TRIM(ADJUSTL(zchar)) 146 146 z2d(:,:) = a_i(:,:,jl) 147 147 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 148 znam = 't_su'//'_htc'// zchar148 znam = 't_su'//'_htc'//TRIM(ADJUSTL(zchar)) 149 149 z2d(:,:) = t_su(:,:,jl) 150 150 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 151 END DO 152 153 DO jl = 1, jpl 154 WRITE(zchar,'(I1)') jl 155 znam = 'tempt_sl1'//'_htc'//zchar 151 znam = 'tempt_sl1'//'_htc'//TRIM(ADJUSTL(zchar)) 156 152 z2d(:,:) = e_s(:,:,1,jl) 157 153 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 158 END DO159 160 DO jl = 1, jpl161 WRITE(zchar,'(I1)') jl162 154 DO jk = 1, nlay_i 163 WRITE(zchar1,'(I 1)') jk164 znam = 'tempt'//'_il'// zchar1//'_htc'//zchar155 WRITE(zchar1,'(I2)') jk 156 znam = 'tempt'//'_il'//TRIM(ADJUSTL(zchar1))//'_htc'//TRIM(ADJUSTL(zchar)) 165 157 z2d(:,:) = e_i(:,:,jk,jl) 166 158 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) … … 177 169 178 170 DO jl = 1, jpl 179 WRITE(zchar,'(I 1)') jl180 znam = 'sxice'//'_htc'// zchar171 WRITE(zchar,'(I2)') jl 172 znam = 'sxice'//'_htc'//TRIM(ADJUSTL(zchar)) 181 173 z2d(:,:) = sxice(:,:,jl) 182 174 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 183 znam = 'syice'//'_htc'// zchar175 znam = 'syice'//'_htc'//TRIM(ADJUSTL(zchar)) 184 176 z2d(:,:) = syice(:,:,jl) 185 177 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 186 znam = 'sxxice'//'_htc'// zchar178 znam = 'sxxice'//'_htc'//TRIM(ADJUSTL(zchar)) 187 179 z2d(:,:) = sxxice(:,:,jl) 188 180 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 189 znam = 'syyice'//'_htc'// zchar181 znam = 'syyice'//'_htc'//TRIM(ADJUSTL(zchar)) 190 182 z2d(:,:) = syyice(:,:,jl) 191 183 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 192 znam = 'sxyice'//'_htc'// zchar184 znam = 'sxyice'//'_htc'//TRIM(ADJUSTL(zchar)) 193 185 z2d(:,:) = sxyice(:,:,jl) 194 186 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 195 znam = 'sxsn'//'_htc'// zchar187 znam = 'sxsn'//'_htc'//TRIM(ADJUSTL(zchar)) 196 188 z2d(:,:) = sxsn(:,:,jl) 197 189 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 198 znam = 'sysn'//'_htc'// zchar190 znam = 'sysn'//'_htc'//TRIM(ADJUSTL(zchar)) 199 191 z2d(:,:) = sysn(:,:,jl) 200 192 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 201 znam = 'sxxsn'//'_htc'// zchar193 znam = 'sxxsn'//'_htc'//TRIM(ADJUSTL(zchar)) 202 194 z2d(:,:) = sxxsn(:,:,jl) 203 195 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 204 znam = 'syysn'//'_htc'// zchar196 znam = 'syysn'//'_htc'//TRIM(ADJUSTL(zchar)) 205 197 z2d(:,:) = syysn(:,:,jl) 206 198 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 207 znam = 'sxysn'//'_htc'// zchar199 znam = 'sxysn'//'_htc'//TRIM(ADJUSTL(zchar)) 208 200 z2d(:,:) = sxysn(:,:,jl) 209 201 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 210 znam = 'sxa'//'_htc'// zchar202 znam = 'sxa'//'_htc'//TRIM(ADJUSTL(zchar)) 211 203 z2d(:,:) = sxa(:,:,jl) 212 204 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 213 znam = 'sya'//'_htc'// zchar205 znam = 'sya'//'_htc'//TRIM(ADJUSTL(zchar)) 214 206 z2d(:,:) = sya(:,:,jl) 215 207 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 216 znam = 'sxxa'//'_htc'// zchar208 znam = 'sxxa'//'_htc'//TRIM(ADJUSTL(zchar)) 217 209 z2d(:,:) = sxxa(:,:,jl) 218 210 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 219 znam = 'syya'//'_htc'// zchar211 znam = 'syya'//'_htc'//TRIM(ADJUSTL(zchar)) 220 212 z2d(:,:) = syya(:,:,jl) 221 213 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 222 znam = 'sxya'//'_htc'// zchar214 znam = 'sxya'//'_htc'//TRIM(ADJUSTL(zchar)) 223 215 z2d(:,:) = sxya(:,:,jl) 224 216 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 225 znam = 'sxc0'//'_htc'// zchar217 znam = 'sxc0'//'_htc'//TRIM(ADJUSTL(zchar)) 226 218 z2d(:,:) = sxc0(:,:,jl) 227 219 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 228 znam = 'syc0'//'_htc'// zchar220 znam = 'syc0'//'_htc'//TRIM(ADJUSTL(zchar)) 229 221 z2d(:,:) = syc0(:,:,jl) 230 222 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 231 znam = 'sxxc0'//'_htc'// zchar223 znam = 'sxxc0'//'_htc'//TRIM(ADJUSTL(zchar)) 232 224 z2d(:,:) = sxxc0(:,:,jl) 233 225 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 234 znam = 'syyc0'//'_htc'// zchar226 znam = 'syyc0'//'_htc'//TRIM(ADJUSTL(zchar)) 235 227 z2d(:,:) = syyc0(:,:,jl) 236 228 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 237 znam = 'sxyc0'//'_htc'// zchar229 znam = 'sxyc0'//'_htc'//TRIM(ADJUSTL(zchar)) 238 230 z2d(:,:) = sxyc0(:,:,jl) 239 231 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 240 znam = 'sxsal'//'_htc'// zchar232 znam = 'sxsal'//'_htc'//TRIM(ADJUSTL(zchar)) 241 233 z2d(:,:) = sxsal(:,:,jl) 242 234 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 243 znam = 'sysal'//'_htc'// zchar235 znam = 'sysal'//'_htc'//TRIM(ADJUSTL(zchar)) 244 236 z2d(:,:) = sysal(:,:,jl) 245 237 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 246 znam = 'sxxsal'//'_htc'// zchar238 znam = 'sxxsal'//'_htc'//TRIM(ADJUSTL(zchar)) 247 239 z2d(:,:) = sxxsal(:,:,jl) 248 240 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 249 znam = 'syysal'//'_htc'// zchar241 znam = 'syysal'//'_htc'//TRIM(ADJUSTL(zchar)) 250 242 z2d(:,:) = syysal(:,:,jl) 251 243 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 252 znam = 'sxysal'//'_htc'// zchar244 znam = 'sxysal'//'_htc'//TRIM(ADJUSTL(zchar)) 253 245 z2d(:,:) = sxysal(:,:,jl) 254 246 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 255 znam = 'sxage'//'_htc'// zchar247 znam = 'sxage'//'_htc'//TRIM(ADJUSTL(zchar)) 256 248 z2d(:,:) = sxage(:,:,jl) 257 249 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 258 znam = 'syage'//'_htc'// zchar250 znam = 'syage'//'_htc'//TRIM(ADJUSTL(zchar)) 259 251 z2d(:,:) = syage(:,:,jl) 260 252 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 261 znam = 'sxxage'//'_htc'// zchar253 znam = 'sxxage'//'_htc'//TRIM(ADJUSTL(zchar)) 262 254 z2d(:,:) = sxxage(:,:,jl) 263 255 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 264 znam = 'syyage'//'_htc'// zchar256 znam = 'syyage'//'_htc'//TRIM(ADJUSTL(zchar)) 265 257 z2d(:,:) = syyage(:,:,jl) 266 258 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 267 znam = 'sxyage'//'_htc'// zchar259 znam = 'sxyage'//'_htc'//TRIM(ADJUSTL(zchar)) 268 260 z2d(:,:) = sxyage(:,:,jl) 269 261 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) … … 277 269 278 270 DO jl = 1, jpl 279 WRITE(zchar,'(I 1)') jl271 WRITE(zchar,'(I2)') jl 280 272 DO jk = 1, nlay_i 281 WRITE(zchar1,'(I 1)') jk282 znam = 'sxe'//'_il'// zchar1//'_htc'//zchar273 WRITE(zchar1,'(I2)') jk 274 znam = 'sxe'//'_il'//TRIM(ADJUSTL(zchar1))//'_htc'//TRIM(ADJUSTL(zchar)) 283 275 z2d(:,:) = sxe(:,:,jk,jl) 284 276 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 285 znam = 'sye'//'_il'// zchar1//'_htc'//zchar277 znam = 'sye'//'_il'//TRIM(ADJUSTL(zchar1))//'_htc'//TRIM(ADJUSTL(zchar)) 286 278 z2d(:,:) = sye(:,:,jk,jl) 287 279 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 288 znam = 'sxxe'//'_il'// zchar1//'_htc'//zchar280 znam = 'sxxe'//'_il'//TRIM(ADJUSTL(zchar1))//'_htc'//TRIM(ADJUSTL(zchar)) 289 281 z2d(:,:) = sxxe(:,:,jk,jl) 290 282 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 291 znam = 'syye'//'_il'// zchar1//'_htc'//zchar283 znam = 'syye'//'_il'//TRIM(ADJUSTL(zchar1))//'_htc'//TRIM(ADJUSTL(zchar)) 292 284 z2d(:,:) = syye(:,:,jk,jl) 293 285 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 294 znam = 'sxye'//'_il'// zchar1//'_htc'//zchar286 znam = 'sxye'//'_il'//TRIM(ADJUSTL(zchar1))//'_htc'//TRIM(ADJUSTL(zchar)) 295 287 z2d(:,:) = sxye(:,:,jk,jl) 296 288 CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) … … 318 310 REAL(wp), POINTER, DIMENSION(:,:) :: z2d 319 311 CHARACTER(len=15) :: znam 320 CHARACTER(len= 1) :: zchar, zchar1312 CHARACTER(len=2) :: zchar, zchar1 321 313 INTEGER :: jlibalt = jprstlib 322 314 LOGICAL :: llok … … 356 348 & ' control of time parameter nrstdt' ) 357 349 350 ! Prognostic variables 358 351 DO jl = 1, jpl 359 WRITE(zchar,'(I 1)') jl360 znam = 'v_i'//'_htc'// zchar352 WRITE(zchar,'(I2)') jl 353 znam = 'v_i'//'_htc'//TRIM(ADJUSTL(zchar)) 361 354 CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 362 355 v_i(:,:,jl) = z2d(:,:) 363 znam = 'v_s'//'_htc'// zchar356 znam = 'v_s'//'_htc'//TRIM(ADJUSTL(zchar)) 364 357 CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 365 358 v_s(:,:,jl) = z2d(:,:) 366 znam = 'smv_i'//'_htc'// zchar359 znam = 'smv_i'//'_htc'//TRIM(ADJUSTL(zchar)) 367 360 CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 368 361 smv_i(:,:,jl) = z2d(:,:) 369 znam = 'oa_i'//'_htc'// zchar362 znam = 'oa_i'//'_htc'//TRIM(ADJUSTL(zchar)) 370 363 CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 371 364 oa_i(:,:,jl) = z2d(:,:) 372 znam = 'a_i'//'_htc'// zchar365 znam = 'a_i'//'_htc'//TRIM(ADJUSTL(zchar)) 373 366 CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 374 367 a_i(:,:,jl) = z2d(:,:) 375 znam = 't_su'//'_htc'// zchar368 znam = 't_su'//'_htc'//TRIM(ADJUSTL(zchar)) 376 369 CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 377 370 t_su(:,:,jl) = z2d(:,:) 378 END DO 379 380 DO jl = 1, jpl 381 WRITE(zchar,'(I1)') jl 382 znam = 'tempt_sl1'//'_htc'//zchar 371 znam = 'tempt_sl1'//'_htc'//TRIM(ADJUSTL(zchar)) 383 372 CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 384 373 e_s(:,:,1,jl) = z2d(:,:) 385 END DO386 387 DO jl = 1, jpl388 WRITE(zchar,'(I1)') jl389 374 DO jk = 1, nlay_i 390 WRITE(zchar1,'(I 1)') jk391 znam = 'tempt'//'_il'// zchar1//'_htc'//zchar375 WRITE(zchar1,'(I2)') jk 376 znam = 'tempt'//'_il'//TRIM(ADJUSTL(zchar1))//'_htc'//TRIM(ADJUSTL(zchar)) 392 377 CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 393 378 e_i(:,:,jk,jl) = z2d(:,:) … … 404 389 405 390 DO jl = 1, jpl 406 WRITE(zchar,'(I 1)') jl407 znam = 'sxice'//'_htc'// zchar391 WRITE(zchar,'(I2)') jl 392 znam = 'sxice'//'_htc'//TRIM(ADJUSTL(zchar)) 408 393 CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 409 394 sxice(:,:,jl) = z2d(:,:) 410 znam = 'syice'//'_htc'// zchar395 znam = 'syice'//'_htc'//TRIM(ADJUSTL(zchar)) 411 396 CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 412 397 syice(:,:,jl) = z2d(:,:) 413 znam = 'sxxice'//'_htc'// zchar398 znam = 'sxxice'//'_htc'//TRIM(ADJUSTL(zchar)) 414 399 CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 415 400 sxxice(:,:,jl) = z2d(:,:) 416 znam = 'syyice'//'_htc'// zchar401 znam = 'syyice'//'_htc'//TRIM(ADJUSTL(zchar)) 417 402 CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 418 403 syyice(:,:,jl) = z2d(:,:) 419 znam = 'sxyice'//'_htc'// zchar404 znam = 'sxyice'//'_htc'//TRIM(ADJUSTL(zchar)) 420 405 CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 421 406 sxyice(:,:,jl) = z2d(:,:) 422 znam = 'sxsn'//'_htc'// zchar407 znam = 'sxsn'//'_htc'//TRIM(ADJUSTL(zchar)) 423 408 CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 424 409 sxsn(:,:,jl) = z2d(:,:) 425 znam = 'sysn'//'_htc'// zchar410 znam = 'sysn'//'_htc'//TRIM(ADJUSTL(zchar)) 426 411 CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 427 412 sysn(:,:,jl) = z2d(:,:) 428 znam = 'sxxsn'//'_htc'// zchar413 znam = 'sxxsn'//'_htc'//TRIM(ADJUSTL(zchar)) 429 414 CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 430 415 sxxsn(:,:,jl) = z2d(:,:) 431 znam = 'syysn'//'_htc'// zchar416 znam = 'syysn'//'_htc'//TRIM(ADJUSTL(zchar)) 432 417 CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 433 418 syysn(:,:,jl) = z2d(:,:) 434 znam = 'sxysn'//'_htc'// zchar419 znam = 'sxysn'//'_htc'//TRIM(ADJUSTL(zchar)) 435 420 CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 436 421 sxysn(:,:,jl) = z2d(:,:) 437 znam = 'sxa'//'_htc'// zchar422 znam = 'sxa'//'_htc'//TRIM(ADJUSTL(zchar)) 438 423 CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 439 424 sxa(:,:,jl) = z2d(:,:) 440 znam = 'sya'//'_htc'// zchar425 znam = 'sya'//'_htc'//TRIM(ADJUSTL(zchar)) 441 426 CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 442 427 sya(:,:,jl) = z2d(:,:) 443 znam = 'sxxa'//'_htc'// zchar428 znam = 'sxxa'//'_htc'//TRIM(ADJUSTL(zchar)) 444 429 CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 445 430 sxxa(:,:,jl) = z2d(:,:) 446 znam = 'syya'//'_htc'// zchar431 znam = 'syya'//'_htc'//TRIM(ADJUSTL(zchar)) 447 432 CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 448 433 syya(:,:,jl) = z2d(:,:) 449 znam = 'sxya'//'_htc'// zchar434 znam = 'sxya'//'_htc'//TRIM(ADJUSTL(zchar)) 450 435 CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 451 436 sxya(:,:,jl) = z2d(:,:) 452 znam = 'sxc0'//'_htc'// zchar437 znam = 'sxc0'//'_htc'//TRIM(ADJUSTL(zchar)) 453 438 CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 454 439 sxc0(:,:,jl) = z2d(:,:) 455 znam = 'syc0'//'_htc'// zchar440 znam = 'syc0'//'_htc'//TRIM(ADJUSTL(zchar)) 456 441 CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 457 442 syc0(:,:,jl) = z2d(:,:) 458 znam = 'sxxc0'//'_htc'// zchar443 znam = 'sxxc0'//'_htc'//TRIM(ADJUSTL(zchar)) 459 444 CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 460 445 sxxc0(:,:,jl) = z2d(:,:) 461 znam = 'syyc0'//'_htc'// zchar446 znam = 'syyc0'//'_htc'//TRIM(ADJUSTL(zchar)) 462 447 CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 463 448 syyc0(:,:,jl) = z2d(:,:) 464 znam = 'sxyc0'//'_htc'// zchar449 znam = 'sxyc0'//'_htc'//TRIM(ADJUSTL(zchar)) 465 450 CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 466 451 sxyc0(:,:,jl) = z2d(:,:) 467 znam = 'sxsal'//'_htc'// zchar452 znam = 'sxsal'//'_htc'//TRIM(ADJUSTL(zchar)) 468 453 CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 469 454 sxsal(:,:,jl) = z2d(:,:) 470 znam = 'sysal'//'_htc'// zchar455 znam = 'sysal'//'_htc'//TRIM(ADJUSTL(zchar)) 471 456 CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 472 457 sysal(:,:,jl) = z2d(:,:) 473 znam = 'sxxsal'//'_htc'// zchar458 znam = 'sxxsal'//'_htc'//TRIM(ADJUSTL(zchar)) 474 459 CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 475 460 sxxsal(:,:,jl) = z2d(:,:) 476 znam = 'syysal'//'_htc'// zchar461 znam = 'syysal'//'_htc'//TRIM(ADJUSTL(zchar)) 477 462 CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 478 463 syysal(:,:,jl) = z2d(:,:) 479 znam = 'sxysal'//'_htc'// zchar464 znam = 'sxysal'//'_htc'//TRIM(ADJUSTL(zchar)) 480 465 CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 481 466 sxysal(:,:,jl) = z2d(:,:) 482 znam = 'sxage'//'_htc'// zchar467 znam = 'sxage'//'_htc'//TRIM(ADJUSTL(zchar)) 483 468 CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 484 469 sxage(:,:,jl) = z2d(:,:) 485 znam = 'syage'//'_htc'// zchar470 znam = 'syage'//'_htc'//TRIM(ADJUSTL(zchar)) 486 471 CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 487 472 syage(:,:,jl) = z2d(:,:) 488 znam = 'sxxage'//'_htc'// zchar473 znam = 'sxxage'//'_htc'//TRIM(ADJUSTL(zchar)) 489 474 CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 490 475 sxxage(:,:,jl) = z2d(:,:) 491 znam = 'syyage'//'_htc'// zchar476 znam = 'syyage'//'_htc'//TRIM(ADJUSTL(zchar)) 492 477 CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 493 478 syyage(:,:,jl) = z2d(:,:) 494 znam = 'sxyage'//'_htc'// zchar479 znam = 'sxyage'//'_htc'//TRIM(ADJUSTL(zchar)) 495 480 CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 496 481 sxyage(:,:,jl)= z2d(:,:) … … 504 489 505 490 DO jl = 1, jpl 506 WRITE(zchar,'(I 1)') jl491 WRITE(zchar,'(I2)') jl 507 492 DO jk = 1, nlay_i 508 WRITE(zchar1,'(I 1)') jk509 znam = 'sxe'//'_il'// zchar1//'_htc'//zchar493 WRITE(zchar1,'(I2)') jk 494 znam = 'sxe'//'_il'//TRIM(ADJUSTL(zchar1))//'_htc'//TRIM(ADJUSTL(zchar)) 510 495 CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 511 496 sxe(:,:,jk,jl) = z2d(:,:) 512 znam = 'sye'//'_il'// zchar1//'_htc'//zchar497 znam = 'sye'//'_il'//TRIM(ADJUSTL(zchar1))//'_htc'//TRIM(ADJUSTL(zchar)) 513 498 CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 514 499 sye(:,:,jk,jl) = z2d(:,:) 515 znam = 'sxxe'//'_il'// zchar1//'_htc'//zchar500 znam = 'sxxe'//'_il'//TRIM(ADJUSTL(zchar1))//'_htc'//TRIM(ADJUSTL(zchar)) 516 501 CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 517 502 sxxe(:,:,jk,jl) = z2d(:,:) 518 znam = 'syye'//'_il'// zchar1//'_htc'//zchar503 znam = 'syye'//'_il'//TRIM(ADJUSTL(zchar1))//'_htc'//TRIM(ADJUSTL(zchar)) 519 504 CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 520 505 syye(:,:,jk,jl) = z2d(:,:) 521 znam = 'sxye'//'_il'// zchar1//'_htc'//zchar506 znam = 'sxye'//'_il'//TRIM(ADJUSTL(zchar1))//'_htc'//TRIM(ADJUSTL(zchar)) 522 507 CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 523 508 sxye(:,:,jk,jl) = z2d(:,:) -
branches/UKMO/dev_r5518_GO6_package_lim3/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90
r6498 r8540 60 60 !!---------------------------------------------------------------------- 61 61 !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) 62 !! $Id $62 !! $Id: limsbc.F90 6963 2016-09-30 12:40:04Z clem $ 63 63 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 64 64 !!---------------------------------------------------------------------- … … 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) 114 IF( iom_use('qsr_oce') ) CALL iom_put( "qsr_oce" , qsr_oce(:,:) * pfrld(:,:) ) ! solar flux at ocean surface 115 IF( iom_use('qns_oce') ) CALL iom_put( "qns_oce" , qns_oce(:,:) * pfrld(:,:) + qemp_oce(:,:) ) ! non-solar flux at ocean surface 116 IF( iom_use('qsr_ice') ) CALL iom_put( "qsr_ice" , SUM( qsr_ice(:,:,:) * a_i_b(:,:,:), dim=3 ) ) ! solar flux at ice surface 117 IF( iom_use('qns_ice') ) CALL iom_put( "qns_ice" , SUM( qns_ice(:,:,:) * a_i_b(:,:,:), dim=3 ) + qemp_ice(:,:) ) ! non-solar flux at ice surface 118 IF( iom_use('qtr_ice') ) CALL iom_put( "qtr_ice" , SUM( ftr_ice(:,:,:) * a_i_b(:,:,:), dim=3 ) ) ! solar flux transmitted thru ice 119 IF( iom_use('qt_oce' ) ) CALL iom_put( "qt_oce" , ( qsr_oce(:,:) + qns_oce(:,:) ) * pfrld(:,:) + qemp_oce(:,:) ) 120 IF( iom_use('qt_ice' ) ) CALL iom_put( "qt_ice" , SUM( ( qns_ice(:,:,:) + qsr_ice(:,:,:) ) & 121 & * a_i_b(:,:,:), dim=3 ) + qemp_ice(:,:) ) 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 112 ! make call for albedo output before it is modified 128 113 CALL wrk_alloc( jpi,jpj, zalb ) 129 114 -
branches/UKMO/dev_r5518_GO6_package_lim3/NEMOGCM/NEMO/LIM_SRC_3/limtab.F90
r6486 r8540 21 21 !!---------------------------------------------------------------------- 22 22 !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2010) 23 !! $Id $23 !! $Id: limtab.F90 4161 2013-11-07 10:01:27Z cetlod $ 24 24 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 25 25 !!---------------------------------------------------------------------- -
branches/UKMO/dev_r5518_GO6_package_lim3/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90
r6498 r8540 56 56 !!---------------------------------------------------------------------- 57 57 !! NEMO/LIM3 3.3 , UCL - NEMO Consortium (2010) 58 !! $Id $58 !! $Id: limthd.F90 7607 2017-01-25 15:37:31Z cetlod $ 59 59 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 60 60 !!---------------------------------------------------------------------- … … 613 613 CALL tab_1d_2d( nbpb, qns_ice(:,:,jl), npb, qns_ice_1d(1:nbpb) , jpi, jpj) 614 614 CALL tab_1d_2d( nbpb, ftr_ice(:,:,jl), npb, ftr_ice_1d(1:nbpb) , jpi, jpj ) 615 !616 615 END SELECT 617 616 … … 633 632 INTEGER :: ios ! Local integer output status for namelist read 634 633 NAMELIST/namicethd/ rn_hnewice, ln_frazil, rn_maxfrazb, rn_vfrazb, rn_Cfrazb, & 635 & rn_himin, rn_betas, rn_kappa_i, nn_conv_dif, rn_terr_dif, nn_ice_thcon, &636 & nn_monocat, ln_it_qnsice634 & rn_himin, rn_betas, rn_kappa_i, nn_conv_dif, rn_terr_dif, nn_ice_thcon, & 635 & rn_cdsn, nn_monocat, ln_it_qnsice 637 636 !!------------------------------------------------------------------- 638 637 ! … … 673 672 WRITE(numout,*)' maximal err. on T for heat diffusion computation rn_terr_dif = ', rn_terr_dif 674 673 WRITE(numout,*)' switch for comp. of thermal conductivity in the ice nn_ice_thcon = ', nn_ice_thcon 674 WRITE(numout,*)' thermal conductivity of the snow rn_cdsn = ', rn_cdsn 675 675 WRITE(numout,*)' check heat conservation in the ice/snow con_i = ', con_i 676 676 WRITE(numout,*)' virtual ITD mono-category parameterizations (1) or not nn_monocat = ', nn_monocat -
branches/UKMO/dev_r5518_GO6_package_lim3/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90
r7993 r8540 38 38 !!---------------------------------------------------------------------- 39 39 !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2010) 40 !! $Id $40 !! $Id: limthd_dh.F90 6469 2016-04-13 15:15:13Z clem $ 41 41 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 42 42 !!---------------------------------------------------------------------- -
branches/UKMO/dev_r5518_GO6_package_lim3/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90
r6486 r8540 32 32 !!---------------------------------------------------------------------- 33 33 !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) 34 !! $Id $34 !! $Id: limthd_dif.F90 7607 2017-01-25 15:37:31Z cetlod $ 35 35 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 36 36 !!---------------------------------------------------------------------- … … 376 376 377 377 ! Effective thickness he (zhe) 378 zfac = 1._wp / ( r cdsn + zkimean )379 zratio_s = r cdsn * zfac378 zfac = 1._wp / ( rn_cdsn + zkimean ) 379 zratio_s = rn_cdsn * zfac 380 380 zratio_i = zkimean * zfac 381 381 zhe = zratio_s * ht_i_1d(ji) + zratio_i * ht_s_1d(ji) … … 400 400 DO ji = kideb, kiut 401 401 zfac = 1. / MAX( epsi10 , zh_s(ji) ) 402 zkappa_s(ji,0) = zghe(ji) * r cdsn * zfac403 zkappa_s(ji,nlay_s) = zghe(ji) * r cdsn * zfac402 zkappa_s(ji,0) = zghe(ji) * rn_cdsn * zfac 403 zkappa_s(ji,nlay_s) = zghe(ji) * rn_cdsn * zfac 404 404 END DO 405 405 406 406 DO jk = 1, nlay_s-1 407 407 DO ji = kideb , kiut 408 zkappa_s(ji,jk) = zghe(ji) * 2.0 * r cdsn / MAX( epsi10, 2.0 * zh_s(ji) )408 zkappa_s(ji,jk) = zghe(ji) * 2.0 * rn_cdsn / MAX( epsi10, 2.0 * zh_s(ji) ) 409 409 END DO 410 410 END DO … … 422 422 zkappa_i(ji,0) = zghe(ji) * ztcond_i(ji,0) * zfac 423 423 zkappa_i(ji,nlay_i) = zghe(ji) * ztcond_i(ji,nlay_i) * zfac 424 zkappa_s(ji,nlay_s) = zghe(ji) * zghe(ji) * 2.0 * r cdsn * ztcond_i(ji,0) / &425 & MAX( epsi10, ( zghe(ji) * ztcond_i(ji,0) * zh_s(ji) + zghe(ji) * r cdsn * zh_i(ji) ) )424 zkappa_s(ji,nlay_s) = zghe(ji) * zghe(ji) * 2.0 * rn_cdsn * ztcond_i(ji,0) / & 425 & MAX( epsi10, ( zghe(ji) * ztcond_i(ji,0) * zh_s(ji) + zghe(ji) * rn_cdsn * zh_i(ji) ) ) 426 426 zkappa_i(ji,0) = zkappa_s(ji,nlay_s) * isnow(ji) + zkappa_i(ji,0) * ( 1._wp - isnow(ji) ) 427 427 END DO … … 697 697 & ( isnow(ji) * t_s_1d(ji,1) + ( 1._wp - isnow(ji) ) * t_i_1d(ji,1) ) ) / zdiagbis(ji,numeqmin(ji)) 698 698 END DO 699 699 700 ! 700 701 !-------------------------------------------------------------------------- -
branches/UKMO/dev_r5518_GO6_package_lim3/NEMOGCM/NEMO/LIM_SRC_3/limthd_ent.F90
r6486 r8540 39 39 !!---------------------------------------------------------------------- 40 40 !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) 41 !! $Id $41 !! $Id: limthd_ent.F90 5134 2015-03-09 17:27:34Z clem $ 42 42 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 43 43 !!---------------------------------------------------------------------- -
branches/UKMO/dev_r5518_GO6_package_lim3/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90
r6498 r8540 40 40 !!---------------------------------------------------------------------- 41 41 !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) 42 !! $Id $42 !! $Id: limthd_lac.F90 6316 2016-02-15 13:35:37Z cetlod $ 43 43 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 44 44 !!---------------------------------------------------------------------- -
branches/UKMO/dev_r5518_GO6_package_lim3/NEMOGCM/NEMO/LIM_SRC_3/limthd_sal.F90
r7993 r8540 33 33 !!---------------------------------------------------------------------- 34 34 !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) 35 !! $Id $35 !! $Id: limthd_sal.F90 6469 2016-04-13 15:15:13Z clem $ 36 36 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 37 37 !!---------------------------------------------------------------------- -
branches/UKMO/dev_r5518_GO6_package_lim3/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90
r7993 r8540 44 44 !!---------------------------------------------------------------------- 45 45 !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) 46 !! $Id $46 !! $Id: limtrp.F90 6476 2016-04-15 12:50:29Z mcastril $ 47 47 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 48 48 !!---------------------------------------------------------------------- -
branches/UKMO/dev_r5518_GO6_package_lim3/NEMOGCM/NEMO/LIM_SRC_3/limupdate1.F90
r6498 r8540 39 39 !!---------------------------------------------------------------------- 40 40 !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) 41 !! $Id $41 !! $Id: limupdate1.F90 7814 2017-03-20 16:21:42Z clem $ 42 42 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 43 43 !!---------------------------------------------------------------------- … … 62 62 63 63 IF( kt == nit000 .AND. lwp ) THEN 64 WRITE(numout,*) ' lim_update1 ' 65 WRITE(numout,*) ' ~~~~~~~~~~~ ' 64 WRITE(numout,*)'' 65 WRITE(numout,*)' lim_update1 ' 66 WRITE(numout,*)' ~~~~~~~~~~~ ' 66 67 ENDIF 67 68 -
branches/UKMO/dev_r5518_GO6_package_lim3/NEMOGCM/NEMO/LIM_SRC_3/limupdate2.F90
r6498 r8540 41 41 !!---------------------------------------------------------------------- 42 42 !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) 43 !! $Id $43 !! $Id: limupdate2.F90 7814 2017-03-20 16:21:42Z clem $ 44 44 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 45 45 !!---------------------------------------------------------------------- … … 62 62 63 63 IF( kt == nit000 .AND. lwp ) THEN 64 WRITE(numout,*) ' lim_update2 ' 65 WRITE(numout,*) ' ~~~~~~~~~~~ ' 64 WRITE(numout,*)'' 65 WRITE(numout,*)' lim_update2 ' 66 WRITE(numout,*)' ~~~~~~~~~~~ ' 66 67 ENDIF 67 68 -
branches/UKMO/dev_r5518_GO6_package_lim3/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90
r7993 r8540 27 27 !! - et_i(jpi,jpj) !total ice thermal content 28 28 !! - smt_i(jpi,jpj) !mean ice salinity 29 !! - ot_i(jpi,jpj) !average ice age29 !! - tm_i (jpi,jpj) !mean ice temperature 30 30 !!====================================================================== 31 31 !! History : - ! 2006-01 (M. Vancoppenolle) Original code … … 54 54 PUBLIC lim_var_eqv2glo 55 55 PUBLIC lim_var_salprof 56 PUBLIC lim_var_icetm57 56 PUBLIC lim_var_bv 58 57 PUBLIC lim_var_salprof1d … … 62 61 !!---------------------------------------------------------------------- 63 62 !! NEMO/LIM3 3.5 , UCL - NEMO Consortium (2011) 64 !! $Id $63 !! $Id: limvar.F90 7814 2017-03-20 16:21:42Z clem $ 65 64 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 66 65 !!---------------------------------------------------------------------- … … 89 88 ! Compute variables 90 89 !-------------------- 91 vt_i (:,:) = 0._wp 92 vt_s (:,:) = 0._wp 93 at_i (:,:) = 0._wp 94 ato_i(:,:) = 1._wp 95 ! 96 DO jl = 1, jpl 97 DO jj = 1, jpj 98 DO ji = 1, jpi 99 ! 100 vt_i(ji,jj) = vt_i(ji,jj) + v_i(ji,jj,jl) ! ice volume 101 vt_s(ji,jj) = vt_s(ji,jj) + v_s(ji,jj,jl) ! snow volume 102 at_i(ji,jj) = at_i(ji,jj) + a_i(ji,jj,jl) ! ice concentration 103 ! 104 rswitch = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi10 ) ) 105 icethi(ji,jj) = vt_i(ji,jj) / MAX( at_i(ji,jj) , epsi10 ) * rswitch ! ice thickness 106 END DO 107 END DO 108 END DO 109 90 ! integrated values 91 vt_i (:,:) = SUM( v_i, dim=3 ) 92 vt_s (:,:) = SUM( v_s, dim=3 ) 93 at_i (:,:) = SUM( a_i, dim=3 ) 94 et_s(:,:) = SUM( SUM( e_s(:,:,:,:), dim=4 ), dim=3 ) 95 et_i(:,:) = SUM( SUM( e_i(:,:,:,:), dim=4 ), dim=3 ) 96 ! 110 97 DO jj = 1, jpj 111 98 DO ji = 1, jpi … … 115 102 116 103 IF( kn > 1 ) THEN 117 et_s (:,:) = 0._wp 118 ot_i (:,:) = 0._wp 104 ! 105 ! mean ice/snow thickness 106 DO jj = 1, jpj 107 DO ji = 1, jpi 108 rswitch = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi10 ) ) 109 htm_i(ji,jj) = vt_i(ji,jj) / MAX( at_i(ji,jj) , epsi10 ) * rswitch 110 htm_s(ji,jj) = vt_s(ji,jj) / MAX( at_i(ji,jj) , epsi10 ) * rswitch 111 ENDDO 112 ENDDO 113 114 ! mean temperature (K), salinity and age 119 115 smt_i(:,:) = 0._wp 120 et_i (:,:) = 0._wp 121 ! 116 tm_i(:,:) = 0._wp 117 tm_su(:,:) = 0._wp 118 om_i (:,:) = 0._wp 122 119 DO jl = 1, jpl 120 123 121 DO jj = 1, jpj 124 122 DO ji = 1, jpi 125 et_s(ji,jj) = et_s(ji,jj) + e_s(ji,jj,1,jl) ! snow heat content 126 rswitch = MAX( 0._wp , SIGN( 1._wp , vt_i(ji,jj) - epsi20 ) ) 127 smt_i(ji,jj) = smt_i(ji,jj) + smv_i(ji,jj,jl) / MAX( vt_i(ji,jj) , epsi20 ) * rswitch ! ice salinity 128 rswitch = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi20 ) ) 129 ot_i(ji,jj) = ot_i(ji,jj) + oa_i(ji,jj,jl) / MAX( at_i(ji,jj) , epsi20 ) * rswitch ! ice age 130 END DO 131 END DO 132 END DO 133 ! 134 DO jl = 1, jpl 123 rswitch = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi10 ) ) 124 tm_su(ji,jj) = tm_su(ji,jj) + rswitch * ( t_su(ji,jj,jl) - rt0 ) * a_i(ji,jj,jl) / MAX( at_i(ji,jj) , epsi10 ) 125 om_i (ji,jj) = om_i (ji,jj) + rswitch * oa_i(ji,jj,jl) / MAX( at_i(ji,jj) , epsi10 ) 126 END DO 127 END DO 128 135 129 DO jk = 1, nlay_i 136 et_i(:,:) = et_i(:,:) + e_i(:,:,jk,jl) ! ice heat content 137 END DO 138 END DO 130 DO jj = 1, jpj 131 DO ji = 1, jpi 132 rswitch = MAX( 0._wp , SIGN( 1._wp , vt_i(ji,jj) - epsi10 ) ) 133 tm_i(ji,jj) = tm_i(ji,jj) + r1_nlay_i * rswitch * ( t_i(ji,jj,jk,jl) - rt0 ) * v_i(ji,jj,jl) & 134 & / MAX( vt_i(ji,jj) , epsi10 ) 135 smt_i(ji,jj) = smt_i(ji,jj) + r1_nlay_i * rswitch * s_i(ji,jj,jk,jl) * v_i(ji,jj,jl) & 136 & / MAX( vt_i(ji,jj) , epsi10 ) 137 END DO 138 END DO 139 END DO 140 END DO 141 tm_i = tm_i + rt0 142 tm_su = tm_su + rt0 139 143 ! 140 144 ENDIF … … 246 250 ! Mean temperature 247 251 !------------------- 248 vt_i (:,:) = 0._wp249 DO jl = 1, jpl250 vt_i(:,:) = vt_i(:,:) + v_i(:,:,jl)251 END DO252 ! integrated values 253 vt_i (:,:) = SUM( v_i, dim=3 ) 254 vt_s (:,:) = SUM( v_s, dim=3 ) 255 at_i (:,:) = SUM( a_i, dim=3 ) 252 256 253 257 tm_i(:,:) = 0._wp … … 398 402 399 403 400 SUBROUTINE lim_var_icetm 401 !!------------------------------------------------------------------ 402 !! *** ROUTINE lim_var_icetm *** 403 !! 404 !! ** Purpose : computes mean sea ice temperature 404 SUBROUTINE lim_var_bv 405 !!------------------------------------------------------------------ 406 !! *** ROUTINE lim_var_bv *** 407 !! 408 !! ** Purpose : computes mean brine volume (%) in sea ice 409 !! 410 !! ** Method : e = - 0.054 * S (ppt) / T (C) 411 !! 412 !! References : Vancoppenolle et al., JGR, 2007 405 413 !!------------------------------------------------------------------ 406 414 INTEGER :: ji, jj, jk, jl ! dummy loop indices 407 415 !!------------------------------------------------------------------ 408 409 ! Mean sea ice temperature 410 vt_i (:,:) = 0._wp 411 DO jl = 1, jpl 412 vt_i(:,:) = vt_i(:,:) + v_i(:,:,jl) 413 END DO 414 415 tm_i(:,:) = 0._wp 416 ! 417 bvm_i(:,:) = 0._wp 418 bv_i (:,:,:) = 0._wp 416 419 DO jl = 1, jpl 417 420 DO jk = 1, nlay_i 418 421 DO jj = 1, jpj 419 422 DO ji = 1, jpi 420 rswitch = MAX( 0._wp , SIGN( 1._wp , vt_i(ji,jj) - epsi10 ) ) 421 tm_i(ji,jj) = tm_i(ji,jj) + r1_nlay_i * rswitch * ( t_i(ji,jj,jk,jl) - rt0 ) * v_i(ji,jj,jl) & 422 & / MAX( vt_i(ji,jj) , epsi10 ) 423 END DO 424 END DO 425 END DO 426 END DO 427 tm_i = tm_i + rt0 428 429 END SUBROUTINE lim_var_icetm 430 431 432 SUBROUTINE lim_var_bv 433 !!------------------------------------------------------------------ 434 !! *** ROUTINE lim_var_bv *** 435 !! 436 !! ** Purpose : computes mean brine volume (%) in sea ice 437 !! 438 !! ** Method : e = - 0.054 * S (ppt) / T (C) 439 !! 440 !! References : Vancoppenolle et al., JGR, 2007 441 !!------------------------------------------------------------------ 442 INTEGER :: ji, jj, jk, jl ! dummy loop indices 443 REAL(wp) :: zbvi ! local scalars 444 !!------------------------------------------------------------------ 445 ! 446 vt_i (:,:) = 0._wp 447 DO jl = 1, jpl 448 vt_i(:,:) = vt_i(:,:) + v_i(:,:,jl) 449 END DO 450 451 bv_i(:,:) = 0._wp 452 DO jl = 1, jpl 453 DO jk = 1, nlay_i 454 DO jj = 1, jpj 455 DO ji = 1, jpi 456 rswitch = ( 1._wp - MAX( 0._wp , SIGN( 1._wp , (t_i(ji,jj,jk,jl) - rt0) + epsi10 ) ) ) 457 zbvi = - rswitch * tmut * s_i(ji,jj,jk,jl) / MIN( t_i(ji,jj,jk,jl) - rt0, - epsi10 ) & 458 & * v_i(ji,jj,jl) * r1_nlay_i 459 rswitch = ( 1._wp - MAX( 0._wp , SIGN( 1._wp , - vt_i(ji,jj) + epsi20 ) ) ) 460 bv_i(ji,jj) = bv_i(ji,jj) + rswitch * zbvi / MAX( vt_i(ji,jj) , epsi20 ) 461 END DO 423 rswitch = ( 1._wp - MAX( 0._wp , SIGN( 1._wp , (t_i(ji,jj,jk,jl) - rt0) + epsi10 ) ) ) 424 bv_i(ji,jj,jl) = bv_i(ji,jj,jl) - rswitch * tmut * s_i(ji,jj,jk,jl) * r1_nlay_i & 425 & / MIN( t_i(ji,jj,jk,jl) - rt0, - epsi10 ) 426 END DO 427 END DO 428 END DO 429 430 DO jj = 1, jpj 431 DO ji = 1, jpi 432 rswitch = MAX( 0._wp , SIGN( 1._wp , vt_i(ji,jj) - epsi10 ) ) 433 bvm_i(ji,jj) = bvm_i(ji,jj) + rswitch * bv_i(ji,jj,jl) * v_i(ji,jj,jl) / MAX( vt_i(ji,jj), epsi10 ) 462 434 END DO 463 435 END DO … … 683 655 INTEGER :: ji, jk, jl ! dummy loop indices 684 656 INTEGER :: ijpij, i_fill, jl0 685 REAL(wp) :: zarg, zV, zconv, zdh 657 REAL(wp) :: zarg, zV, zconv, zdh, zdv 686 658 REAL(wp), DIMENSION(:), INTENT(in) :: zhti, zhts, zai ! input ice/snow variables 687 659 REAL(wp), DIMENSION(:,:), INTENT(inout) :: zht_i, zht_s, za_i ! output ice/snow variables … … 704 676 IF( zhti(ji) > 0._wp ) THEN 705 677 706 ! initialisation of tests 707 itest(:) = 0 678 ! find which category (jl0) the input ice thickness falls into 679 jl0 = jpl 680 DO jl = 1, jpl 681 IF ( ( zhti(ji) >= hi_max(jl-1) ) .AND. ( zhti(ji) < hi_max(jl) ) ) THEN 682 jl0 = jl 683 CYCLE 684 ENDIF 685 END DO 686 687 ! initialisation of tests 688 itest(:) = 0 708 689 709 i_fill = jpl + 1 !==================================== 710 DO WHILE ( ( SUM( itest(:) ) /= 4 ) .AND. ( i_fill >= 2 ) ) ! iterative loop on i_fill categories 711 ! iteration !==================================== 712 i_fill = i_fill - 1 690 i_fill = jpl + 1 !==================================== 691 DO WHILE ( ( SUM( itest(:) ) /= 4 ) .AND. ( i_fill >= 2 ) ) ! iterative loop on i_fill categories 692 ! iteration !==================================== 693 i_fill = i_fill - 1 694 695 ! initialisation of ice variables for each try 696 zht_i(ji,1:jpl) = 0._wp 697 za_i (ji,1:jpl) = 0._wp 698 itest(:) = 0 699 700 ! *** case very thin ice: fill only category 1 701 IF ( i_fill == 1 ) THEN 702 zht_i(ji,1) = zhti(ji) 703 za_i (ji,1) = zai (ji) 704 705 ! *** case ice is thicker: fill categories >1 706 ELSE 707 708 ! Fill ice thicknesses in the (i_fill-1) cat by hmean 709 DO jl = 1, i_fill - 1 710 zht_i(ji,jl) = hi_mean(jl) 711 END DO 712 713 ! Concentrations in the (i_fill-1) categories 714 za_i(ji,jl0) = zai(ji) / SQRT(REAL(jpl)) 715 DO jl = 1, i_fill - 1 716 IF ( jl /= jl0 ) THEN 717 zarg = ( zht_i(ji,jl) - zhti(ji) ) / ( zhti(ji) * 0.5_wp ) 718 za_i(ji,jl) = za_i (ji,jl0) * EXP(-zarg**2) 719 ENDIF 720 END DO 721 722 ! Concentration in the last (i_fill) category 723 za_i(ji,i_fill) = zai(ji) - SUM( za_i(ji,1:i_fill-1) ) 724 725 ! Ice thickness in the last (i_fill) category 726 zV = SUM( za_i(ji,1:i_fill-1) * zht_i(ji,1:i_fill-1) ) 727 zht_i(ji,i_fill) = ( zhti(ji) * zai(ji) - zV ) / MAX( za_i(ji,i_fill), epsi10 ) 728 729 ! clem: correction if concentration of upper cat is greater than lower cat 730 ! (it should be a gaussian around jl0 but sometimes it is not) 731 IF ( jl0 /= jpl ) THEN 732 DO jl = jpl, jl0+1, -1 733 IF ( za_i(ji,jl) > za_i(ji,jl-1) ) THEN 734 zdv = zht_i(ji,jl) * za_i(ji,jl) 735 zht_i(ji,jl ) = 0._wp 736 za_i (ji,jl ) = 0._wp 737 za_i (ji,1:jl-1) = za_i(ji,1:jl-1) + zdv / MAX( REAL(jl-1) * zhti(ji), epsi10 ) 738 END IF 739 ENDDO 740 ENDIF 741 742 ENDIF ! case ice is thick or thin 713 743 714 ! initialisation of ice variables for each try 715 zht_i(ji,1:jpl) = 0._wp 716 za_i (ji,1:jpl) = 0._wp 744 !--------------------- 745 ! Compatibility tests 746 !--------------------- 747 ! Test 1: area conservation 748 zconv = ABS( zai(ji) - SUM( za_i(ji,1:jpl) ) ) 749 IF ( zconv < epsi06 ) itest(1) = 1 717 750 718 ! *** case very thin ice: fill only category 1 719 IF ( i_fill == 1 ) THEN 720 zht_i(ji,1) = zhti(ji) 721 za_i (ji,1) = zai (ji) 722 723 ! *** case ice is thicker: fill categories >1 724 ELSE 725 726 ! Fill ice thicknesses except the last one (i_fill) by hmean 727 DO jl = 1, i_fill - 1 728 zht_i(ji,jl) = hi_mean(jl) 729 END DO 751 ! Test 2: volume conservation 752 zconv = ABS( zhti(ji)*zai(ji) - SUM( za_i(ji,1:jpl)*zht_i(ji,1:jpl) ) ) 753 IF ( zconv < epsi06 ) itest(2) = 1 730 754 731 ! find which category (jl0) the input ice thickness falls into 732 jl0 = i_fill 755 ! Test 3: thickness of the last category is in-bounds ? 756 IF ( zht_i(ji,i_fill) >= hi_max(i_fill-1) ) itest(3) = 1 757 758 ! Test 4: positivity of ice concentrations 759 itest(4) = 1 733 760 DO jl = 1, i_fill 734 IF ( ( zhti(ji) >= hi_max(jl-1) ) .AND. ( zhti(ji) < hi_max(jl) ) ) THEN 735 jl0 = jl 736 CYCLE 737 ENDIF 738 END DO 739 740 ! Concentrations in the (i_fill-1) categories 741 za_i(ji,jl0) = zai(ji) / SQRT(REAL(jpl)) 742 DO jl = 1, i_fill - 1 743 IF ( jl == jl0 ) CYCLE 744 zarg = ( zht_i(ji,jl) - zhti(ji) ) / ( zhti(ji) * 0.5_wp ) 745 za_i(ji,jl) = za_i (ji,jl0) * EXP(-zarg**2) 746 END DO 747 748 ! Concentration in the last (i_fill) category 749 za_i(ji,i_fill) = zai(ji) - SUM( za_i(ji,1:i_fill-1) ) 750 751 ! Ice thickness in the last (i_fill) category 752 zV = SUM( za_i(ji,1:i_fill-1) * zht_i(ji,1:i_fill-1) ) 753 zht_i(ji,i_fill) = ( zhti(ji) * zai(ji) - zV ) / za_i(ji,i_fill) 754 755 ENDIF ! case ice is thick or thin 756 757 !--------------------- 758 ! Compatibility tests 759 !--------------------- 760 ! Test 1: area conservation 761 zconv = ABS( zai(ji) - SUM( za_i(ji,1:jpl) ) ) 762 IF ( zconv < epsi06 ) itest(1) = 1 763 764 ! Test 2: volume conservation 765 zconv = ABS( zhti(ji)*zai(ji) - SUM( za_i(ji,1:jpl)*zht_i(ji,1:jpl) ) ) 766 IF ( zconv < epsi06 ) itest(2) = 1 767 768 ! Test 3: thickness of the last category is in-bounds ? 769 IF ( zht_i(ji,i_fill) >= hi_max(i_fill-1) ) itest(3) = 1 770 771 ! Test 4: positivity of ice concentrations 772 itest(4) = 1 773 DO jl = 1, i_fill 774 IF ( za_i(ji,jl) < 0._wp ) itest(4) = 0 775 END DO 776 !============================ 777 END DO ! end iteration on categories 778 !============================ 761 IF ( za_i(ji,jl) < 0._wp ) itest(4) = 0 762 END DO 763 ! !============================ 764 END DO ! end iteration on categories 765 ! !============================ 779 766 ENDIF ! if zhti > 0 780 767 END DO ! i loop -
branches/UKMO/dev_r5518_GO6_package_lim3/NEMOGCM/NEMO/LIM_SRC_3/limwri.F90
r6498 r8540 17 17 USE sbc_oce ! Surface boundary condition: ocean fields 18 18 USE sbc_ice ! Surface boundary condition: ice fields 19 USE dom_ice20 19 USE ice 21 20 USE limvar … … 36 35 !!---------------------------------------------------------------------- 37 36 !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) 38 !! $Id $37 !! $Id: limwri.F90 6965 2016-09-30 13:34:19Z clem $ 39 38 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 40 39 !!---------------------------------------------------------------------- 41 40 CONTAINS 42 43 #if defined key_dimgout44 # include "limwri_dimg.h90"45 #else46 41 47 42 SUBROUTINE lim_wri( kindic ) … … 59 54 INTEGER :: ji, jj, jk, jl ! dummy loop indices 60 55 REAL(wp) :: z1_365 61 REAL(wp) :: z tmp62 REAL(wp), POINTER, DIMENSION(:,:,:) :: z oi, zei, zt_i, zt_s63 REAL(wp), POINTER, DIMENSION(:,:) :: z2d, z 2da, z2db, zswi ! 2D workspace56 REAL(wp) :: z2da, z2db, ztmp 57 REAL(wp), POINTER, DIMENSION(:,:,:) :: zswi2 58 REAL(wp), POINTER, DIMENSION(:,:) :: z2d, zswi ! 2D workspace 64 59 !!------------------------------------------------------------------- 65 60 66 61 IF( nn_timing == 1 ) CALL timing_start('limwri') 67 62 68 CALL wrk_alloc( jpi, jpj, jpl, z oi, zei, zt_i, zt_s)69 CALL wrk_alloc( jpi, jpj , z2d, z 2da, z2db, zswi )63 CALL wrk_alloc( jpi, jpj, jpl, zswi2 ) 64 CALL wrk_alloc( jpi, jpj , z2d, zswi ) 70 65 71 66 !----------------------------- … … 74 69 z1_365 = 1._wp / 365._wp 75 70 76 CALL lim_var_icetm ! mean sea ice temperature77 78 CALL lim_var_bv ! brine volume 79 80 DO jj = 1, jpj ! presence indicator of ice71 ! brine volume 72 CALL lim_var_bv 73 74 ! tresholds for outputs 75 DO jj = 1, jpj 81 76 DO ji = 1, jpi 82 77 zswi(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06 ) ) 83 78 END DO 84 79 END DO 85 ! 86 ! 87 ! 88 IF ( iom_use( "icethic_cea" ) ) THEN ! mean ice thickness 89 DO jj = 1, jpj 80 DO jl = 1, jpl 81 DO jj = 1, jpj 90 82 DO ji = 1, jpi 91 z 2d(ji,jj) = vt_i(ji,jj) / MAX( at_i(ji,jj), epsi06 ) * zswi(ji,jj)83 zswi2(ji,jj,jl) = MAX( 0._wp , SIGN( 1._wp , a_i(ji,jj,jl) - epsi06 ) ) 92 84 END DO 93 85 END DO 94 CALL iom_put( "icethic_cea" , z2d ) 95 ENDIF 96 97 IF ( iom_use( "snowthic_cea" ) ) THEN ! snow thickness = mean snow thickness over the cell 98 DO jj = 1, jpj 99 DO ji = 1, jpi 100 z2d(ji,jj) = vt_s(ji,jj) / MAX( at_i(ji,jj), epsi06 ) * zswi(ji,jj) 101 END DO 102 END DO 103 CALL iom_put( "snowthic_cea" , z2d ) 104 ENDIF 86 END DO 105 87 ! 88 ! fluxes 89 ! pfrld is the lead fraction at the previous time step (actually between TRP and THD) 90 IF( iom_use('qsr_oce') ) CALL iom_put( "qsr_oce" , qsr_oce(:,:) * pfrld(:,:) ) ! solar flux at ocean surface 91 IF( iom_use('qns_oce') ) CALL iom_put( "qns_oce" , qns_oce(:,:) * pfrld(:,:) + qemp_oce(:,:) ) ! non-solar flux at ocean surface 92 IF( iom_use('qsr_ice') ) CALL iom_put( "qsr_ice" , SUM( qsr_ice(:,:,:) * a_i_b(:,:,:), dim=3 ) ) ! solar flux at ice surface 93 IF( iom_use('qns_ice') ) CALL iom_put( "qns_ice" , SUM( qns_ice(:,:,:) * a_i_b(:,:,:), dim=3 ) + qemp_ice(:,:) ) ! non-solar flux at ice surface 94 IF( iom_use('qtr_ice') ) CALL iom_put( "qtr_ice" , SUM( ftr_ice(:,:,:) * a_i_b(:,:,:), dim=3 ) ) ! solar flux transmitted thru ice 95 IF( iom_use('qt_oce' ) ) CALL iom_put( "qt_oce" , ( qsr_oce(:,:) + qns_oce(:,:) ) * pfrld(:,:) + qemp_oce(:,:) ) 96 IF( iom_use('qt_ice' ) ) CALL iom_put( "qt_ice" , SUM( ( qns_ice(:,:,:) + qsr_ice(:,:,:) ) & 97 & * a_i_b(:,:,:),dim=3 ) + qemp_ice(:,:) ) 98 IF( iom_use('qemp_oce') ) CALL iom_put( "qemp_oce" , qemp_oce(:,:) ) 99 IF( iom_use('qemp_ice') ) CALL iom_put( "qemp_ice" , qemp_ice(:,:) ) 100 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) 101 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) 102 103 ! velocity 106 104 IF ( iom_use( "uice_ipa" ) .OR. iom_use( "vice_ipa" ) .OR. iom_use( "icevel" ) ) THEN 107 105 DO jj = 2 , jpjm1 108 106 DO ji = 2 , jpim1 109 z2da(ji,jj) = ( u_ice(ji,jj) * umask(ji,jj,1) + u_ice(ji-1,jj) * umask(ji-1,jj,1) ) * 0.5_wp 110 z2db(ji,jj) = ( v_ice(ji,jj) * vmask(ji,jj,1) + v_ice(ji,jj-1) * vmask(ji,jj-1,1) ) * 0.5_wp 107 z2da = ( u_ice(ji,jj) * umask(ji,jj,1) + u_ice(ji-1,jj) * umask(ji-1,jj,1) ) * 0.5_wp 108 z2db = ( v_ice(ji,jj) * vmask(ji,jj,1) + v_ice(ji,jj-1) * vmask(ji,jj-1,1) ) * 0.5_wp 109 z2d(ji,jj) = SQRT( z2da * z2da + z2db * z2db ) 111 110 END DO 112 111 END DO 113 CALL lbc_lnk( z2da, 'T', -1. ) 114 CALL lbc_lnk( z2db, 'T', -1. ) 115 CALL iom_put( "uice_ipa" , z2da ) ! ice velocity u component 116 CALL iom_put( "vice_ipa" , z2db ) ! ice velocity v component 117 DO jj = 1, jpj 118 DO ji = 1, jpi 119 z2d(ji,jj) = SQRT( z2da(ji,jj) * z2da(ji,jj) + z2db(ji,jj) * z2db(ji,jj) ) 120 END DO 121 END DO 122 CALL iom_put( "icevel" , z2d ) ! ice velocity module 112 CALL lbc_lnk( z2d, 'T', 1. ) 113 CALL iom_put( "uice_ipa" , u_ice ) ! ice velocity u component 114 CALL iom_put( "vice_ipa" , v_ice ) ! ice velocity v component 115 CALL iom_put( "icevel" , z2d ) ! ice velocity module 123 116 ENDIF 124 117 ! 125 IF ( iom_use( "miceage" ) ) THEN 126 z2d(:,:) = 0.e0 127 DO jl = 1, jpl 128 DO jj = 1, jpj 129 DO ji = 1, jpi 130 rswitch = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - 0.1 ) ) 131 z2d(ji,jj) = z2d(ji,jj) + rswitch * oa_i(ji,jj,jl) / MAX( at_i(ji,jj), 0.1 ) 132 END DO 133 END DO 134 END DO 135 CALL iom_put( "miceage" , z2d * z1_365 ) ! mean ice age 136 ENDIF 137 138 IF ( iom_use( "micet" ) ) THEN 139 DO jj = 1, jpj 140 DO ji = 1, jpi 141 z2d(ji,jj) = ( tm_i(ji,jj) - rt0 ) * zswi(ji,jj) 142 END DO 143 END DO 144 CALL iom_put( "micet" , z2d ) ! mean ice temperature 145 ENDIF 118 IF ( iom_use( "miceage" ) ) CALL iom_put( "miceage" , om_i * zswi * z1_365 ) ! mean ice age 119 IF ( iom_use( "icethic_cea" ) ) CALL iom_put( "icethic_cea" , htm_i * zswi ) ! ice thickness mean 120 IF ( iom_use( "snowthic_cea" ) ) CALL iom_put( "snowthic_cea", htm_s * zswi ) ! snow thickness mean 121 IF ( iom_use( "micet" ) ) CALL iom_put( "micet" , ( tm_i - rt0 ) * zswi ) ! ice mean temperature 122 IF ( iom_use( "icest" ) ) CALL iom_put( "icest" , ( tm_su - rt0 ) * zswi ) ! ice surface temperature 123 IF ( iom_use( "icecolf" ) ) CALL iom_put( "icecolf" , hicol ) ! frazil ice collection thickness 146 124 ! 147 IF ( iom_use( "icest" ) ) THEN148 z2d(:,:) = 0.e0149 DO jl = 1, jpl150 DO jj = 1, jpj151 DO ji = 1, jpi152 z2d(ji,jj) = z2d(ji,jj) + zswi(ji,jj) * ( t_su(ji,jj,jl) - rt0 ) * a_i(ji,jj,jl) / MAX( at_i(ji,jj) , epsi06 )153 END DO154 END DO155 END DO156 CALL iom_put( "icest" , z2d ) ! ice surface temperature157 ENDIF158 159 IF ( iom_use( "icecolf" ) ) CALL iom_put( "icecolf", hicol ) ! frazil ice collection thickness160 161 125 CALL iom_put( "isst" , sst_m ) ! sea surface temperature 162 126 CALL iom_put( "isss" , sss_m ) ! sea surface salinity 163 CALL iom_put( "iceconc" , at_i 164 CALL iom_put( "icevolu" , vt_i 165 CALL iom_put( "icehc" , et_i 166 CALL iom_put( "isnowhc" , et_s 167 CALL iom_put( "ibrinv" , bv _i * 100._wp) ! brine volume127 CALL iom_put( "iceconc" , at_i * zswi ) ! ice concentration 128 CALL iom_put( "icevolu" , vt_i * zswi ) ! ice volume = mean ice thickness over the cell 129 CALL iom_put( "icehc" , et_i * zswi ) ! ice total heat content 130 CALL iom_put( "isnowhc" , et_s * zswi ) ! snow total heat content 131 CALL iom_put( "ibrinv" , bvm_i * zswi * 100. ) ! brine volume 168 132 CALL iom_put( "utau_ice" , utau_ice ) ! wind stress over ice along i-axis at I-point 169 133 CALL iom_put( "vtau_ice" , vtau_ice ) ! wind stress over ice along j-axis at I-point 170 134 CALL iom_put( "snowpre" , sprecip * 86400. ) ! snow precipitation 171 CALL iom_put( "micesalt" , smt_i 172 173 CALL iom_put( "icestr" , strength * 0.001 )! ice strength174 CALL iom_put( "idive" , divu_i * 1.0e8 ) 175 CALL iom_put( "ishear" , shear_i * 1.0e8 ) 176 CALL iom_put( "snowvol" , vt_s 135 CALL iom_put( "micesalt" , smt_i * zswi ) ! mean ice salinity 136 137 CALL iom_put( "icestr" , strength * zswi ) ! ice strength 138 CALL iom_put( "idive" , divu_i * 1.0e8 ) ! divergence 139 CALL iom_put( "ishear" , shear_i * 1.0e8 ) ! shear 140 CALL iom_put( "snowvol" , vt_s * zswi ) ! snow volume 177 141 178 142 CALL iom_put( "icetrp" , diag_trp_vi * rday ) ! ice volume transport … … 183 147 184 148 CALL iom_put( "sfxbog" , sfx_bog * rday ) ! salt flux from bottom growth 185 CALL iom_put( "sfxbom" , sfx_bom * rday ) ! salt flux from bottom melt 186 CALL iom_put( "sfxsum" , sfx_sum * rday ) ! salt flux from surface melt 149 CALL iom_put( "sfxbom" , sfx_bom * rday ) ! salt flux from bottom melting 150 CALL iom_put( "sfxsum" , sfx_sum * rday ) ! salt flux from surface melting 187 151 CALL iom_put( "sfxsni" , sfx_sni * rday ) ! salt flux from snow ice formation 188 152 CALL iom_put( "sfxopw" , sfx_opw * rday ) ! salt flux from open water formation 189 153 CALL iom_put( "sfxdyn" , sfx_dyn * rday ) ! salt flux from ridging rafting 190 CALL iom_put( "sfxres" , sfx_res * rday ) ! salt flux from residual154 CALL iom_put( "sfxres" , sfx_res * rday ) ! salt flux from limupdate (resultant) 191 155 CALL iom_put( "sfxbri" , sfx_bri * rday ) ! salt flux from brines 192 156 CALL iom_put( "sfxsub" , sfx_sub * rday ) ! salt flux from sublimation … … 202 166 CALL iom_put( "vfxbom" , wfx_bom * ztmp ) ! bottom melt 203 167 CALL iom_put( "vfxice" , wfx_ice * ztmp ) ! total ice growth/melt 168 169 IF ( iom_use( "vfxthin" ) ) THEN ! ice production for open water + thin ice (<20cm) => comparable to observations 170 WHERE( htm_i(:,:) < 0.2 .AND. htm_i(:,:) > 0. ) ; z2d = wfx_bog 171 ELSEWHERE ; z2d = 0._wp 172 END WHERE 173 CALL iom_put( "vfxthin", ( wfx_opw + z2d ) * ztmp ) 174 ENDIF 175 176 ztmp = rday / rhosn 177 CALL iom_put( "vfxspr" , wfx_spr * ztmp ) ! precip (snow) 204 178 CALL iom_put( "vfxsnw" , wfx_snw * ztmp ) ! total snw growth/melt 205 CALL iom_put( "vfxsub" , wfx_sub * ztmp ) ! sublimation (snow )206 CALL iom_put( "vfxs pr" , wfx_spr * ztmp ) ! precip (snow)207 179 CALL iom_put( "vfxsub" , wfx_sub * ztmp ) ! sublimation (snow/ice) 180 CALL iom_put( "vfxsub_err" , wfx_err_sub * ztmp ) ! "excess" of sublimation sent to ocean 181 208 182 CALL iom_put( "afxtot" , afx_tot * rday ) ! concentration tendency (total) 209 183 CALL iom_put( "afxdyn" , afx_dyn * rday ) ! concentration tendency (dynamics) … … 225 199 CALL iom_put ('hfxdif' , hfx_dif(:,:) ) ! 226 200 CALL iom_put ('hfxopw' , hfx_opw(:,:) ) ! 227 CALL iom_put ('hfxtur' , fhtur(:,:) * SUM( a_i_b(:,:,:), dim=3) ) ! turbulent heat flux at ice base201 CALL iom_put ('hfxtur' , fhtur(:,:) * SUM( a_i_b(:,:,:), dim=3 ) ) ! turbulent heat flux at ice base 228 202 CALL iom_put ('hfxdhc' , diag_heat(:,:) ) ! Heat content variation in snow and ice 229 203 CALL iom_put ('hfxspr' , hfx_spr(:,:) ) ! Heat content of snow precip 230 204 231 232 IF ( iom_use( "vfxthin" ) ) THEN ! ice production for open water + thin ice (<20cm) => comparable to observations233 DO jj = 1, jpj234 DO ji = 1, jpi235 z2d(ji,jj) = vt_i(ji,jj) / MAX( at_i(ji,jj), epsi06 ) * zswi(ji,jj) ! mean ice thickness236 END DO237 END DO238 WHERE( z2d(:,:) < 0.2 .AND. z2d(:,:) > 0. ) ; z2da = wfx_bog239 ELSEWHERE ; z2da = 0._wp240 END WHERE241 CALL iom_put( "vfxthin", ( wfx_opw + z2da ) * ztmp )242 ENDIF243 205 244 206 !-------------------------------- 245 207 ! Output values for each category 246 208 !-------------------------------- 247 CALL iom_put( "iceconc_cat" , a_i ) ! area for categories 248 CALL iom_put( "icethic_cat" , ht_i ) ! thickness for categories 249 CALL iom_put( "snowthic_cat" , ht_s ) ! snow depth for categories 250 CALL iom_put( "salinity_cat" , sm_i ) ! salinity for categories 251 209 IF ( iom_use( "iceconc_cat" ) ) CALL iom_put( "iceconc_cat" , a_i * zswi2 ) ! area for categories 210 IF ( iom_use( "icethic_cat" ) ) CALL iom_put( "icethic_cat" , ht_i * zswi2 ) ! thickness for categories 211 IF ( iom_use( "snowthic_cat" ) ) CALL iom_put( "snowthic_cat" , ht_s * zswi2 ) ! snow depth for categories 212 IF ( iom_use( "salinity_cat" ) ) CALL iom_put( "salinity_cat" , sm_i * zswi2 ) ! salinity for categories 252 213 ! ice temperature 253 IF ( iom_use( "icetemp_cat" ) ) THEN 254 zt_i(:,:,:) = SUM( t_i(:,:,:,:), dim=3 ) * r1_nlay_i 255 CALL iom_put( "icetemp_cat" , zt_i - rt0 ) 256 ENDIF 257 214 IF ( iom_use( "icetemp_cat" ) ) CALL iom_put( "icetemp_cat", ( SUM( t_i(:,:,:,:), dim=3 ) * r1_nlay_i - rt0 ) * zswi2 ) 258 215 ! snow temperature 259 IF ( iom_use( "snwtemp_cat" ) ) THEN 260 zt_s(:,:,:) = SUM( t_s(:,:,:,:), dim=3 ) * r1_nlay_s 261 CALL iom_put( "snwtemp_cat" , zt_s - rt0 ) 262 ENDIF 263 264 ! Compute ice age 265 IF ( iom_use( "iceage_cat" ) ) THEN 266 DO jl = 1, jpl 267 DO jj = 1, jpj 268 DO ji = 1, jpi 269 rswitch = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - 0.1 ) ) 270 rswitch = rswitch * MAX( 0._wp , SIGN( 1._wp , a_i(ji,jj,jl) - 0.1 ) ) 271 zoi(ji,jj,jl) = oa_i(ji,jj,jl) / MAX( a_i(ji,jj,jl) , 0.1 ) * rswitch 272 END DO 273 END DO 274 END DO 275 CALL iom_put( "iceage_cat" , zoi * z1_365 ) ! ice age for categories 276 ENDIF 277 278 ! Compute brine volume 279 IF ( iom_use( "brinevol_cat" ) ) THEN 280 zei(:,:,:) = 0._wp 281 DO jl = 1, jpl 282 DO jk = 1, nlay_i 283 DO jj = 1, jpj 284 DO ji = 1, jpi 285 rswitch = MAX( 0._wp , SIGN( 1._wp , a_i(ji,jj,jl) - epsi06 ) ) 286 zei(ji,jj,jl) = zei(ji,jj,jl) + 100.0 * & 287 ( - tmut * s_i(ji,jj,jk,jl) / MIN( ( t_i(ji,jj,jk,jl) - rt0 ), - epsi06 ) ) * & 288 rswitch * r1_nlay_i 289 END DO 290 END DO 291 END DO 292 END DO 293 CALL iom_put( "brinevol_cat" , zei ) ! brine volume for categories 294 ENDIF 216 IF ( iom_use( "snwtemp_cat" ) ) CALL iom_put( "snwtemp_cat", ( SUM( t_s(:,:,:,:), dim=3 ) * r1_nlay_s - rt0 ) * zswi2 ) 217 ! ice age 218 IF ( iom_use( "iceage_cat" ) ) CALL iom_put( "iceage_cat" , o_i * zswi2 * z1_365 ) 219 ! brine volume 220 IF ( iom_use( "brinevol_cat" ) ) CALL iom_put( "brinevol_cat", bv_i * 100. * zswi2 ) 295 221 296 222 ! ! Create an output files (output.lim.abort.nc) if S < 0 or u > 20 m/s … … 298 224 ! not yet implemented 299 225 300 CALL wrk_dealloc( jpi, jpj, jpl, z oi, zei, zt_i, zt_s)301 CALL wrk_dealloc( jpi, jpj , z2d, zswi , z2da, z2db)226 CALL wrk_dealloc( jpi, jpj, jpl, zswi2 ) 227 CALL wrk_dealloc( jpi, jpj , z2d, zswi ) 302 228 303 229 IF( nn_timing == 1 ) CALL timing_stop('limwri') 304 230 305 231 END SUBROUTINE lim_wri 306 #endif307 232 308 233 … … 319 244 !! 4.0 ! 2013-06 (C. Rousset) 320 245 !!---------------------------------------------------------------------- 321 INTEGER, INTENT( in ) :: kt ! ocean time-step index) 322 INTEGER, INTENT( in ) :: kid , kh_i 246 INTEGER, INTENT( in ) :: kt ! ocean time-step index) 247 INTEGER, INTENT( in ) :: kid , kh_i 248 INTEGER :: nz_i, jl 249 REAL(wp), DIMENSION(jpl) :: jcat 323 250 !!---------------------------------------------------------------------- 324 325 CALL histdef( kid, "iicethic", "Ice thickness" , "m" , & 326 & jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 327 CALL histdef( kid, "iiceconc", "Ice concentration" , "%" , & 328 & jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 329 CALL histdef( kid, "iicetemp", "Ice temperature" , "C" , & 330 & jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 331 CALL histdef( kid, "iicevelu", "i-Ice speed (I-point)" , "m/s" , & 332 & jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 333 CALL histdef( kid, "iicevelv", "j-Ice speed (I-point)" , "m/s" , & 334 & jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 335 CALL histdef( kid, "iicestru", "i-Wind stress over ice (I-pt)", "Pa", & 336 & jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 337 CALL histdef( kid, "iicestrv", "j-Wind stress over ice (I-pt)", "Pa", & 338 & jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 339 CALL histdef( kid, "iicesflx", "Solar flux over ocean" , "w/m2" , & 340 & jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 341 CALL histdef( kid, "iicenflx", "Non-solar flux over ocean" , "w/m2" , & 251 DO jl = 1, jpl 252 jcat(jl) = REAL(jl) 253 ENDDO 254 255 CALL histvert( kid, "ncatice", "Ice Categories","", jpl, jcat, nz_i, "up") 256 257 CALL histdef( kid, "sithic", "Ice thickness" , "m" , & 258 & jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 259 CALL histdef( kid, "siconc", "Ice concentration" , "%" , & 260 & jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 261 CALL histdef( kid, "sitemp", "Ice temperature" , "C" , & 262 & jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 263 CALL histdef( kid, "sivelu", "i-Ice speed " , "m/s" , & 264 & jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 265 CALL histdef( kid, "sivelv", "j-Ice speed " , "m/s" , & 266 & jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 267 CALL histdef( kid, "sistru", "i-Wind stress over ice " , "Pa" , & 268 & jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 269 CALL histdef( kid, "sistrv", "j-Wind stress over ice " , "Pa" , & 270 & jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 271 CALL histdef( kid, "sisflx", "Solar flux over ocean" , "w/m2" , & 272 & jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 273 CALL histdef( kid, "sinflx", "Non-solar flux over ocean" , "w/m2" , & 342 274 & jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 343 275 CALL histdef( kid, "isnowpre", "Snow precipitation" , "kg/m2/s", & 344 276 & jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 345 CALL histdef( kid, "iicesali", "Ice salinity" , "PSU" , & 346 & jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 347 CALL histdef( kid, "iicevolu", "Ice volume" , "m" , & 348 & jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 349 CALL histdef( kid, "iicedive", "Ice divergence" , "10-8s-1", & 350 & jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 351 CALL histdef( kid, "iicebopr", "Ice bottom production" , "m/s" , & 352 & jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 353 CALL histdef( kid, "iicedypr", "Ice dynamic production" , "m/s" , & 354 & jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 355 CALL histdef( kid, "iicelapr", "Ice open water prod" , "m/s" , & 277 CALL histdef( kid, "sisali", "Ice salinity" , "PSU" , & 278 & jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 279 CALL histdef( kid, "sivolu", "Ice volume" , "m" , & 280 & jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 281 CALL histdef( kid, "sidive", "Ice divergence" , "10-8s-1", & 282 & jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 283 284 CALL histdef( kid, "vfxbog", "Ice bottom production" , "m/s" , & 285 & jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 286 CALL histdef( kid, "vfxdyn", "Ice dynamic production" , "m/s" , & 287 & jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 288 CALL histdef( kid, "vfxopw", "Ice open water prod" , "m/s" , & 356 289 & jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 357 CALL histdef( kid, "iicesipr", "Snow ice production " , "m/s" , & 358 & jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 359 CALL histdef( kid, "iicerepr", "Ice prod from limupdate" , "m/s" , & 360 & jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 361 CALL histdef( kid, "iicebome", "Ice bottom melt" , "m/s" , & 362 & jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 363 CALL histdef( kid, "iicesume", "Ice surface melt" , "m/s" , & 364 & jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 365 CALL histdef( kid, "iisfxdyn", "Salt flux from dynmics" , "" , & 366 & jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 367 CALL histdef( kid, "iisfxres", "Salt flux from limupdate", "" , & 368 & jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 290 CALL histdef( kid, "vfxsni", "Snow ice production " , "m/s" , & 291 & jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 292 CALL histdef( kid, "vfxres", "Ice prod from limupdate" , "m/s" , & 293 & jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 294 CALL histdef( kid, "vfxbom", "Ice bottom melt" , "m/s" , & 295 & jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 296 CALL histdef( kid, "vfxsum", "Ice surface melt" , "m/s" , & 297 & jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 298 299 CALL histdef( kid, "sithicat", "Ice thickness" , "m" , & 300 & jpi, jpj, kh_i, jpl, 1, jpl, nz_i, 32, "inst(x)", rdt, rdt ) 301 CALL histdef( kid, "siconcat", "Ice concentration" , "%" , & 302 & jpi, jpj, kh_i, jpl, 1, jpl, nz_i, 32, "inst(x)", rdt, rdt ) 303 CALL histdef( kid, "sisalcat", "Ice salinity" , "" , & 304 & jpi, jpj, kh_i, jpl, 1, jpl, nz_i, 32, "inst(x)", rdt, rdt ) 305 CALL histdef( kid, "sitemcat", "Ice temperature" , "C" , & 306 & jpi, jpj, kh_i, jpl, 1, jpl, nz_i, 32, "inst(x)", rdt, rdt ) 307 CALL histdef( kid, "snthicat", "Snw thickness" , "m" , & 308 & jpi, jpj, kh_i, jpl, 1, jpl, nz_i, 32, "inst(x)", rdt, rdt ) 309 CALL histdef( kid, "sntemcat", "Snw temperature" , "C" , & 310 & jpi, jpj, kh_i, jpl, 1, jpl, nz_i, 32, "inst(x)", rdt, rdt ) 369 311 370 312 CALL histend( kid, snc4set ) ! end of the file definition 371 313 372 CALL histwrite( kid, " iicethic", kt, icethi, jpi*jpj, (/1/) )373 CALL histwrite( kid, " iiceconc", kt, at_i , jpi*jpj, (/1/) )374 CALL histwrite( kid, " iicetemp", kt, tm_i - rt0 , jpi*jpj, (/1/) )375 CALL histwrite( kid, " iicevelu", kt, u_ice , jpi*jpj, (/1/) )376 CALL histwrite( kid, " iicevelv", kt, v_ice , jpi*jpj, (/1/) )377 CALL histwrite( kid, " iicestru", kt, utau_ice , jpi*jpj, (/1/) )378 CALL histwrite( kid, " iicestrv", kt, vtau_ice , jpi*jpj, (/1/) )379 CALL histwrite( kid, " iicesflx", kt, qsr , jpi*jpj, (/1/) )380 CALL histwrite( kid, " iicenflx", kt, qns , jpi*jpj, (/1/) )314 CALL histwrite( kid, "sithic", kt, htm_i , jpi*jpj, (/1/) ) 315 CALL histwrite( kid, "siconc", kt, at_i , jpi*jpj, (/1/) ) 316 CALL histwrite( kid, "sitemp", kt, tm_i - rt0 , jpi*jpj, (/1/) ) 317 CALL histwrite( kid, "sivelu", kt, u_ice , jpi*jpj, (/1/) ) 318 CALL histwrite( kid, "sivelv", kt, v_ice , jpi*jpj, (/1/) ) 319 CALL histwrite( kid, "sistru", kt, utau_ice , jpi*jpj, (/1/) ) 320 CALL histwrite( kid, "sistrv", kt, vtau_ice , jpi*jpj, (/1/) ) 321 CALL histwrite( kid, "sisflx", kt, qsr , jpi*jpj, (/1/) ) 322 CALL histwrite( kid, "sinflx", kt, qns , jpi*jpj, (/1/) ) 381 323 CALL histwrite( kid, "isnowpre", kt, sprecip , jpi*jpj, (/1/) ) 382 CALL histwrite( kid, "iicesali", kt, smt_i , jpi*jpj, (/1/) ) 383 CALL histwrite( kid, "iicevolu", kt, vt_i , jpi*jpj, (/1/) ) 384 CALL histwrite( kid, "iicedive", kt, divu_i*1.0e8 , jpi*jpj, (/1/) ) 385 386 CALL histwrite( kid, "iicebopr", kt, wfx_bog , jpi*jpj, (/1/) ) 387 CALL histwrite( kid, "iicedypr", kt, wfx_dyn , jpi*jpj, (/1/) ) 388 CALL histwrite( kid, "iicelapr", kt, wfx_opw , jpi*jpj, (/1/) ) 389 CALL histwrite( kid, "iicesipr", kt, wfx_sni , jpi*jpj, (/1/) ) 390 CALL histwrite( kid, "iicerepr", kt, wfx_res , jpi*jpj, (/1/) ) 391 CALL histwrite( kid, "iicebome", kt, wfx_bom , jpi*jpj, (/1/) ) 392 CALL histwrite( kid, "iicesume", kt, wfx_sum , jpi*jpj, (/1/) ) 393 CALL histwrite( kid, "iisfxdyn", kt, sfx_dyn , jpi*jpj, (/1/) ) 394 CALL histwrite( kid, "iisfxres", kt, sfx_res , jpi*jpj, (/1/) ) 324 CALL histwrite( kid, "sisali", kt, smt_i , jpi*jpj, (/1/) ) 325 CALL histwrite( kid, "sivolu", kt, vt_i , jpi*jpj, (/1/) ) 326 CALL histwrite( kid, "sidive", kt, divu_i*1.0e8 , jpi*jpj, (/1/) ) 327 328 CALL histwrite( kid, "vfxbog", kt, wfx_bog , jpi*jpj, (/1/) ) 329 CALL histwrite( kid, "vfxdyn", kt, wfx_dyn , jpi*jpj, (/1/) ) 330 CALL histwrite( kid, "vfxopw", kt, wfx_opw , jpi*jpj, (/1/) ) 331 CALL histwrite( kid, "vfxsni", kt, wfx_sni , jpi*jpj, (/1/) ) 332 CALL histwrite( kid, "vfxres", kt, wfx_res , jpi*jpj, (/1/) ) 333 CALL histwrite( kid, "vfxbom", kt, wfx_bom , jpi*jpj, (/1/) ) 334 CALL histwrite( kid, "vfxsum", kt, wfx_sum , jpi*jpj, (/1/) ) 335 336 CALL histwrite( kid, "sithicat", kt, ht_i , jpi*jpj*jpl, (/1/) ) 337 CALL histwrite( kid, "siconcat", kt, a_i , jpi*jpj*jpl, (/1/) ) 338 CALL histwrite( kid, "sisalcat", kt, sm_i , jpi*jpj*jpl, (/1/) ) 339 CALL histwrite( kid, "sitemcat", kt, tm_i - rt0 , jpi*jpj*jpl, (/1/) ) 340 CALL histwrite( kid, "snthicat", kt, ht_s , jpi*jpj*jpl, (/1/) ) 341 CALL histwrite( kid, "sntemcat", kt, tm_su - rt0 , jpi*jpj*jpl, (/1/) ) 395 342 396 343 ! Close the file -
branches/UKMO/dev_r5518_GO6_package_lim3/NEMOGCM/NEMO/LIM_SRC_3/limwri_dimg.h90
r6486 r8540 2 2 !!---------------------------------------------------------------------- 3 3 !! NEMO/LIM3 3.3 , UCL - NEMO Consortium (2010) 4 !! $Id $4 !! $Id: limwri_dimg.h90 5123 2015-03-04 16:06:03Z clem $ 5 5 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 6 6 !!---------------------------------------------------------------------- -
branches/UKMO/dev_r5518_GO6_package_lim3/NEMOGCM/NEMO/LIM_SRC_3/thd_ice.F90
r6498 r8540 130 130 !!---------------------------------------------------------------------- 131 131 !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) 132 !! $Id $132 !! $Id: thd_ice.F90 6399 2016-03-22 17:17:23Z clem $ 133 133 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 134 134 !!---------------------------------------------------------------------- -
branches/UKMO/dev_r5518_GO6_package_lim3/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90
r6498 r8540 78 78 !!---------------------------------------------------------------------- 79 79 !! NEMO/OPA 4.0 , UCL NEMO Consortium (2011) 80 !! $Id $80 !! $Id: sbcice_lim.F90 7778 2017-03-09 14:19:31Z clem $ 81 81 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 82 82 !!---------------------------------------------------------------------- … … 229 229 CALL lim_sbc_flx( kt ) ! Update surface ocean mass, heat and salt fluxes 230 230 ! 231 IF(ln_limdiaout) CALL lim_diahsb 231 IF(ln_limdiaout) CALL lim_diahsb( kt ) ! Diagnostics and outputs 232 232 ! 233 233 CALL lim_wri( 1 ) ! Ice outputs … … 310 310 numit = nit000 - 1 311 311 ENDIF 312 CALL lim_var_agg( 1)312 CALL lim_var_agg(2) 313 313 CALL lim_var_glo2eqv 314 314 ! 315 315 CALL lim_sbc_init ! ice surface boundary condition 316 ! 317 IF( ln_limdiaout) CALL lim_diahsb_init ! initialization for diags 316 318 ! 317 319 fr_i(:,:) = at_i(:,:) ! initialisation of sea-ice fraction … … 648 650 CONTAINS 649 651 SUBROUTINE sbc_ice_lim ( kt, kblk ) ! Dummy routine 652 INTEGER, INTENT(in) :: kt, kblk 650 653 WRITE(*,*) 'sbc_ice_lim: You should not have seen this print! error?', kt, kblk 651 654 END SUBROUTINE sbc_ice_lim
Note: See TracChangeset
for help on using the changeset viewer.