Changeset 13741
- Timestamp:
- 2020-11-06T14:50:17+01:00 (4 years ago)
- Location:
- NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling
- Files:
-
- 69 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/cfgs/AGRIF_DEMO/EXPREF/1_namelist_cfg
r13553 r13741 358 358 !! !! 359 359 !! namtrd dynamics and/or tracer trends (default: OFF) 360 !! namptr Poleward Transport Diagnostics (default: OFF)361 360 !! namhsb Heat and salt budgets (default: OFF) 362 361 !! namdiu Cool skin and warm layer models (default: OFF) -
NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/cfgs/AGRIF_DEMO/EXPREF/2_namelist_cfg
r13216 r13741 307 307 !! !! 308 308 !! namtrd dynamics and/or tracer trends (default: OFF) 309 !! namptr Poleward Transport Diagnostics (default: OFF)310 309 !! namhsb Heat and salt budgets (default: OFF) 311 310 !! namdiu Cool skin and warm layer models (default: OFF) -
NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/cfgs/AGRIF_DEMO/EXPREF/3_namelist_cfg
r13216 r13741 307 307 !! !! 308 308 !! namtrd dynamics and/or tracer trends (default: OFF) 309 !! namptr Poleward Transport Diagnostics (default: OFF)310 309 !! namhsb Heat and salt budgets (default: OFF) 311 310 !! namdiu Cool skin and warm layer models (default: OFF) -
NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/cfgs/AGRIF_DEMO/EXPREF/file_def_nemo-ice.xml
r9896 r13741 52 52 <field field_ref="normstr" name="normstr" /> 53 53 <field field_ref="sheastr" name="sheastr" /> 54 <field field_ref="isig1" name="isig1" /> 55 <field field_ref="isig2" name="isig2" /> 56 <field field_ref="isig3" name="isig3" /> 54 <field field_ref="sig1_pnorm" name="sig1_pnorm"/> 55 <field field_ref="sig2_pnorm" name="sig2_pnorm"/> 57 56 58 57 <!-- heat fluxes --> -
NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/cfgs/AGRIF_DEMO/EXPREF/namelist_cfg
r13553 r13741 359 359 !! !! 360 360 !! namtrd dynamics and/or tracer trends (default: OFF) 361 !! namptr Poleward Transport Diagnostics (default: OFF)362 361 !! namhsb Heat and salt budgets (default: OFF) 363 362 !! namdiu Cool skin and warm layer models (default: OFF) -
NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/cfgs/AMM12/EXPREF/namelist_cfg
r13553 r13741 345 345 !! !! 346 346 !! namtrd dynamics and/or tracer trends (default: OFF) 347 !! namptr Poleward Transport Diagnostics (default: OFF)348 347 !! namhsb Heat and salt budgets (default: OFF) 349 348 !! namdiu Cool skin and warm layer models (default: OFF) -
NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/cfgs/C1D_PAPA/EXPREF/namelist_cfg
r13553 r13741 414 414 !! !! 415 415 !! namtrd dynamics and/or tracer trends (default: OFF) 416 !! namptr Poleward Transport Diagnostics (default: OFF)417 416 !! namhsb Heat and salt budgets (default: OFF) 418 417 !! namdiu Cool skin and warm layer models (default: OFF) … … 429 428 / 430 429 !----------------------------------------------------------------------- 431 &namptr ! Poleward Transport Diagnostic (default: OFF)432 430 !----------------------------------------------------------------------- 433 431 / -
NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/cfgs/GYRE_BFM/EXPREF/namelist_cfg
r13553 r13741 223 223 !! !! 224 224 !! namtrd dynamics and/or tracer trends (default: OFF) 225 !! namptr Poleward Transport Diagnostics (default: OFF)226 225 !! namhsb Heat and salt budgets (default: OFF) 227 226 !! namdiu Cool skin and warm layer models (default: OFF) -
NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/cfgs/GYRE_PISCES/EXPREF/namelist_cfg
r13553 r13741 217 217 !! !! 218 218 !! namtrd dynamics and/or tracer trends (default: OFF) 219 !! namptr Poleward Transport Diagnostics (default: OFF)220 219 !! namhsb Heat and salt budgets (default: OFF) 221 220 !! namdiu Cool skin and warm layer models (default: OFF) -
NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/cfgs/ORCA2_ICE_ABL/EXPREF/file_def_nemo-ice.xml
r11945 r13741 53 53 <field field_ref="normstr" name="normstr" /> 54 54 <field field_ref="sheastr" name="sheastr" /> 55 <field field_ref="isig1" name="isig1" /> 56 <field field_ref="isig2" name="isig2" /> 57 <field field_ref="isig3" name="isig3" /> 55 <field field_ref="sig1_pnorm" name="sig1_pnorm"/> 56 <field field_ref="sig2_pnorm" name="sig2_pnorm"/> 58 57 59 58 <!-- heat fluxes --> -
NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/cfgs/ORCA2_ICE_ABL/EXPREF/namelist_cfg
r13553 r13741 401 401 !! !! 402 402 !! namtrd dynamics and/or tracer trends (default: OFF) 403 !! namptr Poleward Transport Diagnostics (default: OFF)404 403 !! namhsb Heat and salt budgets (default: OFF) 405 404 !! namdiu Cool skin and warm layer models (default: OFF) -
NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/cfgs/ORCA2_ICE_PISCES/EXPREF/file_def_nemo-ice.xml
r12377 r13741 53 53 <field field_ref="normstr" name="normstr" /> 54 54 <field field_ref="sheastr" name="sheastr" /> 55 <field field_ref="isig1" name="isig1" /> 56 <field field_ref="isig2" name="isig2" /> 57 <field field_ref="isig3" name="isig3" /> 55 <field field_ref="sig1_pnorm" name="sig1_pnorm"/> 56 <field field_ref="sig2_pnorm" name="sig2_pnorm"/> 58 57 59 58 <!-- heat fluxes --> -
NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/cfgs/ORCA2_ICE_PISCES/EXPREF/namelist_cfg
r13553 r13741 400 400 !! !! 401 401 !! namtrd dynamics and/or tracer trends (default: OFF) 402 !! namptr Poleward Transport Diagnostics (default: OFF)403 402 !! namhsb Heat and salt budgets (default: OFF) 404 403 !! namdiu Cool skin and warm layer models (default: OFF) -
NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/cfgs/ORCA2_OFF_PISCES/EXPREF/namelist_cfg
r13553 r13741 199 199 &namdrg ! top/bottom drag coefficient (default: NO selection) 200 200 !----------------------------------------------------------------------- 201 ln_drg_OFF = .true. ! must select one drag coefficient to true even if we don't use it 201 202 / 202 203 !----------------------------------------------------------------------- … … 347 348 &namzdf ! vertical physics manager (default: NO selection) 348 349 !----------------------------------------------------------------------- 350 ln_zdfcst = .true. ! must select one vertical physics to true even if we don't use it 349 351 / 350 352 !----------------------------------------------------------------------- … … 372 374 !! !! 373 375 !! namtrd dynamics and/or tracer trends (default: OFF) 374 !! namptr Poleward Transport Diagnostics (default: OFF)375 376 !! namhsb Heat and salt budgets (default: OFF) 376 377 !! namdiu Cool skin and warm layer models (default: OFF) … … 387 388 / 388 389 !----------------------------------------------------------------------- 389 &namptr ! Poleward Transport Diagnostic (default: OFF)390 390 !----------------------------------------------------------------------- 391 391 / -
NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/cfgs/ORCA2_OFF_TRC/EXPREF/namelist_cfg
r13553 r13741 370 370 !! !! 371 371 !! namtrd dynamics and/or tracer trends (default: OFF) 372 !! namptr Poleward Transport Diagnostics (default: OFF)373 372 !! namhsb Heat and salt budgets (default: OFF) 374 373 !! namdiu Cool skin and warm layer models (default: OFF) … … 385 384 / 386 385 !----------------------------------------------------------------------- 387 &namptr ! Poleward Transport Diagnostic (default: OFF)388 386 !----------------------------------------------------------------------- 389 387 / -
NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/cfgs/ORCA2_SAS_ICE/EXPREF/file_def_nemo-ice.xml
r9896 r13741 52 52 <field field_ref="normstr" name="normstr" /> 53 53 <field field_ref="sheastr" name="sheastr" /> 54 <field field_ref="isig1" name="isig1" /> 55 <field field_ref="isig2" name="isig2" /> 56 <field field_ref="isig3" name="isig3" /> 54 <field field_ref="sig1_pnorm" name="sig1_pnorm"/> 55 <field field_ref="sig2_pnorm" name="sig2_pnorm"/> 57 56 58 57 <!-- heat fluxes --> -
NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/cfgs/ORCA2_SAS_ICE/EXPREF/namelist_cfg
r13553 r13741 171 171 !! !! 172 172 !! namtrd dynamics and/or tracer trends (default: OFF) 173 !! namptr Poleward Transport Diagnostics (default: OFF)174 173 !! namhsb Heat and salt budgets (default: OFF) 175 174 !! namdiu Cool skin and warm layer models (default: OFF) -
NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/cfgs/SHARED/field_def_nemo-ice.xml
r13553 r13741 75 75 <field id="utau_bi" long_name="X-component of ocean bottom stress on sea ice -landfast" standard_name="ocean_bottom_upward_x_stress" unit="N/m2" /> 76 76 <field id="vtau_bi" long_name="Y-component of ocean bottom stress on sea ice -landfast" standard_name="ocean_bottom_upward_y_stress" unit="N/m2" /> 77 <field id="isig1" long_name="1st principal stress component for EVP rhg" unit="" /> 78 <field id="isig2" long_name="2nd principal stress component for EVP rhg" unit="" /> 79 <field id="isig3" long_name="convergence measure for EVP rheology (must be around 1)" unit="" /> 77 <field id="sig1_pnorm" long_name="P-normalized 1st principal stress component" unit="" /> 78 <field id="sig2_pnorm" long_name="P-normalized 2nd principal stress component" unit="" /> 80 79 <field id="normstr" long_name="Average normal stress in sea ice" standard_name="average_normal_stress" unit="N/m" /> 81 80 <field id="sheastr" long_name="Maximum shear stress in sea ice" standard_name="maximum_shear_stress" unit="N/m" /> … … 165 164 166 165 <!-- diags --> 167 <field id="hfxdhc" long_name="Heat content variation in snow and ice (neg = ice cooling)" unit="W/m2" /> 166 <field id="hfxdhc" long_name="Heat content variation in snow and ice (neg = ice cooling)" unit="W/m2" /> 167 <!-- available if ln_icediachk=T --> 168 <field id="icedrift_mass" long_name="Ice mass drift (conservation check)" unit="kg/m2/s" /> 169 <field id="icedrift_salt" long_name="Ice salt drift (conservation check)" unit="kg/m2/s" /> 170 <field id="icedrift_heat" long_name="Ice heat drift (conservation check)" unit="W/m2" /> 168 171 169 172 <!-- sbcssm variables --> … … 400 403 <field field_ref="normstr" name="normstr" /> 401 404 <field field_ref="sheastr" name="sheastr" /> 402 <field field_ref="isig1" name="isig1" /> 403 <field field_ref="isig2" name="isig2" /> 404 <field field_ref="isig3" name="isig3" /> 405 <field field_ref="sig1_pnorm" name="sig1_pnorm"/> 406 <field field_ref="sig2_pnorm" name="sig2_pnorm"/> 405 407 406 408 <!-- heat fluxes --> … … 459 461 <field field_ref="vfxthin" name="vfxthin" /> 460 462 461 <!-- diag error for negative ice volume after advection -->462 <field field_ref="iceneg_pres" name="sineg_pres" />463 <field field_ref="iceneg_volu" name="sineg_volu" />464 <field field_ref="iceneg_hfx" name="sineg_hfx" />465 463 </field_group> 466 464 -
NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/cfgs/SHARED/namelist_ref
r13553 r13741 466 466 !----------------------------------------------------------------------- 467 467 ln_rnf_mouth = .false. ! specific treatment at rivers mouths 468 rn_hrnf = 15.e0 ! depth over which enhanced vertical mixing is used (ln_rnf_mouth=T)469 rn_avt_rnf = 1.e-3 ! value of the additional vertical mixing coef. [m2/s] (ln_rnf_mouth=T)470 rn_rfact = 1.e0 ! multiplicative factor for runoff468 rn_hrnf = 15.e0 ! depth over which enhanced vertical mixing is used (ln_rnf_mouth=T) 469 rn_avt_rnf = 1.e-3 ! value of the additional vertical mixing coef. [m2/s] (ln_rnf_mouth=T) 470 rn_rfact = 1.e0 ! multiplicative factor for runoff 471 471 ln_rnf_depth = .false. ! read in depth information for runoff 472 ln_rnf_tem = .false. ! read in temperature information for runoff 473 ln_rnf_sal = .false. ! read in salinity information for runoff 474 ln_rnf_depth_ini = .false. ! compute depth at initialisation from runoff file 475 rn_rnf_max = 5.735e-4 ! max value of the runoff climatologie over global domain ( ln_rnf_depth_ini = .true ) 476 rn_dep_max = 150. ! depth over which runoffs is spread ( ln_rnf_depth_ini = .true ) 477 nn_rnf_depth_file = 0 ! create (=1) a runoff depth file or not (=0) 478 479 cn_dir = './' ! root directory for the runoff data location 472 ln_rnf_tem = .false. ! read in temperature information for runoff 473 ln_rnf_sal = .false. ! read in salinity information for runoff 474 ln_rnf_icb = .false. ! read iceberg flux 475 ln_rnf_depth_ini = .false. ! compute depth at initialisation from runoff file 476 rn_rnf_max = 5.735e-4 ! max value of the runoff climatologie over global domain ( ln_rnf_depth_ini = .true ) 477 rn_dep_max = 150. ! depth over which runoffs is spread ( ln_rnf_depth_ini = .true ) 478 nn_rnf_depth_file = 0 ! create (=1) a runoff depth file or not (=0) 479 480 cn_dir = './' ! root directory for the runoff data location 480 481 !___________!_________________________!___________________!___________!_____________!________!___________!__________________!__________!_______________! 481 482 ! ! file name ! frequency (hours) ! variable ! time interp.! clim ! 'yearly'/ ! weights filename ! rotation ! land/sea mask ! 482 483 ! ! ! (if <0 months) ! name ! (logical) ! (T/F) ! 'monthly' ! ! pairing ! filename ! 483 484 sn_rnf = 'runoff_core_monthly' , -1. , 'sorunoff', .true. , .true. , 'yearly' , '' , '' , '' 484 sn_cnf = 'runoff_core_monthly' , 0. , 'socoefr0', .false. , .true. , 'yearly' , '' , '' , ''485 sn_cnf = 'runoff_core_monthly' , -12. , 'socoefr0', .false. , .true. , 'yearly' , '' , '' , '' 485 486 sn_s_rnf = 'runoffs' , 24. , 'rosaline', .true. , .true. , 'yearly' , '' , '' , '' 486 487 sn_t_rnf = 'runoffs' , 24. , 'rotemper', .true. , .true. , 'yearly' , '' , '' , '' 487 sn_dep_rnf = 'runoffs' , 0. , 'rodepth' , .false. , .true. , 'yearly' , '' , '' , '' 488 sn_i_rnf = 'NOT USED' , 24. , 'xxxxxxxx', .true. , .true. , 'yearly' , '' , '' , '' 489 sn_dep_rnf = 'runoffs' , -12. , 'rodepth' , .false. , .true. , 'yearly' , '' , '' , '' 488 490 / 489 491 !----------------------------------------------------------------------- … … 1056 1058 ln_dynrnf = .false. ! runoffs option enabled (T) or not (F) 1057 1059 ln_dynrnf_depth = .false. ! runoffs is spread in vertical (T) or not (F) 1058 ! fwbcorr = 3.786e-06! annual global mean of empmr for ssh correction1060 fwbcorr = 0.0 ! annual global mean of empmr for ssh correction 1059 1061 1060 1062 cn_dir = './' ! root directory for the ocean data location … … 1240 1242 !! !! 1241 1243 !! namtrd dynamics and/or tracer trends (default: OFF) 1242 !! namptr Poleward Transport Diagnostics (default: OFF)1243 1244 !! namhsb Heat and salt budgets (default: OFF) 1244 1245 !! namdiu Cool skin and warm layer models (default: OFF) -
NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/cfgs/SPITZ12/EXPREF/file_def_nemo-ice.xml
r11536 r13741 53 53 <field field_ref="normstr" name="normstr" /> 54 54 <field field_ref="sheastr" name="sheastr" /> 55 <field field_ref="isig1" name="isig1" /> 56 <field field_ref="isig2" name="isig2" /> 57 <field field_ref="isig3" name="isig3" /> 55 <field field_ref="sig1_pnorm" name="sig1_pnorm"/> 56 <field field_ref="sig2_pnorm" name="sig2_pnorm"/> 58 57 59 58 <!-- heat fluxes --> -
NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/cfgs/SPITZ12/EXPREF/namelist_cfg
r13553 r13741 359 359 !! !! 360 360 !! namtrd dynamics and/or tracer trends (default: OFF) 361 !! namptr Poleward Transport Diagnostics (default: OFF)362 361 !! namhsb Heat and salt budgets (default: OFF) 363 362 !! namdiu Cool skin and warm layer models (default: OFF) -
NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/cfgs/WED025/EXPREF/file_def_nemo-ice.xml
r12905 r13741 50 50 <field field_ref="normstr" name="normstr" /> 51 51 <field field_ref="sheastr" name="sheastr" /> 52 <field field_ref="isig1" name="isig1" /> 53 <field field_ref="isig2" name="isig2" /> 54 <field field_ref="isig3" name="isig3" /> 52 <field field_ref="sig1_pnorm" name="sig1_pnorm"/> 53 <field field_ref="sig2_pnorm" name="sig2_pnorm"/> 55 54 56 55 <!-- heat fluxes --> -
NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/cfgs/WED025/EXPREF/namelist_cfg
r13553 r13741 567 567 !! !! 568 568 !! namtrd dynamics and/or tracer trends (default: OFF) 569 !! namptr Poleward Transport Diagnostics (default: OFF)570 569 !! namhsb Heat and salt budgets (default: OFF) 571 570 !! namdiu Cool skin and warm layer models (default: OFF) … … 584 583 / 585 584 !----------------------------------------------------------------------- 586 &namptr ! Poleward Transport Diagnostic (default: OFF)587 585 !----------------------------------------------------------------------- 588 586 / -
NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/src/ICE/ice.F90
r13553 r13741 391 391 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: diag_vice !: ice volume variation [m/s] 392 392 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: diag_vsnw !: snw volume variation [m/s] 393 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: diag_aice !: ice conc. variation [s-1] 394 ! 395 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: diag_adv_mass !: advection of mass (kg/m2/s) 396 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: diag_adv_salt !: advection of salt (g/m2/s) 397 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: diag_adv_heat !: advection of heat (W/m2) 393 398 ! 394 399 !!---------------------------------------------------------------------- … … 493 498 ! * Ice diagnostics 494 499 ii = ii + 1 495 ALLOCATE( diag_trp_vi(jpi,jpj) , diag_trp_vs (jpi,jpj) , diag_trp_ei(jpi,jpj), & 496 & diag_trp_es(jpi,jpj) , diag_trp_sv (jpi,jpj) , diag_heat (jpi,jpj), & 497 & diag_sice (jpi,jpj) , diag_vice (jpi,jpj) , diag_vsnw (jpi,jpj), STAT=ierr(ii) ) 500 ALLOCATE( diag_trp_vi(jpi,jpj) , diag_trp_vs (jpi,jpj) , diag_trp_ei(jpi,jpj), & 501 & diag_trp_es(jpi,jpj) , diag_trp_sv (jpi,jpj) , diag_heat (jpi,jpj), & 502 & diag_sice (jpi,jpj) , diag_vice (jpi,jpj) , diag_vsnw (jpi,jpj), diag_aice(jpi,jpj), & 503 & diag_adv_mass(jpi,jpj), diag_adv_salt(jpi,jpj), diag_adv_heat(jpi,jpj), STAT=ierr(ii) ) 498 504 499 505 ! * Ice conservation -
NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/src/ICE/ice1d.F90
r13553 r13741 145 145 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: sst_1d 146 146 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: sss_1d 147 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: frq_m_1d 147 148 148 149 ! convergence check … … 225 226 ! 226 227 ii = ii + 1 227 ALLOCATE( sst_1d(jpij) , sss_1d(jpij) , STAT=ierr(ii) )228 ALLOCATE( sst_1d(jpij) , sss_1d(jpij) , frq_m_1d(jpij) , STAT=ierr(ii) ) 228 229 ! 229 230 ii = ii + 1 -
NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/src/ICE/icecor.F90
r13553 r13741 55 55 INTEGER :: ji, jj, jk, jl ! dummy loop indices 56 56 REAL(wp) :: zsal, zzc 57 REAL(wp), DIMENSION(jpi,jpj) :: zafx ! concentration trends diag58 57 !!---------------------------------------------------------------------- 59 58 ! controls … … 95 94 zsal = sv_i(ji,jj,jl) 96 95 sv_i(ji,jj,jl) = MIN( MAX( rn_simin*v_i(ji,jj,jl) , sv_i(ji,jj,jl) ) , rn_simax*v_i(ji,jj,jl) ) 97 sfx_res(ji,jj) = sfx_res(ji,jj) - ( sv_i(ji,jj,jl) - zsal ) * zzc ! associated salt flux 96 IF( kn /= 0 ) & ! no ice-ocean exchanges if kn=0 (for bdy for instance) otherwise conservation diags will fail 97 & sfx_res(ji,jj) = sfx_res(ji,jj) - ( sv_i(ji,jj,jl) - zsal ) * zzc ! associated salt flux 98 98 END_2D 99 99 END DO 100 100 ENDIF 101 ! !-----------------------------------------------------102 CALL ice_var_zapsmall ! Zap small values !103 ! !-----------------------------------------------------104 101 102 IF( kn /= 0 ) THEN ! no zapsmall if kn=0 (for bdy for instance) because we do not want ice-ocean exchanges (wfx,sfx,hfx) 103 ! otherwise conservation diags will fail 104 ! !----------------------------------------------------- 105 CALL ice_var_zapsmall ! Zap small values ! 106 ! !----------------------------------------------------- 107 ENDIF 105 108 ! !----------------------------------------------------- 106 109 IF( kn == 2 ) THEN ! Ice drift case: Corrections to avoid wrong values ! … … 115 118 CALL lbc_lnk_multi( 'icecor', u_ice, 'U', -1.0_wp, v_ice, 'V', -1.0_wp ) 116 119 ENDIF 117 118 ! !-----------------------------------------------------119 SELECT CASE( kn ) ! Diagnostics !120 ! !-----------------------------------------------------121 CASE( 1 ) !--- dyn trend diagnostics122 !123 IF( ln_icediachk .OR. iom_use('hfxdhc') ) THEN124 diag_heat(:,:) = - SUM(SUM( e_i (:,:,1:nlay_i,:) - e_i_b (:,:,1:nlay_i,:), dim=4 ), dim=3 ) * r1_Dt_ice & ! W.m-2125 & - SUM(SUM( e_s (:,:,1:nlay_s,:) - e_s_b (:,:,1:nlay_s,:), dim=4 ), dim=3 ) * r1_Dt_ice126 diag_sice(:,:) = SUM( sv_i(:,:,:) - sv_i_b(:,:,:) , dim=3 ) * r1_Dt_ice * rhoi127 diag_vice(:,:) = SUM( v_i (:,:,:) - v_i_b (:,:,:) , dim=3 ) * r1_Dt_ice * rhoi128 diag_vsnw(:,:) = SUM( v_s (:,:,:) - v_s_b (:,:,:) , dim=3 ) * r1_Dt_ice * rhos129 ENDIF130 ! ! concentration tendency (dynamics)131 IF( iom_use('afxdyn') .OR. iom_use('afxthd') .OR. iom_use('afxtot') ) THEN132 zafx(:,:) = SUM( a_i(:,:,:) - a_i_b(:,:,:), dim=3 ) * r1_Dt_ice133 CALL iom_put( 'afxdyn' , zafx )134 ENDIF135 !136 CASE( 2 ) !--- thermo trend diagnostics & ice aging137 !138 oa_i(:,:,:) = oa_i(:,:,:) + a_i(:,:,:) * rDt_ice ! ice natural aging incrementation139 !140 IF( ln_icediachk .OR. iom_use('hfxdhc') ) THEN141 diag_heat(:,:) = diag_heat(:,:) &142 & - SUM(SUM( e_i (:,:,1:nlay_i,:) - e_i_b (:,:,1:nlay_i,:), dim=4 ), dim=3 ) * r1_Dt_ice &143 & - SUM(SUM( e_s (:,:,1:nlay_s,:) - e_s_b (:,:,1:nlay_s,:), dim=4 ), dim=3 ) * r1_Dt_ice144 diag_sice(:,:) = diag_sice(:,:) &145 & + SUM( sv_i(:,:,:) - sv_i_b(:,:,:) , dim=3 ) * r1_Dt_ice * rhoi146 diag_vice(:,:) = diag_vice(:,:) &147 & + SUM( v_i (:,:,:) - v_i_b (:,:,:) , dim=3 ) * r1_Dt_ice * rhoi148 diag_vsnw(:,:) = diag_vsnw(:,:) &149 & + SUM( v_s (:,:,:) - v_s_b (:,:,:) , dim=3 ) * r1_Dt_ice * rhos150 CALL iom_put ( 'hfxdhc' , diag_heat )151 ENDIF152 ! ! concentration tendency (total + thermo)153 IF( iom_use('afxdyn') .OR. iom_use('afxthd') .OR. iom_use('afxtot') ) THEN154 zafx(:,:) = zafx(:,:) + SUM( a_i(:,:,:) - a_i_b(:,:,:), dim=3 ) * r1_Dt_ice155 CALL iom_put( 'afxthd' , SUM( a_i(:,:,:) - a_i_b(:,:,:), dim=3 ) * r1_Dt_ice )156 CALL iom_put( 'afxtot' , zafx )157 ENDIF158 !159 END SELECT160 120 ! 161 121 ! controls -
NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/src/ICE/icectl.F90
r13553 r13741 43 43 PUBLIC ice_prt 44 44 PUBLIC ice_prt3D 45 PUBLIC ice_drift_wri 46 PUBLIC ice_drift_init 45 47 46 48 ! thresold rates for conservation … … 49 51 REAL(wp), PARAMETER :: zchk_s = 2.5e-6 ! g/m2/s <=> 1e-6 m of ice per hour spuriously gained/lost (considering s=10g/kg) 50 52 REAL(wp), PARAMETER :: zchk_t = 7.5e-2 ! W/m2 <=> 1e-6 m of ice per hour spuriously gained/lost (considering Lf=3e5J/kg) 53 54 ! for drift outputs 55 CHARACTER(LEN=50) :: clname="icedrift_diagnostics.ascii" ! ascii filename 56 INTEGER :: numicedrift ! outfile unit 57 REAL(wp) :: rdiag_icemass, rdiag_icesalt, rdiag_iceheat 58 REAL(wp) :: rdiag_adv_icemass, rdiag_adv_icesalt, rdiag_adv_iceheat 51 59 52 60 !! * Substitutions … … 132 140 133 141 ! -- advection scheme is conservative? -- ! 134 zvtrp = glob_sum( 'icectl', ( diag_trp_vi * rhoi + diag_trp_vs * rhos ) * e1e2t ) ! must be close to 0 (only for Prather)135 zetrp = glob_sum( 'icectl', ( diag_trp_ei + diag_trp_es ) * e1e2t ) ! must be close to 0 (only for Prather)142 zvtrp = glob_sum( 'icectl', diag_adv_mass * e1e2t ) 143 zetrp = glob_sum( 'icectl', diag_adv_heat * e1e2t ) 136 144 137 145 ! ice area (+epsi10 to set a threshold > 0 when there is no ice) … … 156 164 & WRITE(numout,*) cd_routine,' : violation a_i > amax = ',zdiag_amax 157 165 ! check if advection scheme is conservative 158 ! only check for Prather because Ultimate-Macho uses corrective fluxes (wfx etc) 159 ! so the formulation for conservation is different (and not coded) 160 ! it does not mean UM is not conservative (it is checked with above prints) => update (09/2019): same for Prather now 161 !IF( ln_adv_Pra .AND. ABS(zvtrp) > zchk_m * rn_icechk_glo * zarea .AND. cd_routine == 'icedyn_adv' ) & 162 ! & WRITE(numout,*) cd_routine,' : violation adv scheme [kg] = ',zvtrp * rDt_ice 166 IF( ABS(zvtrp) > zchk_m * rn_icechk_glo * zarea .AND. cd_routine == 'icedyn_adv' ) & 167 & WRITE(numout,*) cd_routine,' : violation adv scheme [kg] = ',zvtrp * rdt_ice 168 IF( ABS(zetrp) > zchk_t * rn_icechk_glo * zarea .AND. cd_routine == 'icedyn_adv' ) & 169 & WRITE(numout,*) cd_routine,' : violation adv scheme [J] = ',zetrp * rdt_ice 163 170 ENDIF 164 171 ! … … 186 193 ! water flux 187 194 ! -- mass diag -- ! 188 zdiag_mass = glob_sum( 'icectl', ( wfx_ice + wfx_snw + wfx_spr + wfx_sub + diag_vice + diag_vsnw ) * e1e2t ) 195 zdiag_mass = glob_sum( 'icectl', ( wfx_ice + wfx_snw + wfx_spr + wfx_sub & 196 & + diag_vice + diag_vsnw - diag_adv_mass ) * e1e2t ) 189 197 190 198 ! -- salt diag -- ! 191 zdiag_salt = glob_sum( 'icectl', ( sfx + diag_sice ) * e1e2t )199 zdiag_salt = glob_sum( 'icectl', ( sfx + diag_sice - diag_adv_salt ) * e1e2t ) 192 200 193 201 ! -- heat diag -- ! 194 ! clem: not the good formulation 195 !!zdiag_heat = glob_sum( 'icectl', ( qt_oce_ai - qt_atm_oi + diag_heat + hfx_thd + hfx_dyn + hfx_res + hfx_sub + hfx_spr & 196 !! & ) * e1e2t ) 202 zdiag_heat = glob_sum( 'icectl', ( qt_oce_ai - qt_atm_oi + diag_heat - diag_adv_heat ) * e1e2t ) 203 ! equivalent to this: 204 !!zdiag_heat = glob_sum( 'icectl', ( -diag_heat + hfx_sum + hfx_bom + hfx_bog + hfx_dif + hfx_opw + hfx_snw & 205 !! & - hfx_thd - hfx_dyn - hfx_res - hfx_sub - hfx_spr & 206 !! & ) * e1e2t ) 197 207 198 208 ! ice area (+epsi10 to set a threshold > 0 when there is no ice) … … 204 214 IF( ABS(zdiag_salt) > zchk_s * rn_icechk_glo * zarea ) & 205 215 & WRITE(numout,*) cd_routine,' : violation salt cons. [g] = ',zdiag_salt * rDt_ice 206 !!IF( ABS(zdiag_heat) > zchk_t * rn_icechk_glo * zarea ) WRITE(numout,*) cd_routine,' : violation heat cons. [J] = ',zdiag_heat * rDt_ice 216 IF( ABS(zdiag_heat) > zchk_t * rn_icechk_glo * zarea ) & 217 & WRITE(numout,*) cd_routine,' : violation heat cons. [J] = ',zdiag_heat * rDt_ice 207 218 ENDIF 208 219 ! … … 725 736 726 737 END SUBROUTINE ice_prt3D 738 739 740 SUBROUTINE ice_drift_wri( kt ) 741 !!------------------------------------------------------------------- 742 !! *** ROUTINE ice_drift_wri *** 743 !! 744 !! ** Purpose : conservation of mass, salt and heat 745 !! write the drift in a ascii file at each time step 746 !! and the total run drifts 747 !!------------------------------------------------------------------- 748 INTEGER, INTENT(in) :: kt ! ice time-step index 749 ! 750 INTEGER :: ji, jj 751 REAL(wp) :: zdiag_mass, zdiag_salt, zdiag_heat, zdiag_adv_mass, zdiag_adv_salt, zdiag_adv_heat 752 REAL(wp), DIMENSION(jpi,jpj) :: zdiag_mass2D, zdiag_salt2D, zdiag_heat2D 753 !!------------------------------------------------------------------- 754 ! 755 IF( kt == nit000 .AND. lwp ) THEN 756 WRITE(numout,*) 757 WRITE(numout,*) 'ice_drift_wri: sea-ice drifts' 758 WRITE(numout,*) '~~~~~~~~~~~~~' 759 ENDIF 760 ! 761 ! 2D budgets (must be close to 0) 762 IF( iom_use('icedrift_mass') .OR. iom_use('icedrift_salt') .OR. iom_use('icedrift_heat') ) THEN 763 DO_2D( 1, 1, 1, 1 ) 764 zdiag_mass2D(ji,jj) = wfx_ice(ji,jj) + wfx_snw(ji,jj) + wfx_spr(ji,jj) + wfx_sub(ji,jj) & 765 & + diag_vice(ji,jj) + diag_vsnw(ji,jj) - diag_adv_mass(ji,jj) 766 zdiag_salt2D(ji,jj) = sfx(ji,jj) + diag_sice(ji,jj) - diag_adv_salt(ji,jj) 767 zdiag_heat2D(ji,jj) = qt_oce_ai(ji,jj) - qt_atm_oi(ji,jj) + diag_heat(ji,jj) - diag_adv_heat(ji,jj) 768 END_2D 769 ! 770 ! write outputs 771 CALL iom_put( 'icedrift_mass', zdiag_mass2D ) 772 CALL iom_put( 'icedrift_salt', zdiag_salt2D ) 773 CALL iom_put( 'icedrift_heat', zdiag_heat2D ) 774 ENDIF 775 776 ! -- mass diag -- ! 777 zdiag_mass = glob_sum( 'icectl', ( wfx_ice + wfx_snw + wfx_spr + wfx_sub & 778 & + diag_vice + diag_vsnw - diag_adv_mass ) * e1e2t ) * rdt_ice 779 zdiag_adv_mass = glob_sum( 'icectl', diag_adv_mass * e1e2t ) * rDt_ice 780 781 ! -- salt diag -- ! 782 zdiag_salt = glob_sum( 'icectl', ( sfx + diag_sice - diag_adv_salt ) * e1e2t ) * rdt_ice * 1.e-3 783 zdiag_adv_salt = glob_sum( 'icectl', diag_adv_salt * e1e2t ) * rDt_ice * 1.e-3 784 785 ! -- heat diag -- ! 786 zdiag_heat = glob_sum( 'icectl', ( qt_oce_ai - qt_atm_oi + diag_heat - diag_adv_heat ) * e1e2t ) 787 zdiag_adv_heat = glob_sum( 'icectl', diag_adv_heat * e1e2t ) 788 789 ! ! write out to file 790 IF( lwp ) THEN 791 ! check global drift (must be close to 0) 792 WRITE(numicedrift,FMT='(2x,i6,3x,a19,4x,f25.5)') kt, 'mass drift [kg]', zdiag_mass 793 WRITE(numicedrift,FMT='(11x, a19,4x,f25.5)') 'salt drift [kg]', zdiag_salt 794 WRITE(numicedrift,FMT='(11x, a19,4x,f25.5)') 'heat drift [W] ', zdiag_heat 795 ! check drift from advection scheme (can be /=0 with bdy but not sure why) 796 WRITE(numicedrift,FMT='(11x, a19,4x,f25.5)') 'mass drift adv [kg]', zdiag_adv_mass 797 WRITE(numicedrift,FMT='(11x, a19,4x,f25.5)') 'salt drift adv [kg]', zdiag_adv_salt 798 WRITE(numicedrift,FMT='(11x, a19,4x,f25.5)') 'heat drift adv [W] ', zdiag_adv_heat 799 ENDIF 800 ! ! drifts 801 rdiag_icemass = rdiag_icemass + zdiag_mass 802 rdiag_icesalt = rdiag_icesalt + zdiag_salt 803 rdiag_iceheat = rdiag_iceheat + zdiag_heat 804 rdiag_adv_icemass = rdiag_adv_icemass + zdiag_adv_mass 805 rdiag_adv_icesalt = rdiag_adv_icesalt + zdiag_adv_salt 806 rdiag_adv_iceheat = rdiag_adv_iceheat + zdiag_adv_heat 807 ! 808 ! ! output drifts and close ascii file 809 IF( kt == nitend - nn_fsbc + 1 .AND. lwp ) THEN 810 ! to ascii file 811 WRITE(numicedrift,*) '******************************************' 812 WRITE(numicedrift,FMT='(3x,a23,6x,E10.2)') 'Run mass drift [kg]', rdiag_icemass 813 WRITE(numicedrift,FMT='(3x,a23,6x,E10.2)') 'Run mass drift adv [kg]', rdiag_adv_icemass 814 WRITE(numicedrift,*) '******************************************' 815 WRITE(numicedrift,FMT='(3x,a23,6x,E10.2)') 'Run salt drift [kg]', rdiag_icesalt 816 WRITE(numicedrift,FMT='(3x,a23,6x,E10.2)') 'Run salt drift adv [kg]', rdiag_adv_icesalt 817 WRITE(numicedrift,*) '******************************************' 818 WRITE(numicedrift,FMT='(3x,a23,6x,E10.2)') 'Run heat drift [W] ', rdiag_iceheat 819 WRITE(numicedrift,FMT='(3x,a23,6x,E10.2)') 'Run heat drift adv [W] ', rdiag_adv_iceheat 820 CLOSE( numicedrift ) 821 ! 822 ! to ocean output 823 WRITE(numout,*) 824 WRITE(numout,*) 'ice_drift_wri: ice drifts information for the run ' 825 WRITE(numout,*) '~~~~~~~~~~~~~' 826 ! check global drift (must be close to 0) 827 WRITE(numout,*) ' sea-ice mass drift [kg] = ', rdiag_icemass 828 WRITE(numout,*) ' sea-ice salt drift [kg] = ', rdiag_icesalt 829 WRITE(numout,*) ' sea-ice heat drift [W] = ', rdiag_iceheat 830 ! check drift from advection scheme (can be /=0 with bdy but not sure why) 831 WRITE(numout,*) ' sea-ice mass drift adv [kg] = ', rdiag_adv_icemass 832 WRITE(numout,*) ' sea-ice salt drift adv [kg] = ', rdiag_adv_icesalt 833 WRITE(numout,*) ' sea-ice heat drift adv [W] = ', rdiag_adv_iceheat 834 ENDIF 835 ! 836 END SUBROUTINE ice_drift_wri 837 838 SUBROUTINE ice_drift_init 839 !!---------------------------------------------------------------------- 840 !! *** ROUTINE ice_drift_init *** 841 !! 842 !! ** Purpose : create output file, initialise arrays 843 !!---------------------------------------------------------------------- 844 ! 845 IF( .NOT.ln_icediachk ) RETURN ! exit 846 ! 847 IF(lwp) THEN 848 WRITE(numout,*) 849 WRITE(numout,*) 'ice_drift_init: Output ice drifts to ',TRIM(clname), ' file' 850 WRITE(numout,*) '~~~~~~~~~~~~~' 851 WRITE(numout,*) 852 ! 853 ! create output ascii file 854 CALL ctl_opn( numicedrift, clname, 'UNKNOWN', 'FORMATTED', 'SEQUENTIAL', 1, numout, lwp, narea ) 855 WRITE(numicedrift,*) 'Timestep Drifts' 856 WRITE(numicedrift,*) '******************************************' 857 ENDIF 858 ! 859 rdiag_icemass = 0._wp 860 rdiag_icesalt = 0._wp 861 rdiag_iceheat = 0._wp 862 rdiag_adv_icemass = 0._wp 863 rdiag_adv_icesalt = 0._wp 864 rdiag_adv_iceheat = 0._wp 865 ! 866 END SUBROUTINE ice_drift_init 727 867 728 868 #else -
NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/src/ICE/icedyn_adv_pra.F90
r13553 r13741 88 88 INTEGER :: ji, jj, jk, jl, jt ! dummy loop indices 89 89 INTEGER :: icycle ! number of sub-timestep for the advection 90 REAL(wp) :: zdt 90 REAL(wp) :: zdt, z1_dt ! - - 91 91 REAL(wp), DIMENSION(1) :: zcflprv, zcflnow ! for global communication 92 92 REAL(wp), DIMENSION(jpi,jpj) :: zati1, zati2 … … 100 100 REAL(wp), DIMENSION(jpi,jpj,nlay_s,jpl) :: z0es 101 101 REAL(wp), DIMENSION(jpi,jpj,nlay_i,jpl) :: z0ei 102 !! diagnostics 103 REAL(wp), DIMENSION(jpi,jpj) :: zdiag_adv_mass, zdiag_adv_salt, zdiag_adv_heat 102 104 !!---------------------------------------------------------------------- 103 105 ! … … 109 111 ELSEWHERE ; zs_i(:,:,:) = 0._wp 110 112 END WHERE 111 DO jl = 1, jpl 112 DO_2D( 0, 0, 0, 0 ) 113 zhip_max(ji,jj,jl) = MAX( epsi20, ph_ip(ji,jj,jl), ph_ip(ji+1,jj ,jl), ph_ip(ji ,jj+1,jl), & 114 & ph_ip(ji-1,jj ,jl), ph_ip(ji ,jj-1,jl), & 115 & ph_ip(ji+1,jj+1,jl), ph_ip(ji-1,jj-1,jl), & 116 & ph_ip(ji+1,jj-1,jl), ph_ip(ji-1,jj+1,jl) ) 117 zhi_max (ji,jj,jl) = MAX( epsi20, ph_i (ji,jj,jl), ph_i (ji+1,jj ,jl), ph_i (ji ,jj+1,jl), & 118 & ph_i (ji-1,jj ,jl), ph_i (ji ,jj-1,jl), & 119 & ph_i (ji+1,jj+1,jl), ph_i (ji-1,jj-1,jl), & 120 & ph_i (ji+1,jj-1,jl), ph_i (ji-1,jj+1,jl) ) 121 zhs_max (ji,jj,jl) = MAX( epsi20, ph_s (ji,jj,jl), ph_s (ji+1,jj ,jl), ph_s (ji ,jj+1,jl), & 122 & ph_s (ji-1,jj ,jl), ph_s (ji ,jj-1,jl), & 123 & ph_s (ji+1,jj+1,jl), ph_s (ji-1,jj-1,jl), & 124 & ph_s (ji+1,jj-1,jl), ph_s (ji-1,jj+1,jl) ) 125 zsi_max (ji,jj,jl) = MAX( epsi20, zs_i (ji,jj,jl), zs_i (ji+1,jj ,jl), zs_i (ji ,jj+1,jl), & 126 & zs_i (ji-1,jj ,jl), zs_i (ji ,jj-1,jl), & 127 & zs_i (ji+1,jj+1,jl), zs_i (ji-1,jj-1,jl), & 128 & zs_i (ji+1,jj-1,jl), zs_i (ji-1,jj+1,jl) ) 129 END_2D 130 END DO 113 CALL icemax3D( ph_i , zhi_max ) 114 CALL icemax3D( ph_s , zhs_max ) 115 CALL icemax3D( ph_ip, zhip_max) 116 CALL icemax3D( zs_i , zsi_max ) 131 117 CALL lbc_lnk_multi( 'icedyn_adv_pra', zhi_max, 'T', 1._wp, zhs_max, 'T', 1._wp, zhip_max, 'T', 1._wp, zsi_max, 'T', 1._wp ) 132 118 ! … … 141 127 ELSEWHERE ; ze_s(:,:,jk,:) = 0._wp 142 128 END WHERE 143 END DO 144 DO jl = 1, jpl 145 DO_3D( 0, 0, 0, 0, 1, nlay_i ) 146 zei_max(ji,jj,jk,jl) = MAX( epsi20, ze_i(ji,jj,jk,jl), ze_i(ji+1,jj ,jk,jl), ze_i(ji ,jj+1,jk,jl), & 147 & ze_i(ji-1,jj ,jk,jl), ze_i(ji ,jj-1,jk,jl), & 148 & ze_i(ji+1,jj+1,jk,jl), ze_i(ji-1,jj-1,jk,jl), & 149 & ze_i(ji+1,jj-1,jk,jl), ze_i(ji-1,jj+1,jk,jl) ) 150 END_3D 151 END DO 152 DO jl = 1, jpl 153 DO_3D( 0, 0, 0, 0, 1, nlay_s ) 154 zes_max(ji,jj,jk,jl) = MAX( epsi20, ze_s(ji,jj,jk,jl), ze_s(ji+1,jj ,jk,jl), ze_s(ji ,jj+1,jk,jl), & 155 & ze_s(ji-1,jj ,jk,jl), ze_s(ji ,jj-1,jk,jl), & 156 & ze_s(ji+1,jj+1,jk,jl), ze_s(ji-1,jj-1,jk,jl), & 157 & ze_s(ji+1,jj-1,jk,jl), ze_s(ji-1,jj+1,jk,jl) ) 158 END_3D 159 END DO 160 CALL lbc_lnk( 'icedyn_adv_pra', zei_max, 'T', 1. ) 161 CALL lbc_lnk( 'icedyn_adv_pra', zes_max, 'T', 1. ) 129 END DO 130 CALL icemax4D( ze_i , zei_max ) 131 CALL icemax4D( ze_s , zes_max ) 132 CALL lbc_lnk( 'icedyn_adv_pra', zei_max, 'T', 1._wp ) 133 CALL lbc_lnk( 'icedyn_adv_pra', zes_max, 'T', 1._wp ) 162 134 ! 163 135 ! … … 175 147 ENDIF 176 148 zdt = rDt_ice / REAL(icycle) 149 z1_dt = 1._wp / zdt 177 150 178 151 ! --- transport --- ! … … 181 154 182 155 DO jt = 1, icycle 156 157 ! diagnostics 158 zdiag_adv_mass(:,:) = SUM( pv_i(:,:,:) , dim=3 ) * rhoi + SUM( pv_s(:,:,:) , dim=3 ) * rhos 159 zdiag_adv_salt(:,:) = SUM( psv_i(:,:,:) , dim=3 ) * rhoi 160 zdiag_adv_heat(:,:) = - SUM(SUM( pe_i(:,:,1:nlay_i,:) , dim=4 ), dim=3 ) & 161 & - SUM(SUM( pe_s(:,:,1:nlay_s,:) , dim=4 ), dim=3 ) 183 162 184 163 ! record at_i before advection (for open water) … … 279 258 CALL adv_x( zdt , zudy , 0._wp , zarea , z0vl , sxvl , sxxvl , syvl , syyvl , sxyvl ) 280 259 ENDIF 281 ENDIF 282 ! 260 ENDIF 261 ! 262 ENDIF 263 264 ! --- Lateral boundary conditions --- ! 265 ! caution: for gradients (sx and sy) the sign changes 266 CALL lbc_lnk_multi( 'icedyn_adv_pra', z0ice , 'T', 1._wp, sxice , 'T', -1._wp, syice , 'T', -1._wp & ! ice volume 267 & , sxxice, 'T', 1._wp, syyice, 'T', 1._wp, sxyice, 'T', 1._wp & 268 & , z0snw , 'T', 1._wp, sxsn , 'T', -1._wp, sysn , 'T', -1._wp & ! snw volume 269 & , sxxsn , 'T', 1._wp, syysn , 'T', 1._wp, sxysn , 'T', 1._wp ) 270 CALL lbc_lnk_multi( 'icedyn_adv_pra', z0smi , 'T', 1._wp, sxsal , 'T', -1._wp, sysal , 'T', -1._wp & ! ice salinity 271 & , sxxsal, 'T', 1._wp, syysal, 'T', 1._wp, sxysal, 'T', 1._wp & 272 & , z0ai , 'T', 1._wp, sxa , 'T', -1._wp, sya , 'T', -1._wp & ! ice concentration 273 & , sxxa , 'T', 1._wp, syya , 'T', 1._wp, sxya , 'T', 1._wp ) 274 CALL lbc_lnk_multi( 'icedyn_adv_pra', z0oi , 'T', 1._wp, sxage , 'T', -1._wp, syage , 'T', -1._wp & ! ice age 275 & , sxxage, 'T', 1._wp, syyage, 'T', 1._wp, sxyage, 'T', 1._wp ) 276 CALL lbc_lnk_multi( 'icedyn_adv_pra', z0es , 'T', 1._wp, sxc0 , 'T', -1._wp, syc0 , 'T', -1._wp & ! snw enthalpy 277 & , sxxc0 , 'T', 1._wp, syyc0 , 'T', 1._wp, sxyc0 , 'T', 1._wp ) 278 CALL lbc_lnk_multi( 'icedyn_adv_pra', z0ei , 'T', 1._wp, sxe , 'T', -1._wp, sye , 'T', -1._wp & ! ice enthalpy 279 & , sxxe , 'T', 1._wp, syye , 'T', 1._wp, sxye , 'T', 1._wp ) 280 IF ( ln_pnd_LEV ) THEN 281 CALL lbc_lnk_multi( 'icedyn_adv_pra', z0ap , 'T', 1._wp, sxap , 'T', -1._wp, syap , 'T', -1._wp & ! melt pond fraction 282 & , sxxap, 'T', 1._wp, syyap, 'T', 1._wp, sxyap, 'T', 1._wp & 283 & , z0vp , 'T', 1._wp, sxvp , 'T', -1._wp, syvp , 'T', -1._wp & ! melt pond volume 284 & , sxxvp, 'T', 1._wp, syyvp, 'T', 1._wp, sxyvp, 'T', 1._wp ) 285 IF ( ln_pnd_lids ) THEN 286 CALL lbc_lnk_multi( 'icedyn_adv_pra', z0vl ,'T', 1._wp, sxvl ,'T', -1._wp, syvl ,'T', -1._wp & ! melt pond lid volume 287 & , sxxvl,'T', 1._wp, syyvl,'T', 1._wp, sxyvl,'T', 1._wp ) 288 ENDIF 283 289 ENDIF 284 290 … … 312 318 END_2D 313 319 CALL lbc_lnk( 'icedyn_adv_pra', pato_i, 'T', 1.0_wp ) 320 ! 321 ! --- diagnostics --- ! 322 diag_adv_mass(:,:) = diag_adv_mass(:,:) + ( SUM( pv_i(:,:,:) , dim=3 ) * rhoi + SUM( pv_s(:,:,:) , dim=3 ) * rhos & 323 & - zdiag_adv_mass(:,:) ) * z1_dt 324 diag_adv_salt(:,:) = diag_adv_salt(:,:) + ( SUM( psv_i(:,:,:) , dim=3 ) * rhoi & 325 & - zdiag_adv_salt(:,:) ) * z1_dt 326 diag_adv_heat(:,:) = diag_adv_heat(:,:) + ( - SUM(SUM( pe_i(:,:,1:nlay_i,:) , dim=4 ), dim=3 ) & 327 & - SUM(SUM( pe_s(:,:,1:nlay_s,:) , dim=4 ), dim=3 ) & 328 & - zdiag_adv_heat(:,:) ) * z1_dt 314 329 ! 315 330 ! --- Ensure non-negative fields --- ! … … 350 365 !! 351 366 INTEGER :: ji, jj, jl, jcat ! dummy loop indices 367 INTEGER :: jj0 ! dummy loop indices 352 368 REAL(wp) :: zs1max, zslpmax, ztemp ! local scalars 353 369 REAL(wp) :: zs1new, zalf , zalfq , zbt ! - - 354 370 REAL(wp) :: zs2new, zalf1, zalf1q, zbt1 ! - - 371 REAL(wp) :: zpsm, zps0 372 REAL(wp) :: zpsx, zpsy, zpsxx, zpsyy, zpsxy 355 373 REAL(wp), DIMENSION(jpi,jpj) :: zf0 , zfx , zfy , zbet ! 2D workspace 356 374 REAL(wp), DIMENSION(jpi,jpj) :: zfm , zfxx , zfyy , zfxy ! - - 357 375 REAL(wp), DIMENSION(jpi,jpj) :: zalg, zalg1, zalg1q ! - - 358 376 !----------------------------------------------------------------------- 377 ! in order to avoid lbc_lnk (communications): 378 ! jj loop must be 1:jpj if adv_x is called first 379 ! and 2:jpj-1 if adv_x is called second 380 jj0 = NINT(pcrh) 359 381 ! 360 382 jcat = SIZE( ps0 , 3 ) ! size of input arrays … … 363 385 ! 364 386 ! Limitation of moments. 365 DO_2D( 0, 0, 1, 1 ) 366 ! Initialize volumes of boxes (=area if adv_x first called, =psm otherwise) 367 psm (ji,jj,jl) = MAX( pcrh * e1e2t(ji,jj) + ( 1.0 - pcrh ) * psm(ji,jj,jl) , epsi20 ) 368 ! 369 zslpmax = MAX( 0._wp, ps0(ji,jj,jl) ) 370 zs1max = 1.5 * zslpmax 371 zs1new = MIN( zs1max, MAX( -zs1max, psx(ji,jj,jl) ) ) 372 zs2new = MIN( 2.0 * zslpmax - 0.3334 * ABS( zs1new ), & 373 & MAX( ABS( zs1new ) - zslpmax, psxx(ji,jj,jl) ) ) 374 rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zslpmax) ) ) * tmask(ji,jj,1) ! Case of empty boxes & Apply mask 375 376 ps0 (ji,jj,jl) = zslpmax 377 psx (ji,jj,jl) = zs1new * rswitch 378 psxx(ji,jj,jl) = zs2new * rswitch 379 psy (ji,jj,jl) = psy (ji,jj,jl) * rswitch 380 psyy(ji,jj,jl) = psyy(ji,jj,jl) * rswitch 381 psxy(ji,jj,jl) = MIN( zslpmax, MAX( -zslpmax, psxy(ji,jj,jl) ) ) * rswitch 382 END_2D 383 384 ! Calculate fluxes and moments between boxes i<-->i+1 385 DO_2D( 0, 0, 1, 1 ) ! Flux from i to i+1 WHEN u GT 0 386 zbet(ji,jj) = MAX( 0._wp, SIGN( 1._wp, put(ji,jj) ) ) 387 zalf = MAX( 0._wp, put(ji,jj) ) * pdt / psm(ji,jj,jl) 388 zalfq = zalf * zalf 389 zalf1 = 1.0 - zalf 390 zalf1q = zalf1 * zalf1 391 ! 392 zfm (ji,jj) = zalf * psm (ji,jj,jl) 393 zf0 (ji,jj) = zalf * ( ps0 (ji,jj,jl) + zalf1 * ( psx(ji,jj,jl) + (zalf1 - zalf) * psxx(ji,jj,jl) ) ) 394 zfx (ji,jj) = zalfq * ( psx (ji,jj,jl) + 3.0 * zalf1 * psxx(ji,jj,jl) ) 395 zfxx(ji,jj) = zalf * psxx(ji,jj,jl) * zalfq 396 zfy (ji,jj) = zalf * ( psy (ji,jj,jl) + zalf1 * psxy(ji,jj,jl) ) 397 zfxy(ji,jj) = zalfq * psxy(ji,jj,jl) 398 zfyy(ji,jj) = zalf * psyy(ji,jj,jl) 399 400 ! Readjust moments remaining in the box. 401 psm (ji,jj,jl) = psm (ji,jj,jl) - zfm(ji,jj) 402 ps0 (ji,jj,jl) = ps0 (ji,jj,jl) - zf0(ji,jj) 403 psx (ji,jj,jl) = zalf1q * ( psx(ji,jj,jl) - 3.0 * zalf * psxx(ji,jj,jl) ) 404 psxx(ji,jj,jl) = zalf1 * zalf1q * psxx(ji,jj,jl) 405 psy (ji,jj,jl) = psy (ji,jj,jl) - zfy(ji,jj) 406 psyy(ji,jj,jl) = psyy(ji,jj,jl) - zfyy(ji,jj) 407 psxy(ji,jj,jl) = zalf1q * psxy(ji,jj,jl) 408 END_2D 409 410 DO_2D( 0, 0, 1, 0 ) ! Flux from i+1 to i when u LT 0. 411 zalf = MAX( 0._wp, -put(ji,jj) ) * pdt / psm(ji+1,jj,jl) 412 zalg (ji,jj) = zalf 413 zalfq = zalf * zalf 414 zalf1 = 1.0 - zalf 415 zalg1 (ji,jj) = zalf1 416 zalf1q = zalf1 * zalf1 417 zalg1q(ji,jj) = zalf1q 418 ! 419 zfm (ji,jj) = zfm (ji,jj) + zalf * psm (ji+1,jj,jl) 420 zf0 (ji,jj) = zf0 (ji,jj) + zalf * ( ps0 (ji+1,jj,jl) & 421 & - zalf1 * ( psx(ji+1,jj,jl) - (zalf1 - zalf ) * psxx(ji+1,jj,jl) ) ) 422 zfx (ji,jj) = zfx (ji,jj) + zalfq * ( psx (ji+1,jj,jl) - 3.0 * zalf1 * psxx(ji+1,jj,jl) ) 423 zfxx (ji,jj) = zfxx(ji,jj) + zalf * psxx(ji+1,jj,jl) * zalfq 424 zfy (ji,jj) = zfy (ji,jj) + zalf * ( psy (ji+1,jj,jl) - zalf1 * psxy(ji+1,jj,jl) ) 425 zfxy (ji,jj) = zfxy(ji,jj) + zalfq * psxy(ji+1,jj,jl) 426 zfyy (ji,jj) = zfyy(ji,jj) + zalf * psyy(ji+1,jj,jl) 427 END_2D 428 429 DO_2D( 0, 0, 0, 0 ) ! Readjust moments remaining in the box. 430 zbt = zbet(ji-1,jj) 431 zbt1 = 1.0 - zbet(ji-1,jj) 432 ! 433 psm (ji,jj,jl) = zbt * psm(ji,jj,jl) + zbt1 * ( psm(ji,jj,jl) - zfm(ji-1,jj) ) 434 ps0 (ji,jj,jl) = zbt * ps0(ji,jj,jl) + zbt1 * ( ps0(ji,jj,jl) - zf0(ji-1,jj) ) 435 psx (ji,jj,jl) = zalg1q(ji-1,jj) * ( psx(ji,jj,jl) + 3.0 * zalg(ji-1,jj) * psxx(ji,jj,jl) ) 436 psxx(ji,jj,jl) = zalg1 (ji-1,jj) * zalg1q(ji-1,jj) * psxx(ji,jj,jl) 437 psy (ji,jj,jl) = zbt * psy (ji,jj,jl) + zbt1 * ( psy (ji,jj,jl) - zfy (ji-1,jj) ) 438 psyy(ji,jj,jl) = zbt * psyy(ji,jj,jl) + zbt1 * ( psyy(ji,jj,jl) - zfyy(ji-1,jj) ) 439 psxy(ji,jj,jl) = zalg1q(ji-1,jj) * psxy(ji,jj,jl) 440 END_2D 441 442 ! Put the temporary moments into appropriate neighboring boxes. 443 DO_2D( 0, 0, 0, 0 ) ! Flux from i to i+1 IF u GT 0. 444 zbt = zbet(ji-1,jj) 445 zbt1 = 1.0 - zbet(ji-1,jj) 446 psm(ji,jj,jl) = zbt * ( psm(ji,jj,jl) + zfm(ji-1,jj) ) + zbt1 * psm(ji,jj,jl) 447 zalf = zbt * zfm(ji-1,jj) / psm(ji,jj,jl) 448 zalf1 = 1.0 - zalf 449 ztemp = zalf * ps0(ji,jj,jl) - zalf1 * zf0(ji-1,jj) 450 ! 451 ps0 (ji,jj,jl) = zbt * ( ps0(ji,jj,jl) + zf0(ji-1,jj) ) + zbt1 * ps0(ji,jj,jl) 452 psx (ji,jj,jl) = zbt * ( zalf * zfx(ji-1,jj) + zalf1 * psx(ji,jj,jl) + 3.0 * ztemp ) + zbt1 * psx(ji,jj,jl) 453 psxx(ji,jj,jl) = zbt * ( zalf * zalf * zfxx(ji-1,jj) + zalf1 * zalf1 * psxx(ji,jj,jl) & 454 & + 5.0 * ( zalf * zalf1 * ( psx (ji,jj,jl) - zfx(ji-1,jj) ) - ( zalf1 - zalf ) * ztemp ) ) & 455 & + zbt1 * psxx(ji,jj,jl) 456 psxy(ji,jj,jl) = zbt * ( zalf * zfxy(ji-1,jj) + zalf1 * psxy(ji,jj,jl) & 457 & + 3.0 * (- zalf1*zfy(ji-1,jj) + zalf * psy(ji,jj,jl) ) ) & 458 & + zbt1 * psxy(ji,jj,jl) 459 psy (ji,jj,jl) = zbt * ( psy (ji,jj,jl) + zfy (ji-1,jj) ) + zbt1 * psy (ji,jj,jl) 460 psyy(ji,jj,jl) = zbt * ( psyy(ji,jj,jl) + zfyy(ji-1,jj) ) + zbt1 * psyy(ji,jj,jl) 461 END_2D 462 463 DO_2D( 0, 0, 0, 0 ) ! Flux from i+1 to i IF u LT 0. 464 zbt = zbet(ji,jj) 465 zbt1 = 1.0 - zbet(ji,jj) 466 psm(ji,jj,jl) = zbt * psm(ji,jj,jl) + zbt1 * ( psm(ji,jj,jl) + zfm(ji,jj) ) 467 zalf = zbt1 * zfm(ji,jj) / psm(ji,jj,jl) 468 zalf1 = 1.0 - zalf 469 ztemp = - zalf * ps0(ji,jj,jl) + zalf1 * zf0(ji,jj) 470 ! 471 ps0 (ji,jj,jl) = zbt * ps0 (ji,jj,jl) + zbt1 * ( ps0(ji,jj,jl) + zf0(ji,jj) ) 472 psx (ji,jj,jl) = zbt * psx (ji,jj,jl) + zbt1 * ( zalf * zfx(ji,jj) + zalf1 * psx(ji,jj,jl) + 3.0 * ztemp ) 473 psxx(ji,jj,jl) = zbt * psxx(ji,jj,jl) + zbt1 * ( zalf * zalf * zfxx(ji,jj) + zalf1 * zalf1 * psxx(ji,jj,jl) & 474 & + 5.0 * ( zalf * zalf1 * ( - psx(ji,jj,jl) + zfx(ji,jj) ) & 475 & + ( zalf1 - zalf ) * ztemp ) ) 476 psxy(ji,jj,jl) = zbt * psxy(ji,jj,jl) + zbt1 * ( zalf * zfxy(ji,jj) + zalf1 * psxy(ji,jj,jl) & 477 & + 3.0 * ( zalf1 * zfy(ji,jj) - zalf * psy(ji,jj,jl) ) ) 478 psy (ji,jj,jl) = zbt * psy (ji,jj,jl) + zbt1 * ( psy (ji,jj,jl) + zfy (ji,jj) ) 479 psyy(ji,jj,jl) = zbt * psyy(ji,jj,jl) + zbt1 * ( psyy(ji,jj,jl) + zfyy(ji,jj) ) 480 END_2D 481 387 DO jj = Njs0 - jj0, Nje0 + jj0 388 389 DO ji = Nis0 - 1, Nie0 + 1 390 391 zpsm = psm (ji,jj,jl) ! optimization 392 zps0 = ps0 (ji,jj,jl) 393 zpsx = psx (ji,jj,jl) 394 zpsxx = psxx(ji,jj,jl) 395 zpsy = psy (ji,jj,jl) 396 zpsyy = psyy(ji,jj,jl) 397 zpsxy = psxy(ji,jj,jl) 398 399 ! Initialize volumes of boxes (=area if adv_x first called, =psm otherwise) 400 zpsm = MAX( pcrh * e1e2t(ji,jj) + ( 1.0 - pcrh ) * zpsm , epsi20 ) 401 ! 402 zslpmax = MAX( 0._wp, zps0 ) 403 zs1max = 1.5 * zslpmax 404 zs1new = MIN( zs1max, MAX( -zs1max, zpsx ) ) 405 zs2new = MIN( 2.0 * zslpmax - 0.3334 * ABS( zs1new ), MAX( ABS( zs1new ) - zslpmax, zpsxx ) ) 406 rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zslpmax) ) ) * tmask(ji,jj,1) ! Case of empty boxes & Apply mask 407 408 zps0 = zslpmax 409 zpsx = zs1new * rswitch 410 zpsxx = zs2new * rswitch 411 zpsy = zpsy * rswitch 412 zpsyy = zpsyy * rswitch 413 zpsxy = MIN( zslpmax, MAX( -zslpmax, zpsxy ) ) * rswitch 414 415 ! Calculate fluxes and moments between boxes i<-->i+1 416 ! ! Flux from i to i+1 WHEN u GT 0 417 zbet(ji,jj) = MAX( 0._wp, SIGN( 1._wp, put(ji,jj) ) ) 418 zalf = MAX( 0._wp, put(ji,jj) ) * pdt / zpsm 419 zalfq = zalf * zalf 420 zalf1 = 1.0 - zalf 421 zalf1q = zalf1 * zalf1 422 ! 423 zfm (ji,jj) = zalf * zpsm 424 zf0 (ji,jj) = zalf * ( zps0 + zalf1 * ( zpsx + (zalf1 - zalf) * zpsxx ) ) 425 zfx (ji,jj) = zalfq * ( zpsx + 3.0 * zalf1 * zpsxx ) 426 zfxx(ji,jj) = zalf * zpsxx * zalfq 427 zfy (ji,jj) = zalf * ( zpsy + zalf1 * zpsxy ) 428 zfxy(ji,jj) = zalfq * zpsxy 429 zfyy(ji,jj) = zalf * zpsyy 430 431 ! ! Readjust moments remaining in the box. 432 zpsm = zpsm - zfm(ji,jj) 433 zps0 = zps0 - zf0(ji,jj) 434 zpsx = zalf1q * ( zpsx - 3.0 * zalf * zpsxx ) 435 zpsxx = zalf1 * zalf1q * zpsxx 436 zpsy = zpsy - zfy (ji,jj) 437 zpsyy = zpsyy - zfyy(ji,jj) 438 zpsxy = zalf1q * zpsxy 439 ! 440 psm (ji,jj,jl) = zpsm ! optimization 441 ps0 (ji,jj,jl) = zps0 442 psx (ji,jj,jl) = zpsx 443 psxx(ji,jj,jl) = zpsxx 444 psy (ji,jj,jl) = zpsy 445 psyy(ji,jj,jl) = zpsyy 446 psxy(ji,jj,jl) = zpsxy 447 ! 448 END DO 449 450 DO ji = Nis0 - 1, Nie0 451 ! ! Flux from i+1 to i when u LT 0. 452 zalf = MAX( 0._wp, -put(ji,jj) ) * pdt / psm(ji+1,jj,jl) 453 zalg (ji,jj) = zalf 454 zalfq = zalf * zalf 455 zalf1 = 1.0 - zalf 456 zalg1 (ji,jj) = zalf1 457 zalf1q = zalf1 * zalf1 458 zalg1q(ji,jj) = zalf1q 459 ! 460 zfm (ji,jj) = zfm (ji,jj) + zalf * psm (ji+1,jj,jl) 461 zf0 (ji,jj) = zf0 (ji,jj) + zalf * ( ps0 (ji+1,jj,jl) & 462 & - zalf1 * ( psx(ji+1,jj,jl) - (zalf1 - zalf ) * psxx(ji+1,jj,jl) ) ) 463 zfx (ji,jj) = zfx (ji,jj) + zalfq * ( psx (ji+1,jj,jl) - 3.0 * zalf1 * psxx(ji+1,jj,jl) ) 464 zfxx (ji,jj) = zfxx(ji,jj) + zalf * psxx(ji+1,jj,jl) * zalfq 465 zfy (ji,jj) = zfy (ji,jj) + zalf * ( psy (ji+1,jj,jl) - zalf1 * psxy(ji+1,jj,jl) ) 466 zfxy (ji,jj) = zfxy(ji,jj) + zalfq * psxy(ji+1,jj,jl) 467 zfyy (ji,jj) = zfyy(ji,jj) + zalf * psyy(ji+1,jj,jl) 468 END DO 469 470 DO ji = Nis0, Nie0 471 ! 472 zpsm = psm (ji,jj,jl) ! optimization 473 zps0 = ps0 (ji,jj,jl) 474 zpsx = psx (ji,jj,jl) 475 zpsxx = psxx(ji,jj,jl) 476 zpsy = psy (ji,jj,jl) 477 zpsyy = psyy(ji,jj,jl) 478 zpsxy = psxy(ji,jj,jl) 479 ! ! Readjust moments remaining in the box. 480 zbt = zbet(ji-1,jj) 481 zbt1 = 1.0 - zbet(ji-1,jj) 482 ! 483 zpsm = zbt * zpsm + zbt1 * ( zpsm - zfm(ji-1,jj) ) 484 zps0 = zbt * zps0 + zbt1 * ( zps0 - zf0(ji-1,jj) ) 485 zpsx = zalg1q(ji-1,jj) * ( zpsx + 3.0 * zalg(ji-1,jj) * zpsxx ) 486 zpsxx = zalg1 (ji-1,jj) * zalg1q(ji-1,jj) * zpsxx 487 zpsy = zbt * zpsy + zbt1 * ( zpsy - zfy (ji-1,jj) ) 488 zpsyy = zbt * zpsyy + zbt1 * ( zpsyy - zfyy(ji-1,jj) ) 489 zpsxy = zalg1q(ji-1,jj) * zpsxy 490 491 ! Put the temporary moments into appropriate neighboring boxes. 492 ! ! Flux from i to i+1 IF u GT 0. 493 zbt = zbet(ji-1,jj) 494 zbt1 = 1.0 - zbet(ji-1,jj) 495 zpsm = zbt * ( zpsm + zfm(ji-1,jj) ) + zbt1 * zpsm 496 zalf = zbt * zfm(ji-1,jj) / zpsm 497 zalf1 = 1.0 - zalf 498 ztemp = zalf * zps0 - zalf1 * zf0(ji-1,jj) 499 ! 500 zps0 = zbt * ( zps0 + zf0(ji-1,jj) ) + zbt1 * zps0 501 zpsx = zbt * ( zalf * zfx(ji-1,jj) + zalf1 * zpsx + 3.0 * ztemp ) + zbt1 * zpsx 502 zpsxx = zbt * ( zalf * zalf * zfxx(ji-1,jj) + zalf1 * zalf1 * zpsxx & 503 & + 5.0 * ( zalf * zalf1 * ( zpsx - zfx(ji-1,jj) ) - ( zalf1 - zalf ) * ztemp ) ) & 504 & + zbt1 * zpsxx 505 zpsxy = zbt * ( zalf * zfxy(ji-1,jj) + zalf1 * zpsxy & 506 & + 3.0 * (- zalf1*zfy(ji-1,jj) + zalf * zpsy ) ) & 507 & + zbt1 * zpsxy 508 zpsy = zbt * ( zpsy + zfy (ji-1,jj) ) + zbt1 * zpsy 509 zpsyy = zbt * ( zpsyy + zfyy(ji-1,jj) ) + zbt1 * zpsyy 510 511 ! ! Flux from i+1 to i IF u LT 0. 512 zbt = zbet(ji,jj) 513 zbt1 = 1.0 - zbet(ji,jj) 514 zpsm = zbt * zpsm + zbt1 * ( zpsm + zfm(ji,jj) ) 515 zalf = zbt1 * zfm(ji,jj) / zpsm 516 zalf1 = 1.0 - zalf 517 ztemp = - zalf * zps0 + zalf1 * zf0(ji,jj) 518 ! 519 zps0 = zbt * zps0 + zbt1 * ( zps0 + zf0(ji,jj) ) 520 zpsx = zbt * zpsx + zbt1 * ( zalf * zfx(ji,jj) + zalf1 * zpsx + 3.0 * ztemp ) 521 zpsxx = zbt * zpsxx + zbt1 * ( zalf * zalf * zfxx(ji,jj) + zalf1 * zalf1 * zpsxx & 522 & + 5.0 * ( zalf * zalf1 * ( - zpsx + zfx(ji,jj) ) & 523 & + ( zalf1 - zalf ) * ztemp ) ) 524 zpsxy = zbt * zpsxy + zbt1 * ( zalf * zfxy(ji,jj) + zalf1 * zpsxy & 525 & + 3.0 * ( zalf1 * zfy(ji,jj) - zalf * zpsy ) ) 526 zpsy = zbt * zpsy + zbt1 * ( zpsy + zfy (ji,jj) ) 527 zpsyy = zbt * zpsyy + zbt1 * ( zpsyy + zfyy(ji,jj) ) 528 ! 529 psm (ji,jj,jl) = zpsm ! optimization 530 ps0 (ji,jj,jl) = zps0 531 psx (ji,jj,jl) = zpsx 532 psxx(ji,jj,jl) = zpsxx 533 psy (ji,jj,jl) = zpsy 534 psyy(ji,jj,jl) = zpsyy 535 psxy(ji,jj,jl) = zpsxy 536 END DO 537 ! 538 END DO 539 ! 482 540 END DO 483 484 !-- Lateral boundary conditions 485 CALL lbc_lnk_multi( 'icedyn_adv_pra', psm(:,:,1:jcat) , 'T', 1.0_wp, ps0 , 'T', 1.0_wp & 486 & , psx , 'T', -1.0_wp, psy , 'T', -1.0_wp & ! caution gradient ==> the sign changes 487 & , psxx , 'T', 1.0_wp, psyy, 'T', 1.0_wp , psxy, 'T', 1.0_wp ) 488 ! 541 ! 489 542 END SUBROUTINE adv_x 490 543 … … 507 560 !! 508 561 INTEGER :: ji, jj, jl, jcat ! dummy loop indices 562 INTEGER :: ji0 ! dummy loop indices 509 563 REAL(wp) :: zs1max, zslpmax, ztemp ! temporary scalars 510 564 REAL(wp) :: zs1new, zalf , zalfq , zbt ! - - 511 565 REAL(wp) :: zs2new, zalf1, zalf1q, zbt1 ! - - 566 REAL(wp) :: zpsm, zps0 567 REAL(wp) :: zpsx, zpsy, zpsxx, zpsyy, zpsxy 512 568 REAL(wp), DIMENSION(jpi,jpj) :: zf0, zfx , zfy , zbet ! 2D workspace 513 569 REAL(wp), DIMENSION(jpi,jpj) :: zfm, zfxx, zfyy, zfxy ! - - 514 570 REAL(wp), DIMENSION(jpi,jpj) :: zalg, zalg1, zalg1q ! - - 515 571 !--------------------------------------------------------------------- 572 ! in order to avoid lbc_lnk (communications): 573 ! ji loop must be 1:jpi if adv_y is called first 574 ! and 2:jpi-1 if adv_y is called second 575 ji0 = NINT(pcrh) 516 576 ! 517 577 jcat = SIZE( ps0 , 3 ) ! size of input arrays … … 520 580 ! 521 581 ! Limitation of moments. 522 DO_2D( 1, 1, 0, 0 ) 523 ! Initialize volumes of boxes (=area if adv_x first called, =psm otherwise) 524 psm(ji,jj,jl) = MAX( pcrh * e1e2t(ji,jj) + ( 1.0 - pcrh ) * psm(ji,jj,jl) , epsi20 ) 525 ! 526 zslpmax = MAX( 0._wp, ps0(ji,jj,jl) ) 582 DO_2D( 1, 1, ji0, ji0 ) 583 ! 584 zpsm = psm (ji,jj,jl) ! optimization 585 zps0 = ps0 (ji,jj,jl) 586 zpsx = psx (ji,jj,jl) 587 zpsxx = psxx(ji,jj,jl) 588 zpsy = psy (ji,jj,jl) 589 zpsyy = psyy(ji,jj,jl) 590 zpsxy = psxy(ji,jj,jl) 591 ! 592 ! Initialize volumes of boxes (=area if adv_y first called, =psm otherwise) 593 zpsm = MAX( pcrh * e1e2t(ji,jj) + ( 1.0 - pcrh ) * zpsm , epsi20 ) 594 ! 595 zslpmax = MAX( 0._wp, zps0 ) 527 596 zs1max = 1.5 * zslpmax 528 zs1new = MIN( zs1max, MAX( -zs1max, psy(ji,jj,jl) ) ) 529 zs2new = MIN( ( 2.0 * zslpmax - 0.3334 * ABS( zs1new ) ), & 530 & MAX( ABS( zs1new )-zslpmax, psyy(ji,jj,jl) ) ) 597 zs1new = MIN( zs1max, MAX( -zs1max, zpsy ) ) 598 zs2new = MIN( ( 2.0 * zslpmax - 0.3334 * ABS( zs1new ) ), MAX( ABS( zs1new )-zslpmax, zpsyy ) ) 531 599 rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zslpmax) ) ) * tmask(ji,jj,1) ! Case of empty boxes & Apply mask 532 600 ! 533 ps0 (ji,jj,jl) = zslpmax 534 psx (ji,jj,jl) = psx (ji,jj,jl) * rswitch 535 psxx(ji,jj,jl) = psxx(ji,jj,jl) * rswitch 536 psy (ji,jj,jl) = zs1new * rswitch 537 psyy(ji,jj,jl) = zs2new * rswitch 538 psxy(ji,jj,jl) = MIN( zslpmax, MAX( -zslpmax, psxy(ji,jj,jl) ) ) * rswitch 539 END_2D 540 541 ! Calculate fluxes and moments between boxes j<-->j+1 542 DO_2D( 1, 1, 0, 0 ) ! Flux from j to j+1 WHEN v GT 0 601 zps0 = zslpmax 602 zpsx = zpsx * rswitch 603 zpsxx = zpsxx * rswitch 604 zpsy = zs1new * rswitch 605 zpsyy = zs2new * rswitch 606 zpsxy = MIN( zslpmax, MAX( -zslpmax, zpsxy ) ) * rswitch 607 608 ! Calculate fluxes and moments between boxes j<-->j+1 609 ! ! Flux from j to j+1 WHEN v GT 0 543 610 zbet(ji,jj) = MAX( 0._wp, SIGN( 1._wp, pvt(ji,jj) ) ) 544 zalf = MAX( 0._wp, pvt(ji,jj) ) * pdt / psm(ji,jj,jl)611 zalf = MAX( 0._wp, pvt(ji,jj) ) * pdt / zpsm 545 612 zalfq = zalf * zalf 546 613 zalf1 = 1.0 - zalf 547 614 zalf1q = zalf1 * zalf1 548 615 ! 549 zfm (ji,jj) = zalf * psm(ji,jj,jl) 550 zf0 (ji,jj) = zalf * ( ps0(ji,jj,jl) + zalf1 * ( psy(ji,jj,jl) + (zalf1-zalf) * psyy(ji,jj,jl) ) ) 551 zfy (ji,jj) = zalfq *( psy(ji,jj,jl) + 3.0*zalf1*psyy(ji,jj,jl) ) 552 zfyy(ji,jj) = zalf * zalfq * psyy(ji,jj,jl) 553 zfx (ji,jj) = zalf * ( psx(ji,jj,jl) + zalf1 * psxy(ji,jj,jl) ) 554 zfxy(ji,jj) = zalfq * psxy(ji,jj,jl) 555 zfxx(ji,jj) = zalf * psxx(ji,jj,jl) 556 ! 557 ! Readjust moments remaining in the box. 558 psm (ji,jj,jl) = psm (ji,jj,jl) - zfm(ji,jj) 559 ps0 (ji,jj,jl) = ps0 (ji,jj,jl) - zf0(ji,jj) 560 psy (ji,jj,jl) = zalf1q * ( psy(ji,jj,jl) -3.0 * zalf * psyy(ji,jj,jl) ) 561 psyy(ji,jj,jl) = zalf1 * zalf1q * psyy(ji,jj,jl) 562 psx (ji,jj,jl) = psx (ji,jj,jl) - zfx(ji,jj) 563 psxx(ji,jj,jl) = psxx(ji,jj,jl) - zfxx(ji,jj) 564 psxy(ji,jj,jl) = zalf1q * psxy(ji,jj,jl) 616 zfm (ji,jj) = zalf * zpsm 617 zf0 (ji,jj) = zalf * ( zps0 + zalf1 * ( zpsy + (zalf1-zalf) * zpsyy ) ) 618 zfy (ji,jj) = zalfq *( zpsy + 3.0*zalf1*zpsyy ) 619 zfyy(ji,jj) = zalf * zalfq * zpsyy 620 zfx (ji,jj) = zalf * ( zpsx + zalf1 * zpsxy ) 621 zfxy(ji,jj) = zalfq * zpsxy 622 zfxx(ji,jj) = zalf * zpsxx 623 ! 624 ! ! Readjust moments remaining in the box. 625 zpsm = zpsm - zfm(ji,jj) 626 zps0 = zps0 - zf0(ji,jj) 627 zpsy = zalf1q * ( zpsy -3.0 * zalf * zpsyy ) 628 zpsyy = zalf1 * zalf1q * zpsyy 629 zpsx = zpsx - zfx(ji,jj) 630 zpsxx = zpsxx - zfxx(ji,jj) 631 zpsxy = zalf1q * zpsxy 632 ! 633 psm (ji,jj,jl) = zpsm ! optimization 634 ps0 (ji,jj,jl) = zps0 635 psx (ji,jj,jl) = zpsx 636 psxx(ji,jj,jl) = zpsxx 637 psy (ji,jj,jl) = zpsy 638 psyy(ji,jj,jl) = zpsyy 639 psxy(ji,jj,jl) = zpsxy 565 640 END_2D 566 641 ! 567 DO_2D( 1, 0, 0, 0 ) ! Flux from j+1 to j when v LT 0. 642 DO_2D( 1, 0, ji0, ji0 ) 643 ! ! Flux from j+1 to j when v LT 0. 568 644 zalf = MAX( 0._wp, -pvt(ji,jj) ) * pdt / psm(ji,jj+1,jl) 569 645 zalg (ji,jj) = zalf … … 584 660 END_2D 585 661 586 ! Readjust moments remaining in the box.587 DO_2D( 0, 0, 0, 0 )662 DO_2D( 0, 0, ji0, ji0 ) 663 ! ! Readjust moments remaining in the box. 588 664 zbt = zbet(ji,jj-1) 589 665 zbt1 = ( 1.0 - zbet(ji,jj-1) ) 590 666 ! 591 psm (ji,jj,jl) = zbt * psm(ji,jj,jl) + zbt1 * ( psm(ji,jj,jl) - zfm(ji,jj-1) ) 592 ps0 (ji,jj,jl) = zbt * ps0(ji,jj,jl) + zbt1 * ( ps0(ji,jj,jl) - zf0(ji,jj-1) ) 593 psy (ji,jj,jl) = zalg1q(ji,jj-1) * ( psy(ji,jj,jl) + 3.0 * zalg(ji,jj-1) * psyy(ji,jj,jl) ) 594 psyy(ji,jj,jl) = zalg1 (ji,jj-1) * zalg1q(ji,jj-1) * psyy(ji,jj,jl) 595 psx (ji,jj,jl) = zbt * psx (ji,jj,jl) + zbt1 * ( psx (ji,jj,jl) - zfx (ji,jj-1) ) 596 psxx(ji,jj,jl) = zbt * psxx(ji,jj,jl) + zbt1 * ( psxx(ji,jj,jl) - zfxx(ji,jj-1) ) 597 psxy(ji,jj,jl) = zalg1q(ji,jj-1) * psxy(ji,jj,jl) 667 zpsm = psm (ji,jj,jl) ! optimization 668 zps0 = ps0 (ji,jj,jl) 669 zpsx = psx (ji,jj,jl) 670 zpsxx = psxx(ji,jj,jl) 671 zpsy = psy (ji,jj,jl) 672 zpsyy = psyy(ji,jj,jl) 673 zpsxy = psxy(ji,jj,jl) 674 ! 675 zpsm = zbt * zpsm + zbt1 * ( zpsm - zfm(ji,jj-1) ) 676 zps0 = zbt * zps0 + zbt1 * ( zps0 - zf0(ji,jj-1) ) 677 zpsy = zalg1q(ji,jj-1) * ( zpsy + 3.0 * zalg(ji,jj-1) * zpsyy ) 678 zpsyy = zalg1 (ji,jj-1) * zalg1q(ji,jj-1) * zpsyy 679 zpsx = zbt * zpsx + zbt1 * ( zpsx - zfx (ji,jj-1) ) 680 zpsxx = zbt * zpsxx + zbt1 * ( zpsxx - zfxx(ji,jj-1) ) 681 zpsxy = zalg1q(ji,jj-1) * zpsxy 682 683 ! Put the temporary moments into appropriate neighboring boxes. 684 ! ! Flux from j to j+1 IF v GT 0. 685 zbt = zbet(ji,jj-1) 686 zbt1 = 1.0 - zbet(ji,jj-1) 687 zpsm = zbt * ( zpsm + zfm(ji,jj-1) ) + zbt1 * zpsm 688 zalf = zbt * zfm(ji,jj-1) / zpsm 689 zalf1 = 1.0 - zalf 690 ztemp = zalf * zps0 - zalf1 * zf0(ji,jj-1) 691 ! 692 zps0 = zbt * ( zps0 + zf0(ji,jj-1) ) + zbt1 * zps0 693 zpsy = zbt * ( zalf * zfy(ji,jj-1) + zalf1 * zpsy + 3.0 * ztemp ) & 694 & + zbt1 * zpsy 695 zpsyy = zbt * ( zalf * zalf * zfyy(ji,jj-1) + zalf1 * zalf1 * zpsyy & 696 & + 5.0 * ( zalf * zalf1 * ( zpsy - zfy(ji,jj-1) ) - ( zalf1 - zalf ) * ztemp ) ) & 697 & + zbt1 * zpsyy 698 zpsxy = zbt * ( zalf * zfxy(ji,jj-1) + zalf1 * zpsxy & 699 & + 3.0 * (- zalf1 * zfx(ji,jj-1) + zalf * zpsx ) ) & 700 & + zbt1 * zpsxy 701 zpsx = zbt * ( zpsx + zfx (ji,jj-1) ) + zbt1 * zpsx 702 zpsxx = zbt * ( zpsxx + zfxx(ji,jj-1) ) + zbt1 * zpsxx 703 704 ! ! Flux from j+1 to j IF v LT 0. 705 zbt = zbet(ji,jj) 706 zbt1 = 1.0 - zbet(ji,jj) 707 zpsm = zbt * zpsm + zbt1 * ( zpsm + zfm(ji,jj) ) 708 zalf = zbt1 * zfm(ji,jj) / zpsm 709 zalf1 = 1.0 - zalf 710 ztemp = - zalf * zps0 + zalf1 * zf0(ji,jj) 711 ! 712 zps0 = zbt * zps0 + zbt1 * ( zps0 + zf0(ji,jj) ) 713 zpsy = zbt * zpsy + zbt1 * ( zalf * zfy(ji,jj) + zalf1 * zpsy + 3.0 * ztemp ) 714 zpsyy = zbt * zpsyy + zbt1 * ( zalf * zalf * zfyy(ji,jj) + zalf1 * zalf1 * zpsyy & 715 & + 5.0 * ( zalf * zalf1 * ( - zpsy + zfy(ji,jj) ) & 716 & + ( zalf1 - zalf ) * ztemp ) ) 717 zpsxy = zbt * zpsxy + zbt1 * ( zalf * zfxy(ji,jj) + zalf1 * zpsxy & 718 & + 3.0 * ( zalf1 * zfx(ji,jj) - zalf * zpsx ) ) 719 zpsx = zbt * zpsx + zbt1 * ( zpsx + zfx (ji,jj) ) 720 zpsxx = zbt * zpsxx + zbt1 * ( zpsxx + zfxx(ji,jj) ) 721 ! 722 psm (ji,jj,jl) = zpsm ! optimization 723 ps0 (ji,jj,jl) = zps0 724 psx (ji,jj,jl) = zpsx 725 psxx(ji,jj,jl) = zpsxx 726 psy (ji,jj,jl) = zpsy 727 psyy(ji,jj,jl) = zpsyy 728 psxy(ji,jj,jl) = zpsxy 598 729 END_2D 599 600 ! Put the temporary moments into appropriate neighboring boxes. 601 DO_2D( 0, 0, 0, 0 ) ! Flux from j to j+1 IF v GT 0. 602 zbt = zbet(ji,jj-1) 603 zbt1 = 1.0 - zbet(ji,jj-1) 604 psm(ji,jj,jl) = zbt * ( psm(ji,jj,jl) + zfm(ji,jj-1) ) + zbt1 * psm(ji,jj,jl) 605 zalf = zbt * zfm(ji,jj-1) / psm(ji,jj,jl) 606 zalf1 = 1.0 - zalf 607 ztemp = zalf * ps0(ji,jj,jl) - zalf1 * zf0(ji,jj-1) 608 ! 609 ps0(ji,jj,jl) = zbt * ( ps0(ji,jj,jl) + zf0(ji,jj-1) ) + zbt1 * ps0(ji,jj,jl) 610 psy(ji,jj,jl) = zbt * ( zalf * zfy(ji,jj-1) + zalf1 * psy(ji,jj,jl) + 3.0 * ztemp ) & 611 & + zbt1 * psy(ji,jj,jl) 612 psyy(ji,jj,jl) = zbt * ( zalf * zalf * zfyy(ji,jj-1) + zalf1 * zalf1 * psyy(ji,jj,jl) & 613 & + 5.0 * ( zalf * zalf1 * ( psy(ji,jj,jl) - zfy(ji,jj-1) ) - ( zalf1 - zalf ) * ztemp ) ) & 614 & + zbt1 * psyy(ji,jj,jl) 615 psxy(ji,jj,jl) = zbt * ( zalf * zfxy(ji,jj-1) + zalf1 * psxy(ji,jj,jl) & 616 & + 3.0 * (- zalf1 * zfx(ji,jj-1) + zalf * psx(ji,jj,jl) ) ) & 617 & + zbt1 * psxy(ji,jj,jl) 618 psx (ji,jj,jl) = zbt * ( psx (ji,jj,jl) + zfx (ji,jj-1) ) + zbt1 * psx (ji,jj,jl) 619 psxx(ji,jj,jl) = zbt * ( psxx(ji,jj,jl) + zfxx(ji,jj-1) ) + zbt1 * psxx(ji,jj,jl) 620 END_2D 621 622 DO_2D( 0, 0, 0, 0 ) ! Flux from j+1 to j IF v LT 0. 623 zbt = zbet(ji,jj) 624 zbt1 = 1.0 - zbet(ji,jj) 625 psm(ji,jj,jl) = zbt * psm(ji,jj,jl) + zbt1 * ( psm(ji,jj,jl) + zfm(ji,jj) ) 626 zalf = zbt1 * zfm(ji,jj) / psm(ji,jj,jl) 627 zalf1 = 1.0 - zalf 628 ztemp = - zalf * ps0(ji,jj,jl) + zalf1 * zf0(ji,jj) 629 ! 630 ps0 (ji,jj,jl) = zbt * ps0 (ji,jj,jl) + zbt1 * ( ps0(ji,jj,jl) + zf0(ji,jj) ) 631 psy (ji,jj,jl) = zbt * psy (ji,jj,jl) + zbt1 * ( zalf * zfy(ji,jj) + zalf1 * psy(ji,jj,jl) + 3.0 * ztemp ) 632 psyy(ji,jj,jl) = zbt * psyy(ji,jj,jl) + zbt1 * ( zalf * zalf * zfyy(ji,jj) + zalf1 * zalf1 * psyy(ji,jj,jl) & 633 & + 5.0 * ( zalf * zalf1 * ( - psy(ji,jj,jl) + zfy(ji,jj) ) & 634 & + ( zalf1 - zalf ) * ztemp ) ) 635 psxy(ji,jj,jl) = zbt * psxy(ji,jj,jl) + zbt1 * ( zalf * zfxy(ji,jj) + zalf1 * psxy(ji,jj,jl) & 636 & + 3.0 * ( zalf1 * zfx(ji,jj) - zalf * psx(ji,jj,jl) ) ) 637 psx (ji,jj,jl) = zbt * psx (ji,jj,jl) + zbt1 * ( psx (ji,jj,jl) + zfx (ji,jj) ) 638 psxx(ji,jj,jl) = zbt * psxx(ji,jj,jl) + zbt1 * ( psxx(ji,jj,jl) + zfxx(ji,jj) ) 639 END_2D 640 730 ! 641 731 END DO 642 643 !-- Lateral boundary conditions644 CALL lbc_lnk_multi( 'icedyn_adv_pra', psm(:,:,1:jcat) , 'T', 1.0_wp, ps0 , 'T', 1.0_wp &645 & , psx , 'T', -1.0_wp, psy , 'T', -1.0_wp & ! caution gradient ==> the sign changes646 & , psxx , 'T', 1.0_wp, psyy, 'T', 1.0_wp , psxy, 'T', 1.0_wp )647 732 ! 648 733 END SUBROUTINE adv_y … … 872 957 ! 873 958 ! ! ice thickness 874 CALL iom_get( numrir, jpdom_auto, 'sxice' , sxice )875 CALL iom_get( numrir, jpdom_auto, 'syice' , syice )959 CALL iom_get( numrir, jpdom_auto, 'sxice' , sxice , psgn = -1._wp ) 960 CALL iom_get( numrir, jpdom_auto, 'syice' , syice , psgn = -1._wp ) 876 961 CALL iom_get( numrir, jpdom_auto, 'sxxice', sxxice ) 877 962 CALL iom_get( numrir, jpdom_auto, 'syyice', syyice ) 878 963 CALL iom_get( numrir, jpdom_auto, 'sxyice', sxyice ) 879 964 ! ! snow thickness 880 CALL iom_get( numrir, jpdom_auto, 'sxsn' , sxsn )881 CALL iom_get( numrir, jpdom_auto, 'sysn' , sysn )965 CALL iom_get( numrir, jpdom_auto, 'sxsn' , sxsn , psgn = -1._wp ) 966 CALL iom_get( numrir, jpdom_auto, 'sysn' , sysn , psgn = -1._wp ) 882 967 CALL iom_get( numrir, jpdom_auto, 'sxxsn' , sxxsn ) 883 968 CALL iom_get( numrir, jpdom_auto, 'syysn' , syysn ) 884 969 CALL iom_get( numrir, jpdom_auto, 'sxysn' , sxysn ) 885 970 ! ! ice concentration 886 CALL iom_get( numrir, jpdom_auto, 'sxa' , sxa )887 CALL iom_get( numrir, jpdom_auto, 'sya' , sya )971 CALL iom_get( numrir, jpdom_auto, 'sxa' , sxa , psgn = -1._wp ) 972 CALL iom_get( numrir, jpdom_auto, 'sya' , sya , psgn = -1._wp ) 888 973 CALL iom_get( numrir, jpdom_auto, 'sxxa' , sxxa ) 889 974 CALL iom_get( numrir, jpdom_auto, 'syya' , syya ) 890 975 CALL iom_get( numrir, jpdom_auto, 'sxya' , sxya ) 891 976 ! ! ice salinity 892 CALL iom_get( numrir, jpdom_auto, 'sxsal' , sxsal )893 CALL iom_get( numrir, jpdom_auto, 'sysal' , sysal )977 CALL iom_get( numrir, jpdom_auto, 'sxsal' , sxsal , psgn = -1._wp ) 978 CALL iom_get( numrir, jpdom_auto, 'sysal' , sysal , psgn = -1._wp ) 894 979 CALL iom_get( numrir, jpdom_auto, 'sxxsal', sxxsal ) 895 980 CALL iom_get( numrir, jpdom_auto, 'syysal', syysal ) 896 981 CALL iom_get( numrir, jpdom_auto, 'sxysal', sxysal ) 897 982 ! ! ice age 898 CALL iom_get( numrir, jpdom_auto, 'sxage' , sxage )899 CALL iom_get( numrir, jpdom_auto, 'syage' , syage )983 CALL iom_get( numrir, jpdom_auto, 'sxage' , sxage , psgn = -1._wp ) 984 CALL iom_get( numrir, jpdom_auto, 'syage' , syage , psgn = -1._wp ) 900 985 CALL iom_get( numrir, jpdom_auto, 'sxxage', sxxage ) 901 986 CALL iom_get( numrir, jpdom_auto, 'syyage', syyage ) … … 904 989 DO jk = 1, nlay_s 905 990 WRITE(zchar1,'(I2.2)') jk 906 znam = 'sxc0'//'_l'//zchar1 ; CALL iom_get( numrir, jpdom_auto, znam , z3d ) ; sxc0 (:,:,jk,:) = z3d(:,:,:)907 znam = 'syc0'//'_l'//zchar1 ; CALL iom_get( numrir, jpdom_auto, znam , z3d ) ; syc0 (:,:,jk,:) = z3d(:,:,:)991 znam = 'sxc0'//'_l'//zchar1 ; CALL iom_get( numrir, jpdom_auto, znam , z3d, psgn = -1._wp ) ; sxc0 (:,:,jk,:) = z3d(:,:,:) 992 znam = 'syc0'//'_l'//zchar1 ; CALL iom_get( numrir, jpdom_auto, znam , z3d, psgn = -1._wp ) ; syc0 (:,:,jk,:) = z3d(:,:,:) 908 993 znam = 'sxxc0'//'_l'//zchar1 ; CALL iom_get( numrir, jpdom_auto, znam , z3d ) ; sxxc0(:,:,jk,:) = z3d(:,:,:) 909 994 znam = 'syyc0'//'_l'//zchar1 ; CALL iom_get( numrir, jpdom_auto, znam , z3d ) ; syyc0(:,:,jk,:) = z3d(:,:,:) … … 913 998 DO jk = 1, nlay_i 914 999 WRITE(zchar1,'(I2.2)') jk 915 znam = 'sxe'//'_l'//zchar1 ; CALL iom_get( numrir, jpdom_auto, znam , z3d ) ; sxe (:,:,jk,:) = z3d(:,:,:)916 znam = 'sye'//'_l'//zchar1 ; CALL iom_get( numrir, jpdom_auto, znam , z3d ) ; sye (:,:,jk,:) = z3d(:,:,:)1000 znam = 'sxe'//'_l'//zchar1 ; CALL iom_get( numrir, jpdom_auto, znam , z3d, psgn = -1._wp ) ; sxe (:,:,jk,:) = z3d(:,:,:) 1001 znam = 'sye'//'_l'//zchar1 ; CALL iom_get( numrir, jpdom_auto, znam , z3d, psgn = -1._wp ) ; sye (:,:,jk,:) = z3d(:,:,:) 917 1002 znam = 'sxxe'//'_l'//zchar1 ; CALL iom_get( numrir, jpdom_auto, znam , z3d ) ; sxxe(:,:,jk,:) = z3d(:,:,:) 918 1003 znam = 'syye'//'_l'//zchar1 ; CALL iom_get( numrir, jpdom_auto, znam , z3d ) ; syye(:,:,jk,:) = z3d(:,:,:) … … 921 1006 ! 922 1007 IF( ln_pnd_LEV ) THEN ! melt pond fraction 923 IF( iom_varid( numr or, 'sxap', ldstop = .FALSE. ) > 0 ) THEN924 CALL iom_get( numrir, jpdom_auto, 'sxap' , sxap )925 CALL iom_get( numrir, jpdom_auto, 'syap' , syap )1008 IF( iom_varid( numrir, 'sxap', ldstop = .FALSE. ) > 0 ) THEN 1009 CALL iom_get( numrir, jpdom_auto, 'sxap' , sxap , psgn = -1._wp ) 1010 CALL iom_get( numrir, jpdom_auto, 'syap' , syap , psgn = -1._wp ) 926 1011 CALL iom_get( numrir, jpdom_auto, 'sxxap', sxxap ) 927 1012 CALL iom_get( numrir, jpdom_auto, 'syyap', syyap ) 928 1013 CALL iom_get( numrir, jpdom_auto, 'sxyap', sxyap ) 929 1014 ! ! melt pond volume 930 CALL iom_get( numrir, jpdom_auto, 'sxvp' , sxvp )931 CALL iom_get( numrir, jpdom_auto, 'syvp' , syvp )1015 CALL iom_get( numrir, jpdom_auto, 'sxvp' , sxvp , psgn = -1._wp ) 1016 CALL iom_get( numrir, jpdom_auto, 'syvp' , syvp , psgn = -1._wp ) 932 1017 CALL iom_get( numrir, jpdom_auto, 'sxxvp', sxxvp ) 933 1018 CALL iom_get( numrir, jpdom_auto, 'syyvp', syyvp ) … … 939 1024 ! 940 1025 IF ( ln_pnd_lids ) THEN ! melt pond lid volume 941 IF( iom_varid( numr or, 'sxvl', ldstop = .FALSE. ) > 0 ) THEN942 CALL iom_get( numrir, jpdom_auto, 'sxvl' , sxvl )943 CALL iom_get( numrir, jpdom_auto, 'syvl' , syvl )1026 IF( iom_varid( numrir, 'sxvl', ldstop = .FALSE. ) > 0 ) THEN 1027 CALL iom_get( numrir, jpdom_auto, 'sxvl' , sxvl , psgn = -1._wp ) 1028 CALL iom_get( numrir, jpdom_auto, 'syvl' , syvl , psgn = -1._wp ) 944 1029 CALL iom_get( numrir, jpdom_auto, 'sxxvl', sxxvl ) 945 1030 CALL iom_get( numrir, jpdom_auto, 'syyvl', syyvl ) … … 1056 1141 END SUBROUTINE adv_pra_rst 1057 1142 1143 SUBROUTINE icemax3D( pice , pmax ) 1144 !!--------------------------------------------------------------------- 1145 !! *** ROUTINE icemax3D *** 1146 !! ** Purpose : compute the max of the 9 points around 1147 !!---------------------------------------------------------------------- 1148 REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: pice ! input 1149 REAL(wp), DIMENSION(:,:,:) , INTENT(out) :: pmax ! output 1150 REAL(wp), DIMENSION(2:jpim1,jpj) :: zmax ! temporary array 1151 INTEGER :: ji, jj, jl ! dummy loop indices 1152 !!---------------------------------------------------------------------- 1153 DO jl = 1, jpl 1154 DO jj = Njs0-1, Nje0+1 1155 DO ji = Nis0, Nie0 1156 zmax(ji,jj) = MAX( epsi20, pice(ji,jj,jl), pice(ji-1,jj,jl), pice(ji+1,jj,jl) ) 1157 END DO 1158 END DO 1159 DO jj = Njs0, Nje0 1160 DO ji = Nis0, Nie0 1161 pmax(ji,jj,jl) = MAX( epsi20, zmax(ji,jj), zmax(ji,jj-1), zmax(ji,jj+1) ) 1162 END DO 1163 END DO 1164 END DO 1165 END SUBROUTINE icemax3D 1166 1167 SUBROUTINE icemax4D( pice , pmax ) 1168 !!--------------------------------------------------------------------- 1169 !! *** ROUTINE icemax4D *** 1170 !! ** Purpose : compute the max of the 9 points around 1171 !!---------------------------------------------------------------------- 1172 REAL(wp), DIMENSION(:,:,:,:) , INTENT(in ) :: pice ! input 1173 REAL(wp), DIMENSION(:,:,:,:) , INTENT(out) :: pmax ! output 1174 REAL(wp), DIMENSION(2:jpim1,jpj) :: zmax ! temporary array 1175 INTEGER :: jlay, ji, jj, jk, jl ! dummy loop indices 1176 !!---------------------------------------------------------------------- 1177 jlay = SIZE( pice , 3 ) ! size of input arrays 1178 DO jl = 1, jpl 1179 DO jk = 1, jlay 1180 DO jj = Njs0-1, Nje0+1 1181 DO ji = Nis0, Nie0 1182 zmax(ji,jj) = MAX( epsi20, pice(ji,jj,jk,jl), pice(ji-1,jj,jk,jl), pice(ji+1,jj,jk,jl) ) 1183 END DO 1184 END DO 1185 DO jj = Njs0, Nje0 1186 DO ji = Nis0, Nie0 1187 pmax(ji,jj,jk,jl) = MAX( epsi20, zmax(ji,jj), zmax(ji,jj-1), zmax(ji,jj+1) ) 1188 END DO 1189 END DO 1190 END DO 1191 END DO 1192 END SUBROUTINE icemax4D 1193 1058 1194 #else 1059 1195 !!---------------------------------------------------------------------- -
NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/src/ICE/icedyn_adv_umx.F90
r13553 r13741 92 92 INTEGER :: icycle ! number of sub-timestep for the advection 93 93 REAL(wp) :: zamsk ! 1 if advection of concentration, 0 if advection of other tracers 94 REAL(wp) :: zdt, z vi_cen94 REAL(wp) :: zdt, z1_dt, zvi_cen 95 95 REAL(wp), DIMENSION(1) :: zcflprv, zcflnow ! for global communication 96 96 REAL(wp), DIMENSION(jpi,jpj) :: zudy, zvdx, zcu_box, zcv_box … … 104 104 ! 105 105 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zuv_ho, zvv_ho, zuv_ups, zvv_ups, z1_vi, z1_vs 106 !! diagnostics 107 REAL(wp), DIMENSION(jpi,jpj) :: zdiag_adv_mass, zdiag_adv_salt, zdiag_adv_heat 106 108 !!---------------------------------------------------------------------- 107 109 ! … … 113 115 ELSEWHERE ; zs_i(:,:,:) = 0._wp 114 116 END WHERE 115 DO jl = 1, jpl 116 DO_2D( 0, 0, 0, 0 ) 117 zhip_max(ji,jj,jl) = MAX( epsi20, ph_ip(ji,jj,jl), ph_ip(ji+1,jj ,jl), ph_ip(ji ,jj+1,jl), & 118 & ph_ip(ji-1,jj ,jl), ph_ip(ji ,jj-1,jl), & 119 & ph_ip(ji+1,jj+1,jl), ph_ip(ji-1,jj-1,jl), & 120 & ph_ip(ji+1,jj-1,jl), ph_ip(ji-1,jj+1,jl) ) 121 zhi_max (ji,jj,jl) = MAX( epsi20, ph_i (ji,jj,jl), ph_i (ji+1,jj ,jl), ph_i (ji ,jj+1,jl), & 122 & ph_i (ji-1,jj ,jl), ph_i (ji ,jj-1,jl), & 123 & ph_i (ji+1,jj+1,jl), ph_i (ji-1,jj-1,jl), & 124 & ph_i (ji+1,jj-1,jl), ph_i (ji-1,jj+1,jl) ) 125 zhs_max (ji,jj,jl) = MAX( epsi20, ph_s (ji,jj,jl), ph_s (ji+1,jj ,jl), ph_s (ji ,jj+1,jl), & 126 & ph_s (ji-1,jj ,jl), ph_s (ji ,jj-1,jl), & 127 & ph_s (ji+1,jj+1,jl), ph_s (ji-1,jj-1,jl), & 128 & ph_s (ji+1,jj-1,jl), ph_s (ji-1,jj+1,jl) ) 129 zsi_max (ji,jj,jl) = MAX( epsi20, zs_i (ji,jj,jl), zs_i (ji+1,jj ,jl), zs_i (ji ,jj+1,jl), & 130 & zs_i (ji-1,jj ,jl), zs_i (ji ,jj-1,jl), & 131 & zs_i (ji+1,jj+1,jl), zs_i (ji-1,jj-1,jl), & 132 & zs_i (ji+1,jj-1,jl), zs_i (ji-1,jj+1,jl) ) 133 END_2D 134 END DO 117 CALL icemax3D( ph_i , zhi_max ) 118 CALL icemax3D( ph_s , zhs_max ) 119 CALL icemax3D( ph_ip, zhip_max) 120 CALL icemax3D( zs_i , zsi_max ) 135 121 CALL lbc_lnk_multi( 'icedyn_adv_umx', zhi_max, 'T', 1._wp, zhs_max, 'T', 1._wp, zhip_max, 'T', 1._wp, zsi_max, 'T', 1._wp ) 136 122 ! … … 145 131 ELSEWHERE ; ze_s(:,:,jk,:) = 0._wp 146 132 END WHERE 147 END DO 148 DO jl = 1, jpl 149 DO_3D( 0, 0, 0, 0, 1, nlay_i ) 150 zei_max(ji,jj,jk,jl) = MAX( epsi20, ze_i(ji,jj,jk,jl), ze_i(ji+1,jj ,jk,jl), ze_i(ji ,jj+1,jk,jl), & 151 & ze_i(ji-1,jj ,jk,jl), ze_i(ji ,jj-1,jk,jl), & 152 & ze_i(ji+1,jj+1,jk,jl), ze_i(ji-1,jj-1,jk,jl), & 153 & ze_i(ji+1,jj-1,jk,jl), ze_i(ji-1,jj+1,jk,jl) ) 154 END_3D 155 END DO 156 DO jl = 1, jpl 157 DO_3D( 0, 0, 0, 0, 1, nlay_s ) 158 zes_max(ji,jj,jk,jl) = MAX( epsi20, ze_s(ji,jj,jk,jl), ze_s(ji+1,jj ,jk,jl), ze_s(ji ,jj+1,jk,jl), & 159 & ze_s(ji-1,jj ,jk,jl), ze_s(ji ,jj-1,jk,jl), & 160 & ze_s(ji+1,jj+1,jk,jl), ze_s(ji-1,jj-1,jk,jl), & 161 & ze_s(ji+1,jj-1,jk,jl), ze_s(ji-1,jj+1,jk,jl) ) 162 END_3D 163 END DO 164 CALL lbc_lnk( 'icedyn_adv_pra', zei_max, 'T', 1. ) 165 CALL lbc_lnk( 'icedyn_adv_pra', zes_max, 'T', 1. ) 133 END DO 134 CALL icemax4D( ze_i , zei_max ) 135 CALL icemax4D( ze_s , zes_max ) 136 CALL lbc_lnk( 'icedyn_adv_umx', zei_max, 'T', 1._wp ) 137 CALL lbc_lnk( 'icedyn_adv_umx', zes_max, 'T', 1._wp ) 166 138 ! 167 139 ! … … 179 151 ENDIF 180 152 zdt = rDt_ice / REAL(icycle) 153 z1_dt = 1._wp / zdt 181 154 182 155 ! --- transport --- ! … … 207 180 !---------------! 208 181 DO jt = 1, icycle 182 183 ! diagnostics 184 zdiag_adv_mass(:,:) = SUM( pv_i(:,:,:) , dim=3 ) * rhoi + SUM( pv_s(:,:,:) , dim=3 ) * rhos 185 zdiag_adv_salt(:,:) = SUM( psv_i(:,:,:) , dim=3 ) * rhoi 186 zdiag_adv_heat(:,:) = - SUM(SUM( pe_i(:,:,1:nlay_i,:) , dim=4 ), dim=3 ) & 187 & - SUM(SUM( pe_s(:,:,1:nlay_s,:) , dim=4 ), dim=3 ) 209 188 210 189 ! record at_i before advection (for open water) … … 377 356 ENDIF 378 357 ENDIF 358 359 ! --- Lateral boundary conditions --- ! 360 IF ( ln_pnd_LEV .AND. ln_pnd_lids ) THEN 361 CALL lbc_lnk_multi( 'icedyn_adv_umx', pa_i,'T',1._wp, pv_i,'T',1._wp, pv_s,'T',1._wp, psv_i,'T',1._wp, poa_i,'T',1._wp & 362 & , pa_ip,'T',1._wp, pv_ip,'T',1._wp, pv_il,'T',1._wp ) 363 ELSEIF( ln_pnd_LEV .AND. .NOT.ln_pnd_lids ) THEN 364 CALL lbc_lnk_multi( 'icedyn_adv_umx', pa_i,'T',1._wp, pv_i,'T',1._wp, pv_s,'T',1._wp, psv_i,'T',1._wp, poa_i,'T',1._wp & 365 & , pa_ip,'T',1._wp, pv_ip,'T',1._wp ) 366 ELSE 367 CALL lbc_lnk_multi( 'icedyn_adv_umx', pa_i,'T',1._wp, pv_i,'T',1._wp, pv_s,'T',1._wp, psv_i,'T',1._wp, poa_i,'T',1._wp ) 368 ENDIF 369 CALL lbc_lnk( 'icedyn_adv_umx', pe_i, 'T', 1._wp ) 370 CALL lbc_lnk( 'icedyn_adv_umx', pe_s, 'T', 1._wp ) 379 371 ! 380 372 !== Open water area ==! … … 384 376 & - ( zudy(ji,jj) - zudy(ji-1,jj) + zvdx(ji,jj) - zvdx(ji,jj-1) ) * r1_e1e2t(ji,jj) * zdt 385 377 END_2D 386 CALL lbc_lnk( 'icedyn_adv_umx', pato_i, 'T', 1.0_wp ) 387 ! 378 CALL lbc_lnk( 'icedyn_adv_umx', pato_i, 'T', 1._wp ) 379 ! 380 ! --- diagnostics --- ! 381 diag_adv_mass(:,:) = diag_adv_mass(:,:) + ( SUM( pv_i(:,:,:) , dim=3 ) * rhoi + SUM( pv_s(:,:,:) , dim=3 ) * rhos & 382 & - zdiag_adv_mass(:,:) ) * z1_dt 383 diag_adv_salt(:,:) = diag_adv_salt(:,:) + ( SUM( psv_i(:,:,:) , dim=3 ) * rhoi & 384 & - zdiag_adv_salt(:,:) ) * z1_dt 385 diag_adv_heat(:,:) = diag_adv_heat(:,:) + ( - SUM(SUM( pe_i(:,:,1:nlay_i,:) , dim=4 ), dim=3 ) & 386 & - SUM(SUM( pe_s(:,:,1:nlay_s,:) , dim=4 ), dim=3 ) & 387 & - zdiag_adv_heat(:,:) ) * z1_dt 388 388 ! 389 389 ! --- Ensure non-negative fields and in-bound thicknesses --- ! … … 445 445 !! work on H (and not V). It is partly related to the multi-category approach 446 446 !! Therefore, after advection we limit the thickness to the largest value of the 9-points around (only if ice 447 !! concentration is small). Since we do not limit S and T, large values can occur at the edge but it does not really matter 448 !! since sv_i and e_i are still good. 447 !! concentration is small). We also limit S and T. 449 448 !!---------------------------------------------------------------------- 450 449 REAL(wp) , INTENT(in ) :: pamsk ! advection of concentration (1) or other tracers (0) … … 490 489 IF( pamsk == 0._wp ) THEN 491 490 DO jl = 1, jpl 492 DO_2D( 1, 0, 1, 0 )491 DO_2D( 0, 0, 1, 0 ) 493 492 IF( ABS( pu(ji,jj) ) > epsi10 ) THEN 494 493 zfu_ho (ji,jj,jl) = zfu_ho (ji,jj,jl) * puc (ji,jj,jl) / pu(ji,jj) … … 499 498 ENDIF 500 499 ! 500 END_2D 501 DO_2D( 1, 0, 0, 0 ) 501 502 IF( ABS( pv(ji,jj) ) > epsi10 ) THEN 502 503 zfv_ho (ji,jj,jl) = zfv_ho (ji,jj,jl) * pvc (ji,jj,jl) / pv(ji,jj) … … 533 534 IF( PRESENT( pua_ho ) ) THEN 534 535 DO jl = 1, jpl 535 DO_2D( 1, 0, 1, 0 ) 536 pua_ho (ji,jj,jl) = zfu_ho (ji,jj,jl) ; pva_ho (ji,jj,jl) = zfv_ho (ji,jj,jl) 537 pua_ups(ji,jj,jl) = zfu_ups(ji,jj,jl) ; pva_ups(ji,jj,jl) = zfv_ups(ji,jj,jl) 536 DO_2D( 0, 0, 1, 0 ) 537 pua_ho (ji,jj,jl) = zfu_ho (ji,jj,jl) 538 pua_ups(ji,jj,jl) = zfu_ups(ji,jj,jl) 539 END_2D 540 DO_2D( 1, 0, 0, 0 ) 541 pva_ho (ji,jj,jl) = zfv_ho (ji,jj,jl) 542 pva_ups(ji,jj,jl) = zfv_ups(ji,jj,jl) 538 543 END_2D 539 544 END DO … … 549 554 END_2D 550 555 END DO 551 CALL lbc_lnk( 'icedyn_adv_umx', ptc, 'T', 1.0_wp )552 556 ! 553 557 END SUBROUTINE adv_umx … … 588 592 ! 589 593 DO jl = 1, jpl !-- flux in x-direction 590 DO_2D( 1, 0, 1, 0 )594 DO_2D( 1, 1, 1, 0 ) 591 595 pfu_ups(ji,jj,jl) = MAX( pu(ji,jj), 0._wp ) * pt(ji,jj,jl) + MIN( pu(ji,jj), 0._wp ) * pt(ji+1,jj,jl) 592 596 END_2D … … 594 598 ! 595 599 DO jl = 1, jpl !-- first guess of tracer from u-flux 596 DO_2D( 0, 0, 0, 0 )600 DO_2D( 1, 1, 0, 0 ) 597 601 ztra = - ( pfu_ups(ji,jj,jl) - pfu_ups(ji-1,jj,jl) ) & 598 602 & + ( pu (ji,jj ) - pu (ji-1,jj ) ) * pt(ji,jj,jl) * (1.-pamsk) … … 601 605 END_2D 602 606 END DO 603 CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1.0_wp )604 607 ! 605 608 DO jl = 1, jpl !-- flux in y-direction 606 DO_2D( 1, 0, 1, 0 )609 DO_2D( 1, 0, 0, 0 ) 607 610 pfv_ups(ji,jj,jl) = MAX( pv(ji,jj), 0._wp ) * zpt(ji,jj,jl) + MIN( pv(ji,jj), 0._wp ) * zpt(ji,jj+1,jl) 608 611 END_2D … … 612 615 ! 613 616 DO jl = 1, jpl !-- flux in y-direction 614 DO_2D( 1, 0, 1, 0)617 DO_2D( 1, 0, 1, 1 ) 615 618 pfv_ups(ji,jj,jl) = MAX( pv(ji,jj), 0._wp ) * pt(ji,jj,jl) + MIN( pv(ji,jj), 0._wp ) * pt(ji,jj+1,jl) 616 619 END_2D … … 618 621 ! 619 622 DO jl = 1, jpl !-- first guess of tracer from v-flux 620 DO_2D( 0, 0, 0, 0)623 DO_2D( 0, 0, 1, 1 ) 621 624 ztra = - ( pfv_ups(ji,jj,jl) - pfv_ups(ji,jj-1,jl) ) & 622 625 & + ( pv (ji,jj ) - pv (ji,jj-1 ) ) * pt(ji,jj,jl) * (1.-pamsk) … … 625 628 END_2D 626 629 END DO 627 CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1.0_wp )628 630 ! 629 631 DO jl = 1, jpl !-- flux in x-direction 630 DO_2D( 1, 0, 1, 0 )632 DO_2D( 0, 0, 1, 0 ) 631 633 pfu_ups(ji,jj,jl) = MAX( pu(ji,jj), 0._wp ) * zpt(ji,jj,jl) + MIN( pu(ji,jj), 0._wp ) * zpt(ji+1,jj,jl) 632 634 END_2D … … 677 679 ! 678 680 DO jl = 1, jpl 679 DO_2D( 1, 0, 1, 0 )681 DO_2D( 1, 1, 1, 0 ) 680 682 pfu_ho(ji,jj,jl) = 0.5_wp * pu(ji,jj) * ( pt(ji,jj,jl) + pt(ji+1,jj ,jl) ) 683 END_2D 684 DO_2D( 1, 0, 1, 1 ) 681 685 pfv_ho(ji,jj,jl) = 0.5_wp * pv(ji,jj) * ( pt(ji,jj,jl) + pt(ji ,jj+1,jl) ) 682 686 END_2D … … 695 699 ! 696 700 DO jl = 1, jpl !-- flux in x-direction 697 DO_2D( 1, 0, 1, 0 )701 DO_2D( 1, 1, 1, 0 ) 698 702 pfu_ho(ji,jj,jl) = 0.5_wp * pu(ji,jj) * ( pt(ji,jj,jl) + pt(ji+1,jj,jl) ) 699 703 END_2D … … 702 706 703 707 DO jl = 1, jpl !-- first guess of tracer from u-flux 704 DO_2D( 0, 0, 0, 0 )708 DO_2D( 1, 1, 0, 0 ) 705 709 ztra = - ( pfu_ho(ji,jj,jl) - pfu_ho(ji-1,jj,jl) ) & 706 710 & + ( pu (ji,jj ) - pu (ji-1,jj ) ) * pt(ji,jj,jl) * (1.-pamsk) … … 709 713 END_2D 710 714 END DO 711 CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1.0_wp )712 715 713 716 DO jl = 1, jpl !-- flux in y-direction 714 DO_2D( 1, 0, 1, 0 )717 DO_2D( 1, 0, 0, 0 ) 715 718 pfv_ho(ji,jj,jl) = 0.5_wp * pv(ji,jj) * ( zpt(ji,jj,jl) + zpt(ji,jj+1,jl) ) 716 719 END_2D … … 721 724 ! 722 725 DO jl = 1, jpl !-- flux in y-direction 723 DO_2D( 1, 0, 1, 0)726 DO_2D( 1, 0, 1, 1 ) 724 727 pfv_ho(ji,jj,jl) = 0.5_wp * pv(ji,jj) * ( pt(ji,jj,jl) + pt(ji,jj+1,jl) ) 725 728 END_2D … … 728 731 ! 729 732 DO jl = 1, jpl !-- first guess of tracer from v-flux 730 DO_2D( 0, 0, 0, 0)733 DO_2D( 0, 0, 1, 1 ) 731 734 ztra = - ( pfv_ho(ji,jj,jl) - pfv_ho(ji,jj-1,jl) ) & 732 735 & + ( pv (ji,jj ) - pv (ji,jj-1 ) ) * pt(ji,jj,jl) * (1.-pamsk) … … 735 738 END_2D 736 739 END DO 737 CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1.0_wp )738 740 ! 739 741 DO jl = 1, jpl !-- flux in x-direction 740 DO_2D( 1, 0, 1, 0 )742 DO_2D( 0, 0, 1, 0 ) 741 743 pfu_ho(ji,jj,jl) = 0.5_wp * pu(ji,jj) * ( zpt(ji,jj,jl) + zpt(ji+1,jj,jl) ) 742 744 END_2D … … 895 897 ! 896 898 DO jl = 1, jpl 897 DO_2D( 1, 0, 1, 0 )899 DO_2D( 0, 0, 1, 0 ) 898 900 pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * ( pt(ji+1,jj,jl) + pt(ji,jj,jl) & 899 901 & - SIGN( 1._wp, pu(ji,jj) ) * ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) ) … … 904 906 ! 905 907 DO jl = 1, jpl 906 DO_2D( 1, 0, 1, 0 )908 DO_2D( 0, 0, 1, 0 ) 907 909 zcu = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 908 910 pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * ( pt(ji+1,jj,jl) + pt(ji,jj,jl) & … … 914 916 ! 915 917 DO jl = 1, jpl 916 DO_2D( 1, 0, 1, 0 )918 DO_2D( 0, 0, 1, 0 ) 917 919 zcu = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 918 920 zdx2 = e1u(ji,jj) * e1u(ji,jj) … … 928 930 ! 929 931 DO jl = 1, jpl 930 DO_2D( 1, 0, 1, 0 )932 DO_2D( 0, 0, 1, 0 ) 931 933 zcu = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 932 934 zdx2 = e1u(ji,jj) * e1u(ji,jj) … … 942 944 ! 943 945 DO jl = 1, jpl 944 DO_2D( 1, 0, 1, 0 )946 DO_2D( 0, 0, 1, 0 ) 945 947 zcu = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 946 948 zdx2 = e1u(ji,jj) * e1u(ji,jj) … … 963 965 IF( ll_neg ) THEN 964 966 DO jl = 1, jpl 965 DO_2D( 1, 0, 1, 0 )967 DO_2D( 0, 0, 1, 0 ) 966 968 IF( pt_u(ji,jj,jl) < 0._wp .OR. ( imsk_small(ji,jj,jl) == 0 .AND. pamsk == 0. ) ) THEN 967 969 pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * ( pt(ji+1,jj,jl) + pt(ji,jj,jl) & … … 973 975 ! !-- High order flux in i-direction --! 974 976 DO jl = 1, jpl 975 DO_2D( 1, 0, 1, 0 )977 DO_2D( 0, 0, 1, 0 ) 976 978 pfu_ho(ji,jj,jl) = pu(ji,jj) * pt_u(ji,jj,jl) 977 979 END_2D … … 1031 1033 CASE( 1 ) !== 1st order central TIM ==! (Eq. 21) 1032 1034 DO jl = 1, jpl 1033 DO_2D( 1, 0, 1, 0 )1035 DO_2D( 1, 0, 0, 0 ) 1034 1036 pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * ( pt(ji,jj+1,jl) + pt(ji,jj,jl) & 1035 1037 & - SIGN( 1._wp, pv(ji,jj) ) * ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) ) … … 1039 1041 CASE( 2 ) !== 2nd order central TIM ==! (Eq. 23) 1040 1042 DO jl = 1, jpl 1041 DO_2D( 1, 0, 1, 0 )1043 DO_2D( 1, 0, 0, 0 ) 1042 1044 zcv = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 1043 1045 pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * ( pt(ji,jj+1,jl) + pt(ji,jj,jl) & … … 1048 1050 CASE( 3 ) !== 3rd order central TIM ==! (Eq. 24) 1049 1051 DO jl = 1, jpl 1050 DO_2D( 1, 0, 1, 0 )1052 DO_2D( 1, 0, 0, 0 ) 1051 1053 zcv = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 1052 1054 zdy2 = e2v(ji,jj) * e2v(ji,jj) … … 1061 1063 CASE( 4 ) !== 4th order central TIM ==! (Eq. 27) 1062 1064 DO jl = 1, jpl 1063 DO_2D( 1, 0, 1, 0 )1065 DO_2D( 1, 0, 0, 0 ) 1064 1066 zcv = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 1065 1067 zdy2 = e2v(ji,jj) * e2v(ji,jj) … … 1074 1076 CASE( 5 ) !== 5th order central TIM ==! (Eq. 29) 1075 1077 DO jl = 1, jpl 1076 DO_2D( 1, 0, 1, 0 )1078 DO_2D( 1, 0, 0, 0 ) 1077 1079 zcv = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 1078 1080 zdy2 = e2v(ji,jj) * e2v(ji,jj) … … 1095 1097 IF( ll_neg ) THEN 1096 1098 DO jl = 1, jpl 1097 DO_2D( 1, 0, 1, 0 )1099 DO_2D( 1, 0, 0, 0 ) 1098 1100 IF( pt_v(ji,jj,jl) < 0._wp .OR. ( jmsk_small(ji,jj,jl) == 0 .AND. pamsk == 0. ) ) THEN 1099 1101 pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * ( ( pt(ji,jj+1,jl) + pt(ji,jj,jl) ) & … … 1105 1107 ! !-- High order flux in j-direction --! 1106 1108 DO jl = 1, jpl 1107 DO_2D( 1, 0, 1, 0 )1109 DO_2D( 1, 0, 0, 0 ) 1108 1110 pfv_ho(ji,jj,jl) = pv(ji,jj) * pt_v(ji,jj,jl) 1109 1111 END_2D … … 1141 1143 ! -------------------------------------------------- 1142 1144 DO jl = 1, jpl 1143 DO_2D( 1, 0, 1, 0 )1145 DO_2D( 0, 0, 1, 0 ) 1144 1146 pfu_ho(ji,jj,jl) = pfu_ho(ji,jj,jl) - pfu_ups(ji,jj,jl) 1147 END_2D 1148 DO_2D( 1, 0, 0, 0 ) 1145 1149 pfv_ho(ji,jj,jl) = pfv_ho(ji,jj,jl) - pfv_ups(ji,jj,jl) 1146 1150 END_2D … … 1248 1252 ! --------------------------------- 1249 1253 DO jl = 1, jpl 1250 DO_2D( 1, 0, 1, 0 )1254 DO_2D( 0, 0, 1, 0 ) 1251 1255 zau = MIN( 1._wp , zbetdo(ji,jj,jl) , zbetup(ji+1,jj,jl) ) 1252 1256 zbu = MIN( 1._wp , zbetup(ji,jj,jl) , zbetdo(ji+1,jj,jl) ) … … 1259 1263 END_2D 1260 1264 1261 DO_2D( 1, 0, 1, 0 )1265 DO_2D( 1, 0, 0, 0 ) 1262 1266 zav = MIN( 1._wp , zbetdo(ji,jj,jl) , zbetup(ji,jj+1,jl) ) 1263 1267 zbv = MIN( 1._wp , zbetup(ji,jj,jl) , zbetdo(ji,jj+1,jl) ) … … 1616 1620 END SUBROUTINE Hsnow 1617 1621 1622 SUBROUTINE icemax3D( pice , pmax ) 1623 !!--------------------------------------------------------------------- 1624 !! *** ROUTINE icemax3D *** 1625 !! ** Purpose : compute the max of the 9 points around 1626 !!---------------------------------------------------------------------- 1627 REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: pice ! input 1628 REAL(wp), DIMENSION(:,:,:) , INTENT(out) :: pmax ! output 1629 REAL(wp), DIMENSION(2:jpim1,jpj) :: zmax ! temporary array 1630 INTEGER :: ji, jj, jl ! dummy loop indices 1631 !!---------------------------------------------------------------------- 1632 DO jl = 1, jpl 1633 DO jj = Njs0-1, Nje0+1 1634 DO ji = Nis0, Nie0 1635 zmax(ji,jj) = MAX( epsi20, pice(ji,jj,jl), pice(ji-1,jj,jl), pice(ji+1,jj,jl) ) 1636 END DO 1637 END DO 1638 DO jj = Njs0, Nje0 1639 DO ji = Nis0, Nie0 1640 pmax(ji,jj,jl) = MAX( epsi20, zmax(ji,jj), zmax(ji,jj-1), zmax(ji,jj+1) ) 1641 END DO 1642 END DO 1643 END DO 1644 END SUBROUTINE icemax3D 1645 1646 SUBROUTINE icemax4D( pice , pmax ) 1647 !!--------------------------------------------------------------------- 1648 !! *** ROUTINE icemax4D *** 1649 !! ** Purpose : compute the max of the 9 points around 1650 !!---------------------------------------------------------------------- 1651 REAL(wp), DIMENSION(:,:,:,:) , INTENT(in ) :: pice ! input 1652 REAL(wp), DIMENSION(:,:,:,:) , INTENT(out) :: pmax ! output 1653 REAL(wp), DIMENSION(2:jpim1,jpj) :: zmax ! temporary array 1654 INTEGER :: jlay, ji, jj, jk, jl ! dummy loop indices 1655 !!---------------------------------------------------------------------- 1656 jlay = SIZE( pice , 3 ) ! size of input arrays 1657 DO jl = 1, jpl 1658 DO jk = 1, jlay 1659 DO jj = Njs0-1, Nje0+1 1660 DO ji = Nis0, Nie0 1661 zmax(ji,jj) = MAX( epsi20, pice(ji,jj,jk,jl), pice(ji-1,jj,jk,jl), pice(ji+1,jj,jk,jl) ) 1662 END DO 1663 END DO 1664 DO jj = Njs0, Nje0 1665 DO ji = Nis0, Nie0 1666 pmax(ji,jj,jk,jl) = MAX( epsi20, zmax(ji,jj), zmax(ji,jj-1), zmax(ji,jj+1) ) 1667 END DO 1668 END DO 1669 END DO 1670 END DO 1671 END SUBROUTINE icemax4D 1618 1672 1619 1673 #else -
NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/src/ICE/icedyn_rdgrft.F90
r13553 r13741 349 349 ELSEIF( zGsum(ji,jl-1) < rn_gstar ) THEN 350 350 apartf(ji,jl) = z1_gstar * ( rn_gstar - zGsum(ji,jl-1) ) * & 351 & ( 2._wp - ( zGsum(ji,jl-1) + rn_gstar 351 & ( 2._wp - ( zGsum(ji,jl-1) + rn_gstar ) * z1_gstar ) 352 352 ELSE 353 353 apartf(ji,jl) = 0._wp -
NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/src/ICE/icedyn_rhg_evp.F90
r13553 r13741 148 148 ! 149 149 REAL(wp), DIMENSION(jpi,jpj) :: zds ! shear 150 REAL(wp), DIMENSION(jpi,jpj) :: zten_i ! tension 150 151 REAL(wp), DIMENSION(jpi,jpj) :: zs1, zs2, zs12 ! stress tensor components 151 152 REAL(wp), DIMENSION(jpi,jpj) :: zsshdyn ! array used for the calculation of ice surface slope: … … 170 171 REAL(wp), DIMENSION(jpi,jpj) :: zu_ice, zv_ice 171 172 !! --- diags 172 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zsig1, zsig2, zsig3 173 REAL(wp) :: zsig1, zsig2, zsig12, zfac, z1_strength 174 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zsig_I, zsig_II, zsig1_p, zsig2_p 173 175 !! --- SIMIP diags 174 176 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_xmtrp_ice ! X-component of ice mass transport (kg/s) … … 726 728 & ) * r1_e1e2t(ji,jj) 727 729 zdt2 = zdt * zdt 730 731 zten_i(ji,jj) = zdt 728 732 729 733 ! shear**2 at T points (doc eq. A16) … … 741 745 742 746 ! delta at T points 743 z delta(ji,jj) = SQRT( pdivu_i(ji,jj) * pdivu_i(ji,jj) + ( zdt2 + zds2 ) * z1_ecc2 )744 rswitch = 1._wp - MAX( 0._wp, SIGN( 1._wp, -z delta(ji,jj)) ) ! 0 if delta=0745 pdelta_i(ji,jj) = z delta(ji,jj) + rn_creepl * rswitch747 zfac = SQRT( pdivu_i(ji,jj) * pdivu_i(ji,jj) + ( zdt2 + zds2 ) * z1_ecc2 ) ! delta 748 rswitch = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zfac ) ) ! 0 if delta=0 749 pdelta_i(ji,jj) = zfac + rn_creepl * rswitch ! delta+creepl 746 750 747 751 END_2D 748 CALL lbc_lnk_multi( 'icedyn_rhg_evp', pshear_i, 'T', 1.0_wp, pdivu_i, 'T', 1.0_wp, pdelta_i, 'T', 1.0_wp ) 752 CALL lbc_lnk_multi( 'icedyn_rhg_evp', pshear_i, 'T', 1._wp, pdivu_i, 'T', 1._wp, pdelta_i, 'T', 1._wp, zten_i, 'T', 1._wp, & 753 & zs1 , 'T', 1._wp, zs2 , 'T', 1._wp, zs12 , 'F', 1._wp ) 749 754 750 755 ! --- Store the stress tensor for the next time step --- ! 751 CALL lbc_lnk_multi( 'icedyn_rhg_evp', zs1, 'T', 1.0_wp, zs2, 'T', 1.0_wp, zs12, 'F', 1.0_wp )752 756 pstress1_i (:,:) = zs1 (:,:) 753 757 pstress2_i (:,:) = zs2 (:,:) … … 778 782 IF( iom_use('icestr') ) CALL iom_put( 'icestr' , strength * zmsk00 ) ! strength 779 783 780 ! --- stress tensor--- !781 IF( iom_use(' isig1') .OR. iom_use('isig2') .OR. iom_use('isig3') .OR. iom_use('normstr') .OR. iom_use('sheastr') ) THEN782 ! 783 ALLOCATE( zsig 1(jpi,jpj) , zsig2(jpi,jpj) , zsig3(jpi,jpj) )784 ! --- Stress tensor invariants (SIMIP diags) --- ! 785 IF( iom_use('normstr') .OR. iom_use('sheastr') ) THEN 786 ! 787 ALLOCATE( zsig_I(jpi,jpj) , zsig_II(jpi,jpj) ) 784 788 ! 785 DO_2D( 0, 0, 0, 0 ) 786 zdum1 = ( zmsk00(ji-1,jj) * pstress12_i(ji-1,jj) + zmsk00(ji ,jj-1) * pstress12_i(ji ,jj-1) + & ! stress12_i at T-point 787 & zmsk00(ji ,jj) * pstress12_i(ji ,jj) + zmsk00(ji-1,jj-1) * pstress12_i(ji-1,jj-1) ) & 788 & / MAX( 1._wp, zmsk00(ji-1,jj) + zmsk00(ji,jj-1) + zmsk00(ji,jj) + zmsk00(ji-1,jj-1) ) 789 790 zshear = SQRT( pstress2_i(ji,jj) * pstress2_i(ji,jj) + 4._wp * zdum1 * zdum1 ) ! shear stress 791 792 zdum2 = zmsk00(ji,jj) / MAX( 1._wp, strength(ji,jj) ) 793 794 !! zsig1(ji,jj) = 0.5_wp * zdum2 * ( pstress1_i(ji,jj) + zshear ) ! principal stress (y-direction, see Hunke & Dukowicz 2002) 795 !! zsig2(ji,jj) = 0.5_wp * zdum2 * ( pstress1_i(ji,jj) - zshear ) ! principal stress (x-direction, see Hunke & Dukowicz 2002) 796 !! zsig3(ji,jj) = zdum2**2 * ( ( pstress1_i(ji,jj) + strength(ji,jj) )**2 + ( rn_ecc * zshear )**2 ) ! quadratic relation linking compressive stress to shear stress 797 !! ! (scheme converges if this value is ~1, see Bouillon et al 2009 (eq. 11)) 798 zsig1(ji,jj) = 0.5_wp * zdum2 * ( pstress1_i(ji,jj) ) ! compressive stress, see Bouillon et al. 2015 799 zsig2(ji,jj) = 0.5_wp * zdum2 * ( zshear ) ! shear stress 800 zsig3(ji,jj) = zdum2**2 * ( ( pstress1_i(ji,jj) + strength(ji,jj) )**2 + ( rn_ecc * zshear )**2 ) 801 END_2D 802 CALL lbc_lnk_multi( 'icedyn_rhg_evp', zsig1, 'T', 1.0_wp, zsig2, 'T', 1.0_wp, zsig3, 'T', 1.0_wp ) 803 ! 804 CALL iom_put( 'isig1' , zsig1 ) 805 CALL iom_put( 'isig2' , zsig2 ) 806 CALL iom_put( 'isig3' , zsig3 ) 807 ! 808 ! Stress tensor invariants (normal and shear stress N/m) 809 IF( iom_use('normstr') ) CALL iom_put( 'normstr' , ( zs1(:,:) + zs2(:,:) ) * zmsk00(:,:) ) ! Normal stress 810 IF( iom_use('sheastr') ) CALL iom_put( 'sheastr' , SQRT( ( zs1(:,:) - zs2(:,:) )**2 + 4*zs12(:,:)**2 ) * zmsk00(:,:) ) ! Shear stress 811 812 DEALLOCATE( zsig1 , zsig2 , zsig3 ) 789 DO_2D( 1, 1, 1, 1 ) 790 791 ! Ice stresses 792 ! sigma1, sigma2, sigma12 are some useful recombination of the stresses (Hunke and Dukowicz MWR 2002, Bouillon et al., OM2013) 793 ! These are NOT stress tensor components, neither stress invariants, neither stress principal components 794 ! I know, this can be confusing... 795 zfac = strength(ji,jj) / ( pdelta_i(ji,jj) + rn_creepl ) 796 zsig1 = zfac * ( pdivu_i(ji,jj) - pdelta_i(ji,jj) ) 797 zsig2 = zfac * z1_ecc2 * zten_i(ji,jj) 798 zsig12 = zfac * z1_ecc2 * pshear_i(ji,jj) 799 800 ! Stress invariants (sigma_I, sigma_II, Coon 1974, Feltham 2008) 801 zsig_I (ji,jj) = zsig1 * 0.5_wp ! 1st stress invariant, aka average normal stress, aka negative pressure 802 zsig_II(ji,jj) = SQRT ( MAX( 0._wp, zsig2 * zsig2 * 0.25_wp + zsig12 ) ) ! 2nd '' '', aka maximum shear stress 803 804 END_2D 805 ! 806 ! Stress tensor invariants (normal and shear stress N/m) - SIMIP diags - definitions following Coon (1974) and Feltham (2008) 807 IF( iom_use('normstr') ) CALL iom_put( 'normstr', zsig_I (:,:) * zmsk00(:,:) ) ! Normal stress 808 IF( iom_use('sheastr') ) CALL iom_put( 'sheastr', zsig_II(:,:) * zmsk00(:,:) ) ! Maximum shear stress 809 810 DEALLOCATE ( zsig_I, zsig_II ) 811 812 ENDIF 813 814 ! --- Normalized stress tensor principal components --- ! 815 ! This are used to plot the normalized yield curve, see Lemieux & Dupont, 2020 816 ! Recommendation 1 : we use ice strength, not replacement pressure 817 ! Recommendation 2 : need to use deformations at PREVIOUS iterate for viscosities 818 IF( iom_use('sig1_pnorm') .OR. iom_use('sig2_pnorm') ) THEN 819 ! 820 ALLOCATE( zsig1_p(jpi,jpj) , zsig2_p(jpi,jpj) , zsig_I(jpi,jpj) , zsig_II(jpi,jpj) ) 821 ! 822 DO_2D( 1, 1, 1, 1 ) 823 824 ! Ice stresses computed with **viscosities** (delta, p/delta) at **previous** iterates 825 ! and **deformations** at current iterates 826 ! following Lemieux & Dupont (2020) 827 zfac = zp_delt(ji,jj) 828 zsig1 = zfac * ( pdivu_i(ji,jj) - ( zdelta(ji,jj) + rn_creepl ) ) 829 zsig2 = zfac * z1_ecc2 * zten_i(ji,jj) 830 zsig12 = zfac * z1_ecc2 * pshear_i(ji,jj) 831 832 ! Stress invariants (sigma_I, sigma_II, Coon 1974, Feltham 2008), T-point 833 zsig_I(ji,jj) = zsig1 * 0.5_wp ! 1st stress invariant, aka average normal stress, aka negative pressure 834 zsig_II(ji,jj) = SQRT ( MAX( 0._wp, zsig2 * zsig2 * 0.25_wp + zsig12 ) ) ! 2nd '' '', aka maximum shear stress 835 836 ! Normalized principal stresses (used to display the ellipse) 837 z1_strength = 1._wp / MAX( 1._wp, strength(ji,jj) ) 838 zsig1_p(ji,jj) = ( zsig_I(ji,jj) + zsig_II(ji,jj) ) * z1_strength 839 zsig2_p(ji,jj) = ( zsig_I(ji,jj) - zsig_II(ji,jj) ) * z1_strength 840 END_2D 841 ! 842 CALL iom_put( 'sig1_pnorm' , zsig1_p ) 843 CALL iom_put( 'sig2_pnorm' , zsig2_p ) 844 845 DEALLOCATE( zsig1_p , zsig2_p , zsig_I, zsig_II ) 846 813 847 ENDIF 814 848 … … 946 980 istatus = NF90_PUT_VAR( ncvgid, nvarid, (/zresm/), (/it/), (/1/) ) 947 981 ! close file 948 IF( kt == nitend ) istatus = NF90_CLOSE(ncvgid)982 IF( kt == nitend - nn_fsbc + 1 ) istatus = NF90_CLOSE(ncvgid) 949 983 ENDIF 950 984 -
NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/src/ICE/iceitd.F90
r13553 r13741 627 627 END_2D 628 628 ! 629 !!clem CALL tab_2d_1d( npti, nptidx(1:npti), h_i_1d(1:npti), h_i(:,:,jl) ) 630 CALL tab_2d_1d( npti, nptidx(1:npti), a_i_1d(1:npti), a_i(:,:,jl) )631 CALL tab_2d_1d( npti, nptidx(1:npti), v_i_1d(1:npti), v_i(:,:,jl) )632 !633 DO ji = 1, npti634 jdonor(ji,jl) = jl635 ! how much of a_i you send in cat sup is somewhat arbitrary636 !!clem: these do not work properly after a restart (I do not know why) => not sure it is still true 637 !! zdaice(ji,jl) = a_i_1d(ji) * ( h_i_1d(ji) - hi_max(jl) + epsi10 ) / h_i_1d(ji) 638 !! zdvice(ji,jl) = v_i_1d(ji) - ( a_i_1d(ji) - zdaice(ji,jl) ) * ( hi_max(jl) - epsi10 ) 639 !!clem: these do not work properly after a restart (I do not know why) => not sure it is still true 640 !! zdaice(ji,jl) = a_i_1d(ji) 641 !! zdvice(ji,jl) = v_i_1d(ji)642 !!clem: these are from UCL and work ok 643 zdaice(ji,jl) = a_i_1d(ji) * 0.5_wp644 zdvice(ji,jl) = v_i_1d(ji) - zdaice(ji,jl) * ( hi_max(jl) + hi_max(jl-1)) * 0.5_wp645 END DO646 !647 IF( npti > 0 ) THEN629 IF( npti > 0 ) THEN 630 !!clem CALL tab_2d_1d( npti, nptidx(1:npti), h_i_1d(1:npti), h_i(:,:,jl) ) 631 CALL tab_2d_1d( npti, nptidx(1:npti), a_i_1d(1:npti), a_i(:,:,jl) ) 632 CALL tab_2d_1d( npti, nptidx(1:npti), v_i_1d(1:npti), v_i(:,:,jl) ) 633 ! 634 DO ji = 1, npti 635 jdonor(ji,jl) = jl 636 ! how much of a_i you send in cat sup is somewhat arbitrary 637 !!clem: these do not work properly after a restart (I do not know why) => not sure it is still true 638 !! zdaice(ji,jl) = a_i_1d(ji) * ( h_i_1d(ji) - hi_max(jl) + epsi10 ) / h_i_1d(ji) 639 !! zdvice(ji,jl) = v_i_1d(ji) - ( a_i_1d(ji) - zdaice(ji,jl) ) * ( hi_max(jl) - epsi10 ) 640 !!clem: these do not work properly after a restart (I do not know why) => not sure it is still true 641 !! zdaice(ji,jl) = a_i_1d(ji) 642 !! zdvice(ji,jl) = v_i_1d(ji) 643 !!clem: these are from UCL and work ok 644 zdaice(ji,jl) = a_i_1d(ji) * 0.5_wp 645 zdvice(ji,jl) = v_i_1d(ji) - zdaice(ji,jl) * ( hi_max(jl) + hi_max(jl-1) ) * 0.5_wp 646 END DO 647 ! 648 648 CALL itd_shiftice( jdonor(1:npti,:), zdaice(1:npti,:), zdvice(1:npti,:) ) ! Shift jl=>jl+1 649 649 ! Reset shift parameters … … 666 666 END_2D 667 667 ! 668 CALL tab_2d_1d( npti, nptidx(1:npti), a_i_1d(1:npti), a_i(:,:,jl+1) ) ! jl+1 is ok669 CALL tab_2d_1d( npti, nptidx(1:npti), v_i_1d(1:npti), v_i(:,:,jl+1) ) ! jl+1 is ok670 DO ji = 1, npti671 jdonor(ji,jl) = jl + 1672 zdaice(ji,jl) = a_i_1d(ji)673 zdvice(ji,jl) = v_i_1d(ji)674 END DO675 !676 668 IF( npti > 0 ) THEN 669 CALL tab_2d_1d( npti, nptidx(1:npti), a_i_1d(1:npti), a_i(:,:,jl+1) ) ! jl+1 is ok 670 CALL tab_2d_1d( npti, nptidx(1:npti), v_i_1d(1:npti), v_i(:,:,jl+1) ) ! jl+1 is ok 671 DO ji = 1, npti 672 jdonor(ji,jl) = jl + 1 673 zdaice(ji,jl) = a_i_1d(ji) 674 zdvice(ji,jl) = v_i_1d(ji) 675 END DO 676 ! 677 677 CALL itd_shiftice( jdonor(1:npti,:), zdaice(1:npti,:), zdvice(1:npti,:) ) ! Shift jl+1=>jl 678 678 ! Reset shift parameters -
NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/src/ICE/icestp.F90
r13553 r13741 55 55 USE icedyn ! sea-ice: dynamics 56 56 USE icethd ! sea-ice: thermodynamics 57 USE icecor ! sea-ice: corrections58 57 USE iceupdate ! sea-ice: sea surface boundary condition update 59 58 USE icedia ! sea-ice: budget diagnostics … … 86 85 PUBLIC ice_init ! called by sbcmod.F90 87 86 87 !! * Substitutions 88 # include "do_loop_substitute.h90" 88 89 !!---------------------------------------------------------------------- 89 90 !! NEMO/ICE 4.0 , NEMO Consortium (2018) … … 160 161 IF( ln_icedyn .AND. .NOT.lk_c1d ) & 161 162 & CALL ice_dyn( kt, Kmm ) ! -- Ice dynamics 163 ! 164 CALL diag_trends( 1 ) ! record dyn trends 162 165 ! 163 166 ! !== lateral boundary conditions ==! … … 188 191 IF( ln_icethd ) CALL ice_thd( kt ) ! -- Ice thermodynamics 189 192 ! 190 CALL ice_cor( kt , 2 ) ! -- Corrections 191 ! 193 CALL diag_trends( 2 ) ! record thermo trends 192 194 CALL ice_var_glo2eqv ! necessary calls (at least for coupling) 193 195 CALL ice_var_agg( 2 ) ! necessary calls (at least for coupling) … … 196 198 ! 197 199 IF( ln_icediahsb ) CALL ice_dia( kt ) ! -- Diagnostics outputs 200 ! 201 IF( ln_icediachk ) CALL ice_drift_wri( kt ) ! -- Diagnostics outputs for conservation 198 202 ! 199 203 CALL ice_wri( kt ) ! -- Ice outputs … … 263 267 CALL ice_thd_init ! set ice thermodynics parameters (clem: important to call it first for melt ponds) 264 268 ! 265 ! ! Initial sea-ice state 266 CALL ice_istate_init 269 CALL ice_sbc_init ! set ice-ocean and ice-atm. coupling parameters 270 ! 271 CALL ice_istate_init ! Initial sea-ice state 267 272 IF ( ln_rstart .OR. nn_iceini_file == 2 ) THEN 268 273 CALL ice_rst_read( Kbb, Kmm, Kaa ) ! start from a restart file … … 273 278 CALL ice_var_agg(1) 274 279 ! 275 CALL ice_sbc_init ! set ice-ocean and ice-atm. coupling parameters276 !277 280 CALL ice_dyn_init ! set ice dynamics parameters 278 281 ! … … 282 285 ! 283 286 CALL ice_dia_init ! initialization for diags 287 ! 288 CALL ice_drift_init ! initialization for diags of conservation 284 289 ! 285 290 fr_i (:,:) = at_i(:,:) ! initialisation of sea-ice fraction … … 342 347 ENDIF 343 348 ! 344 IF( ln_bdy .AND. ln_icediachk ) CALL ctl_warn('par_init: online conservation check does not work with BDY')345 !346 349 rDt_ice = REAL(nn_fsbc) * rn_Dt !--- sea-ice timestep and its inverse 347 350 r1_Dt_ice = 1._wp / rDt_ice … … 393 396 !! of the time step 394 397 !!---------------------------------------------------------------------- 395 INTEGER :: ji, jj ! dummy loop index 396 !!---------------------------------------------------------------------- 397 sfx (:,:) = 0._wp ; 398 sfx_bri(:,:) = 0._wp ; sfx_lam(:,:) = 0._wp 399 sfx_sni(:,:) = 0._wp ; sfx_opw(:,:) = 0._wp 400 sfx_bog(:,:) = 0._wp ; sfx_dyn(:,:) = 0._wp 401 sfx_bom(:,:) = 0._wp ; sfx_sum(:,:) = 0._wp 402 sfx_res(:,:) = 0._wp ; sfx_sub(:,:) = 0._wp 403 ! 404 wfx_snw(:,:) = 0._wp ; wfx_ice(:,:) = 0._wp 405 wfx_sni(:,:) = 0._wp ; wfx_opw(:,:) = 0._wp 406 wfx_bog(:,:) = 0._wp ; wfx_dyn(:,:) = 0._wp 407 wfx_bom(:,:) = 0._wp ; wfx_sum(:,:) = 0._wp 408 wfx_res(:,:) = 0._wp ; wfx_sub(:,:) = 0._wp 409 wfx_spr(:,:) = 0._wp ; wfx_lam(:,:) = 0._wp 410 wfx_snw_dyn(:,:) = 0._wp ; wfx_snw_sum(:,:) = 0._wp 411 wfx_snw_sub(:,:) = 0._wp ; wfx_ice_sub(:,:) = 0._wp 412 wfx_snw_sni(:,:) = 0._wp 413 wfx_pnd(:,:) = 0._wp 414 415 hfx_thd(:,:) = 0._wp ; 416 hfx_snw(:,:) = 0._wp ; hfx_opw(:,:) = 0._wp 417 hfx_bog(:,:) = 0._wp ; hfx_dyn(:,:) = 0._wp 418 hfx_bom(:,:) = 0._wp ; hfx_sum(:,:) = 0._wp 419 hfx_res(:,:) = 0._wp ; hfx_sub(:,:) = 0._wp 420 hfx_spr(:,:) = 0._wp ; hfx_dif(:,:) = 0._wp 421 hfx_err_dif(:,:) = 0._wp 422 wfx_err_sub(:,:) = 0._wp 423 ! 424 diag_heat(:,:) = 0._wp ; diag_sice(:,:) = 0._wp 425 diag_vice(:,:) = 0._wp ; diag_vsnw(:,:) = 0._wp 426 427 ! SIMIP diagnostics 428 qcn_ice_bot(:,:,:) = 0._wp ; qcn_ice_top(:,:,:) = 0._wp ! conductive fluxes 429 t_si (:,:,:) = rt0 ! temp at the ice-snow interface 430 431 tau_icebfr (:,:) = 0._wp ! landfast ice param only (clem: important to keep the init here) 432 cnd_ice (:,:,:) = 0._wp ! initialisation: effective conductivity at the top of ice/snow (ln_cndflx=T) 433 qcn_ice (:,:,:) = 0._wp ! initialisation: conductive flux (ln_cndflx=T & ln_cndemule=T) 434 qtr_ice_bot(:,:,:) = 0._wp ! initialization: part of solar radiation transmitted through the ice needed at least for outputs 435 qsb_ice_bot(:,:) = 0._wp ! (needed if ln_icethd=F) 436 ! 437 ! for control checks (ln_icediachk) 438 diag_trp_vi(:,:) = 0._wp ; diag_trp_vs(:,:) = 0._wp 439 diag_trp_ei(:,:) = 0._wp ; diag_trp_es(:,:) = 0._wp 440 diag_trp_sv(:,:) = 0._wp 441 398 INTEGER :: ji, jj, jl ! dummy loop index 399 !!---------------------------------------------------------------------- 400 401 DO_2D( 1, 1, 1, 1 ) 402 sfx (ji,jj) = 0._wp ; 403 sfx_bri(ji,jj) = 0._wp ; sfx_lam(ji,jj) = 0._wp 404 sfx_sni(ji,jj) = 0._wp ; sfx_opw(ji,jj) = 0._wp 405 sfx_bog(ji,jj) = 0._wp ; sfx_dyn(ji,jj) = 0._wp 406 sfx_bom(ji,jj) = 0._wp ; sfx_sum(ji,jj) = 0._wp 407 sfx_res(ji,jj) = 0._wp ; sfx_sub(ji,jj) = 0._wp 408 ! 409 wfx_snw(ji,jj) = 0._wp ; wfx_ice(ji,jj) = 0._wp 410 wfx_sni(ji,jj) = 0._wp ; wfx_opw(ji,jj) = 0._wp 411 wfx_bog(ji,jj) = 0._wp ; wfx_dyn(ji,jj) = 0._wp 412 wfx_bom(ji,jj) = 0._wp ; wfx_sum(ji,jj) = 0._wp 413 wfx_res(ji,jj) = 0._wp ; wfx_sub(ji,jj) = 0._wp 414 wfx_spr(ji,jj) = 0._wp ; wfx_lam(ji,jj) = 0._wp 415 wfx_snw_dyn(ji,jj) = 0._wp ; wfx_snw_sum(ji,jj) = 0._wp 416 wfx_snw_sub(ji,jj) = 0._wp ; wfx_ice_sub(ji,jj) = 0._wp 417 wfx_snw_sni(ji,jj) = 0._wp 418 wfx_pnd(ji,jj) = 0._wp 419 420 hfx_thd(ji,jj) = 0._wp ; 421 hfx_snw(ji,jj) = 0._wp ; hfx_opw(ji,jj) = 0._wp 422 hfx_bog(ji,jj) = 0._wp ; hfx_dyn(ji,jj) = 0._wp 423 hfx_bom(ji,jj) = 0._wp ; hfx_sum(ji,jj) = 0._wp 424 hfx_res(ji,jj) = 0._wp ; hfx_sub(ji,jj) = 0._wp 425 hfx_spr(ji,jj) = 0._wp ; hfx_dif(ji,jj) = 0._wp 426 hfx_err_dif(ji,jj) = 0._wp 427 wfx_err_sub(ji,jj) = 0._wp 428 ! 429 diag_heat(ji,jj) = 0._wp ; diag_sice(ji,jj) = 0._wp 430 diag_vice(ji,jj) = 0._wp ; diag_vsnw(ji,jj) = 0._wp 431 432 tau_icebfr (ji,jj) = 0._wp ! landfast ice param only (clem: important to keep the init here) 433 qsb_ice_bot(ji,jj) = 0._wp ! (needed if ln_icethd=F) 434 435 fhld(ji,jj) = 0._wp ! needed if ln_icethd=F 436 437 ! for control checks (ln_icediachk) 438 diag_trp_vi(ji,jj) = 0._wp ; diag_trp_vs(ji,jj) = 0._wp 439 diag_trp_ei(ji,jj) = 0._wp ; diag_trp_es(ji,jj) = 0._wp 440 diag_trp_sv(ji,jj) = 0._wp 441 ! 442 diag_adv_mass(ji,jj) = 0._wp 443 diag_adv_salt(ji,jj) = 0._wp 444 diag_adv_heat(ji,jj) = 0._wp 445 END_2D 446 447 DO jl = 1, jpl 448 DO_2D( 1, 1, 1, 1 ) 449 ! SIMIP diagnostics 450 t_si (ji,jj,jl) = rt0 ! temp at the ice-snow interface 451 qcn_ice_bot(ji,jj,jl) = 0._wp 452 qcn_ice_top(ji,jj,jl) = 0._wp ! conductive fluxes 453 cnd_ice (ji,jj,jl) = 0._wp ! effective conductivity at the top of ice/snow (ln_cndflx=T) 454 qcn_ice (ji,jj,jl) = 0._wp ! conductive flux (ln_cndflx=T & ln_cndemule=T) 455 qtr_ice_bot(ji,jj,jl) = 0._wp ! part of solar radiation transmitted through the ice needed at least for outputs 456 END_2D 457 ENDDO 458 442 459 END SUBROUTINE diag_set0 460 461 462 SUBROUTINE diag_trends( kn ) 463 !!---------------------------------------------------------------------- 464 !! *** ROUTINE diag_trends *** 465 !! 466 !! ** purpose : diagnostics of the trends. Used for conservation purposes 467 !! and outputs 468 !!---------------------------------------------------------------------- 469 INTEGER, INTENT(in) :: kn ! 1 = after dyn ; 2 = after thermo 470 !!---------------------------------------------------------------------- 471 ! 472 ! --- trends of heat, salt, mass (used for conservation controls) 473 IF( ln_icediachk .OR. iom_use('hfxdhc') ) THEN 474 ! 475 diag_heat(:,:) = diag_heat(:,:) & 476 & - SUM(SUM( e_i (:,:,1:nlay_i,:) - e_i_b (:,:,1:nlay_i,:), dim=4 ), dim=3 ) * r1_Dt_ice & 477 & - SUM(SUM( e_s (:,:,1:nlay_s,:) - e_s_b (:,:,1:nlay_s,:), dim=4 ), dim=3 ) * r1_Dt_ice 478 diag_sice(:,:) = diag_sice(:,:) & 479 & + SUM( sv_i(:,:,:) - sv_i_b(:,:,:) , dim=3 ) * r1_Dt_ice * rhoi 480 diag_vice(:,:) = diag_vice(:,:) & 481 & + SUM( v_i (:,:,:) - v_i_b (:,:,:) , dim=3 ) * r1_Dt_ice * rhoi 482 diag_vsnw(:,:) = diag_vsnw(:,:) & 483 & + SUM( v_s (:,:,:) - v_s_b (:,:,:) , dim=3 ) * r1_Dt_ice * rhos 484 ! 485 IF( kn == 2 ) CALL iom_put ( 'hfxdhc' , diag_heat ) ! output of heat trend 486 ! 487 ENDIF 488 ! 489 ! --- trends of concentration (used for simip outputs) 490 IF( iom_use('afxdyn') .OR. iom_use('afxthd') .OR. iom_use('afxtot') ) THEN 491 ! 492 diag_aice(:,:) = diag_aice(:,:) + SUM( a_i(:,:,:) - a_i_b(:,:,:), dim=3 ) * r1_Dt_ice 493 ! 494 IF( kn == 1 ) CALL iom_put( 'afxdyn' , diag_aice ) ! dyn trend 495 IF( kn == 2 ) CALL iom_put( 'afxthd' , SUM( a_i(:,:,:) - a_i_b(:,:,:), dim=3 ) * r1_Dt_ice ) ! thermo trend 496 IF( kn == 2 ) CALL iom_put( 'afxtot' , diag_aice ) ! total trend 497 ! 498 ENDIF 499 ! 500 END SUBROUTINE diag_trends 443 501 444 502 #else -
NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/src/ICE/icethd.F90
r13553 r13741 18 18 USE ice ! sea-ice: variables 19 19 !!gm list trop longue ==>>> why not passage en argument d'appel ? 20 USE sbc_oce , ONLY : sss_m, sst_m, e3t_m, utau, vtau, ssu_m, ssv_m, frq_m, qns_tot, qsr_tot,sprecip, ln_cpl20 USE sbc_oce , ONLY : sss_m, sst_m, e3t_m, utau, vtau, ssu_m, ssv_m, frq_m, sprecip, ln_cpl 21 21 USE sbc_ice , ONLY : qsr_oce, qns_oce, qemp_oce, qsr_ice, qns_ice, dqns_ice, evap_ice, qprec_ice, qevap_ice, & 22 22 & qml_ice, qcn_ice, qtr_ice_top … … 30 30 USE icethd_pnd ! sea-ice: melt ponds 31 31 USE iceitd ! sea-ice: remapping thickness distribution 32 USE icecor ! sea-ice: corrections 32 33 USE icetab ! sea-ice: 1D <==> 2D transformation 33 34 USE icevar ! sea-ice: operations … … 52 53 LOGICAL :: ln_icedO ! activate ice growth in open-water (T) or not (F) 53 54 LOGICAL :: ln_icedS ! activate gravity drainage and flushing (T) or not (F) 54 LOGICAL :: ln_leadhfx ! 55 LOGICAL :: ln_leadhfx ! heat in the leads is used to melt sea-ice before warming the ocean 55 56 56 57 !! for convergence tests … … 91 92 ! 92 93 INTEGER :: ji, jj, jk, jl ! dummy loop indices 93 REAL(wp) :: zfric_u, zqld, zqfr, zqfr_neg 94 REAL(wp), PARAMETER :: zfric_umin = 0._wp 95 REAL(wp), PARAMETER :: zch = 0.0057_wp 96 REAL(wp), DIMENSION(jpi,jpj) :: zu_io, zv_io, zfric ! ice-ocean velocity (m/s) and frictional velocity (m2/s2)94 REAL(wp) :: zfric_u, zqld, zqfr, zqfr_neg, zqfr_pos 95 REAL(wp), PARAMETER :: zfric_umin = 0._wp ! lower bound for the friction velocity (cice value=5.e-04) 96 REAL(wp), PARAMETER :: zch = 0.0057_wp ! heat transfer coefficient 97 REAL(wp), DIMENSION(jpi,jpj) :: zu_io, zv_io, zfric, zvel ! ice-ocean velocity (m/s) and frictional velocity (m2/s2) 97 98 ! 98 99 !!------------------------------------------------------------------- … … 124 125 & ( zu_io(ji,jj) * zu_io(ji,jj) + zu_io(ji-1,jj) * zu_io(ji-1,jj) & 125 126 & + zv_io(ji,jj) * zv_io(ji,jj) + zv_io(ji,jj-1) * zv_io(ji,jj-1) ) ) * tmask(ji,jj,1) 127 zvel(ji,jj) = 0.5_wp * SQRT( ( u_ice(ji-1,jj) + u_ice(ji,jj) ) * ( u_ice(ji-1,jj) + u_ice(ji,jj) ) + & 128 & ( v_ice(ji,jj-1) + v_ice(ji,jj) ) * ( v_ice(ji,jj-1) + v_ice(ji,jj) ) ) 126 129 END_2D 127 130 ELSE ! if no ice dynamics => transmit directly the atmospheric stress to the ocean … … 130 133 & ( utau(ji,jj) * utau(ji,jj) + utau(ji-1,jj) * utau(ji-1,jj) & 131 134 & + vtau(ji,jj) * vtau(ji,jj) + vtau(ji,jj-1) * vtau(ji,jj-1) ) ) * tmask(ji,jj,1) 135 zvel(ji,jj) = 0._wp 132 136 END_2D 133 137 ENDIF 134 CALL lbc_lnk ( 'icethd', zfric, 'T',1.0_wp )138 CALL lbc_lnk_multi( 'icethd', zfric, 'T', 1.0_wp, zvel, 'T', 1.0_wp ) 135 139 ! 136 140 !--------------------------------------------------------------------! … … 140 144 rswitch = tmask(ji,jj,1) * MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi10 ) ) ! 0 if no ice 141 145 ! 142 ! ! solar irradiance transmission at the mixed layer bottom and used in the lead heat budget143 ! ! practically no "direct lateral ablation"144 !145 ! ! net downward heat flux from the ice to the ocean, expressed as a function of ocean146 ! ! temperature and turbulent mixing (McPhee, 1992)147 !148 146 ! --- Energy received in the lead from atm-oce exchanges, zqld is defined everywhere (J.m-2) --- ! 149 147 zqld = tmask(ji,jj,1) * rDt_ice * & … … 151 149 & ( 1._wp - at_i_b(ji,jj) ) * qns_oce(ji,jj) + qemp_oce(ji,jj) ) 152 150 153 ! --- Energy needed to bring ocean surface layer until its freezing (mostly<0 but >0 if supercooling, J.m-2) --- ! 151 ! --- Energy needed to bring ocean surface layer until its freezing, zqfr is defined everywhere (J.m-2) --- ! 152 ! (mostly<0 but >0 if supercooling) 154 153 zqfr = rho0 * rcp * e3t_m(ji,jj) * ( t_bo(ji,jj) - ( sst_m(ji,jj) + rt0 ) ) * tmask(ji,jj,1) ! both < 0 (t_bo < sst) and > 0 (t_bo > sst) 155 154 zqfr_neg = MIN( zqfr , 0._wp ) ! only < 0 156 157 ! --- Sensible ocean-to-ice heat flux (mostly>0 but <0 if supercooling, W/m2) 155 zqfr_pos = MAX( zqfr , 0._wp ) ! only > 0 156 157 ! --- Sensible ocean-to-ice heat flux (W/m2) --- ! 158 ! (mostly>0 but <0 if supercooling) 158 159 zfric_u = MAX( SQRT( zfric(ji,jj) ), zfric_umin ) 159 qsb_ice_bot(ji,jj) = rswitch * rho0 * rcp * zch * zfric_u * ( ( sst_m(ji,jj) + rt0 ) - t_bo(ji,jj) ) ! W.m-2 160 161 qsb_ice_bot(ji,jj) = rswitch * MIN( qsb_ice_bot(ji,jj), - zqfr_neg * r1_Dt_ice / MAX( at_i(ji,jj), epsi10 ) ) 160 qsb_ice_bot(ji,jj) = rswitch * rho0 * rcp * zch * zfric_u * ( ( sst_m(ji,jj) + rt0 ) - t_bo(ji,jj) ) 161 162 162 ! upper bound for qsb_ice_bot: the heat retrieved from the ocean must be smaller than the heat necessary to reach 163 163 ! the freezing point, so that we do not have SST < T_freeze 164 ! This implies: - ( qsb_ice_bot(ji,jj) * at_i(ji,jj) * rtdice ) - zqfr >= 0 165 166 !-- Energy Budget of the leads (J.m-2), source of ice growth in open water. Must be < 0 to form ice 167 qlead(ji,jj) = MIN( 0._wp , zqld - ( qsb_ice_bot(ji,jj) * at_i(ji,jj) * rDt_ice ) - zqfr ) 168 169 ! If there is ice and leads are warming => transfer energy from the lead budget and use it for bottom melting 170 ! If the grid cell is fully covered by ice (no leads) => transfer energy from the lead budget to the ice bottom budget 171 IF( ( zqld >= 0._wp .AND. at_i(ji,jj) > 0._wp ) .OR. at_i(ji,jj) >= (1._wp - epsi10) ) THEN 172 IF( ln_leadhfx ) THEN ; fhld(ji,jj) = rswitch * zqld * r1_Dt_ice / MAX( at_i(ji,jj), epsi10 ) ! divided by at_i since this is (re)multiplied by a_i in icethd_dh.F90 173 ELSE ; fhld(ji,jj) = 0._wp 174 ENDIF 164 ! This implies: qsb_ice_bot(ji,jj) * at_i(ji,jj) * rtdice <= - zqfr_neg 165 ! The following formulation is ok for both normal conditions and supercooling 166 qsb_ice_bot(ji,jj) = rswitch * MIN( qsb_ice_bot(ji,jj), - zqfr_neg * r1_Dt_ice / MAX( at_i(ji,jj), epsi10 ) ) 167 168 ! --- Energy Budget of the leads (qlead, J.m-2) --- ! 169 ! qlead is the energy received from the atm. in the leads. 170 ! If warming (zqld >= 0), then the energy in the leads is used to melt ice (bottom melting) => fhld (W/m2) 171 ! If cooling (zqld < 0), then the energy in the leads is used to grow ice in open water => qlead (J.m-2) 172 IF( zqld >= 0._wp .AND. at_i(ji,jj) > 0._wp ) THEN 173 ! upper bound for fhld: fhld should be equal to zqld 174 ! but we have to make sure that this heat will not make the sst drop below the freezing point 175 ! so the max heat that can be pulled out of the ocean is zqld - qsb - zqfr_pos 176 ! The following formulation is ok for both normal conditions and supercooling 177 fhld (ji,jj) = rswitch * MAX( 0._wp, ( zqld - zqfr_pos ) * r1_Dt_ice / MAX( at_i(ji,jj), epsi10 ) & ! divided by at_i since this is (re)multiplied by a_i in icethd_dh.F90 178 & - qsb_ice_bot(ji,jj) ) 175 179 qlead(ji,jj) = 0._wp 176 180 ELSE 177 181 fhld (ji,jj) = 0._wp 182 ! upper bound for qlead: qlead should be equal to zqld 183 ! but before using this heat for ice formation, we suppose that the ocean cools down till the freezing point. 184 ! The energy for this cooling down is zqfr. Also some heat will be removed from the ocean from turbulent fluxes (qsb) 185 ! and freezing point is reached if zqfr = zqld - qsb*a/dt 186 ! so the max heat that can be pulled out of the ocean is zqld - qsb - zqfr 187 ! The following formulation is ok for both normal conditions and supercooling 188 qlead(ji,jj) = MIN( 0._wp , zqld - ( qsb_ice_bot(ji,jj) * at_i(ji,jj) * rDt_ice ) - zqfr ) 178 189 ENDIF 179 190 ! 180 ! Net heat flux on top of the ice-ocean [W.m-2] 181 ! --------------------------------------------- 182 qt_atm_oi(ji,jj) = qns_tot(ji,jj) + qsr_tot(ji,jj) 191 ! If ice is landfast and ice concentration reaches its max 192 ! => stop ice formation in open water 193 IF( zvel(ji,jj) <= 5.e-04_wp .AND. at_i(ji,jj) >= rn_amax_2d(ji,jj)-epsi06 ) qlead(ji,jj) = 0._wp 194 ! 195 ! If the grid cell is almost fully covered by ice (no leads) 196 ! => stop ice formation in open water 197 IF( at_i(ji,jj) >= (1._wp - epsi10) ) qlead(ji,jj) = 0._wp 198 ! 199 ! If ln_leadhfx is false 200 ! => do not use energy of the leads to melt sea-ice 201 IF( .NOT.ln_leadhfx ) fhld(ji,jj) = 0._wp 202 ! 183 203 END_2D 184 204 … … 191 211 ENDIF 192 212 193 ! ---------------------------------------------------------------------194 ! Net heat flux on top of the ocean after ice thermo (1st step) [W.m-2]195 ! ---------------------------------------------------------------------196 ! First step here : non solar + precip - qlead - qsensible197 ! Second step in icethd_dh : heat remaining if total melt (zq_rema)198 ! Third step in iceupdate.F90 : heat from ice-ocean mass exchange (zf_mass) + solar199 qt_oce_ai(:,:) = ( 1._wp - at_i_b(:,:) ) * qns_oce(:,:) + qemp_oce(:,:) & ! Non solar heat flux received by the ocean200 & - qlead(:,:) * r1_Dt_ice & ! heat flux taken from the ocean where there is open water ice formation201 & - at_i (:,:) * qsb_ice_bot(:,:) & ! heat flux taken by sensible flux202 & - at_i (:,:) * fhld (:,:) ! heat flux taken during bottom growth/melt203 ! ! (fhld should be 0 while bott growth)204 213 !-------------------------------------------------------------------------------------------! 205 214 ! Thermodynamic computation (only on grid points covered by ice) => loop over ice categories … … 255 264 IF( ln_icedO ) CALL ice_thd_do ! --- Frazil ice growth in leads --- ! 256 265 ! 266 CALL ice_cor( kt , 2 ) ! --- Corrections --- ! 267 ! 268 oa_i(:,:,:) = oa_i(:,:,:) + a_i(:,:,:) * rdt_ice ! ice natural aging incrementation 269 ! 257 270 ! convergence tests 258 271 IF( ln_zdf_chkcvg ) THEN … … 418 431 CALL tab_2d_1d( npti, nptidx(1:npti), hfx_res_1d (1:npti), hfx_res ) 419 432 CALL tab_2d_1d( npti, nptidx(1:npti), hfx_err_dif_1d(1:npti), hfx_err_dif ) 420 CALL tab_2d_1d( npti, nptidx(1:npti), qt_oce_ai_1d (1:npti), qt_oce_ai )421 433 ! 422 434 ! ocean surface fields 423 435 CALL tab_2d_1d( npti, nptidx(1:npti), sst_1d(1:npti), sst_m ) 424 436 CALL tab_2d_1d( npti, nptidx(1:npti), sss_1d(1:npti), sss_m ) 437 CALL tab_2d_1d( npti, nptidx(1:npti), frq_m_1d(1:npti), frq_m ) 425 438 ! 426 439 ! to update ice age … … 510 523 CALL tab_1d_2d( npti, nptidx(1:npti), hfx_res_1d (1:npti), hfx_res ) 511 524 CALL tab_1d_2d( npti, nptidx(1:npti), hfx_err_dif_1d(1:npti), hfx_err_dif ) 512 CALL tab_1d_2d( npti, nptidx(1:npti), qt_oce_ai_1d (1:npti), qt_oce_ai )513 525 ! 514 526 CALL tab_1d_2d( npti, nptidx(1:npti), qns_ice_1d (1:npti), qns_ice (:,:,kl) ) -
NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/src/ICE/icethd_dh.F90
r13553 r13741 139 139 ! 140 140 DO ji = 1, npti 141 zf_tt(ji) = qcn_ice_bot_1d(ji) + qsb_ice_bot_1d(ji) + fhld_1d(ji) 141 zf_tt(ji) = qcn_ice_bot_1d(ji) + qsb_ice_bot_1d(ji) + fhld_1d(ji) + qtr_ice_bot_1d(ji) * frq_m_1d(ji) 142 142 zq_bot(ji) = MAX( 0._wp, zf_tt(ji) * rDt_ice ) 143 143 END DO … … 556 556 ! 557 557 ! Remaining heat flux (W.m-2) is sent to the ocean heat budget 558 qt_oce_ai_1d(ji) = qt_oce_ai_1d(ji) + ( zq_rema(ji) * a_i_1d(ji) ) * r1_Dt_ice558 !!hfx_res_1d(ji) = hfx_res_1d(ji) + ( zq_rema(ji) * a_i_1d(ji) ) * r1_Dt_ice 559 559 560 560 IF( ln_icectl .AND. zq_rema(ji) < 0. .AND. lwp ) WRITE(numout,*) 'ALERTE zq_rema <0 = ', zq_rema(ji) -
NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/src/ICE/icethd_do.F90
r13553 r13741 131 131 132 132 ! Default new ice thickness 133 WHERE( qlead(:,:) < 0._wp .AND. tau_icebfr(:,:) == 0._wp ) ; ht_i_new(:,:) = rn_hinew ! if cooling and no landfast 134 ELSEWHERE ; ht_i_new(:,:) = 0._wp 133 WHERE( qlead(:,:) < 0._wp ) ! cooling 134 ht_i_new(:,:) = rn_hinew 135 ELSEWHERE 136 ht_i_new(:,:) = 0._wp 135 137 END WHERE 136 138 … … 146 148 ! 147 149 DO_2D( 0, 0, 0, 0 ) 148 IF ( qlead(ji,jj) < 0._wp .AND. tau_icebfr(ji,jj) == 0._wp ) THEN ! activated if cooling and no landfast150 IF ( qlead(ji,jj) < 0._wp ) THEN ! cooling 149 151 ! -- Wind stress -- ! 150 152 ztaux = ( utau_ice(ji-1,jj ) * umask(ji-1,jj ,1) & … … 198 200 ! 2) Compute thickness, salinity, enthalpy, age, area and volume of new ice 199 201 !------------------------------------------------------------------------------! 200 ! This occurs if open water energy budget is negative (cooling) and there is no landfast ice202 ! it occurs if cooling 201 203 202 204 ! Identify grid points where new ice forms 203 205 npti = 0 ; nptidx(:) = 0 204 206 DO_2D( 1, 1, 1, 1 ) 205 IF ( qlead(ji,jj) < 0._wp .AND. tau_icebfr(ji,jj) == 0._wp) THEN207 IF ( qlead(ji,jj) < 0._wp ) THEN 206 208 npti = npti + 1 207 209 nptidx( npti ) = (jj - 1) * jpi + ji -
NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/src/ICE/iceupdate.F90
r13553 r13741 24 24 USE traqsr ! add penetration of solar flux in the calculation of heat budget 25 25 USE icectl ! sea-ice: control prints 26 USE bdy_oce , ONLY : ln_bdy27 26 USE zdfdrg , ONLY : ln_drgice_imp 28 27 ! … … 92 91 ! 93 92 INTEGER :: ji, jj, jl, jk ! dummy loop indices 94 REAL(wp) :: zqmass ! Heat flux associated with mass exchange ice->ocean (W.m-2)95 93 REAL(wp) :: zqsr ! New solar flux received by the ocean 96 94 REAL(wp), DIMENSION(jpi,jpj) :: z2d ! 2D workspace … … 103 101 WRITE(numout,*)'~~~~~~~~~~~~~~' 104 102 ENDIF 103 104 ! Net heat flux on top of the ice-ocean (W.m-2) 105 !---------------------------------------------- 106 qt_atm_oi(:,:) = qns_tot(:,:) + qsr_tot(:,:) 105 107 106 108 ! --- case we bypass ice thermodynamics --- ! … … 115 117 DO_2D( 1, 1, 1, 1 ) 116 118 117 ! Solar heat flux reaching the ocean = zqsr (W.m-2)119 ! Solar heat flux reaching the ocean (max) = zqsr (W.m-2) 118 120 !--------------------------------------------------- 119 121 zqsr = qsr_tot(ji,jj) - SUM( a_i_b(ji,jj,:) * ( qsr_ice(ji,jj,:) - qtr_ice_bot(ji,jj,:) ) ) … … 121 123 ! Total heat flux reaching the ocean = qt_oce_ai (W.m-2) 122 124 !--------------------------------------------------- 123 zqmass = hfx_thd(ji,jj) + hfx_dyn(ji,jj) + hfx_res(ji,jj) ! heat flux from snow is 0 (T=0 degC) 124 qt_oce_ai(ji,jj) = qt_oce_ai(ji,jj) + zqmass + zqsr 125 126 ! Add the residual from heat diffusion equation and sublimation (W.m-2) 127 !---------------------------------------------------------------------- 128 qt_oce_ai(ji,jj) = qt_oce_ai(ji,jj) + hfx_err_dif(ji,jj) + & 129 & ( hfx_sub(ji,jj) - SUM( qevap_ice(ji,jj,:) * a_i_b(ji,jj,:) ) ) 130 125 qt_oce_ai(ji,jj) = qt_atm_oi(ji,jj) - hfx_sum(ji,jj) - hfx_bom(ji,jj) - hfx_bog(ji,jj) & 126 & - hfx_dif(ji,jj) - hfx_opw(ji,jj) - hfx_snw(ji,jj) & 127 & + hfx_thd(ji,jj) + hfx_dyn(ji,jj) + hfx_res(ji,jj) & 128 & + hfx_sub(ji,jj) - SUM( qevap_ice(ji,jj,:) * a_i_b(ji,jj,:) ) + hfx_spr(ji,jj) 129 131 130 ! New qsr and qns used to compute the oceanic heat flux at the next time step 132 131 !---------------------------------------------------------------------------- 133 qsr(ji,jj) = zqsr 132 ! if warming and some ice remains, then we suppose that the whole solar flux has been consumed to melt the ice 133 ! else ( cooling or no ice left ), then we suppose that no solar flux has been consumed 134 ! 135 IF( fhld(ji,jj) > 0._wp .AND. at_i(ji,jj) > 0._wp ) THEN !-- warming and some ice remains 136 ! solar flux transmitted thru the 1st level of the ocean (i.e. not used by sea-ice) 137 qsr(ji,jj) = ( 1._wp - at_i_b(ji,jj) ) * qsr_oce(ji,jj) * ( 1._wp - frq_m(ji,jj) ) & 138 ! + solar flux transmitted thru ice and the 1st ocean level (also not used by sea-ice) 139 & + SUM( a_i_b(ji,jj,:) * qtr_ice_bot(ji,jj,:) ) * ( 1._wp - frq_m(ji,jj) ) 140 ! 141 ELSE !-- cooling or no ice left 142 qsr(ji,jj) = zqsr 143 ENDIF 144 ! 145 ! the non-solar is simply derived from the solar flux 134 146 qns(ji,jj) = qt_oce_ai(ji,jj) - zqsr 135 147 136 148 ! Mass flux at the atm. surface 137 149 !----------------------------------- … … 140 152 ! Mass flux at the ocean surface 141 153 !------------------------------------ 142 ! case of realistic freshwater flux (Tartinville et al., 2001) (presently ACTIVATED) 143 ! ------------------------------------------------------------------------------------- 144 ! The idea of this approach is that the system that we consider is the ICE-OCEAN system 145 ! Thus FW flux = External ( E-P+snow melt) 146 ! Salt flux = Exchanges in the ice-ocean system then converted into FW 147 ! Associated to Ice formation AND Ice melting 148 ! Even if i see Ice melting as a FW and SALT flux 149 ! 150 ! mass flux from ice/ocean 154 ! ice-ocean mass flux 151 155 wfx_ice(ji,jj) = wfx_bog(ji,jj) + wfx_bom(ji,jj) + wfx_sum(ji,jj) + wfx_sni(ji,jj) & 152 156 & + wfx_opw(ji,jj) + wfx_dyn(ji,jj) + wfx_res(ji,jj) + wfx_lam(ji,jj) + wfx_pnd(ji,jj) 153 154 ! add the snow melt water to snow mass flux to the ocean157 158 ! snw-ocean mass flux 155 159 wfx_snw(ji,jj) = wfx_snw_sni(ji,jj) + wfx_snw_dyn(ji,jj) + wfx_snw_sum(ji,jj) 156 157 ! mass flux at the ocean/ice interface 158 fmmflx(ji,jj) = - ( wfx_ice(ji,jj) + wfx_snw(ji,jj) + wfx_err_sub(ji,jj) ) ! F/M mass flux save at least for biogeochemical model 159 emp(ji,jj) = emp_oce(ji,jj) - wfx_ice(ji,jj) - wfx_snw(ji,jj) - wfx_err_sub(ji,jj) ! mass flux + F/M mass flux (always ice/ocean mass exchange) 160 160 161 ! total mass flux at the ocean/ice interface 162 fmmflx(ji,jj) = - wfx_ice(ji,jj) - wfx_snw(ji,jj) - wfx_err_sub(ji,jj) ! ice-ocean mass flux saved at least for biogeochemical model 163 emp (ji,jj) = emp_oce(ji,jj) - wfx_ice(ji,jj) - wfx_snw(ji,jj) - wfx_err_sub(ji,jj) ! atm-ocean + ice-ocean mass flux 161 164 162 165 ! Salt flux at the ocean surface … … 262 265 CALL iom_put ('hfxdif' , hfx_dif ) ! heat flux used for ice temperature change 263 266 CALL iom_put ('hfxsnw' , hfx_snw ) ! heat flux used for snow melt 264 CALL iom_put ('hfxerr' , hfx_err_dif ) ! heat flux error after heat diffusion (included in qt_oce_ai)267 CALL iom_put ('hfxerr' , hfx_err_dif ) ! heat flux error after heat diffusion 265 268 266 269 ! heat fluxes associated with mass exchange (freeze/melt/precip...) … … 279 282 !--------- 280 283 #if ! defined key_agrif 281 IF( ln_icediachk .AND. .NOT. ln_bdy) CALL ice_cons_final('iceupdate') ! conservation284 IF( ln_icediachk ) CALL ice_cons_final('iceupdate') ! conservation 282 285 #endif 283 IF( ln_icectl 284 IF( sn_cfctl%l_prtctl 285 IF( ln_timing 286 IF( ln_icectl ) CALL ice_prt (kt, iiceprt, jiceprt, 3, 'Final state ice_update') ! prints 287 IF( sn_cfctl%l_prtctl ) CALL ice_prt3D ('iceupdate') ! prints 288 IF( ln_timing ) CALL timing_stop ('ice_update') ! timing 286 289 ! 287 290 END SUBROUTINE ice_update_flx -
NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/src/OCE/BDY/bdyice.F90
r13553 r13741 61 61 !!---------------------------------------------------------------------- 62 62 ! controls 63 IF( ln_timing ) CALL timing_start('bdy_ice_thd') ! timing 64 IF( ln_icediachk ) CALL ice_cons_hsm(0,'bdy_ice_thd', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation 65 IF( ln_icediachk ) CALL ice_cons2D (0,'bdy_ice_thd', diag_v, diag_s, diag_t, diag_fv, diag_fs, diag_ft) ! conservation 63 IF( ln_timing ) CALL timing_start('bdy_ice_thd') ! timing 66 64 ! 67 65 CALL ice_var_glo2eqv … … 110 108 ! 111 109 ! controls 112 IF( ln_icectl ) CALL ice_prt ( kt, iiceprt, jiceprt, 1, ' - ice thermo bdy - ' ) ! prints 113 IF( ln_icediachk ) CALL ice_cons_hsm(1,'bdy_ice_thd', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation 114 IF( ln_icediachk ) CALL ice_cons2D (1,'bdy_ice_thd', diag_v, diag_s, diag_t, diag_fv, diag_fs, diag_ft) ! conservation 115 IF( ln_timing ) CALL timing_stop ('bdy_ice_thd') ! timing 110 IF( ln_icectl ) CALL ice_prt ( kt, iiceprt, jiceprt, 1, ' - ice thermo bdy - ' ) ! prints 111 IF( ln_timing ) CALL timing_stop ('bdy_ice_thd') ! timing 116 112 ! 117 113 END SUBROUTINE bdy_ice -
NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/src/OCE/DIA/diaptr.F90
r13552 r13741 42 42 END INTERFACE 43 43 44 PUBLIC ptr_sj ! call by tra_ldf & tra_adv routines45 PUBLIC ptr_sjk !46 PUBLIC dia_ptr_init ! call in memogcm47 44 PUBLIC dia_ptr ! call in step module 48 45 PUBLIC dia_ptr_hst ! called from tra_ldf/tra_adv routines 49 46 50 ! !!** namelist namptr **51 47 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: hstr_adv, hstr_ldf, hstr_eiv !: Heat/Salt TRansports(adv, diff, Bolus.) 52 48 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: hstr_ove, hstr_btr, hstr_vtr !: heat Salt TRansports(overturn, baro, merional) 53 49 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: pvtr_int, pzon_int !: Other zonal integrals 54 50 55 LOGICAL , PUBLIC :: l_diaptr !: tracers trend flag (set from namelist in trdini) 56 INTEGER, PARAMETER, PUBLIC :: nptr = 5 ! (glo, atl, pac, ind, ipc) 57 INTEGER, PARAMETER :: jp_msk = 3 58 INTEGER, PARAMETER :: jp_vtr = 4 51 LOGICAL, PUBLIC :: l_diaptr !: tracers trend flag 52 INTEGER, PARAMETER :: jp_msk = 3 53 INTEGER, PARAMETER :: jp_vtr = 4 59 54 60 55 REAL(wp) :: rc_sv = 1.e-6_wp ! conversion from m3/s to Sverdrup … … 65 60 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: btmsk34 ! mask out Southern Ocean (=0 south of 34°S) 66 61 67 LOGICAL :: ll_init = .TRUE. !: tracers trend flag (set from namelist in trdini)68 62 LOGICAL :: ll_init = .TRUE. !: tracers trend flag 63 69 64 !! * Substitutions 70 65 # include "do_loop_substitute.h90" … … 89 84 IF( ln_timing ) CALL timing_start('dia_ptr') 90 85 91 IF( kt == nit000 .AND. ll_init ) CALL dia_ptr_init 86 IF( kt == nit000 .AND. ll_init ) CALL dia_ptr_init ! -> will define l_diaptr and nbasin 92 87 ! 93 88 IF( l_diaptr ) THEN … … 123 118 ! 124 119 !overturning calculation 125 REAL(wp), DIMENSION(jpj,jpk,nptr) :: sjk, r1_sjk, v_msf ! i-mean i-k-surface and its inverse 126 REAL(wp), DIMENSION(jpj,jpk,nptr) :: zt_jk, zs_jk ! i-mean T and S, j-Stream-Function 127 128 REAL(wp), DIMENSION(jpi,jpj,jpk,nptr) :: z4d1, z4d2 129 REAL(wp), DIMENSION(jpi,jpj,nptr) :: z3dtr ! i-mean T and S, j-Stream-Function 130 !!---------------------------------------------------------------------- 131 120 REAL(wp), DIMENSION(:,:,: ), ALLOCATABLE :: sjk, r1_sjk, v_msf ! i-mean i-k-surface and its inverse 121 REAL(wp), DIMENSION(:,:,: ), ALLOCATABLE :: zt_jk, zs_jk ! i-mean T and S, j-Stream-Function 122 123 REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: z4d1, z4d2 124 REAL(wp), DIMENSION(:,:,: ), ALLOCATABLE :: z3dtr 125 !!---------------------------------------------------------------------- 126 ! 127 ALLOCATE( z3dtr(jpi,jpj,nbasin) ) 128 ! 132 129 IF( PRESENT( pvtr ) ) THEN 133 130 IF( iom_use( 'zomsf' ) ) THEN ! effective MSF 134 DO jn = 1, nptr ! by sub-basins 131 ALLOCATE( z4d1(jpi,jpj,jpk,nbasin) ) 132 ! 133 DO jn = 1, nbasin ! by sub-basins 135 134 z4d1(1,:,:,jn) = pvtr_int(:,:,jp_vtr,jn) ! zonal cumulative effective transport excluding closed seas 136 135 DO jk = jpkm1, 1, -1 … … 142 141 END DO 143 142 CALL iom_put( 'zomsf', z4d1 * rc_sv ) 143 ! 144 DEALLOCATE( z4d1 ) 144 145 ENDIF 145 146 IF( iom_use( 'sopstove' ) .OR. iom_use( 'sophtove' ) ) THEN 146 DO jn = 1, nptr 147 ALLOCATE( sjk(jpj,jpk,nbasin), r1_sjk(jpj,jpk,nbasin), v_msf(jpj,jpk,nbasin), & 148 & zt_jk(jpj,jpk,nbasin), zs_jk(jpj,jpk,nbasin) ) 149 ! 150 DO jn = 1, nbasin 147 151 sjk(:,:,jn) = pvtr_int(:,:,jp_msk,jn) 148 152 r1_sjk(:,:,jn) = 0._wp … … 156 160 ! 157 161 ENDDO 158 DO jn = 1, n ptr162 DO jn = 1, nbasin 159 163 z3dtr(1,:,jn) = hstr_ove(:,jp_tem,jn) * rc_pwatt ! (conversion in PW) 160 164 DO ji = 1, jpi … … 163 167 ENDDO 164 168 CALL iom_put( 'sophtove', z3dtr ) 165 DO jn = 1, n ptr169 DO jn = 1, nbasin 166 170 z3dtr(1,:,jn) = hstr_ove(:,jp_sal,jn) * rc_ggram ! (conversion in Gg) 167 171 DO ji = 1, jpi … … 170 174 ENDDO 171 175 CALL iom_put( 'sopstove', z3dtr ) 176 ! 177 DEALLOCATE( sjk, r1_sjk, v_msf, zt_jk, zs_jk ) 172 178 ENDIF 173 179 174 180 IF( iom_use( 'sopstbtr' ) .OR. iom_use( 'sophtbtr' ) ) THEN 175 181 ! Calculate barotropic heat and salt transport here 176 DO jn = 1, nptr 182 ALLOCATE( sjk(jpj,1,nbasin), r1_sjk(jpj,1,nbasin) ) 183 ! 184 DO jn = 1, nbasin 177 185 sjk(:,1,jn) = SUM( pvtr_int(:,:,jp_msk,jn), 2 ) 178 186 r1_sjk(:,1,jn) = 0._wp … … 186 194 ! 187 195 ENDDO 188 DO jn = 1, n ptr196 DO jn = 1, nbasin 189 197 z3dtr(1,:,jn) = hstr_btr(:,jp_tem,jn) * rc_pwatt ! (conversion in PW) 198 ! TODO: Change these loop indices in the next commit 190 199 DO ji = 1, jpi 191 200 z3dtr(ji,:,jn) = z3dtr(1,:,jn) … … 193 202 ENDDO 194 203 CALL iom_put( 'sophtbtr', z3dtr ) 195 DO jn = 1, n ptr204 DO jn = 1, nbasin 196 205 z3dtr(1,:,jn) = hstr_btr(:,jp_sal,jn) * rc_ggram ! (conversion in Gg) 197 206 DO ji = 1, jpi … … 200 209 ENDDO 201 210 CALL iom_put( 'sopstbtr', z3dtr ) 202 ENDIF 211 ! 212 DEALLOCATE( sjk, r1_sjk ) 213 ENDIF 203 214 ! 204 215 hstr_ove(:,:,:) = 0._wp ! Zero before next timestep … … 207 218 ELSE 208 219 IF( iom_use( 'zotem' ) .OR. iom_use( 'zosal' ) .OR. iom_use( 'zosrf' ) ) THEN ! i-mean i-k-surface 209 ! 210 DO jn = 1, nptr 220 ALLOCATE( z4d1(jpi,jpj,jpk,nbasin), z4d2(jpi,jpj,jpk,nbasin) ) 221 ! 222 DO jn = 1, nbasin 211 223 z4d1(1,:,:,jn) = pzon_int(:,:,jp_msk,jn) 212 224 DO ji = 2, jpi … … 216 228 CALL iom_put( 'zosrf', z4d1 ) 217 229 ! 218 DO jn = 1, n ptr230 DO jn = 1, nbasin 219 231 z4d2(1,:,:,jn) = pzon_int(:,:,jp_tem,jn) / MAX( z4d1(1,:,:,jn), 10.e-15 ) 220 232 DO ji = 2, jpi … … 224 236 CALL iom_put( 'zotem', z4d2 ) 225 237 ! 226 DO jn = 1, n ptr238 DO jn = 1, nbasin 227 239 z4d2(1,:,:,jn) = pzon_int(:,:,jp_sal,jn) / MAX( z4d1(1,:,:,jn), 10.e-15 ) 228 240 DO ji = 2, jpi … … 232 244 CALL iom_put( 'zosal', z4d2 ) 233 245 ! 246 DEALLOCATE( z4d1, z4d2 ) 234 247 ENDIF 235 248 ! … … 237 250 IF( iom_use( 'sophtadv' ) .OR. iom_use( 'sopstadv' ) ) THEN 238 251 ! 239 DO jn = 1, n ptr252 DO jn = 1, nbasin 240 253 z3dtr(1,:,jn) = hstr_adv(:,jp_tem,jn) * rc_pwatt ! (conversion in PW) 241 254 DO ji = 1, jpi … … 244 257 ENDDO 245 258 CALL iom_put( 'sophtadv', z3dtr ) 246 DO jn = 1, n ptr259 DO jn = 1, nbasin 247 260 z3dtr(1,:,jn) = hstr_adv(:,jp_sal,jn) * rc_ggram ! (conversion in Gg) 248 261 DO ji = 1, jpi … … 255 268 IF( iom_use( 'sophtldf' ) .OR. iom_use( 'sopstldf' ) ) THEN 256 269 ! 257 DO jn = 1, n ptr270 DO jn = 1, nbasin 258 271 z3dtr(1,:,jn) = hstr_ldf(:,jp_tem,jn) * rc_pwatt ! (conversion in PW) 259 272 DO ji = 1, jpi … … 262 275 ENDDO 263 276 CALL iom_put( 'sophtldf', z3dtr ) 264 DO jn = 1, n ptr277 DO jn = 1, nbasin 265 278 z3dtr(1,:,jn) = hstr_ldf(:,jp_sal,jn) * rc_ggram ! (conversion in Gg) 266 279 DO ji = 1, jpi … … 273 286 IF( iom_use( 'sophteiv' ) .OR. iom_use( 'sopsteiv' ) ) THEN 274 287 ! 275 DO jn = 1, n ptr288 DO jn = 1, nbasin 276 289 z3dtr(1,:,jn) = hstr_eiv(:,jp_tem,jn) * rc_pwatt ! (conversion in PW) 277 290 DO ji = 1, jpi … … 280 293 ENDDO 281 294 CALL iom_put( 'sophteiv', z3dtr ) 282 DO jn = 1, n ptr295 DO jn = 1, nbasin 283 296 z3dtr(1,:,jn) = hstr_eiv(:,jp_sal,jn) * rc_ggram ! (conversion in Gg) 284 297 DO ji = 1, jpi … … 290 303 ! 291 304 IF( iom_use( 'sopstvtr' ) .OR. iom_use( 'sophtvtr' ) ) THEN 292 DO jn = 1, n ptr305 DO jn = 1, nbasin 293 306 z3dtr(1,:,jn) = hstr_vtr(:,jp_tem,jn) * rc_pwatt ! (conversion in PW) 294 307 DO ji = 1, jpi … … 297 310 ENDDO 298 311 CALL iom_put( 'sophtvtr', z3dtr ) 299 DO jn = 1, n ptr312 DO jn = 1, nbasin 300 313 z3dtr(1,:,jn) = hstr_vtr(:,jp_sal,jn) * rc_ggram ! (conversion in Gg) 301 314 DO ji = 1, jpi … … 322 335 pzon_int(:,:,:,:) = 0._wp 323 336 ENDIF 337 ! 338 DEALLOCATE( z3dtr ) 339 ! 324 340 END SUBROUTINE dia_ptr_iom 325 341 … … 339 355 INTEGER , INTENT(in) :: Kmm ! time level index 340 356 REAL(wp), DIMENSION(ST_2D(nn_hls),jpk), INTENT(in), OPTIONAL :: pvtr ! j-effective transport 341 REAL(wp), DIMENSION( ST_2D(nn_hls),jpk):: zmask ! 3D workspace342 REAL(wp), DIMENSION( ST_2D(nn_hls),jpk,jpts):: zts ! 4D workspace343 REAL(wp), DIMENSION( ST_1Dj(nn_hls),jpk,nptr):: sjk, v_msf ! Zonal sum: i-k surface area, j-effective transport344 REAL(wp), DIMENSION( ST_1Dj(nn_hls),jpk,nptr):: zt_jk, zs_jk ! Zonal sum: i-k surface area * (T, S)357 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: zmask ! 3D workspace 358 REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: zts ! 4D workspace 359 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: sjk, v_msf ! Zonal sum: i-k surface area, j-effective transport 360 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: zt_jk, zs_jk ! Zonal sum: i-k surface area * (T, S) 345 361 REAL(wp) :: zsfc, zvfc ! i-k surface area 346 362 INTEGER :: ji, jj, jk, jn ! dummy loop indices … … 350 366 ! i sum of effective j transport excluding closed seas 351 367 IF( iom_use( 'zomsf' ) .OR. iom_use( 'sopstove' ) .OR. iom_use( 'sophtove' ) ) THEN 352 DO jn = 1, nptr 368 ALLOCATE( v_msf(ST_1Dj(nn_hls),jpk,nbasin) ) 369 370 DO jn = 1, nbasin 353 371 v_msf(:,:,jn) = ptr_sjk( pvtr(:,:,:), btmsk34(:,:,jn) ) 354 372 ENDDO 355 373 356 374 CALL ptr_sum( pvtr_int(:,:,jp_vtr,:), v_msf(:,:,:) ) 375 376 DEALLOCATE( v_msf ) 357 377 ENDIF 358 378 … … 360 380 IF( iom_use( 'sopstove' ) .OR. iom_use( 'sophtove' ) .OR. & 361 381 & iom_use( 'sopstbtr' ) .OR. iom_use( 'sophtbtr' ) ) THEN 382 ALLOCATE( zmask(ST_2D(nn_hls),jpk), zts(ST_2D(nn_hls),jpk,jpts), & 383 & sjk(ST_1Dj(nn_hls),jpk,nbasin), & 384 & zt_jk(ST_1Dj(nn_hls),jpk,nbasin), zs_jk(ST_1Dj(nn_hls),jpk,nbasin) ) 385 362 386 zmask(:,:,:) = 0._wp 363 387 zts(:,:,:,:) = 0._wp … … 370 394 END_3D 371 395 372 DO jn = 1, n ptr396 DO jn = 1, nbasin 373 397 sjk(:,:,jn) = ptr_sjk( zmask(:,:,:) , btmsk(:,:,jn) ) 374 398 zt_jk(:,:,jn) = ptr_sjk( zts(:,:,:,jp_tem), btmsk(:,:,jn) ) … … 379 403 CALL ptr_sum( pvtr_int(:,:,jp_tem,:), zt_jk(:,:,:) ) 380 404 CALL ptr_sum( pvtr_int(:,:,jp_sal,:), zs_jk(:,:,:) ) 405 406 DEALLOCATE( zmask, zts, sjk, zt_jk, zs_jk ) 381 407 ENDIF 382 408 ELSE 383 409 ! i sum of j surface area - temperature/salinity product on T grid 384 410 IF( iom_use( 'zotem' ) .OR. iom_use( 'zosal' ) .OR. iom_use( 'zosrf' ) ) THEN 411 ALLOCATE( zmask(ST_2D(nn_hls),jpk), zts(ST_2D(nn_hls),jpk,jpts), & 412 & sjk(ST_1Dj(nn_hls),jpk,nbasin), & 413 & zt_jk(ST_1Dj(nn_hls),jpk,nbasin), zs_jk(ST_1Dj(nn_hls),jpk,nbasin) ) 414 385 415 zmask(:,:,:) = 0._wp 386 416 zts(:,:,:,:) = 0._wp … … 393 423 END_3D 394 424 395 DO jn = 1, n ptr425 DO jn = 1, nbasin 396 426 sjk(:,:,jn) = ptr_sjk( zmask(:,:,:) , btmsk(:,:,jn) ) 397 427 zt_jk(:,:,jn) = ptr_sjk( zts(:,:,:,jp_tem), btmsk(:,:,jn) ) … … 402 432 CALL ptr_sum( pzon_int(:,:,jp_tem,:), zt_jk(:,:,:) ) 403 433 CALL ptr_sum( pzon_int(:,:,jp_sal,:), zs_jk(:,:,:) ) 434 435 DEALLOCATE( zmask, zts, sjk, zt_jk, zs_jk ) 404 436 ENDIF 405 437 406 438 ! i-k sum of j surface area - temperature/salinity product on V grid 407 439 IF( iom_use( 'sopstvtr' ) .OR. iom_use( 'sophtvtr' ) ) THEN 440 ALLOCATE( zts(ST_2D(nn_hls),jpk,jpts) ) 441 408 442 zts(:,:,:,:) = 0._wp 409 443 … … 416 450 CALL dia_ptr_hst( jp_tem, 'vtr', zts(:,:,:,jp_tem) ) 417 451 CALL dia_ptr_hst( jp_sal, 'vtr', zts(:,:,:,jp_sal) ) 452 453 DEALLOCATE( zts ) 418 454 ENDIF 419 455 ENDIF … … 425 461 !! *** ROUTINE dia_ptr_init *** 426 462 !! 427 !! ** Purpose : Initialization , namelist read463 !! ** Purpose : Initialization 428 464 !!---------------------------------------------------------------------- 429 465 INTEGER :: inum, jn ! local integers … … 432 468 !!---------------------------------------------------------------------- 433 469 434 l_diaptr = .FALSE. 435 IF( iom_use( 'zomsf' ) .OR. iom_use( 'zotem' ) .OR. iom_use( 'zosal' ) .OR. & 436 & iom_use( 'zosrf' ) .OR. iom_use( 'sopstove' ) .OR. iom_use( 'sophtove' ) .OR. & 437 & iom_use( 'sopstbtr' ) .OR. iom_use( 'sophtbtr' ) .OR. iom_use( 'sophtadv' ) .OR. & 438 & iom_use( 'sopstadv' ) .OR. iom_use( 'sophtldf' ) .OR. iom_use( 'sopstldf' ) .OR. & 439 & iom_use( 'sophteiv' ) .OR. iom_use( 'sopsteiv' ) .OR. iom_use( 'sopstvtr' ) .OR. & 440 & iom_use( 'sophtvtr' ) .OR. iom_use( 'uocetr_vsum_cumul' ) ) l_diaptr = .TRUE. 441 470 ! l_diaptr is defined with iom_use 471 ! --> dia_ptr_init must be done after the call to iom_init 472 ! --> cannot be .TRUE. without cpp key: key_iom --> nbasin define by iom_init is initialized 473 l_diaptr = iom_use( 'zomsf' ) .OR. iom_use( 'zotem' ) .OR. iom_use( 'zosal' ) .OR. & 474 & iom_use( 'zosrf' ) .OR. iom_use( 'sopstove' ) .OR. iom_use( 'sophtove' ) .OR. & 475 & iom_use( 'sopstbtr' ) .OR. iom_use( 'sophtbtr' ) .OR. iom_use( 'sophtadv' ) .OR. & 476 & iom_use( 'sopstadv' ) .OR. iom_use( 'sophtldf' ) .OR. iom_use( 'sopstldf' ) .OR. & 477 & iom_use( 'sophteiv' ) .OR. iom_use( 'sopsteiv' ) .OR. iom_use( 'sopstvtr' ) .OR. & 478 & iom_use( 'sophtvtr' ) .OR. iom_use( 'uocetr_vsum_cumul' ) 442 479 443 480 IF(lwp) THEN ! Control print … … 445 482 WRITE(numout,*) 'dia_ptr_init : poleward transport and msf initialization' 446 483 WRITE(numout,*) '~~~~~~~~~~~~' 447 WRITE(numout,*) ' Namelist namptr : set ptr parameters'448 484 WRITE(numout,*) ' Poleward heat & salt transport (T) or not (F) l_diaptr = ', l_diaptr 449 485 ENDIF … … 452 488 ! 453 489 IF( dia_ptr_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'dia_ptr_init : unable to allocate arrays' ) 454 490 ! 455 491 rc_pwatt = rc_pwatt * rho0_rcp ! conversion from K.s-1 to PetaWatt 456 492 rc_ggram = rc_ggram * rho0 ! conversion from m3/s to Gg/s … … 458 494 IF( lk_mpp ) CALL mpp_ini_znl( numout ) ! Define MPI communicator for zonal sum 459 495 460 btmsk(:,:,:) = 0._wp 461 btmsk(:,:,1) = tmask_i(:,:) 462 CALL iom_open( 'subbasins', inum, ldstop = .FALSE. ) 463 CALL iom_get( inum, jpdom_global, 'atlmsk', btmsk(:,:,2) ) ! Atlantic basin 464 CALL iom_get( inum, jpdom_global, 'pacmsk', btmsk(:,:,3) ) ! Pacific basin 465 CALL iom_get( inum, jpdom_global, 'indmsk', btmsk(:,:,4) ) ! Indian basin 466 CALL iom_close( inum ) 467 btmsk(:,:,5) = MAX ( btmsk(:,:,3), btmsk(:,:,4) ) ! Indo-Pacific basin 468 DO jn = 2, nptr 469 btmsk(:,:,jn) = btmsk(:,:,jn) * tmask_i(:,:) ! interior domain only 496 btmsk(:,:,1) = tmask_i(:,:) 497 IF( nbasin == 5 ) THEN ! nbasin has been initialized in iom_init to define the axis "basin" 498 CALL iom_open( 'subbasins', inum ) 499 CALL iom_get( inum, jpdom_global, 'atlmsk', btmsk(:,:,2) ) ! Atlantic basin 500 CALL iom_get( inum, jpdom_global, 'pacmsk', btmsk(:,:,3) ) ! Pacific basin 501 CALL iom_get( inum, jpdom_global, 'indmsk', btmsk(:,:,4) ) ! Indian basin 502 CALL iom_close( inum ) 503 btmsk(:,:,5) = MAX ( btmsk(:,:,3), btmsk(:,:,4) ) ! Indo-Pacific basin 504 ENDIF 505 DO jn = 2, nbasin 506 btmsk(:,:,jn) = btmsk(:,:,jn) * tmask_i(:,:) ! interior domain only 470 507 END DO 471 508 ! JD : modification so that overturning streamfunction is available in Atlantic at 34S to compare with observations … … 476 513 END WHERE 477 514 btmsk34(:,:,1) = btmsk(:,:,1) 478 DO jn = 2, n ptr479 btmsk34(:,:,jn) = btmsk(:,:,jn) * zmsk(:,:) ! interior domain only515 DO jn = 2, nbasin 516 btmsk34(:,:,jn) = btmsk(:,:,jn) * zmsk(:,:) ! interior domain only 480 517 ENDDO 481 518 … … 508 545 CHARACTER(len=3) , INTENT(in) :: cptr ! transport type 'adv'/'ldf'/'eiv' 509 546 REAL(wp), DIMENSION(ST_2D(nn_hls),jpk) , INTENT(in) :: pvflx ! 3D input array of advection/diffusion 510 REAL(wp), DIMENSION(ST_1Dj(nn_hls),n ptr) :: zsj !547 REAL(wp), DIMENSION(ST_1Dj(nn_hls),nbasin) :: zsj ! 511 548 INTEGER :: jn ! 512 549 513 DO jn = 1, n ptr550 DO jn = 1, nbasin 514 551 zsj(:,jn) = ptr_sj( pvflx(:,:,:), btmsk(:,:,jn) ) 515 552 ENDDO … … 536 573 !! *** ROUTINE ptr_sum_2d *** 537 574 !!---------------------------------------------------------------------- 538 !! ** Purpose : Add two 2D arrays with (j,n ptr) dimensions575 !! ** Purpose : Add two 2D arrays with (j,nbasin) dimensions 539 576 !! 540 577 !! ** Method : - phstr = phstr + pva … … 543 580 !! ** Action : phstr 544 581 !!---------------------------------------------------------------------- 545 REAL(wp), DIMENSION(jpj,n ptr) , INTENT(inout) :: phstr !546 REAL(wp), DIMENSION(ST_1Dj(nn_hls),n ptr), INTENT(in) :: pva !582 REAL(wp), DIMENSION(jpj,nbasin) , INTENT(inout) :: phstr ! 583 REAL(wp), DIMENSION(ST_1Dj(nn_hls),nbasin), INTENT(in) :: pva ! 547 584 INTEGER :: jj 548 585 #if defined key_mpp_mpi 549 INTEGER, DIMENSION(1) :: ish1d550 INTEGER, DIMENSION(2) :: ish2d551 REAL(wp), DIMENSION(jpj*n ptr):: zwork586 INTEGER, DIMENSION(1) :: ish1d 587 INTEGER, DIMENSION(2) :: ish2d 588 REAL(wp), DIMENSION(jpj*nbasin) :: zwork 552 589 #endif 553 590 … … 558 595 #if defined key_mpp_mpi 559 596 IF( ntile == 0 .OR. ntile == nijtile ) THEN 560 ish1d(1) = jpj*n ptr561 ish2d(1) = jpj ; ish2d(2) = n ptr597 ish1d(1) = jpj*nbasin 598 ish2d(1) = jpj ; ish2d(2) = nbasin 562 599 zwork(:) = RESHAPE( phstr(:,:), ish1d ) 563 600 CALL mpp_sum( 'diaptr', zwork, ish1d(1), ncomm_znl ) … … 572 609 !! *** ROUTINE ptr_sum_3d *** 573 610 !!---------------------------------------------------------------------- 574 !! ** Purpose : Add two 3D arrays with (j,k,n ptr) dimensions611 !! ** Purpose : Add two 3D arrays with (j,k,nbasin) dimensions 575 612 !! 576 613 !! ** Method : - phstr = phstr + pva … … 579 616 !! ** Action : phstr 580 617 !!---------------------------------------------------------------------- 581 REAL(wp), DIMENSION(jpj,jpk,n ptr) , INTENT(inout) :: phstr !582 REAL(wp), DIMENSION(ST_1Dj(nn_hls),jpk,n ptr), INTENT(in) :: pva !618 REAL(wp), DIMENSION(jpj,jpk,nbasin) , INTENT(inout) :: phstr ! 619 REAL(wp), DIMENSION(ST_1Dj(nn_hls),jpk,nbasin), INTENT(in) :: pva ! 583 620 INTEGER :: jj, jk 584 621 #if defined key_mpp_mpi 585 622 INTEGER, DIMENSION(1) :: ish1d 586 623 INTEGER, DIMENSION(3) :: ish3d 587 REAL(wp), DIMENSION(jpj*jpk*n ptr) :: zwork624 REAL(wp), DIMENSION(jpj*jpk*nbasin) :: zwork 588 625 #endif 589 626 … … 596 633 #if defined key_mpp_mpi 597 634 IF( ntile == 0 .OR. ntile == nijtile ) THEN 598 ish1d(1) = jpj*jpk*n ptr599 ish3d(1) = jpj ; ish3d(2) = jpk ; ish3d(3) = n ptr635 ish1d(1) = jpj*jpk*nbasin 636 ish3d(1) = jpj ; ish3d(2) = jpk ; ish3d(3) = nbasin 600 637 zwork(:) = RESHAPE( phstr(:,:,:), ish1d ) 601 638 CALL mpp_sum( 'diaptr', zwork, ish1d(1), ncomm_znl ) … … 615 652 ierr(:) = 0 616 653 ! 654 ! nbasin has been initialized in iom_init to define the axis "basin" 655 ! 617 656 IF( .NOT. ALLOCATED( btmsk ) ) THEN 618 ALLOCATE( btmsk(jpi,jpj,n ptr) , btmsk34(jpi,jpj,nptr), &619 & hstr_adv(jpj,jpts,n ptr), hstr_eiv(jpj,jpts,nptr), &620 & hstr_ove(jpj,jpts,n ptr), hstr_btr(jpj,jpts,nptr), &621 & hstr_ldf(jpj,jpts,n ptr), hstr_vtr(jpj,jpts,nptr), STAT=ierr(1) )622 ! 623 ALLOCATE( pvtr_int(jpj,jpk,jpts+2,n ptr), &624 & pzon_int(jpj,jpk,jpts+1,n ptr), STAT=ierr(2) )657 ALLOCATE( btmsk(jpi,jpj,nbasin) , btmsk34(jpi,jpj,nbasin), & 658 & hstr_adv(jpj,jpts,nbasin), hstr_eiv(jpj,jpts,nbasin), & 659 & hstr_ove(jpj,jpts,nbasin), hstr_btr(jpj,jpts,nbasin), & 660 & hstr_ldf(jpj,jpts,nbasin), hstr_vtr(jpj,jpts,nbasin), STAT=ierr(1) ) 661 ! 662 ALLOCATE( pvtr_int(jpj,jpk,jpts+2,nbasin), & 663 & pzon_int(jpj,jpk,jpts+1,nbasin), STAT=ierr(2) ) 625 664 ! 626 665 dia_ptr_alloc = MAXVAL( ierr ) -
NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/src/OCE/DIU/diu_bulk.F90
r13295 r13741 22 22 23 23 ! Namelist parameters 24 LOGICAL, PUBLIC :: ln_diurnal 25 LOGICAL, PUBLIC :: ln_diurnal_only 24 LOGICAL, PUBLIC :: ln_diurnal = .false. ! force definition if diurnal_sst_bulk_init is not called 25 LOGICAL, PUBLIC :: ln_diurnal_only = .false. ! force definition if diurnal_sst_bulk_init is not called 26 26 27 27 ! Parameters -
NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/src/OCE/DOM/closea.F90
r13286 r13741 38 38 LOGICAL, PUBLIC :: ln_clo_rnf !: closed sea treated as runoff (update rnf mask) 39 39 40 LOGICAL, PUBLIC :: l_sbc_clo !: T => net evap/precip over closed seas spread outover the globe/river mouth 41 LOGICAL, PUBLIC :: l_clo_rnf !: T => Some closed seas output freshwater (RNF) to specified runoff points. 42 43 INTEGER, PUBLIC :: ncsg !: number of closed seas global mappings (inferred from closea_mask_glo field) 44 INTEGER, PUBLIC :: ncsr !: number of closed seas rnf mappings (inferred from closea_mask_rnf field) 45 INTEGER, PUBLIC :: ncse !: number of closed seas empmr mappings (inferred from closea_mask_emp field) 40 ! WARNING: keep default definitions in the following lines as dom_clo is called only if ln_closea = .true. 41 LOGICAL, PUBLIC :: l_sbc_clo = .FALSE. !: T => net evap/precip over closed seas spread outover the globe/river mouth 42 LOGICAL, PUBLIC :: l_clo_rnf = .FALSE. !: T => Some closed seas output freshwater (RNF) to specified runoff points. 43 44 INTEGER, PUBLIC :: ncsg = 0 !: number of closed seas global mappings (inferred from closea_mask_glo field) 45 INTEGER, PUBLIC :: ncsr = 0 !: number of closed seas rnf mappings (inferred from closea_mask_rnf field) 46 INTEGER, PUBLIC :: ncse = 0 !: number of closed seas empmr mappings (inferred from closea_mask_emp field) 46 47 47 48 INTEGER, PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:) :: mask_opnsea, mask_csundef !: mask defining the open sea and the undefined closed sea -
NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/src/OCE/DOM/daymod.F90
r13286 r13741 82 82 ndt05 = NINT( 0.5 * rn_Dt ) 83 83 84 IF( .NOT. l_offline ) CALL day_rst( nit000, 'READ' ) 85 84 lrst_oce = .NOT. l_offline ! force definition of offline 85 IF( lrst_oce ) CALL day_rst( nit000, 'READ' ) 86 86 87 ! set the calandar from ndastp (read in restart file and namelist) 87 88 nyear = ndastp / 10000 -
NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/src/OCE/DOM/dom_oce.F90
r13514 r13741 230 230 231 231 !!---------------------------------------------------------------------- 232 !! variable defined here to avoid circular dependencies... 233 !! --------------------------------------------------------------------- 234 INTEGER, PUBLIC :: nbasin ! number of basin to be considered in diaprt (glo, atl, pac, ind, ipc) 235 236 !!---------------------------------------------------------------------- 232 237 !! agrif domain 233 238 !!---------------------------------------------------------------------- -
NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/src/OCE/DOM/domain.F90
r13553 r13741 121 121 WRITE(numout,*) ' cn_cfg = ', TRIM( cn_cfg ), ' nn_cfg = ', nn_cfg 122 122 ENDIF 123 nn_wxios = 0124 ln_xios_read = .FALSE.125 123 ! 126 124 ! !== Reference coordinate system ==! -
NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/src/OCE/DYN/divhor.F90
r13553 r13741 78 78 ! 79 79 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) !== Horizontal divergence ==! 80 hdiv(ji,jj,jk) = ( e2u(ji ,jj) * e3u(ji ,jj,jk,Kmm) * uu(ji ,jj,jk,Kmm) &80 hdiv(ji,jj,jk) = ( e2u(ji ,jj) * e3u(ji ,jj,jk,Kmm) * uu(ji ,jj,jk,Kmm) & 81 81 & - e2u(ji-1,jj) * e3u(ji-1,jj,jk,Kmm) * uu(ji-1,jj,jk,Kmm) & 82 82 & + e1v(ji,jj ) * e3v(ji,jj ,jk,Kmm) * vv(ji,jj ,jk,Kmm) & -
NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/src/OCE/DYN/wet_dry.F90
r13553 r13741 57 57 REAL(wp), PUBLIC :: ssh_ref !: height of z=0 with respect to the geoid; 58 58 59 LOGICAL, PUBLIC :: ll_wd !: Wetting/drying activation switch if either ln_wd_il or ln_wd_dl59 LOGICAL, PUBLIC :: ll_wd = .FALSE. !: Wetting/drying activation switch (ln_wd_il or ln_wd_dl) <- default def if wad_init not called 60 60 61 61 PUBLIC wad_init ! initialisation routine called by step.F90 … … 111 111 112 112 r_rn_wdmin1 = 1 / rn_wdmin1 113 ll_wd = .FALSE.114 113 IF( ln_wd_il .OR. ln_wd_dl ) THEN 115 114 ll_wd = .TRUE. -
NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/src/OCE/FLO/flo_oce.F90
r11536 r13741 19 19 !! ---------------- 20 20 LOGICAL, PUBLIC :: ln_floats !: Activate floats or not 21 INTEGER, PUBLIC :: jpnfl 21 INTEGER, PUBLIC :: jpnfl = 0 !: total number of floats during the run 22 22 INTEGER, PUBLIC :: jpnnewflo !: number of floats added in a new run 23 23 INTEGER, PUBLIC :: jpnrstflo !: number of floats for the restart -
NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/src/OCE/ICB/icbtrj.F90
r13062 r13741 35 35 PUBLIC icb_trj_end ! routine called in icbstp.F90 module 36 36 37 INTEGER :: num_traj 37 INTEGER :: num_traj = 0 38 38 INTEGER :: n_dim, m_dim 39 39 INTEGER :: ntrajid -
NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/src/OCE/IOM/iom.F90
r13553 r13741 123 123 REAL(wp), DIMENSION(2,jpkam1) :: za_bnds ! ABL vertical boundaries 124 124 LOGICAL :: ll_closedef = .TRUE. 125 LOGICAL :: ll_exist 125 126 !!---------------------------------------------------------------------- 126 127 ! … … 230 231 CALL iom_set_axis_attr( "ghw_abl", bounds=za_bnds ) 231 232 232 CALL iom_set_axis_attr( "nfloat", (/ (REAL(ji,wp), ji=1,jpnfl) /) )233 CALL iom_set_axis_attr( "nfloat", (/ (REAL(ji,wp), ji=1,jpnfl) /) ) 233 234 # if defined key_si3 234 235 CALL iom_set_axis_attr( "ncatice", (/ (REAL(ji,wp), ji=1,jpl) /) ) … … 243 244 CALL iom_set_axis_attr( "iax_26C", (/ REAL(26,wp) /) ) ! strange syntaxe and idea... 244 245 CALL iom_set_axis_attr( "iax_28C", (/ REAL(28,wp) /) ) ! strange syntaxe and idea... 245 CALL iom_set_axis_attr( "basin" , (/ (REAL(ji,wp), ji=1,5) /) ) 246 ! for diaprt, we need to define an axis which size can be 1 (default) or 5 (if the file subbasins.nc exists) 247 INQUIRE( FILE = 'subbasins.nc', EXIST = ll_exist ) 248 nbasin = 1 + 4 * COUNT( (/ll_exist/) ) 249 CALL iom_set_axis_attr( "basin" , (/ (REAL(ji,wp), ji=1,nbasin) /) ) 246 250 ENDIF 247 251 ! -
NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/src/OCE/IOM/iom_def.F90
r13286 r13741 33 33 INTEGER, PUBLIC :: iom_open_init = 0 !: used to initialize iom_file(:)%nfid to 0 34 34 !XIOS write restart 35 LOGICAL, PUBLIC :: lwxios 36 INTEGER, PUBLIC :: nxioso !: type of restart file when writing using XIOS 1 - single, 2 - multiple35 LOGICAL, PUBLIC :: lwxios = .FALSE. !: write single file restart using XIOS 36 INTEGER, PUBLIC :: nxioso = 0 !: type of restart file when writing using XIOS 1 - single, 2 - multiple 37 37 !XIOS read restart 38 LOGICAL, PUBLIC :: lrxios 38 LOGICAL, PUBLIC :: lrxios = .FALSE. !: read single file restart using XIOS 39 39 LOGICAL, PUBLIC :: lxios_sini = .FALSE. ! is restart in a single file 40 40 LOGICAL, PUBLIC :: lxios_set = .FALSE. -
NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/src/OCE/ISF/isf_oce.F90
r12077 r13741 74 74 ! 75 75 ! 2.1 -------- ice shelf cavity parameter -------------- 76 LOGICAL , PUBLIC :: l_isfoasis 76 LOGICAL , PUBLIC :: l_isfoasis = .FALSE. 77 77 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: risfload !: ice shelf load 78 78 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: fwfisf_oasis -
NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/src/OCE/LBC/lib_mpp.F90
r13553 r13741 511 511 ALLOCATE(todelay(idvar)%y1d(isz)) 512 512 todelay(idvar)%y1d(:) = CMPLX(todelay(idvar)%z1d(:), 0., wp) ! create %y1d, complex variable needed by mpi_sumdd 513 ndelayid(idvar) = MPI_REQUEST_NULL ! initialised request to a valid value 513 514 END IF 514 515 ENDIF … … 518 519 ALLOCATE(todelay(idvar)%z1d(isz), todelay(idvar)%y1d(isz)) ! allocate also %z1d as used for the restart 519 520 CALL mpi_allreduce( y_in(:), todelay(idvar)%y1d(:), isz, MPI_DOUBLE_COMPLEX, mpi_sumdd, ilocalcomm, ierr ) ! get %y1d 520 todelay(idvar)%z1d(:) = REAL(todelay(idvar)%y1d(:), wp) ! define %z1d from %y1d521 ENDIF 522 523 IF( ndelayid(idvar) > 0 )CALL mpp_delay_rcv( idvar ) ! make sure %z1d is received521 ndelayid(idvar) = MPI_REQUEST_NULL 522 ENDIF 523 524 CALL mpp_delay_rcv( idvar ) ! make sure %z1d is received 524 525 525 526 ! send back pout from todelay(idvar)%z1d defined at previous call … … 530 531 IF( ln_timing ) CALL tic_tac( .TRUE., ld_global = .TRUE.) 531 532 CALL mpi_allreduce( y_in(:), todelay(idvar)%y1d(:), isz, MPI_DOUBLE_COMPLEX, mpi_sumdd, ilocalcomm, ierr ) 532 ndelayid(idvar) = 1533 ndelayid(idvar) = MPI_REQUEST_NULL 533 534 IF( ln_timing ) CALL tic_tac(.FALSE., ld_global = .TRUE.) 534 535 # else … … 591 592 DEALLOCATE(todelay(idvar)%z1d) 592 593 ndelayid(idvar) = -1 ! do as if we had no restart 594 ELSE 595 ndelayid(idvar) = MPI_REQUEST_NULL 593 596 END IF 594 597 ENDIF … … 598 601 ALLOCATE(todelay(idvar)%z1d(isz)) 599 602 CALL mpi_allreduce( p_in(:), todelay(idvar)%z1d(:), isz, MPI_DOUBLE_PRECISION, mpi_max, ilocalcomm, ierr ) ! get %z1d 600 ENDIF 601 602 IF( ndelayid(idvar) > 0 ) CALL mpp_delay_rcv( idvar ) ! make sure %z1d is received 603 ndelayid(idvar) = MPI_REQUEST_NULL 604 ENDIF 605 606 CALL mpp_delay_rcv( idvar ) ! make sure %z1d is received 603 607 604 608 ! send back pout from todelay(idvar)%z1d defined at previous call … … 606 610 607 611 ! send p_in into todelay(idvar)%z1d with a non-blocking communication 612 ! (PM) Should we get rid of MPI2 option ? MPI3 was release in 2013. Who is still using MPI2 ? 608 613 # if defined key_mpi2 609 614 IF( ln_timing ) CALL tic_tac( .TRUE., ld_global = .TRUE.) 610 CALL mpi_allreduce( p_in(:), todelay(idvar)%z1d(:), isz, MPI_TYPE, mpi_max, ilocalcomm, ndelayid(idvar),ierr )615 CALL mpi_allreduce( p_in(:), todelay(idvar)%z1d(:), isz, MPI_TYPE, mpi_max, ilocalcomm, ierr ) 611 616 IF( ln_timing ) CALL tic_tac(.FALSE., ld_global = .TRUE.) 612 617 # else … … 631 636 !!---------------------------------------------------------------------- 632 637 #if defined key_mpp_mpi 633 IF( ndelayid(kid) /= -2 ) THEN 634 #if ! defined key_mpi2 635 IF( ln_timing ) CALL tic_tac( .TRUE., ld_global = .TRUE.) 636 CALL mpi_wait( ndelayid(kid), MPI_STATUS_IGNORE, ierr ) ! make sure todelay(kid) is received 637 IF( ln_timing ) CALL tic_tac(.FALSE., ld_global = .TRUE.) 638 #endif 639 IF( ASSOCIATED(todelay(kid)%y1d) ) todelay(kid)%z1d(:) = REAL(todelay(kid)%y1d(:), wp) ! define %z1d from %y1d 640 ndelayid(kid) = -2 ! add flag to know that mpi_wait was already called on kid 641 ENDIF 638 IF( ln_timing ) CALL tic_tac( .TRUE., ld_global = .TRUE.) 639 ! test on ndelayid(kid) useless as mpi_wait return immediatly if the request handle is MPI_REQUEST_NULL 640 CALL mpi_wait( ndelayid(kid), MPI_STATUS_IGNORE, ierr ) ! after this ndelayid(kid) = MPI_REQUEST_NULL 641 IF( ln_timing ) CALL tic_tac( .FALSE., ld_global = .TRUE.) 642 IF( ASSOCIATED(todelay(kid)%y1d) ) todelay(kid)%z1d(:) = REAL(todelay(kid)%y1d(:), wp) ! define %z1d from %y1d 642 643 #endif 643 644 END SUBROUTINE mpp_delay_rcv -
NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/src/OCE/LDF/ldftra.F90
r13553 r13741 246 246 ENDIF 247 247 ! 248 IF( ln_ldfeiv .AND. .NOT.( ln_traldf_iso .OR. ln_traldf_triad ) ) & 249 & CALL ctl_stop( 'ln_ldfeiv=T requires iso-neutral laplacian diffusion' ) 250 IF( ln_isfcav .AND. ln_traldf_triad ) & 251 & CALL ctl_stop( ' ice shelf cavity and traldf_triad not tested' ) 248 IF( ln_isfcav .AND. ln_traldf_triad ) CALL ctl_stop( ' ice shelf cavity and traldf_triad not tested' ) 252 249 ! 253 250 IF( nldf_tra == np_lap_i .OR. nldf_tra == np_lap_it .OR. & … … 541 538 IF( ln_traldf_blp ) CALL ctl_stop( 'ldf_eiv_init: eddy induced velocity ONLY with laplacian diffusivity' ) 542 539 ! 540 IF( .NOT.( ln_traldf_iso .OR. ln_traldf_triad ) ) & 541 & CALL ctl_stop( 'ln_ldfeiv=T requires iso-neutral laplacian diffusion' ) 543 542 ! != allocate the aei arrays 544 543 ALLOCATE( aeiu(jpi,jpj,jpk), aeiv(jpi,jpj,jpk), STAT=ierr ) -
NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/src/OCE/SBC/sbcfwb.F90
r13286 r13741 94 94 snwice_mass_b(:,:) = 0.e0 ! no sea-ice model is being used : no snow+ice mass 95 95 snwice_mass (:,:) = 0.e0 96 snwice_fmass (:,:) = 0.e0 96 97 #endif 97 98 ! -
NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/src/OCE/ZDF/zdfdrg.F90
r13553 r13741 383 383 IF(ll_bot) zmsk_boost(:,:) = zmsk_boost(:,:) * ssmask(:,:) ! x seafloor mask 384 384 ! 385 l_log_not_linssh = .FALSE. ! default definition 385 386 ! 386 387 SELECT CASE( ndrg ) -
NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/src/OCE/ZDF/zdfgls.F90
r13553 r13741 815 815 WRITE(numout,*) ' Ice-ocean roughness (used if nn_z0_ice/=0) rn_hsri = ', rn_hsri 816 816 WRITE(numout,*) 817 WRITE(numout,*) ' Namelist namdrg_top/_bot: used values:'818 WRITE(numout,*) ' top ocean cavity roughness (m) rn_z0(_top) = ', r_z0_top819 WRITE(numout,*) ' Bottom seafloor roughness (m) rn_z0(_bot) = ', r_z0_bot820 WRITE(numout,*)821 817 ENDIF 822 818 -
NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/src/OCE/ZDF/zdfphy.F90
r13553 r13741 337 337 ! 338 338 END SUBROUTINE zdf_phy 339 340 339 341 INTEGER FUNCTION zdf_phy_alloc() 340 342 !!---------------------------------------------------------------------- -
NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/src/OCE/ZDF/zdftke.F90
r13553 r13741 678 678 CALL ctl_stop( 'zdf_tke_init: wrong value for nn_eice, should be 0,1,2, or 3') 679 679 END SELECT 680 IF( .NOT.ln_drg_OFF ) THEN681 WRITE(numout,*)682 WRITE(numout,*) ' Namelist namdrg_top/_bot: used values:'683 WRITE(numout,*) ' top ocean cavity roughness (m) rn_z0(_top)= ', r_z0_top684 WRITE(numout,*) ' Bottom seafloor roughness (m) rn_z0(_bot)= ', r_z0_bot685 ENDIF686 680 WRITE(numout,*) 687 681 WRITE(numout,*) ' ==>>> critical Richardson nb with your parameters ri_cri = ', ri_cri -
NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/src/OCE/nemogcm.F90
r13286 r13741 54 54 USE asminc ! assimilation increments 55 55 USE asmbkg ! writing out state trajectory 56 USE diaptr ! poleward transports (dia_ptr_init routine)57 56 USE diadct ! sections transports (dia_dct_init routine) 58 57 USE diaobs ! Observation diagnostics (dia_obs_init routine) … … 472 471 ! ! Lateral physics 473 472 CALL ldf_tra_init ! Lateral ocean tracer physics 474 CALL ldf_eiv_init ! eddy induced velocity param. 473 CALL ldf_eiv_init ! eddy induced velocity param. must be done after ldf_tra_init 475 474 CALL ldf_dyn_init ! Lateral ocean momentum physics 476 475 … … 510 509 CALL flo_init( Nnn ) ! drifting Floats 511 510 IF( ln_diacfl ) CALL dia_cfl_init ! Initialise CFL diagnostics 512 ! CALL dia_ptr_init ! Poleward TRansports initialization513 511 CALL dia_dct_init ! Sections tranports 514 512 CALL dia_hsb_init( Nnn ) ! heat content, salt content and volume budgets -
NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/src/OCE/stpctl.F90
r13553 r13741 67 67 REAL(wp) :: zzz ! local real 68 68 REAL(wp), DIMENSION(9) :: zmax, zmaxlocal 69 LOGICAL :: ll_wrtstp, ll_colruns, ll_wrtruns 69 LOGICAL :: ll_wrtstp, ll_colruns, ll_wrtruns, ll_0oce 70 70 LOGICAL, DIMENSION(jpi,jpj,jpk) :: llmsk 71 71 CHARACTER(len=20) :: clname … … 125 125 ! 126 126 llmsk(Nis0:Nie0,Njs0:Nje0,1) = ssmask(Nis0:Nie0,Njs0:Nje0) == 1._wp ! define only the inner domain 127 ! 128 ll_0oce = .NOT. ANY( llmsk(:,:,1) ) ! no ocean point in the inner domain? 129 ! 127 130 IF( ll_wd ) THEN 128 131 zmax(1) = MAXVAL( ABS( ssh(:,:,Kmm) + ssh_ref ), mask = llmsk(:,:,1) ) ! ssh max … … 149 152 ENDIF 150 153 zmax(9) = REAL( nstop, wp ) ! stop indicator 154 ! 151 155 ! !== get global extrema ==! 152 156 ! !== done by all processes if writting run.stat ==! 153 157 IF( ll_colruns ) THEN 154 158 zmaxlocal(:) = zmax(:) 155 CALL mpp_max( "stpctl", zmax ) ! max over the global domain 159 CALL mpp_max( "stpctl", zmax ) ! max over the global domain: ok even of ll_0oce = .true. 156 160 nstop = NINT( zmax(9) ) ! update nstop indicator (now sheared among all local domains) 157 ENDIF 161 ELSE 162 ! if no ocean point: MAXVAL returns -HUGE => we must overwrite this value to avoid error handling bellow. 163 IF( ll_0oce ) zmax(1:4) = (/ 0._wp, 0._wp, -1._wp, 1._wp /) ! default "valid" values... 164 ENDIF 165 ! 166 zmax(3) = -zmax(3) ! move back from max(-zz) to min(zz) : easier to manage! 167 zmax(5) = -zmax(5) ! move back from max(-zz) to min(zz) : easier to manage! 168 IF( ll_colruns ) THEN 169 zmaxlocal(3) = -zmaxlocal(3) ! move back from max(-zz) to min(zz) : easier to manage! 170 zmaxlocal(5) = -zmaxlocal(5) ! move back from max(-zz) to min(zz) : easier to manage! 171 ENDIF 172 ! 158 173 ! !== write "run.stat" files ==! 159 174 ! !== done only by 1st subdomain at writting timestep ==! 160 175 IF( ll_wrtruns ) THEN 161 WRITE(numrun,9500) kt, zmax(1), zmax(2), -zmax(3), zmax(4) 162 istatus = NF90_PUT_VAR( nrunid, nvarid(1), (/ zmax(1)/), (/kt/), (/1/) ) 163 istatus = NF90_PUT_VAR( nrunid, nvarid(2), (/ zmax(2)/), (/kt/), (/1/) ) 164 istatus = NF90_PUT_VAR( nrunid, nvarid(3), (/-zmax(3)/), (/kt/), (/1/) ) 165 istatus = NF90_PUT_VAR( nrunid, nvarid(4), (/ zmax(4)/), (/kt/), (/1/) ) 166 istatus = NF90_PUT_VAR( nrunid, nvarid(5), (/-zmax(5)/), (/kt/), (/1/) ) 167 istatus = NF90_PUT_VAR( nrunid, nvarid(6), (/ zmax(6)/), (/kt/), (/1/) ) 168 IF( ln_zad_Aimp ) THEN 169 istatus = NF90_PUT_VAR( nrunid, nvarid(7), (/ zmax(7)/), (/kt/), (/1/) ) 170 istatus = NF90_PUT_VAR( nrunid, nvarid(8), (/ zmax(8)/), (/kt/), (/1/) ) 171 ENDIF 176 WRITE(numrun,9500) kt, zmax(1), zmax(2), zmax(3), zmax(4) 177 DO ji = 1, 6 + 2 * COUNT( (/ln_zad_Aimp/) ) 178 istatus = NF90_PUT_VAR( nrunid, nvarid(ji), (/zmax(ji)/), (/kt/), (/1/) ) 179 END DO 172 180 IF( kt == nitend ) istatus = NF90_CLOSE(nrunid) 173 181 END IF … … 175 183 ! !== done by all processes at every time step ==! 176 184 ! 177 IF( 178 & 179 & zmax(3) >= 0._wp .OR. & ! negative or zero sea surface salinity180 & 181 & 182 & 183 & 185 IF( zmax(1) > 20._wp .OR. & ! too large sea surface height ( > 20 m ) 186 & zmax(2) > 10._wp .OR. & ! too large velocity ( > 10 m/s) 187 & zmax(3) <= 0._wp .OR. & ! negative or zero sea surface salinity 188 & zmax(4) >= 100._wp .OR. & ! too large sea surface salinity ( > 100 ) 189 & zmax(4) < 0._wp .OR. & ! too large sea surface salinity (keep this line for sea-ice) 190 & ISNAN( zmax(1) + zmax(2) + zmax(3) ) .OR. & ! NaN encounter in the tests 191 & ABS( zmax(1) + zmax(2) + zmax(3) ) > HUGE(1._wp) ) THEN ! Infinity encounter in the tests 184 192 ! 185 193 iloc(:,:) = 0 … … 221 229 ! 222 230 WRITE(ctmp1,*) ' stp_ctl: |ssh| > 20 m or |U| > 10 m/s or S <= 0 or S >= 100 or NaN encounter in the tests' 223 CALL wrt_line( ctmp2, kt, '|ssh| max', 224 CALL wrt_line( ctmp3, kt, '|U| max', 225 CALL wrt_line( ctmp4, kt, 'Sal min', -zmax(3), iloc(:,3), iareasum(3), iareamin(3), iareamax(3) )226 CALL wrt_line( ctmp5, kt, 'Sal max', 231 CALL wrt_line( ctmp2, kt, '|ssh| max', zmax(1), iloc(:,1), iareasum(1), iareamin(1), iareamax(1) ) 232 CALL wrt_line( ctmp3, kt, '|U| max', zmax(2), iloc(:,2), iareasum(2), iareamin(2), iareamax(2) ) 233 CALL wrt_line( ctmp4, kt, 'Sal min', zmax(3), iloc(:,3), iareasum(3), iareamin(3), iareamax(3) ) 234 CALL wrt_line( ctmp5, kt, 'Sal max', zmax(4), iloc(:,4), iareasum(4), iareamin(4), iareamax(4) ) 227 235 IF( Agrif_Root() ) THEN 228 236 WRITE(ctmp6,*) ' ===> output of last computed fields in output.abort* files' -
NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/src/OCE/timing.F90
r13553 r13741 424 424 s_timer => s_timer_root 425 425 DO WHILE ( ASSOCIATED( s_timer%next ) ) 426 IF (.NOT. ASSOCIATED(s_timer%next)) EXIT426 IF (.NOT. ASSOCIATED(s_timer%next)) EXIT 427 427 IF ( s_timer%tsum_clock < s_timer%next%tsum_clock ) THEN 428 428 ALLOCATE(s_wrk) … … 432 432 ll_ord = .FALSE. 433 433 CYCLE 434 ENDIF 435 IF( ASSOCIATED(s_timer%next) ) s_timer => s_timer%next436 END DO 434 ENDIF 435 IF( ASSOCIATED(s_timer%next) ) s_timer => s_timer%next 436 END DO 437 437 IF( ll_ord ) EXIT 438 438 END DO … … 447 447 clfmt = '(1x,a,4x,f12.3,6x,f12.3,x,f12.3,2x,f12.3,6x,f7.3,2x,i9)' 448 448 DO WHILE ( ASSOCIATED(s_timer) ) 449 WRITE(numtime,TRIM(clfmt)) s_timer%cname, & 450 & s_timer%tsum_clock,s_timer%tsum_clock*100./t_elaps(2), & 451 & s_timer%tsum_cpu ,s_timer%tsum_cpu*100./t_cpu(2) , & 452 & s_timer%tsum_cpu/s_timer%tsum_clock, s_timer%niter 449 IF( s_timer%tsum_clock > 0._wp ) & 450 WRITE(numtime,TRIM(clfmt)) s_timer%cname, & 451 & s_timer%tsum_clock,s_timer%tsum_clock*100./t_elaps(2), & 452 & s_timer%tsum_cpu ,s_timer%tsum_cpu*100./t_cpu(2) , & 453 & s_timer%tsum_cpu/s_timer%tsum_clock, s_timer%niter 453 454 s_timer => s_timer%next 454 455 END DO … … 613 614 clfmt = '((A),E15.7,2x,f6.2,5x,f12.2,5x,f6.2,5x,f7.2,2x,f12.2,4x,f6.2,2x,f9.2)' 614 615 DO WHILE ( ASSOCIATED(sl_timer_ave) ) 615 WRITE(numtime,TRIM(clfmt)) sl_timer_ave%cname(1:18), & 616 & sl_timer_ave%tsum_clock,sl_timer_ave%tsum_clock*100.*jpnij/tot_etime, & 617 & sl_timer_ave%tsum_cpu ,sl_timer_ave%tsum_cpu*100.*jpnij/tot_ctime , & 618 & sl_timer_ave%tsum_cpu/sl_timer_ave%tsum_clock, & 619 & sl_timer_ave%tmax_clock*100.*jpnij/tot_etime, & 620 & sl_timer_ave%tmin_clock*100.*jpnij/tot_etime, & 621 & sl_timer_ave%niter/REAL(jpnij) 616 IF( sl_timer_ave%tsum_clock > 0. ) & 617 WRITE(numtime,TRIM(clfmt)) sl_timer_ave%cname(1:18), & 618 & sl_timer_ave%tsum_clock,sl_timer_ave%tsum_clock*100.*jpnij/tot_etime, & 619 & sl_timer_ave%tsum_cpu ,sl_timer_ave%tsum_cpu*100.*jpnij/tot_ctime , & 620 & sl_timer_ave%tsum_cpu/sl_timer_ave%tsum_clock, & 621 & sl_timer_ave%tmax_clock*100.*jpnij/tot_etime, & 622 & sl_timer_ave%tmin_clock*100.*jpnij/tot_etime, & 623 & sl_timer_ave%niter/REAL(jpnij) 622 624 sl_timer_ave => sl_timer_ave%next 623 625 END DO -
NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/src/OFF/nemogcm.F90
r13286 r13741 40 40 USE sbcmod ! surface boundary condition (sbc_init routine) 41 41 USE phycst ! physical constant (par_cst routine) 42 USE zdfphy ! vertical physics manager (zdf_phy_init routine) 42 43 USE dtadyn ! Lecture and Interpolation of the dynamical fields 43 44 USE trcini ! Initilization of the passive tracers … … 49 50 USE trcnam ! passive tracer : namelist 50 51 USE trcrst ! passive tracer restart 51 USE diaptr ! Need to initialise this as some variables are used in if statements later52 52 USE sbc_oce , ONLY : ln_rnf 53 53 USE sbcrnf ! surface boundary condition : runoffs … … 73 73 74 74 CHARACTER (len=64) :: cform_aaa="( /, 'AAAAAAAA', / ) " ! flag for output listing 75 #if defined key_mpp_mpi 76 ! need MPI_Wtime 77 INCLUDE 'mpif.h' 78 #endif 75 79 76 80 !!---------------------------------------------------------------------- … … 96 100 !!---------------------------------------------------------------------- 97 101 INTEGER :: istp ! time step index 102 REAL(wp):: zstptiming ! elapsed time for 1 time step 98 103 !!---------------------------------------------------------------------- 99 104 … … 114 119 ! 115 120 DO WHILE ( istp <= nitend .AND. nstop == 0 ) !== OFF time-stepping ==! 121 122 IF( ln_timing ) THEN 123 zstptiming = MPI_Wtime() 124 IF ( istp == ( nit000 + 1 ) ) elapsed_time = zstptiming 125 IF ( istp == nitend ) elapsed_time = zstptiming - elapsed_time 126 ENDIF 116 127 ! 117 128 IF( istp /= nit000 ) CALL day ( istp ) ! Calendar (day was already called at nit000 in day_init) … … 147 158 #endif 148 159 #endif 149 160 CALL stp_ctl ( istp ) ! Time loop: control and print 150 161 istp = istp + 1 162 163 IF( lwp .AND. ln_timing ) WRITE(numtime,*) 'timing step ', istp-1, ' : ', MPI_Wtime() - zstptiming 164 151 165 END DO 152 166 ! … … 333 347 334 348 CALL sbc_init( Nbb, Nnn, Naa ) ! Forcings : surface module 335 CALL bdy_init ! Open boundaries initialisation 349 CALL bdy_init ! Open boundaries initialisation 350 351 CALL zdf_phy_init( Nnn ) ! Vertical physics 336 352 337 353 ! ! Tracer physics 338 354 CALL ldf_tra_init ! Lateral ocean tracer physics 339 CALL ldf_eiv_init ! Eddy induced velocity param 355 CALL ldf_eiv_init ! Eddy induced velocity param. must be done after ldf_tra_init 340 356 CALL tra_ldf_init ! lateral mixing 341 357 IF( l_ldfslp ) CALL ldf_slp_init ! slope of lateral mixing … … 351 367 CALL dta_dyn_init( Nbb, Nnn, Naa ) ! Initialization for the dynamics 352 368 #endif 353 354 369 CALL trc_init( Nbb, Nnn, Naa ) ! Passive tracers initialization 355 CALL dia_ptr_init ! Poleward TRansports initialization356 370 357 371 IF(lwp) WRITE(numout,cform_aaa) ! Flag AAAAAAA -
NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/src/SAS/nemogcm.F90
r13553 r13741 2 2 !!====================================================================== 3 3 !! *** MODULE nemogcm *** 4 !! StandAlone Surface module : surface fluxes + sea-ice + iceberg floats 4 !! StandAlone Surface module : surface fluxes + sea-ice + iceberg floats + ABL 5 5 !!====================================================================== 6 6 !! History : 3.6 ! 2011-11 (S. Alderson, G. Madec) original code … … 58 58 59 59 #if defined key_mpp_mpi 60 ! need MPI_Wtime 60 61 INCLUDE 'mpif.h' 61 62 #endif … … 83 84 !!---------------------------------------------------------------------- 84 85 INTEGER :: istp ! time step index 86 REAL(wp):: zstptiming ! elapsed time for 1 time step 85 87 !!---------------------------------------------------------------------- 86 88 ! … … 93 95 #if defined key_agrif 94 96 Kbb_a = Nbb; Kmm_a = Nnn; Krhs_a = Nrhs ! agrif_oce module copies of time level indices 95 CALL Agrif_Declare_Var ! " " " " " DYN/TRA 97 CALL Agrif_Declare_Var ! " " " " " DYN/TRA 96 98 # if defined key_top 97 99 CALL Agrif_Declare_Var_top ! " " " " " TOP … … 107 109 ! !== time stepping ==! 108 110 ! !-----------------------! 111 ! 112 ! !== set the model time-step ==! 113 ! 109 114 istp = nit000 110 115 ! … … 124 129 END DO 125 130 ! 126 # else131 # else 127 132 ! 128 133 IF( .NOT.ln_diurnal_only ) THEN !== Standard time-stepping ==! 129 134 ! 130 135 DO WHILE( istp <= nitend .AND. nstop == 0 ) 131 #if defined key_mpp_mpi 136 132 137 ncom_stp = istp 133 IF ( istp == ( nit000 + 1 ) ) elapsed_time = MPI_Wtime() 134 IF ( istp == nitend ) elapsed_time = MPI_Wtime() - elapsed_time 135 #endif 138 IF( ln_timing ) THEN 139 zstptiming = MPI_Wtime() 140 IF ( istp == ( nit000 + 1 ) ) elapsed_time = zstptiming 141 IF ( istp == nitend ) elapsed_time = zstptiming - elapsed_time 142 ENDIF 143 136 144 CALL stp ( istp ) 137 145 istp = istp + 1 146 147 IF( lwp .AND. ln_timing ) WRITE(numtime,*) 'timing step ', istp-1, ' : ', MPI_Wtime() - zstptiming 148 138 149 END DO 139 150 ! … … 305 316 WRITE(numout,*) " ) ) \) |`\ \) '. \ ( ( " 306 317 WRITE(numout,*) " ( ( \_/ '-._\ ) ) " 307 WRITE(numout,*) " ) ) jgs `( ( "318 WRITE(numout,*) " ) ) jgs ` ( ( " 308 319 WRITE(numout,*) " ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ " 309 320 WRITE(numout,*) … … 367 378 & CALL prt_ctl_init ! Print control 368 379 380 IF( ln_rstart ) CALL rst_read_open 369 381 CALL day_init ! model calendar (using both namelist and restart infos) 370 IF( ln_rstart ) CALL rst_read_open 371 372 ! ! external forcing 382 373 383 #if defined key_agrif 374 384 uu(:,:,:,:) = 0.0_wp ; vv(:,:,:,:) = 0.0_wp ; ts(:,:,:,:,:) = 0.0_wp ! needed for interp done at initialization phase 375 385 #endif 386 ! ! external forcing 376 387 CALL sbc_init( Nbb, Nnn, Naa ) ! Forcings : surface module 377 388 -
NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/src/SAS/stpctl.F90
r13553 r13741 68 68 REAL(wp) :: zzz ! local real 69 69 REAL(wp), DIMENSION(4) :: zmax, zmaxlocal 70 LOGICAL :: ll_wrtstp, ll_colruns, ll_wrtruns 70 LOGICAL :: ll_wrtstp, ll_colruns, ll_wrtruns, ll_0oce 71 71 LOGICAL, DIMENSION(jpi,jpj) :: llmsk 72 72 CHARACTER(len=20) :: clname … … 124 124 llmsk(:,Nje1: jpj) = .FALSE. 125 125 ! 126 llmsk(Nis0:Nie0,Njs0:Nje0) = tmask(Nis0:Nie0,Njs0:Nje0,1) == 1._wp ! test only the inner domain 127 IF( COUNT( llmsk(:,:) ) > 0 ) THEN ! avoid huge values sent back for land processors... 128 zmax(1) = MAXVAL( vt_i (:,:) , mask = llmsk ) ! max ice thickness 129 zmax(2) = MAXVAL( ABS( u_ice(:,:) ) , mask = llmsk ) ! max ice velocity (zonal only) 130 zmax(3) = MAXVAL( -tm_i (:,:) + rt0, mask = llmsk ) ! min ice temperature (in degC) 131 ELSE 132 IF( ll_colruns ) THEN ! default value: must not be kept when calling mpp_max -> must be as small as possible 133 zmax(1:3) = -HUGE(1._wp) 134 ELSE ! default value: must not give true for any of the tests bellow (-> avoid manipulating HUGE...) 135 zmax(1:3) = 0._wp 136 ENDIF 137 ENDIF 138 zmax(4) = REAL( nstop, wp ) ! stop indicator 126 llmsk(Nis0:Nie0,Njs0:Nje0) = tmask(Nis0:Nie0,Njs0:Nje0,1) == 1._wp ! test only the inner domain 127 ! 128 ll_0oce = .NOT. ANY( llmsk(:,:) ) ! no ocean point in the inner domain? 129 ! 130 zmax(1) = MAXVAL( vt_i (:,:) , mask = llmsk ) ! max ice thickness 131 zmax(2) = MAXVAL( ABS( u_ice(:,:) ) , mask = llmsk ) ! max ice velocity (zonal only) 132 zmax(3) = MAXVAL( -tm_i (:,:) + rt0, mask = llmsk ) ! min ice temperature (in degC) 133 zmax(4) = REAL( nstop, wp ) ! stop indicator 134 ! 139 135 ! !== get global extrema ==! 140 136 ! !== done by all processes if writting run.stat ==! 141 137 IF( ll_colruns ) THEN 142 138 zmaxlocal(:) = zmax(:) 143 CALL mpp_max( "stpctl", zmax ) ! max over the global domain 139 CALL mpp_max( "stpctl", zmax ) ! max over the global domain: ok even of ll_0oce = .true. 144 140 nstop = NINT( zmax(4) ) ! update nstop indicator (now sheared among all local domains) 145 ENDIF 141 ELSE 142 ! if no ocean point: MAXVAL returns -HUGE => we must overwrite this value to avoid error handling bellow. 143 IF( ll_0oce ) zmax(1:3) = 0._wp ! default "valid" values... 144 ENDIF 145 ! 146 zmax(3) = -zmax(3) ! move back from max(-zz) to min(zz) : easier to manage! 147 IF( ll_colruns ) zmaxlocal(3) = -zmaxlocal(3) ! move back from max(-zz) to min(zz) : easier to manage! 148 ! 146 149 ! !== write "run.stat" files ==! 147 150 ! !== done only by 1st subdomain at writting timestep ==! 148 151 IF( ll_wrtruns ) THEN 149 WRITE(numrun,9500) kt, zmax(1), zmax(2), -zmax(3)150 istatus = NF90_PUT_VAR( nrunid, nvarid(1), (/ zmax(1)/), (/kt/), (/1/) )151 istatus = NF90_PUT_VAR( nrunid, nvarid(2), (/ zmax(2)/), (/kt/), (/1/) )152 istatus = NF90_PUT_VAR( nrunid, nvarid(3), (/-zmax(3)/), (/kt/), (/1/) )152 WRITE(numrun,9500) kt, zmax(1), zmax(2), zmax(3) 153 DO ji = 1, 3 154 istatus = NF90_PUT_VAR( nrunid, nvarid(ji), (/zmax(ji)/), (/kt/), (/1/) ) 155 END DO 153 156 IF( kt == nitend ) istatus = NF90_CLOSE(nrunid) 154 157 END IF … … 158 161 IF( zmax(1) > 100._wp .OR. & ! too large ice thickness maximum ( > 100 m) 159 162 & zmax(2) > 10._wp .OR. & ! too large ice velocity ( > 10 m/s) 160 & zmax(3) >101._wp .OR. & ! too cold ice temperature ( < -100 degC)163 & zmax(3) < -101._wp .OR. & ! too cold ice temperature ( < -100 degC) 161 164 & ISNAN( zmax(1) + zmax(2) + zmax(3) ) .OR. & ! NaN encounter in the tests 162 165 & ABS( zmax(1) + zmax(2) + zmax(3) ) > HUGE(1._wp) ) THEN ! Infinity encounter in the tests … … 192 195 ! 193 196 WRITE(ctmp1,*) ' stp_ctl: ice_thick > 100 m or |ice_vel| > 10 m/s or ice_temp < -100 degC or NaN encounter in the tests' 194 CALL wrt_line( ctmp2, kt, 'ice_thick max', 195 CALL wrt_line( ctmp3, kt, '|ice_vel| max', 196 CALL wrt_line( ctmp4, kt, 'ice_temp min', -zmax(3), iloc(:,3), iareasum(3), iareamin(3), iareamax(3) )197 CALL wrt_line( ctmp2, kt, 'ice_thick max', zmax(1), iloc(:,1), iareasum(1), iareamin(1), iareamax(1) ) 198 CALL wrt_line( ctmp3, kt, '|ice_vel| max', zmax(2),