Changeset 7607
- Timestamp:
- 2017-01-25T16:37:31+01:00 (8 years ago)
- Location:
- branches/2015/nemo_v3_6_STABLE/NEMOGCM
- Files:
-
- 21 edited
-
CONFIG/ORCA2_LIM3/EXP00/namelist_ice_cfg (modified) (4 diffs)
-
CONFIG/SHARED/field_def.xml (modified) (1 diff)
-
CONFIG/SHARED/namelist_ice_lim3_ref (modified) (2 diffs)
-
CONFIG/SHARED/namelist_ref (modified) (1 diff)
-
NEMO/LIM_SRC_3/ice.F90 (modified) (4 diffs)
-
NEMO/LIM_SRC_3/limdyn.F90 (modified) (3 diffs)
-
NEMO/LIM_SRC_3/limhdf.F90 (modified) (2 diffs)
-
NEMO/LIM_SRC_3/limthd.F90 (modified) (3 diffs)
-
NEMO/LIM_SRC_3/limthd_dif.F90 (modified) (4 diffs)
-
NEMO/OPA_SRC/DOM/phycst.F90 (modified) (3 diffs)
-
NEMO/OPA_SRC/SBC/albedo.F90 (modified) (6 diffs)
-
NEMO/OPA_SRC/SBC/sbc_oce.F90 (modified) (2 diffs)
-
NEMO/OPA_SRC/SBC/sbccpl.F90 (modified) (12 diffs)
-
NEMO/OPA_SRC/SBC/sbcisf.F90 (modified) (9 diffs)
-
NEMO/OPA_SRC/stpctl.F90 (modified) (1 diff)
-
NEMO/TOP_SRC/PISCES/P4Z/p4zche.F90 (modified) (18 diffs)
-
NEMO/TOP_SRC/PISCES/P4Z/p4zflx.F90 (modified) (3 diffs)
-
NEMO/TOP_SRC/PISCES/P4Z/p4zlys.F90 (modified) (5 diffs)
-
NEMO/TOP_SRC/PISCES/P4Z/p4zsms.F90 (modified) (4 diffs)
-
NEMO/TOP_SRC/PISCES/sms_pisces.F90 (modified) (2 diffs)
-
NEMO/TOP_SRC/trcstp.F90 (modified) (4 diffs)
Legend:
- Unmodified
- Added
- Removed
-
branches/2015/nemo_v3_6_STABLE/NEMOGCM/CONFIG/ORCA2_LIM3/EXP00/namelist_ice_cfg
r4690 r7607 1 1 !!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 2 !! NEMO/LIM- 2 : Ice configuration namelist. Overwrites SHARED/namelist_ice_lim2_ref2 !! NEMO/LIM-3 : Ice configuration namelist. Overwrites SHARED/namelist_ice_lim3_ref 3 3 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 4 4 … … 15 15 !----------------------------------------------------------------------- 16 16 / 17 !------------------------------------------------------------------------------ 18 &namicehdf ! Ice horizontal diffusion 19 !------------------------------------------------------------------------------ 20 nn_ahi0 = 2 ! horizontal diffusivity computation 21 / 17 22 !----------------------------------------------------------------------- 18 23 &namicethd ! ice thermodynamic … … 24 29 / 25 30 !----------------------------------------------------------------------- 26 &namiceitdme ! parameters for mechanical redistribution of ice 31 &namiceitdme ! parameters for mechanical redistribution of ice 27 32 !----------------------------------------------------------------------- 28 33 / … … 31 36 !----------------------------------------------------------------------- 32 37 / 38 !------------------------------------------------------------------------------ 39 &namiceitd ! Ice discretization 40 !------------------------------------------------------------------------------ 41 / -
branches/2015/nemo_v3_6_STABLE/NEMOGCM/CONFIG/SHARED/field_def.xml
r7511 r7607 225 225 <field id="icealb_cea" long_name="Ice albedo (cell average)" standard_name="sea_ice_albedo" unit="1" /> 226 226 <field id="calving_cea" long_name="Calving" standard_name="water_flux_into_sea_water_from_icebergs" unit="kg/m2/s" /> 227 <field id="iceberg_cea" long_name="Iceberg" standard_name="water_flux_into_sea_water_from_icebergs" unit="kg/m2/s" /> 228 <field id="iceshelf_cea" long_name="Iceshelf" standard_name="water_flux_into_sea_water_from_iceshelf" unit="kg/m2/s" /> 227 229 228 230 <!-- available if key_oasis3 + conservative method --> -
branches/2015/nemo_v3_6_STABLE/NEMOGCM/CONFIG/SHARED/namelist_ice_lim3_ref
r6316 r7607 69 69 rn_relast = 0.333 ! ratio of elastic timescale to ice time step: Telast = dt_ice * rn_relast 70 70 ! advised value: 1/3 (rn_nevp=120) or 1/9 (rn_nevp=300) 71 nn_ahi0 = 2 ! horizontal diffusivity computation72 ! 0: use rn_ahi0_ref73 ! 1: use rn_ahi0_ref x mean grid cell length / ( 2deg mean grid cell length )74 ! 2: use rn_ahi0_ref x grid cell length / ( 2deg mean grid cell length )75 rn_ahi0_ref = 350.0 ! horizontal sea ice diffusivity (m2/s)76 ! if nn_ahi0 > 0, rn_ahi0_ref is the reference value at a nominal 2 deg resolution77 71 / 78 72 !------------------------------------------------------------------------------ 79 73 &namicehdf ! Ice horizontal diffusion 80 74 !------------------------------------------------------------------------------ 75 nn_ahi0 = -1 ! horizontal diffusivity computation 76 ! -1: no diffusion (bypass limhdf) 77 ! 0: use rn_ahi0_ref 78 ! 1: use rn_ahi0_ref x mean grid cell length / ( 2deg mean grid cell length ) 79 ! 2: use rn_ahi0_ref x grid cell length / ( 2deg mean grid cell length ) 80 rn_ahi0_ref = 350.0 ! horizontal sea ice diffusivity (m2/s) 81 ! if nn_ahi0 > 0, rn_ahi0_ref is the reference value at a nominal 2 deg resolution 81 82 nn_convfrq = 5 ! convergence check frequency of the Crant-Nicholson scheme (perf. optimization) 82 83 / … … 98 99 ! 0: k = k0 + beta.S/T (Untersteiner, 1964) 99 100 ! 1: k = k0 + beta1.S/T - beta2.T (Pringle et al., 2007) 101 rn_cdsn = 0.31 ! thermal conductivity of the snow (0.31 W/m/K, Maykut and Untersteiner, 1971) 102 ! Obs: 0.1-0.5 (Lecomte et al, JAMES 2013) 100 103 nn_monocat = 0 ! virtual ITD mono-category parameterizations (1, jpl = 1 only) or not (0) 101 104 ! 2: simple piling instead of ridging --- temporary option -
branches/2015/nemo_v3_6_STABLE/NEMOGCM/CONFIG/SHARED/namelist_ref
r6319 r7607 500 500 &namsbc_alb ! albedo parameters 501 501 !----------------------------------------------------------------------- 502 nn_ice_alb = 0 ! parameterization of ice/snow albedo 503 ! 0: Shine & Henderson-Sellers (JGR 1985) 504 ! 1: "home made" based on Brandt et al. (J. Climate 2005) 505 ! and Grenfell & Perovich (JGR 2004) 506 rn_albice = 0.53 ! albedo of bare puddled ice (values from 0.49 to 0.58) 507 ! 0.53 (default) => if nn_ice_alb=0 508 ! 0.50 (default) => if nn_ice_alb=1 502 nn_ice_alb = 1 ! parameterization of ice/snow albedo 503 ! 0: Shine & Henderson-Sellers (JGR 1985), giving clear-sky albedo 504 ! 1: "home made" based on Brandt et al. (JClim 2005) and Grenfell & Perovich (JGR 2004), 505 ! giving cloud-sky albedo 506 rn_alb_sdry = 0.85 ! dry snow albedo : 0.80 (nn_ice_alb = 0); 0.85 (nn_ice_alb = 1); obs 0.85-0.87 (cloud-sky) 507 rn_alb_smlt = 0.75 ! melting snow albedo : 0.65 ( '' ) ; 0.75 ( '' ) ; obs 0.72-0.82 ( '' ) 508 rn_alb_idry = 0.60 ! dry ice albedo : 0.72 ( '' ) ; 0.60 ( '' ) ; obs 0.54-0.65 ( '' ) 509 rn_alb_imlt = 0.50 ! bare puddled ice albedo : 0.53 ( '' ) ; 0.50 ( '' ) ; obs 0.49-0.58 ( '' ) 509 510 / 510 511 !----------------------------------------------------------------------- -
branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/LIM_SRC_3/ice.F90
r6963 r7607 212 212 REAL(wp), PUBLIC :: rn_betas !: coef. for partitioning of snowfall between leads and sea ice 213 213 REAL(wp), PUBLIC :: rn_kappa_i !: coef. for the extinction of radiation Grenfell et al. (2006) [1/m] 214 REAL(wp), PUBLIC :: rn_cdsn !: thermal conductivity of the snow [W/m/K] 214 215 REAL(wp), PUBLIC :: nn_conv_dif !: maximal number of iterations for heat diffusion 215 216 REAL(wp), PUBLIC :: rn_terr_dif !: maximal tolerated error (C) for heat diffusion … … 320 321 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: o_i !: Sea-Ice Age (days) 321 322 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: oa_i !: Sea-Ice Age times ice area (days) 323 322 324 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: bv_i !: brine volume 323 325 … … 406 408 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: diag_vice !: ice volume variation [m/s] 407 409 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: diag_vsnw !: snw volume variation [m/s] 410 408 411 ! 409 412 !!---------------------------------------------------------------------- … … 463 466 & et_i (jpi,jpj) , et_s (jpi,jpj) , tm_i (jpi,jpj) , bvm_i(jpi,jpj) , & 464 467 & smt_i(jpi,jpj) , tm_su(jpi,jpj) , htm_i(jpi,jpj) , htm_s(jpi,jpj) , & 465 & om_i (jpi,jpj) 468 & om_i (jpi,jpj) , STAT=ierr(ii) ) 466 469 ii = ii + 1 467 470 ALLOCATE( t_s(jpi,jpj,nlay_s,jpl) , e_s(jpi,jpj,nlay_s,jpl) , STAT=ierr(ii) ) -
branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/LIM_SRC_3/limdyn.F90
r5123 r7607 244 244 INTEGER :: ios ! Local integer output status for namelist read 245 245 NAMELIST/namicedyn/ nn_icestr, ln_icestr_bvf, rn_pe_rdg, rn_pstar, rn_crhg, rn_cio, rn_creepl, rn_ecc, & 246 & nn_nevp, rn_relast, nn_ahi0, rn_ahi0_ref 247 INTEGER :: ji, jj 248 REAL(wp) :: za00, zd_max 246 & nn_nevp, rn_relast 249 247 !!------------------------------------------------------------------- 250 248 … … 272 270 WRITE(numout,*) ' number of iterations for subcycling nn_nevp = ', nn_nevp 273 271 WRITE(numout,*) ' ratio of elastic timescale over ice time step rn_relast = ', rn_relast 274 WRITE(numout,*) ' horizontal diffusivity calculation nn_ahi0 = ', nn_ahi0275 WRITE(numout,*) ' horizontal diffusivity coeff. (orca2 grid) rn_ahi0_ref = ', rn_ahi0_ref276 272 ENDIF 277 273 ! … … 279 275 rhoco = rau0 * rn_cio 280 276 ! 281 ! Diffusion coefficients282 SELECT CASE( nn_ahi0 )283 284 CASE( 0 )285 ahiu(:,:) = rn_ahi0_ref286 ahiv(:,:) = rn_ahi0_ref287 288 IF(lwp) WRITE(numout,*) ''289 IF(lwp) WRITE(numout,*) ' laplacian operator: ahim constant = rn_ahi0_ref'290 291 CASE( 1 )292 293 zd_max = MAX( MAXVAL( e1t(:,:) ), MAXVAL( e2t(:,:) ) )294 IF( lk_mpp ) CALL mpp_max( zd_max ) ! max over the global domain295 296 ahiu(:,:) = rn_ahi0_ref * zd_max * 1.e-05_wp ! 1.e05 = 100km = max grid space at 60° latitude in orca2297 ! (60° = min latitude for ice cover)298 ahiv(:,:) = rn_ahi0_ref * zd_max * 1.e-05_wp299 300 IF(lwp) WRITE(numout,*) ''301 IF(lwp) WRITE(numout,*) ' laplacian operator: ahim proportional to max of e1 e2 over the domain (', zd_max, ')'302 IF(lwp) WRITE(numout,*) ' value for ahim = ', rn_ahi0_ref * zd_max * 1.e-05_wp303 304 CASE( 2 )305 306 zd_max = MAX( MAXVAL( e1t(:,:) ), MAXVAL( e2t(:,:) ) )307 IF( lk_mpp ) CALL mpp_max( zd_max ) ! max over the global domain308 309 za00 = rn_ahi0_ref * 1.e-05_wp ! 1.e05 = 100km = max grid space at 60° latitude in orca2310 ! (60° = min latitude for ice cover)311 DO jj = 1, jpj312 DO ji = 1, jpi313 ahiu(ji,jj) = za00 * MAX( e1t(ji,jj), e2t(ji,jj) ) * umask(ji,jj,1)314 ahiv(ji,jj) = za00 * MAX( e1f(ji,jj), e2f(ji,jj) ) * vmask(ji,jj,1)315 END DO316 END DO317 !318 IF(lwp) WRITE(numout,*) ''319 IF(lwp) WRITE(numout,*) ' laplacian operator: ahim proportional to e1'320 IF(lwp) WRITE(numout,*) ' maximum grid-spacing = ', zd_max, ' maximum value for ahim = ', za00*zd_max321 322 END SELECT323 324 277 END SUBROUTINE lim_dyn_init 325 278 -
branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/LIM_SRC_3/limhdf.F90
r6476 r7607 233 233 !!------------------------------------------------------------------- 234 234 INTEGER :: ios ! Local integer output status for namelist read 235 NAMELIST/namicehdf/ nn_convfrq 236 !!------------------------------------------------------------------- 237 ! 238 IF(lwp) THEN 239 WRITE(numout,*) 240 WRITE(numout,*) 'lim_hdf : Ice horizontal diffusion' 241 WRITE(numout,*) '~~~~~~~' 242 ENDIF 235 NAMELIST/namicehdf/ nn_ahi0, rn_ahi0_ref, nn_convfrq 236 INTEGER :: ji, jj 237 REAL(wp) :: za00, zd_max 238 !!------------------------------------------------------------------- 243 239 ! 244 240 REWIND( numnam_ice_ref ) ! Namelist namicehdf in reference namelist : Ice horizontal diffusion … … 253 249 IF(lwp) THEN ! control print 254 250 WRITE(numout,*) 255 WRITE(numout,*)' Namelist of ice parameters for ice horizontal diffusion computation ' 256 WRITE(numout,*)' convergence check frequency of the Crant-Nicholson scheme nn_convfrq = ', nn_convfrq 251 WRITE(numout,*) 'lim_hdf_init : Ice horizontal diffusion' 252 WRITE(numout,*) '~~~~~~~~~~~' 253 WRITE(numout,*) ' horizontal diffusivity calculation nn_ahi0 = ', nn_ahi0 254 WRITE(numout,*) ' horizontal diffusivity coeff. (orca2 grid) rn_ahi0_ref = ', rn_ahi0_ref 255 WRITE(numout,*) ' convergence check frequency of the Crant-Nicholson scheme nn_convfrq = ', nn_convfrq 257 256 ENDIF 257 ! 258 ! Diffusion coefficients 259 SELECT CASE( nn_ahi0 ) 260 261 CASE( -1 ) 262 ahiu(:,:) = 0._wp 263 ahiv(:,:) = 0._wp 264 265 IF(lwp) WRITE(numout,*) '' 266 IF(lwp) WRITE(numout,*) ' No sea-ice diffusion applied' 267 268 CASE( 0 ) 269 ahiu(:,:) = rn_ahi0_ref 270 ahiv(:,:) = rn_ahi0_ref 271 272 IF(lwp) WRITE(numout,*) '' 273 IF(lwp) WRITE(numout,*) ' laplacian operator: ahim constant = rn_ahi0_ref' 274 275 CASE( 1 ) 276 277 zd_max = MAX( MAXVAL( e1t(:,:) ), MAXVAL( e2t(:,:) ) ) 278 IF( lk_mpp ) CALL mpp_max( zd_max ) ! max over the global domain 279 280 ahiu(:,:) = rn_ahi0_ref * zd_max * 1.e-05_wp ! 1.e05 = 100km = max grid space at 60deg latitude in orca2 281 ! (60deg = min latitude for ice cover) 282 ahiv(:,:) = rn_ahi0_ref * zd_max * 1.e-05_wp 283 284 IF(lwp) WRITE(numout,*) '' 285 IF(lwp) WRITE(numout,*) ' laplacian operator: ahim proportional to max of e1 e2 over the domain (', zd_max, ')' 286 IF(lwp) WRITE(numout,*) ' value for ahim = ', rn_ahi0_ref * zd_max * 1.e-05_wp 287 288 CASE( 2 ) 289 290 zd_max = MAX( MAXVAL( e1t(:,:) ), MAXVAL( e2t(:,:) ) ) 291 IF( lk_mpp ) CALL mpp_max( zd_max ) ! max over the global domain 292 293 za00 = rn_ahi0_ref * 1.e-05_wp ! 1.e05 = 100km = max grid space at 60deg latitude in orca2 294 ! (60deg = min latitude for ice cover) 295 DO jj = 1, jpj 296 DO ji = 1, jpi 297 ahiu(ji,jj) = za00 * MAX( e1t(ji,jj), e2t(ji,jj) ) * umask(ji,jj,1) 298 ahiv(ji,jj) = za00 * MAX( e1f(ji,jj), e2f(ji,jj) ) * vmask(ji,jj,1) 299 END DO 300 END DO 301 ! 302 IF(lwp) WRITE(numout,*) '' 303 IF(lwp) WRITE(numout,*) ' laplacian operator: ahim proportional to e1' 304 IF(lwp) WRITE(numout,*) ' maximum grid-spacing = ', zd_max, ' maximum value for ahim = ', za00*zd_max 305 306 END SELECT 258 307 ! 259 308 END SUBROUTINE lim_hdf_init -
branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90
r6399 r7607 613 613 CALL tab_1d_2d( nbpb, qns_ice(:,:,jl), npb, qns_ice_1d(1:nbpb) , jpi, jpj) 614 614 CALL tab_1d_2d( nbpb, ftr_ice(:,:,jl), npb, ftr_ice_1d(1:nbpb) , jpi, jpj ) 615 !616 615 END SELECT 617 616 … … 633 632 INTEGER :: ios ! Local integer output status for namelist read 634 633 NAMELIST/namicethd/ rn_hnewice, ln_frazil, rn_maxfrazb, rn_vfrazb, rn_Cfrazb, & 635 & rn_himin, rn_betas, rn_kappa_i, nn_conv_dif, rn_terr_dif, nn_ice_thcon, &636 & nn_monocat, ln_it_qnsice634 & rn_himin, rn_betas, rn_kappa_i, nn_conv_dif, rn_terr_dif, nn_ice_thcon, & 635 & rn_cdsn, nn_monocat, ln_it_qnsice 637 636 !!------------------------------------------------------------------- 638 637 ! … … 673 672 WRITE(numout,*)' maximal err. on T for heat diffusion computation rn_terr_dif = ', rn_terr_dif 674 673 WRITE(numout,*)' switch for comp. of thermal conductivity in the ice nn_ice_thcon = ', nn_ice_thcon 674 WRITE(numout,*)' thermal conductivity of the snow rn_cdsn = ', rn_cdsn 675 675 WRITE(numout,*)' check heat conservation in the ice/snow con_i = ', con_i 676 676 WRITE(numout,*)' virtual ITD mono-category parameterizations (1) or not nn_monocat = ', nn_monocat -
branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90
r5512 r7607 376 376 377 377 ! Effective thickness he (zhe) 378 zfac = 1._wp / ( r cdsn + zkimean )379 zratio_s = r cdsn * zfac378 zfac = 1._wp / ( rn_cdsn + zkimean ) 379 zratio_s = rn_cdsn * zfac 380 380 zratio_i = zkimean * zfac 381 381 zhe = zratio_s * ht_i_1d(ji) + zratio_i * ht_s_1d(ji) … … 400 400 DO ji = kideb, kiut 401 401 zfac = 1. / MAX( epsi10 , zh_s(ji) ) 402 zkappa_s(ji,0) = zghe(ji) * r cdsn * zfac403 zkappa_s(ji,nlay_s) = zghe(ji) * r cdsn * zfac402 zkappa_s(ji,0) = zghe(ji) * rn_cdsn * zfac 403 zkappa_s(ji,nlay_s) = zghe(ji) * rn_cdsn * zfac 404 404 END DO 405 405 406 406 DO jk = 1, nlay_s-1 407 407 DO ji = kideb , kiut 408 zkappa_s(ji,jk) = zghe(ji) * 2.0 * r cdsn / MAX( epsi10, 2.0 * zh_s(ji) )408 zkappa_s(ji,jk) = zghe(ji) * 2.0 * rn_cdsn / MAX( epsi10, 2.0 * zh_s(ji) ) 409 409 END DO 410 410 END DO … … 422 422 zkappa_i(ji,0) = zghe(ji) * ztcond_i(ji,0) * zfac 423 423 zkappa_i(ji,nlay_i) = zghe(ji) * ztcond_i(ji,nlay_i) * zfac 424 zkappa_s(ji,nlay_s) = zghe(ji) * zghe(ji) * 2.0 * r cdsn * ztcond_i(ji,0) / &425 & MAX( epsi10, ( zghe(ji) * ztcond_i(ji,0) * zh_s(ji) + zghe(ji) * r cdsn * zh_i(ji) ) )424 zkappa_s(ji,nlay_s) = zghe(ji) * zghe(ji) * 2.0 * rn_cdsn * ztcond_i(ji,0) / & 425 & MAX( epsi10, ( zghe(ji) * ztcond_i(ji,0) * zh_s(ji) + zghe(ji) * rn_cdsn * zh_i(ji) ) ) 426 426 zkappa_i(ji,0) = zkappa_s(ji,nlay_s) * isnow(ji) + zkappa_i(ji,0) * ( 1._wp - isnow(ji) ) 427 427 END DO … … 697 697 & ( isnow(ji) * t_s_1d(ji,1) + ( 1._wp - isnow(ji) ) * t_i_1d(ji,1) ) ) / zdiagbis(ji,numeqmin(ji)) 698 698 END DO 699 699 700 ! 700 701 !-------------------------------------------------------------------------- -
branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/OPA_SRC/DOM/phycst.F90
r5147 r7607 65 65 #if defined key_lim3 || defined key_cice 66 66 REAL(wp), PUBLIC :: rhoic = 917._wp !: volumic mass of sea ice [kg/m3] 67 REAL(wp), PUBLIC :: rcdic = 2.034396_wp !: thermal conductivity of fresh ice 68 REAL(wp), PUBLIC :: rcdsn = 0.31_wp !: thermal conductivity of snow 69 REAL(wp), PUBLIC :: cpic = 2067.0_wp !: specific heat for ice 67 REAL(wp), PUBLIC :: rcdic = 2.034396_wp !: thermal conductivity of fresh ice [W/m/K] 68 REAL(wp), PUBLIC :: cpic = 2067.0_wp !: specific heat of fresh ice [J/kg/K] 70 69 REAL(wp), PUBLIC :: lsub = 2.834e+6_wp !: pure ice latent heat of sublimation [J/kg] 71 70 REAL(wp), PUBLIC :: lfus = 0.334e+6_wp !: latent heat of fusion of fresh ice [J/kg] … … 83 82 REAL(wp), PUBLIC :: xlic = 300.33e+6_wp !: volumetric latent heat fusion of ice [J/m3] 84 83 REAL(wp), PUBLIC :: xsn = 2.8e+6_wp !: volumetric latent heat of sublimation of snow [J/m3] 84 #endif 85 #if defined key_cice 86 REAL(wp), PUBLIC :: rcdsn = 0.31_wp !: thermal conductivity of snow [W/m/K], now namelist parameter for LIM3 85 87 #endif 86 88 #if defined key_lim3 … … 177 179 IF(lwp) THEN 178 180 WRITE(numout,*) 181 #if defined key_cice 179 182 WRITE(numout,*) ' thermal conductivity of the snow = ', rcdsn , ' J/s/m/K' 180 WRITE(numout,*) ' thermal conductivity of the ice = ', rcdic , ' J/s/m/K' 183 #endif 184 WRITE(numout,*) ' thermal conductivity of pure ice = ', rcdic , ' J/s/m/K' 181 185 WRITE(numout,*) ' fresh ice specific heat = ', cpic , ' J/kg/K' 182 186 WRITE(numout,*) ' latent heat of fusion of fresh ice / snow = ', lfus , ' J/kg' -
branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/OPA_SRC/SBC/albedo.F90
r6323 r7607 39 39 ! !!* namelist namsbc_alb 40 40 INTEGER :: nn_ice_alb 41 REAL(wp) :: rn_alb ice41 REAL(wp) :: rn_alb_sdry, rn_alb_smlt, rn_alb_idry, rn_alb_imlt 42 42 43 43 !!---------------------------------------------------------------------- … … 101 101 IF( albd_init == 0 ) CALL albedo_init ! initialization 102 102 103 ralb_sf = rn_alb_sdry ! dry snow 104 ralb_sm = rn_alb_smlt ! melting snow 105 ralb_if = rn_alb_idry ! bare frozen ice 106 ralb_im = rn_alb_imlt ! bare puddled ice 103 107 104 108 SELECT CASE ( nn_ice_alb ) … … 109 113 CASE( 0 ) 110 114 111 ralb_sf = 0.80 ! dry snow112 ralb_sm = 0.65 ! melting snow113 ralb_if = 0.72 ! bare frozen ice114 ralb_im = rn_albice! bare puddled ice115 115 !ralb_sf = 0.80 ! dry snow 116 !ralb_sm = 0.65 ! melting snow 117 !ralb_if = 0.72 ! bare frozen ice 118 !ralb_im = ... ! bare puddled ice 119 116 120 ! Computation of ice albedo (free of snow) 117 121 WHERE ( ph_snw == 0._wp .AND. pt_ice >= rt0_ice ) ; zalb(:,:,:) = ralb_im … … 163 167 CASE( 1 ) 164 168 165 ralb_im = rn_albice! bare puddled ice169 ! ralb_im = ... ! bare puddled ice 166 170 ! compilation of values from literature 167 168 169 171 ! ralb_sf = 0.85 ! dry snow 172 ! ralb_sm = 0.75 ! melting snow 173 ! ralb_if = 0.60 ! bare frozen ice 170 174 ! Perovich et al 2002 (Sheba) => the only dataset for which all types of ice/snow were retrieved 171 175 ! ralb_sf = 0.85 ! dry snow … … 248 252 !!---------------------------------------------------------------------- 249 253 INTEGER :: ios ! Local integer output status for namelist read 250 NAMELIST/namsbc_alb/ nn_ice_alb, rn_alb ice254 NAMELIST/namsbc_alb/ nn_ice_alb, rn_alb_sdry, rn_alb_smlt, rn_alb_idry , rn_alb_imlt 251 255 !!---------------------------------------------------------------------- 252 256 ! … … 268 272 WRITE(numout,*) ' Namelist namsbc_alb : albedo ' 269 273 WRITE(numout,*) ' choose the albedo parameterization nn_ice_alb = ', nn_ice_alb 270 WRITE(numout,*) ' albedo of bare puddled ice rn_albice = ', rn_albice 274 WRITE(numout,*) ' albedo of dry snow rn_alb_sdry = ', rn_alb_sdry 275 WRITE(numout,*) ' albedo of melting snow rn_alb_smlt = ', rn_alb_smlt 276 WRITE(numout,*) ' albedo of dry ice rn_alb_idry = ', rn_alb_idry 277 WRITE(numout,*) ' albedo of bare puddled ice rn_alb_imlt = ', rn_alb_imlt 271 278 ENDIF 272 279 ! -
branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90
r5407 r7607 113 113 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: rnf , rnf_b !: river runoff [Kg/m2/s] 114 114 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: fwfisf , fwfisf_b !: ice shelf melting [Kg/m2/s] 115 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: fwficb , fwficb_b !: iceberg melting [Kg/m2/s] 115 116 !! 116 117 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: sbc_tsc, sbc_tsc_b !: sbc content trend [K.m/s] jpi,jpj,jpts … … 164 165 ! 165 166 ALLOCATE( fwfisf (jpi,jpj), rnf (jpi,jpj) , sbc_tsc (jpi,jpj,jpts) , qsr_hc (jpi,jpj,jpk) , & 166 & fwfisf_b(jpi,jpj), rnf_b(jpi,jpj) , sbc_tsc_b(jpi,jpj,jpts) , qsr_hc_b(jpi,jpj,jpk) , STAT=ierr(3) ) 167 & fwfisf_b(jpi,jpj), rnf_b(jpi,jpj) , sbc_tsc_b(jpi,jpj,jpts) , qsr_hc_b(jpi,jpj,jpk) , & 168 & fwficb (jpi,jpj), fwficb_b(jpi,jpj), STAT=ierr(3) ) 167 169 ! 168 170 ALLOCATE( tprecip(jpi,jpj) , sprecip(jpi,jpj) , fr_i(jpi,jpj) , & -
branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90
r7494 r7607 43 43 USE eosbn2 44 44 USE sbcrnf , ONLY : l_rnfcpl 45 USE sbcisf , ONLY : l_isfcpl 45 46 #if defined key_cpl_carbon_cycle 46 47 USE p4zflx, ONLY : oce_co2 … … 105 106 INTEGER, PARAMETER :: jpr_e3t1st = 41 ! first T level thickness 106 107 INTEGER, PARAMETER :: jpr_fraqsr = 42 ! fraction of solar net radiation absorbed in the first ocean level 107 INTEGER, PARAMETER :: jprcv = 42 ! total number of fields received 108 INTEGER, PARAMETER :: jpr_isf = 43 109 INTEGER, PARAMETER :: jpr_icb = 44 110 INTEGER, PARAMETER :: jprcv = 44 ! total number of fields received 108 111 109 112 INTEGER, PARAMETER :: jps_fice = 1 ! ice fraction sent to the atmosphere … … 149 152 ! Received from the atmosphere ! 150 153 TYPE(FLD_C) :: sn_rcv_w10m, sn_rcv_taumod, sn_rcv_tau, sn_rcv_dqnsdt, sn_rcv_qsr, sn_rcv_qns, sn_rcv_emp, sn_rcv_rnf 151 TYPE(FLD_C) :: sn_rcv_cal, sn_rcv_iceflx, sn_rcv_co2 154 TYPE(FLD_C) :: sn_rcv_cal, sn_rcv_iceflx, sn_rcv_co2, sn_rcv_icb, sn_rcv_isf 152 155 ! Other namelist parameters ! 153 156 INTEGER :: nn_cplmodel ! Maximum number of models to/from which NEMO is potentialy sending/receiving data … … 219 222 & sn_rcv_w10m, sn_rcv_taumod, sn_rcv_tau , sn_rcv_dqnsdt, sn_rcv_qsr, & 220 223 & sn_rcv_qns , sn_rcv_emp , sn_rcv_rnf , sn_rcv_cal , sn_rcv_iceflx, & 221 & sn_rcv_co2 , nn_cplmodel , ln_usecplmask224 & sn_rcv_co2 , sn_rcv_icb , sn_rcv_isf, nn_cplmodel , ln_usecplmask 222 225 !!--------------------------------------------------------------------- 223 226 ! … … 258 261 WRITE(numout,*)' runoffs = ', TRIM(sn_rcv_rnf%cldes ), ' (', TRIM(sn_rcv_rnf%clcat ), ')' 259 262 WRITE(numout,*)' calving = ', TRIM(sn_rcv_cal%cldes ), ' (', TRIM(sn_rcv_cal%clcat ), ')' 263 WRITE(numout,*)' iceberg = ', TRIM(sn_rcv_icb%cldes ), ' (', TRIM(sn_rcv_icb%clcat ), ')' 264 WRITE(numout,*)' ice shelf = ', TRIM(sn_rcv_isf%cldes ), ' (', TRIM(sn_rcv_isf%clcat ), ')' 260 265 WRITE(numout,*)' sea ice heat fluxes = ', TRIM(sn_rcv_iceflx%cldes), ' (', TRIM(sn_rcv_iceflx%clcat), ')' 261 266 WRITE(numout,*)' atm co2 = ', TRIM(sn_rcv_co2%cldes ), ' (', TRIM(sn_rcv_co2%clcat ), ')' … … 397 402 END SELECT 398 403 399 ! ! ------------------------- ! 400 ! ! Runoffs & Calving ! 401 ! ! ------------------------- ! 404 405 ! ! ---------------------------------------------------- ! 406 ! ! Runoffs, Calving, Iceberg, Iceshelf cavities ! 407 ! ! ---------------------------------------------------- ! 402 408 srcv(jpr_rnf )%clname = 'O_Runoff' 403 409 IF( TRIM( sn_rcv_rnf%cldes ) == 'coupled' ) THEN … … 409 415 ENDIF 410 416 ! 411 srcv(jpr_cal )%clname = 'OCalving' ; IF( TRIM( sn_rcv_cal%cldes ) == 'coupled' ) srcv(jpr_cal)%laction = .TRUE. 417 srcv(jpr_cal)%clname = 'OCalving' ; IF( TRIM( sn_rcv_cal%cldes) == 'coupled' ) srcv(jpr_cal)%laction = .TRUE. 418 srcv(jpr_isf)%clname = 'OIcshelf' ; IF( TRIM( sn_rcv_isf%cldes) == 'coupled' ) srcv(jpr_isf)%laction = .TRUE. 419 srcv(jpr_icb)%clname = 'OIceberg' ; IF( TRIM( sn_rcv_icb%cldes) == 'coupled' ) srcv(jpr_icb)%laction = .TRUE. 420 421 IF( srcv(jpr_isf)%laction .AND. nn_isf > 0 ) THEN 422 l_isfcpl = .TRUE. ! -> no need to read isf in sbcisf 423 IF(lwp) WRITE(numout,*) 424 IF(lwp) WRITE(numout,*) ' iceshelf received from oasis ' 425 ENDIF 412 426 413 427 ! ! ------------------------- ! … … 1071 1085 ENDIF 1072 1086 ! 1087 ! 1073 1088 ! ! runoffs and calving (added in emp) 1074 IF( srcv(jpr_rnf)%laction ) rnf(:,:) = frcv(jpr_rnf)%z3(:,:,1)1089 IF( srcv(jpr_rnf)%laction ) rnf(:,:) = frcv(jpr_rnf)%z3(:,:,1) 1075 1090 IF( srcv(jpr_cal)%laction ) zemp(:,:) = zemp(:,:) - frcv(jpr_cal)%z3(:,:,1) 1091 1092 IF( srcv(jpr_icb)%laction ) THEN 1093 fwficb(:,:) = frcv(jpr_icb)%z3(:,:,1) 1094 rnf(:,:) = rnf(:,:) + fwficb(:,:) ! iceberg added to runfofs 1095 ENDIF 1096 IF( srcv(jpr_isf)%laction ) fwfisf(:,:) = - frcv(jpr_isf)%z3(:,:,1) ! fresh water flux from the isf (fwfisf <0 mean melting) 1076 1097 1077 1098 IF( ln_mixcpl ) THEN ; emp(:,:) = emp(:,:) * xcplmask(:,:,0) + zemp(:,:) * zmsk(:,:) … … 1091 1112 ENDIF 1092 1113 ENDIF 1114 ! 1115 IF( srcv(jpr_icb)%laction ) zqns(:,:) = zqns(:,:) - frcv(jpr_icb)%z3(:,:,1) * lfus ! remove heat content associated to iceberg melting 1116 ! 1093 1117 IF( ln_mixcpl ) THEN ; qns(:,:) = qns(:,:) * xcplmask(:,:,0) + zqns(:,:) * zmsk(:,:) 1094 1118 ELSE ; qns(:,:) = zqns(:,:) … … 1468 1492 ENDIF 1469 1493 1494 IF( srcv(jpr_icb)%laction ) THEN 1495 fwficb(:,:) = frcv(jpr_icb)%z3(:,:,1) 1496 rnf(:,:) = rnf(:,:) + fwficb(:,:) ! iceberg added to runoffs 1497 CALL iom_put( 'iceberg_cea', frcv(jpr_icb)%z3(:,:,1) ) 1498 ENDIF 1499 IF( srcv(jpr_isf)%laction ) THEN 1500 fwfisf(:,:) = - frcv(jpr_isf)%z3(:,:,1) ! fresh water flux from the isf (fwfisf <0 mean melting) 1501 CALL iom_put( 'iceshelf_cea', frcv(jpr_isf)%z3(:,:,1) ) 1502 ENDIF 1503 1504 1470 1505 IF( ln_mixcpl ) THEN 1471 1506 emp_tot(:,:) = emp_tot(:,:) * xcplmask(:,:,0) + zemp_tot(:,:) * zmsk(:,:) … … 1504 1539 ENDIF 1505 1540 1541 1542 IF( srcv(jpr_icb)%laction ) THEN 1543 fwficb(:,:) = frcv(jpr_icb)%z3(:,:,1) 1544 rnf(:,:) = rnf(:,:) + fwficb(:,:) ! iceberg added to runoffs 1545 CALL iom_put( 'iceberg_cea', frcv(jpr_icb)%z3(:,:,1) ) 1546 ENDIF 1547 IF( srcv(jpr_isf)%laction ) THEN 1548 fwfisf(:,:) = - frcv(jpr_isf)%z3(:,:,1) ! fresh water flux from the isf (fwfisf <0 mean melting) 1549 CALL iom_put( 'iceshelf_cea', frcv(jpr_isf)%z3(:,:,1) ) 1550 ENDIF 1551 1552 1506 1553 IF( ln_mixcpl ) THEN 1507 1554 emp_tot(:,:) = emp_tot(:,:) * xcplmask(:,:,0) + zemp_tot(:,:) * zmsk(:,:) … … 1570 1617 IF( iom_use('hflx_cal_cea') ) CALL iom_put( 'hflx_cal_cea', - frcv(jpr_cal)%z3(:,:,1) * lfus ) ! heat flux from calving 1571 1618 ENDIF 1619 1620 !!chris 1621 !! The heat content associated to the ice shelf in removed in the routine sbcisf.F90 1622 ! 1623 IF( srcv(jpr_icb)%laction ) zqns_tot(:,:) = zqns_tot(:,:) - frcv(jpr_icb)%z3(:,:,1) * lfus ! remove heat content associated to iceberg melting 1624 ! 1625 !! ! 1572 1626 1573 1627 #if defined key_lim3 -
branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/OPA_SRC/SBC/sbcisf.F90
r7494 r7607 55 55 INTEGER, PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:) :: misfkt, misfkb !:Level of ice shelf base 56 56 57 LOGICAL, PUBLIC :: l_isfcpl = .false. ! isf recieved from oasis 58 57 59 58 60 REAL(wp), PUBLIC, SAVE :: rcpi = 2000.0_wp ! phycst ? … … 94 96 REAL(wp), DIMENSION(:,:,:), POINTER :: zfwfisf3d, zqhcisf3d, zqlatisf3d 95 97 REAL(wp), DIMENSION(:,: ), POINTER :: zqhcisf2d 98 REAL(wp) :: zhisf 99 96 100 ! 97 101 !!--------------------------------------------------------------------- … … 142 146 misfkt(:,:) = mikt(:,:) ! same indice for bg03 et cav => used in isfdiv 143 147 ELSE IF ((nn_isf == 3) .OR. (nn_isf == 2)) THEN 144 ALLOCATE( sf_rnfisf(1), STAT=ierror ) 145 ALLOCATE( sf_rnfisf(1)%fnow(jpi,jpj,1), sf_rnfisf(1)%fdta(jpi,jpj,1,2) ) 146 CALL fld_fill( sf_rnfisf, (/ sn_rnfisf /), cn_dirisf, 'sbc_isf_init', 'read fresh water flux isf data', 'namsbc_isf' ) 148 IF( .NOT.l_isfcpl ) THEN 149 ALLOCATE( sf_rnfisf(1), STAT=ierror ) 150 ALLOCATE( sf_rnfisf(1)%fnow(jpi,jpj,1), sf_rnfisf(1)%fdta(jpi,jpj,1,2) ) 151 CALL fld_fill( sf_rnfisf, (/ sn_rnfisf /), cn_dirisf, 'sbc_isf_init', 'read fresh water flux isf data', 'namsbc_isf' ) 152 ENDIF 147 153 148 154 !: read effective lenght (BG03) … … 185 191 186 192 ! load variable used in fldread (use for temporal interpolation of isf fwf forcing) 187 ALLOCATE( sf_fwfisf(1), sf_qisf(1), STAT=ierror ) 188 ALLOCATE( sf_fwfisf(1)%fnow(jpi,jpj,1), sf_fwfisf(1)%fdta(jpi,jpj,1,2) ) 189 ALLOCATE( sf_qisf(1)%fnow(jpi,jpj,1), sf_qisf(1)%fdta(jpi,jpj,1,2) ) 190 CALL fld_fill( sf_fwfisf, (/ sn_fwfisf /), cn_dirisf, 'sbc_isf_init', 'read fresh water flux isf data', 'namsbc_isf' ) 191 !CALL fld_fill( sf_qisf , (/ sn_qisf /), cn_dirisf, 'sbc_isf_init', 'read heat flux isf data' , 'namsbc_isf' ) 193 IF( .NOT.l_isfcpl ) THEN 194 ALLOCATE( sf_fwfisf(1), sf_qisf(1), STAT=ierror ) 195 ALLOCATE( sf_fwfisf(1)%fnow(jpi,jpj,1), sf_fwfisf(1)%fdta(jpi,jpj,1,2) ) 196 ALLOCATE( sf_qisf(1)%fnow(jpi,jpj,1), sf_qisf(1)%fdta(jpi,jpj,1,2) ) 197 CALL fld_fill( sf_fwfisf, (/ sn_fwfisf /), cn_dirisf, 'sbc_isf_init', 'read fresh water flux isf data', 'namsbc_isf' ) 198 !CALL fld_fill( sf_qisf , (/ sn_qisf /), cn_dirisf, 'sbc_isf_init', 'read heat flux isf data' , 'namsbc_isf' ) 199 ENDIF 192 200 END IF 193 201 … … 242 250 CALL iom_put('vtbl',vtbl(:,:)) 243 251 ! compute fwf and heat flux 244 CALL sbc_isf_cav (kt) 252 IF( .NOT.l_isfcpl ) THEN ; CALL sbc_isf_cav (kt) 253 ELSE ; qisf(:,:) = fwfisf(:,:) * lfusisf ! heat flux 254 ENDIF 245 255 246 256 ELSE IF (nn_isf == 2) THEN … … 251 261 ELSE IF (nn_isf == 3) THEN 252 262 ! specified runoff in depth (Mathiot et al., XXXX in preparation) 253 CALL fld_read ( kt, nn_fsbc, sf_rnfisf ) 254 fwfisf(:,:) = - sf_rnfisf(1)%fnow(:,:,1) ! fresh water flux from the isf (fwfisf <0 mean melting) 263 IF( .NOT.l_isfcpl ) THEN 264 CALL fld_read ( kt, nn_fsbc, sf_rnfisf ) 265 fwfisf(:,:) = - sf_rnfisf(1)%fnow(:,:,1) ! fresh water flux from the isf (fwfisf <0 mean melting) 266 ENDIF 255 267 qisf(:,:) = fwfisf(:,:) * lfusisf ! heat flux 256 268 stbl(:,:) = soce … … 258 270 ELSE IF (nn_isf == 4) THEN 259 271 ! specified fwf and heat flux forcing beneath the ice shelf 260 CALL fld_read ( kt, nn_fsbc, sf_fwfisf ) 261 !CALL fld_read ( kt, nn_fsbc, sf_qisf ) 262 fwfisf(:,:) = sf_fwfisf(1)%fnow(:,:,1) ! fwf 272 IF( .NOT.l_isfcpl ) THEN 273 CALL fld_read ( kt, nn_fsbc, sf_fwfisf ) 274 !CALL fld_read ( kt, nn_fsbc, sf_qisf ) 275 fwfisf(:,:) = sf_fwfisf(1)%fnow(:,:,1) ! fwf 276 ENDIF 263 277 qisf(:,:) = fwfisf(:,:) * lfusisf ! heat flux 264 278 !qisf(:,:) = sf_qisf(1)%fnow(:,:,1) ! heat flux … … 288 302 CALL lbc_lnk(qisf(:,:) ,'T',1.) 289 303 290 !============================================================================================================================================= 291 IF ( iom_use('fwfisf3d') .OR. iom_use('qlatisf3d') .OR. iom_use('qhcisf3d') .OR. iom_use('qhcisf')) THEN 304 ! Diagnostics 305 IF( iom_use('fwfisf3d') .OR. iom_use('qlatisf3d') .OR. iom_use('qhcisf3d') .OR. iom_use('qhcisf')) THEN 306 ! 292 307 CALL wrk_alloc( jpi,jpj,jpk, zfwfisf3d, zqhcisf3d, zqlatisf3d ) 293 308 CALL wrk_alloc( jpi,jpj, zqhcisf2d ) 294 309 ! 295 310 zfwfisf3d(:,:,:) = 0.0_wp ! 3d ice shelf melting (kg/m2/s) 296 311 zqhcisf3d(:,:,:) = 0.0_wp ! 3d heat content flux (W/m2) 297 312 zqlatisf3d(:,:,:)= 0.0_wp ! 3d ice shelf melting latent heat flux (W/m2) 298 313 zqhcisf2d(:,:) = fwfisf(:,:) * zt_frz * rcp ! 2d heat content flux (W/m2) 299 314 ! 300 315 DO jj = 1,jpj 301 316 DO ji = 1,jpi … … 303 318 ikb = misfkb(ji,jj) 304 319 DO jk = ikt, ikb - 1 305 zfwfisf3d (ji,jj,jk) = zfwfisf3d (ji,jj,jk) + fwfisf (ji,jj) * r1_hisf_tbl(ji,jj) * fse3t(ji,jj,jk) 306 zqhcisf3d (ji,jj,jk) = zqhcisf3d (ji,jj,jk) + zqhcisf2d(ji,jj) * r1_hisf_tbl(ji,jj) * fse3t(ji,jj,jk) 307 zqlatisf3d(ji,jj,jk) = zqlatisf3d(ji,jj,jk) + qisf (ji,jj) * r1_hisf_tbl(ji,jj) * fse3t(ji,jj,jk) 320 zhisf = r1_hisf_tbl(ji,jj) * fse3t(ji,jj,jk) 321 zfwfisf3d (ji,jj,jk) = zfwfisf3d (ji,jj,jk) + fwfisf(ji,jj) * zhisf 322 zqhcisf3d (ji,jj,jk) = zqhcisf3d (ji,jj,jk) + zqhcisf2d(ji,jj) * zhisf 323 zqlatisf3d(ji,jj,jk) = zqlatisf3d(ji,jj,jk) + qisf(ji,jj) * zhisf 308 324 END DO 309 zfwfisf3d (ji,jj,jk) = zfwfisf3d (ji,jj,jk) + fwfisf (ji,jj) * r1_hisf_tbl(ji,jj) * ralpha(ji,jj) * fse3t(ji,jj,jk) 310 zqhcisf3d (ji,jj,jk) = zqhcisf3d (ji,jj,jk) + zqhcisf2d(ji,jj) * r1_hisf_tbl(ji,jj) * ralpha(ji,jj) * fse3t(ji,jj,jk) 311 zqlatisf3d(ji,jj,jk) = zqlatisf3d(ji,jj,jk) + qisf (ji,jj) * r1_hisf_tbl(ji,jj) * ralpha(ji,jj) * fse3t(ji,jj,jk) 325 jk = ikb 326 zhisf = r1_hisf_tbl(ji,jj) * fse3t(ji,jj,jk) 327 zfwfisf3d (ji,jj,jk) = zfwfisf3d (ji,jj,jk) + fwfisf (ji,jj) * zhisf * ralpha(ji,jj) 328 zqhcisf3d (ji,jj,jk) = zqhcisf3d (ji,jj,jk) + zqhcisf2d(ji,jj) * zhisf * ralpha(ji,jj) 329 zqlatisf3d(ji,jj,jk) = zqlatisf3d(ji,jj,jk) + qisf (ji,jj) * zhisf * ralpha(ji,jj) 312 330 END DO 313 331 END DO 314 315 CALL iom_put( 'fwfisf3d' , zfwfisf3d (:,:,:))316 CALL iom_put( 'qlatisf3d', zqlatisf3d(:,:,:))317 CALL iom_put( 'qhcisf3d' , zqhcisf3d (:,:,:))318 CALL iom_put( 'qhcisf' , zqhcisf2d (:,: ))319 332 ! 333 CALL iom_put( 'fwfisf3d' , zfwfisf3d (:,:,:) ) 334 CALL iom_put( 'qlatisf3d', zqlatisf3d(:,:,:) ) 335 CALL iom_put( 'qhcisf3d' , zqhcisf3d (:,:,:) ) 336 CALL iom_put( 'qhcisf' , zqhcisf2d (:,: ) ) 337 ! 320 338 CALL wrk_dealloc( jpi,jpj,jpk, zfwfisf3d, zqhcisf3d, zqlatisf3d ) 321 339 CALL wrk_dealloc( jpi,jpj, zqhcisf2d ) 340 ! 322 341 END IF 323 !============================================================================================================================================= 324 342 ! 325 343 IF( kt == nit000 ) THEN ! set the forcing field at nit000 - 1 ! 326 344 IF( ln_rstart .AND. & ! Restart: read in restart file -
branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/OPA_SRC/stpctl.F90
r6204 r7607 180 180 ENDIF 181 181 182 9200 FORMAT('it:', i8, ' iter:', i4, ' r: ', e16.10, ' b: ',e16.10)183 9300 FORMAT(' it :', i8, ' ssh2: ', e16.10, ' Umax: ',e16.10,' Smin: ',e16.10)182 9200 FORMAT('it:', i8, ' iter:', i4, ' r: ',d23.16, ' b: ',d23.16 ) 183 9300 FORMAT(' it :', i8, ' ssh2: ', d23.16, ' Umax: ',d23.16,' Smin: ',d23.16) 184 184 ! 185 185 END SUBROUTINE stp_ctl -
branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zche.F90
r6943 r7607 11 11 !! 2.0 ! 2007-12 (C. Ethe, G. Madec) F90 12 12 !! ! 2011-02 (J. Simeon, J.Orr ) update O2 solubility constants 13 !! 3.6 ! 2016-03 (O. Aumont) Change chemistry to MOCSY standards 13 14 !!---------------------------------------------------------------------- 14 15 #if defined key_pisces 15 16 !!---------------------------------------------------------------------- 16 !! 'key_pisces 'PISCES bio-model17 !! 'key_pisces*' PISCES bio-model 17 18 !!---------------------------------------------------------------------- 18 19 !! p4z_che : Sea water chemistry computed following OCMIP protocol … … 22 23 USE sms_pisces ! PISCES Source Minus Sink variables 23 24 USE lib_mpp ! MPP library 25 USE eosbn2, ONLY : nn_eos 24 26 25 27 IMPLICIT NONE 26 28 PRIVATE 27 29 28 PUBLIC p4z_che ! 29 PUBLIC p4z_che_alloc ! 30 31 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: sio3eq ! chemistry of Si 32 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: fekeq ! chemistry of Fe 33 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: chemc ! Solubilities of O2 and CO2 34 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: chemo2 ! Solubilities of O2 and CO2 30 PUBLIC p4z_che ! 31 PUBLIC p4z_che_alloc ! 32 PUBLIC p4z_che_ahini ! 33 PUBLIC p4z_che_solve_hi ! 34 35 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: sio3eq ! chemistry of Si 36 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: fekeq ! chemistry of Fe 37 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: chemc ! Solubilities of O2 and CO2 38 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: chemo2 ! Solubilities of O2 and CO2 39 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: fesol ! solubility of Fe 35 40 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: tempis ! In situ temperature 41 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: salinprac ! Practical salinity 42 43 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: akb3 !: ??? 44 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: akw3 !: ??? 45 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: akf3 !: ??? 46 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: aks3 !: ??? 47 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: ak1p3 !: ??? 48 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: ak2p3 !: ??? 49 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: ak3p3 !: ??? 50 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: aksi3 !: ??? 51 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: borat !: ??? 52 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: fluorid !: ??? 53 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: sulfat !: ??? 54 55 !!* Variable for chemistry of the CO2 cycle 36 56 37 57 REAL(wp), PUBLIC :: atcox = 0.20946 ! units atm 38 58 39 REAL(wp) :: salchl = 1. / 1.80655 ! conversion factor for salinity --> chlorinity (Wooster et al. 1969)40 59 REAL(wp) :: o2atm = 1. / ( 1000. * 0.20946 ) 41 60 42 REAL(wp) :: rgas = 83.14472 ! universal gas constants 43 REAL(wp) :: oxyco = 1. / 22.4144 ! converts from liters of an ideal gas to moles 44 45 REAL(wp) :: bor1 = 0.00023 ! borat constants 46 REAL(wp) :: bor2 = 1. / 10.82 47 48 REAL(wp) :: st1 = 0.14 ! constants for calculate concentrations for sulfate 49 REAL(wp) :: st2 = 1./96.062 ! (Morris & Riley 1966) 50 51 REAL(wp) :: ft1 = 0.000067 ! constants for calculate concentrations for fluorides 52 REAL(wp) :: ft2 = 1./18.9984 ! (Dickson & Riley 1979 ) 53 54 ! ! volumetric solubility constants for o2 in ml/L 55 REAL(wp) :: ox0 = 2.00856 ! from Table 1 for Eq 8 of Garcia and Gordon, 1992. 56 REAL(wp) :: ox1 = 3.22400 ! corrects for moisture and fugacity, but not total atmospheric pressure 57 REAL(wp) :: ox2 = 3.99063 ! Original PISCES code noted this was a solubility, but 58 REAL(wp) :: ox3 = 4.80299 ! was in fact a bunsen coefficient with units L-O2/(Lsw atm-O2) 59 REAL(wp) :: ox4 = 9.78188e-1 ! Hence, need to divide EXP( zoxy ) by 1000, ml-O2 => L-O2 60 REAL(wp) :: ox5 = 1.71069 ! and atcox = 0.20946 to add the 1/atm dimension. 61 REAL(wp) :: ox6 = -6.24097e-3 62 REAL(wp) :: ox7 = -6.93498e-3 63 REAL(wp) :: ox8 = -6.90358e-3 64 REAL(wp) :: ox9 = -4.29155e-3 65 REAL(wp) :: ox10 = -3.11680e-7 61 REAL(wp) :: rgas = 83.14472 ! universal gas constants 62 REAL(wp) :: oxyco = 1. / 22.4144 ! converts from liters of an ideal gas to moles 66 63 67 64 ! ! coeff. for seawater pressure correction : millero 95 68 65 ! ! AGRIF doesn't like the DATA instruction 69 REAL(wp) :: devk11 = -25.5 70 REAL(wp) :: devk12 = -15.82 71 REAL(wp) :: devk13 = -29.48 72 REAL(wp) :: devk14 = -25.60 73 REAL(wp) :: devk15 = -48.76 66 REAL(wp) :: devk10 = -25.5 67 REAL(wp) :: devk11 = -15.82 68 REAL(wp) :: devk12 = -29.48 69 REAL(wp) :: devk13 = -20.02 70 REAL(wp) :: devk14 = -18.03 71 REAL(wp) :: devk15 = -9.78 72 REAL(wp) :: devk16 = -48.76 73 REAL(wp) :: devk17 = -14.51 74 REAL(wp) :: devk18 = -23.12 75 REAL(wp) :: devk19 = -26.57 76 REAL(wp) :: devk110 = -29.48 74 77 ! 75 REAL(wp) :: devk21 = 0.1271 76 REAL(wp) :: devk22 = -0.0219 77 REAL(wp) :: devk23 = 0.1622 78 REAL(wp) :: devk24 = 0.2324 79 REAL(wp) :: devk25 = 0.5304 78 REAL(wp) :: devk20 = 0.1271 79 REAL(wp) :: devk21 = -0.0219 80 REAL(wp) :: devk22 = 0.1622 81 REAL(wp) :: devk23 = 0.1119 82 REAL(wp) :: devk24 = 0.0466 83 REAL(wp) :: devk25 = -0.0090 84 REAL(wp) :: devk26 = 0.5304 85 REAL(wp) :: devk27 = 0.1211 86 REAL(wp) :: devk28 = 0.1758 87 REAL(wp) :: devk29 = 0.2020 88 REAL(wp) :: devk210 = 0.1622 80 89 ! 90 REAL(wp) :: devk30 = 0. 81 91 REAL(wp) :: devk31 = 0. 82 REAL(wp) :: devk32 = 0. 83 REAL(wp) :: devk33 = 2.608E-3 84 REAL(wp) :: devk34 = -3.6246E-3 85 REAL(wp) :: devk35 = 0. 92 REAL(wp) :: devk32 = 2.608E-3 93 REAL(wp) :: devk33 = -1.409e-3 94 REAL(wp) :: devk34 = 0.316e-3 95 REAL(wp) :: devk35 = -0.942e-3 96 REAL(wp) :: devk36 = 0. 97 REAL(wp) :: devk37 = -0.321e-3 98 REAL(wp) :: devk38 = -2.647e-3 99 REAL(wp) :: devk39 = -3.042e-3 100 REAL(wp) :: devk310 = -2.6080e-3 86 101 ! 87 REAL(wp) :: devk41 = -3.08E-3 88 REAL(wp) :: devk42 = 1.13E-3 89 REAL(wp) :: devk43 = -2.84E-3 90 REAL(wp) :: devk44 = -5.13E-3 91 REAL(wp) :: devk45 = -11.76E-3 102 REAL(wp) :: devk40 = -3.08E-3 103 REAL(wp) :: devk41 = 1.13E-3 104 REAL(wp) :: devk42 = -2.84E-3 105 REAL(wp) :: devk43 = -5.13E-3 106 REAL(wp) :: devk44 = -4.53e-3 107 REAL(wp) :: devk45 = -3.91e-3 108 REAL(wp) :: devk46 = -11.76e-3 109 REAL(wp) :: devk47 = -2.67e-3 110 REAL(wp) :: devk48 = -5.15e-3 111 REAL(wp) :: devk49 = -4.08e-3 112 REAL(wp) :: devk410 = -2.84e-3 92 113 ! 93 REAL(wp) :: devk51 = 0.0877E-3 94 REAL(wp) :: devk52 = -0.1475E-3 95 REAL(wp) :: devk53 = 0. 96 REAL(wp) :: devk54 = 0.0794E-3 97 REAL(wp) :: devk55 = 0.3692E-3 114 REAL(wp) :: devk50 = 0.0877E-3 115 REAL(wp) :: devk51 = -0.1475E-3 116 REAL(wp) :: devk52 = 0. 117 REAL(wp) :: devk53 = 0.0794E-3 118 REAL(wp) :: devk54 = 0.09e-3 119 REAL(wp) :: devk55 = 0.054e-3 120 REAL(wp) :: devk56 = 0.3692E-3 121 REAL(wp) :: devk57 = 0.0427e-3 122 REAL(wp) :: devk58 = 0.09e-3 123 REAL(wp) :: devk59 = 0.0714e-3 124 REAL(wp) :: devk510 = 0.0 125 ! 126 ! General parameters 127 REAL(wp), PARAMETER :: pp_rdel_ah_target = 1.E-4_wp 128 REAL(wp), PARAMETER :: pp_ln10 = 2.302585092994045684018_wp 129 130 ! Maximum number of iterations for each method 131 INTEGER, PARAMETER :: jp_maxniter_atgen = 20 132 133 ! Bookkeeping variables for each method 134 ! - SOLVE_AT_GENERAL 135 INTEGER :: niter_atgen = jp_maxniter_atgen 98 136 99 137 !!* Substitution … … 115 153 !!--------------------------------------------------------------------- 116 154 INTEGER :: ji, jj, jk 117 REAL(wp) :: ztkel, zt , zt2, zsal , zsal2 , zbuf1 , zbuf2155 REAL(wp) :: ztkel, ztkel1, zt , zsal , zsal2 , zbuf1 , zbuf2 118 156 REAL(wp) :: ztgg , ztgg2, ztgg3 , ztgg4 , ztgg5 119 157 REAL(wp) :: zpres, ztc , zcl , zcpexp, zoxy , zcpexp2 120 158 REAL(wp) :: zsqrt, ztr , zlogt , zcek1, zc1, zplat 121 REAL(wp) :: zis , zis2 , zsal15, zisqrt, za1 159 REAL(wp) :: zis , zis2 , zsal15, zisqrt, za1, za2 122 160 REAL(wp) :: zckb , zck1 , zck2 , zckw , zak1 , zak2 , zakb , zaksp0, zakw 161 REAL(wp) :: zck1p, zck2p, zck3p, zcksi, zak1p, zak2p, zak3p, zaksi 123 162 REAL(wp) :: zst , zft , zcks , zckf , zaksp1 163 REAL(wp) :: total2free, free2SWS, total2SWS, SWS2total 164 124 165 !!--------------------------------------------------------------------- 125 166 ! 126 167 IF( nn_timing == 1 ) CALL timing_start('p4z_che') 168 ! 169 ! Computation of chemical constants require practical salinity 170 ! Thus, when TEOS08 is used, absolute salinity is converted to 171 ! practical salinity 172 ! ------------------------------------------------------------- 173 IF (nn_eos == -1) THEN 174 salinprac(:,:,:) = tsn(:,:,:,jp_sal) * 35.0 / 35.16504 175 ELSE 176 salinprac(:,:,:) = tsn(:,:,:,jp_sal) 177 ENDIF 178 127 179 ! 128 180 ! Computations of chemical constants require in situ temperature … … 135 187 DO ji = 1, jpi 136 188 zpres = fsdept(ji,jj,jk) / 1000. 137 za1 = 0.04 * ( 1.0 + 0.185 * tsn(ji,jj,jk,jp_tem) + 0.035 * ( tsn(ji,jj,jk,jp_sal) - 35.0) )189 za1 = 0.04 * ( 1.0 + 0.185 * tsn(ji,jj,jk,jp_tem) + 0.035 * (salinprac(ji,jj,jk) - 35.0) ) 138 190 za2 = 0.0075 * ( 1.0 - tsn(ji,jj,jk,jp_tem) / 30.0 ) 139 191 tempis(ji,jj,jk) = tsn(ji,jj,jk,jp_tem) - za1 * zpres + za2 * zpres**2 … … 151 203 ztkel = tempis(ji,jj,1) + 273.15 152 204 zt = ztkel * 0.01 153 zt2 = zt * zt 154 zsal = tsn(ji,jj,1,jp_sal) + ( 1.- tmask(ji,jj,1) ) * 35. 155 zsal2 = zsal * zsal 156 zlogt = LOG( zt ) 205 zsal = salinprac(ji,jj,1) + ( 1.- tmask(ji,jj,1) ) * 35. 157 206 ! ! LN(K0) OF SOLUBILITY OF CO2 (EQ. 12, WEISS, 1980) 158 207 ! ! AND FOR THE ATMOSPHERE FOR NON IDEAL GAS 159 208 zcek1 = 9345.17/ztkel - 60.2409 + 23.3585 * LOG(zt) + zsal*(0.023517 - 0.00023656*ztkel & 160 209 & + 0.0047036e-4*ztkel**2) 161 ! ! SET SOLUBILITIES OF O2 AND CO2 162 chemc(ji,jj,1) = EXP( zcek1 ) * 1.e-6 * rhop(ji,jj,1) / 1000. ! mol/(kg uatm) 210 chemc(ji,jj,1) = EXP( zcek1 ) * 1E-6 * rhop(ji,jj,1) / 1000. ! mol/(L atm) 163 211 chemc(ji,jj,2) = -1636.75 + 12.0408*ztkel - 0.0327957*ztkel**2 + 0.0000316528*ztkel**3 164 212 chemc(ji,jj,3) = 57.7 - 0.118*ztkel … … 176 224 DO ji = 1, jpi 177 225 ztkel = tempis(ji,jj,jk) + 273.15 178 zsal = tsn(ji,jj,jk,jp_sal) + ( 1.- tmask(ji,jj,jk) ) * 35.226 zsal = salinprac(ji,jj,jk) + ( 1.- tmask(ji,jj,jk) ) * 35. 179 227 zsal2 = zsal * zsal 180 228 ztgg = LOG( ( 298.15 - tempis(ji,jj,jk) ) / ztkel ) ! Set the GORDON & GARCIA scaled temperature … … 183 231 ztgg4 = ztgg3 * ztgg 184 232 ztgg5 = ztgg4 * ztgg 185 zoxy = ox0 + ox1 * ztgg + ox2 * ztgg2 + ox3 * ztgg3 + ox4 * ztgg4 + ox5 * ztgg5 & 186 + zsal * ( ox6 + ox7 * ztgg + ox8 * ztgg2 + ox9 * ztgg3 ) + ox10 * zsal2 233 234 zoxy = 2.00856 + 3.22400 * ztgg + 3.99063 * ztgg2 + 4.80299 * ztgg3 & 235 & + 9.78188e-1 * ztgg4 + 1.71069 * ztgg5 + zsal * ( -6.24097e-3 & 236 & - 6.93498e-3 * ztgg - 6.90358e-3 * ztgg2 - 4.29155e-3 * ztgg3 ) & 237 & - 3.11680e-7 * zsal2 187 238 chemo2(ji,jj,jk) = ( EXP( zoxy ) * o2atm ) * oxyco * atcox ! mol/(L atm) 188 239 END DO … … 209 260 ! SET ABSOLUTE TEMPERATURE 210 261 ztkel = tempis(ji,jj,jk) + 273.15 211 zsal = tsn(ji,jj,jk,jp_sal) + ( 1.-tmask(ji,jj,jk) ) * 35.262 zsal = salinprac(ji,jj,jk) + ( 1.-tmask(ji,jj,jk) ) * 35. 212 263 zsqrt = SQRT( zsal ) 213 264 zsal15 = zsqrt * zsal … … 220 271 221 272 ! CHLORINITY (WOOSTER ET AL., 1969) 222 zcl = zsal * salchl273 zcl = zsal / 1.80655 223 274 224 275 ! TOTAL SULFATE CONCENTR. [MOLES/kg soln] 225 zst = st1 * zcl * st2276 zst = 0.14 * zcl /96.062 226 277 227 278 ! TOTAL FLUORIDE CONCENTR. [MOLES/kg soln] 228 zft = ft1 * zcl * ft2279 zft = 0.000067 * zcl /18.9984 229 280 230 281 ! DISSOCIATION CONSTANT FOR SULFATES on free H scale (Dickson 1990) … … 234 285 & - 2698. * ztr * zis**1.5 + 1776.* ztr * zis2 & 235 286 & + LOG(1.0 - 0.001005 * zsal)) 236 !237 aphscale(ji,jj,jk) = ( 1. + zst / zcks )238 287 239 288 ! DISSOCIATION CONSTANT FOR FLUORIDES on free H scale (Dickson and Riley 79) … … 249 298 & * zlogt + 0.053105*zsqrt*ztkel 250 299 251 252 300 ! DISSOCIATION COEFFICIENT FOR CARBONATE ACCORDING TO 253 301 ! MEHRBACH (1973) REFIT BY MILLERO (1995), seawater scale … … 257 305 - 0.01781*zsal + 0.0001122*zsal*zsal) 258 306 259 ! PKW (H2O) (DICKSON AND RILEY, 1979) 260 zckw = -13847.26*ztr + 148.9652 - 23.6521 * zlogt & 261 & + (118.67*ztr - 5.977 + 1.0495 * zlogt) & 262 & * zsqrt - 0.01615 * zsal 307 ! PKW (H2O) (MILLERO, 1995) from composite data 308 zckw = -13847.26 * ztr + 148.9652 - 23.6521 * zlogt + ( 118.67 * ztr & 309 - 5.977 + 1.0495 * zlogt ) * zsqrt - 0.01615 * zsal 310 311 ! CONSTANTS FOR PHOSPHATE (MILLERO, 1995) 312 zck1p = -4576.752*ztr + 115.540 - 18.453*zlogt & 313 & + (-106.736*ztr + 0.69171) * zsqrt & 314 & + (-0.65643*ztr - 0.01844) * zsal 315 316 zck2p = -8814.715*ztr + 172.1033 - 27.927*zlogt & 317 & + (-160.340*ztr + 1.3566)*zsqrt & 318 & + (0.37335*ztr - 0.05778)*zsal 319 320 zck3p = -3070.75*ztr - 18.126 & 321 & + (17.27039*ztr + 2.81197) * zsqrt & 322 & + (-44.99486*ztr - 0.09984) * zsal 323 324 ! CONSTANT FOR SILICATE, MILLERO (1995) 325 zcksi = -8904.2*ztr + 117.400 - 19.334*zlogt & 326 & + (-458.79*ztr + 3.5913) * zisqrt & 327 & + (188.74*ztr - 1.5998) * zis & 328 & + (-12.1652*ztr + 0.07871) * zis2 & 329 & + LOG(1.0 - 0.001005*zsal) 263 330 264 331 ! APPARENT SOLUBILITY PRODUCT K'SP OF CALCITE IN SEAWATER … … 268 335 & - 0.07711*zsal + 0.0041249*zsal15 269 336 337 ! CONVERT FROM DIFFERENT PH SCALES 338 total2free = 1.0/(1.0 + zst/zcks) 339 free2SWS = 1. + zst/zcks + zft/(zckf*total2free) 340 total2SWS = total2free * free2SWS 341 SWS2total = 1.0 / total2SWS 342 270 343 ! K1, K2 OF CARBONIC ACID, KB OF BORIC ACID, KW (H2O) (LIT.?) 271 zak1 = 10**(zck1) 272 zak2 = 10**(zck2) 273 zakb = EXP( zckb )344 zak1 = 10**(zck1) * total2SWS 345 zak2 = 10**(zck2) * total2SWS 346 zakb = EXP( zckb ) * total2SWS 274 347 zakw = EXP( zckw ) 275 348 zaksp1 = 10**(zaksp0) 349 zak1p = exp( zck1p ) 350 zak2p = exp( zck2p ) 351 zak3p = exp( zck3p ) 352 zaksi = exp( zcksi ) 353 zckf = zckf * total2SWS 276 354 277 355 ! FORMULA FOR CPEXP AFTER EDMOND & GIESKES (1970) … … 285 363 ! FORMULA ON P. 1286 IS RIGHT AND CONSISTENT WITH THE 286 364 ! SIGN IN PARTIAL MOLAR VOLUME CHANGE AS SHOWN ON P. 1285)) 287 zcpexp = zpres / (rgas*ztkel)288 zcpexp2 = zpres * z pres/(rgas*ztkel)365 zcpexp = zpres / (rgas*ztkel) 366 zcpexp2 = zpres * zcpexp 289 367 290 368 ! KB OF BORIC ACID, K1,K2 OF CARBONIC ACID PRESSURE … … 292 370 ! (CF. BROECKER ET AL., 1982) 293 371 294 zbuf1 = - ( devk11 + devk21 * ztc + devk31 * ztc * ztc ) 372 zbuf1 = - ( devk10 + devk20 * ztc + devk30 * ztc * ztc ) 373 zbuf2 = 0.5 * ( devk40 + devk50 * ztc ) 374 ak13(ji,jj,jk) = zak1 * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 ) 375 376 zbuf1 = - ( devk11 + devk21 * ztc + devk31 * ztc * ztc ) 295 377 zbuf2 = 0.5 * ( devk41 + devk51 * ztc ) 296 ak 13(ji,jj,jk) = zak1* EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )378 ak23(ji,jj,jk) = zak2 * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 ) 297 379 298 380 zbuf1 = - ( devk12 + devk22 * ztc + devk32 * ztc * ztc ) 299 381 zbuf2 = 0.5 * ( devk42 + devk52 * ztc ) 300 ak 23(ji,jj,jk) = zak2* EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )382 akb3(ji,jj,jk) = zakb * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 ) 301 383 302 384 zbuf1 = - ( devk13 + devk23 * ztc + devk33 * ztc * ztc ) 303 385 zbuf2 = 0.5 * ( devk43 + devk53 * ztc ) 304 ak b3(ji,jj,jk) = zakb* EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )386 akw3(ji,jj,jk) = zakw * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 ) 305 387 306 388 zbuf1 = - ( devk14 + devk24 * ztc + devk34 * ztc * ztc ) 307 389 zbuf2 = 0.5 * ( devk44 + devk54 * ztc ) 308 akw3(ji,jj,jk) = zakw * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 ) 309 390 aks3(ji,jj,jk) = zcks * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 ) 391 392 zbuf1 = - ( devk15 + devk25 * ztc + devk35 * ztc * ztc ) 393 zbuf2 = 0.5 * ( devk45 + devk55 * ztc ) 394 akf3(ji,jj,jk) = zckf * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 ) 395 396 zbuf1 = - ( devk17 + devk27 * ztc + devk37 * ztc * ztc ) 397 zbuf2 = 0.5 * ( devk47 + devk57 * ztc ) 398 ak1p3(ji,jj,jk) = zak1p * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 ) 399 400 zbuf1 = - ( devk18 + devk28 * ztc + devk38 * ztc * ztc ) 401 zbuf2 = 0.5 * ( devk48 + devk58 * ztc ) 402 ak2p3(ji,jj,jk) = zak2p * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 ) 403 404 zbuf1 = - ( devk19 + devk29 * ztc + devk39 * ztc * ztc ) 405 zbuf2 = 0.5 * ( devk49 + devk59 * ztc ) 406 ak3p3(ji,jj,jk) = zak3p * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 ) 407 408 zbuf1 = - ( devk110 + devk210 * ztc + devk310 * ztc * ztc ) 409 zbuf2 = 0.5 * ( devk410 + devk510 * ztc ) 410 aksi3(ji,jj,jk) = zaksi * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 ) 411 412 ! CONVERT FROM DIFFERENT PH SCALES 413 total2free = 1.0/(1.0 + zst/aks3(ji,jj,jk)) 414 free2SWS = 1. + zst/aks3(ji,jj,jk) + zft/akf3(ji,jj,jk) 415 total2SWS = total2free * free2SWS 416 SWS2total = 1.0 / total2SWS 417 418 ! Convert to total scale 419 ak13(ji,jj,jk) = ak13(ji,jj,jk) * SWS2total 420 ak23(ji,jj,jk) = ak23(ji,jj,jk) * SWS2total 421 akb3(ji,jj,jk) = akb3(ji,jj,jk) * SWS2total 422 akw3(ji,jj,jk) = akw3(ji,jj,jk) * SWS2total 423 ak1p3(ji,jj,jk) = ak1p3(ji,jj,jk) * SWS2total 424 ak2p3(ji,jj,jk) = ak2p3(ji,jj,jk) * SWS2total 425 ak3p3(ji,jj,jk) = ak3p3(ji,jj,jk) * SWS2total 426 aksi3(ji,jj,jk) = aksi3(ji,jj,jk) * SWS2total 427 akf3(ji,jj,jk) = akf3(ji,jj,jk) / total2free 310 428 311 429 ! APPARENT SOLUBILITY PRODUCT K'SP OF CALCITE 312 430 ! AS FUNCTION OF PRESSURE FOLLOWING MILLERO 313 431 ! (P. 1285) AND BERNER (1976) 314 zbuf1 = - ( devk1 5 + devk25 * ztc + devk35* ztc * ztc )315 zbuf2 = 0.5 * ( devk4 5 + devk55* ztc )432 zbuf1 = - ( devk16 + devk26 * ztc + devk36 * ztc * ztc ) 433 zbuf2 = 0.5 * ( devk46 + devk56 * ztc ) 316 434 aksp(ji,jj,jk) = zaksp1 * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 ) 317 435 318 ! TOTAL BORATE CONCENTR. [MOLES/L] 319 borat(ji,jj,jk) = bor1 * zcl * bor2 436 ! TOTAL F, S, and BORATE CONCENTR. [MOLES/L] 437 borat(ji,jj,jk) = 0.0002414 * zcl / 10.811 438 sulfat(ji,jj,jk) = zst 439 fluorid(ji,jj,jk) = zft 320 440 321 441 ! Iron and SIO3 saturation concentration from ... 322 442 sio3eq(ji,jj,jk) = EXP( LOG( 10.) * ( 6.44 - 968. / ztkel ) ) * 1.e-6 323 fekeq (ji,jj,jk) = 10**( 17.27 - 1565.7 / ( 273.15 + ztc ) ) 324 443 fekeq (ji,jj,jk) = 10**( 17.27 - 1565.7 / ztkel ) 444 445 ! Liu and Millero (1999) only valid 5 - 50 degC 446 ztkel1 = MAX( 5. , tempis(ji,jj,jk) ) + 273.16 447 fesol(ji,jj,jk,1) = 10**((-13.486) - (0.1856* (zis**0.5)) + (0.3073*zis) + (5254.0/ztkel1)) 448 fesol(ji,jj,jk,2) = 10**(2.517 - (0.885*(zis**0.5)) + (0.2139 * zis) - (1320.0/ztkel1) ) 449 fesol(ji,jj,jk,3) = 10**(0.4511 - (0.3305*(zis**0.5)) - (1996.0/ztkel1) ) 450 fesol(ji,jj,jk,4) = 10**(-0.2965 - (0.7881*(zis**0.5)) - (4086.0/ztkel1) ) 451 fesol(ji,jj,jk,5) = 10**(4.4466 - (0.8505*(zis**0.5)) - (7980.0/ztkel1) ) 325 452 END DO 326 453 END DO … … 331 458 END SUBROUTINE p4z_che 332 459 460 SUBROUTINE p4z_che_ahini( p_hini ) 461 !!--------------------------------------------------------------------- 462 !! *** ROUTINE ahini_for_at *** 463 !! 464 !! Subroutine returns the root for the 2nd order approximation of the 465 !! DIC -- B_T -- A_CB equation for [H+] (reformulated as a cubic 466 !! polynomial) around the local minimum, if it exists. 467 !! Returns * 1E-03_wp if p_alkcb <= 0 468 !! * 1E-10_wp if p_alkcb >= 2*p_dictot + p_bortot 469 !! * 1E-07_wp if 0 < p_alkcb < 2*p_dictot + p_bortot 470 !! and the 2nd order approximation does not have 471 !! a solution 472 !!--------------------------------------------------------------------- 473 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(OUT) :: p_hini 474 INTEGER :: ji, jj, jk 475 REAL(wp) :: zca1, zba1 476 REAL(wp) :: zd, zsqrtd, zhmin 477 REAL(wp) :: za2, za1, za0 478 REAL(wp) :: p_dictot, p_bortot, p_alkcb 479 480 IF( nn_timing == 1 ) CALL timing_start('p4z_che_ahini') 481 ! 482 DO jk = 1, jpk 483 DO jj = 1, jpj 484 DO ji = 1, jpi 485 p_alkcb = trb(ji,jj,jk,jptal) * 1000. / (rhop(ji,jj,jk) + rtrn) 486 p_dictot = trb(ji,jj,jk,jpdic) * 1000. / (rhop(ji,jj,jk) + rtrn) 487 p_bortot = borat(ji,jj,jk) 488 IF (p_alkcb <= 0.) THEN 489 p_hini(ji,jj,jk) = 1.e-3 490 ELSEIF (p_alkcb >= (2.*p_dictot + p_bortot)) THEN 491 p_hini(ji,jj,jk) = 1.e-10_wp 492 ELSE 493 zca1 = p_dictot/( p_alkcb + rtrn ) 494 zba1 = p_bortot/ (p_alkcb + rtrn ) 495 ! Coefficients of the cubic polynomial 496 za2 = aKb3(ji,jj,jk)*(1. - zba1) + ak13(ji,jj,jk)*(1.-zca1) 497 za1 = ak13(ji,jj,jk)*akb3(ji,jj,jk)*(1. - zba1 - zca1) & 498 & + ak13(ji,jj,jk)*ak23(ji,jj,jk)*(1. - (zca1+zca1)) 499 za0 = ak13(ji,jj,jk)*ak23(ji,jj,jk)*akb3(ji,jj,jk)*(1. - zba1 - (zca1+zca1)) 500 ! Taylor expansion around the minimum 501 zd = za2*za2 - 3.*za1 ! Discriminant of the quadratic equation 502 ! for the minimum close to the root 503 504 IF(zd > 0.) THEN ! If the discriminant is positive 505 zsqrtd = SQRT(zd) 506 IF(za2 < 0) THEN 507 zhmin = (-za2 + zsqrtd)/3. 508 ELSE 509 zhmin = -za1/(za2 + zsqrtd) 510 ENDIF 511 p_hini(ji,jj,jk) = zhmin + SQRT(-(za0 + zhmin*(za1 + zhmin*(za2 + zhmin)))/zsqrtd) 512 ELSE 513 p_hini(ji,jj,jk) = 1.e-7 514 ENDIF 515 ! 516 ENDIF 517 END DO 518 END DO 519 END DO 520 ! 521 IF( nn_timing == 1 ) CALL timing_stop('p4z_che_ahini') 522 ! 523 END SUBROUTINE p4z_che_ahini 524 525 !=============================================================================== 526 SUBROUTINE anw_infsup( p_alknw_inf, p_alknw_sup ) 527 528 ! Subroutine returns the lower and upper bounds of "non-water-selfionization" 529 ! contributions to total alkalinity (the infimum and the supremum), i.e 530 ! inf(TA - [OH-] + [H+]) and sup(TA - [OH-] + [H+]) 531 532 ! Argument variables 533 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(OUT) :: p_alknw_inf 534 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(OUT) :: p_alknw_sup 535 536 p_alknw_inf(:,:,:) = -trb(:,:,:,jppo4) * 1000. / (rhop(:,:,:) + rtrn) - sulfat(:,:,:) & 537 & - fluorid(:,:,:) 538 p_alknw_sup(:,:,:) = (2. * trb(:,:,:,jpdic) + 2. * trb(:,:,:,jppo4) + trb(:,:,:,jpsil) ) & 539 & * 1000. / (rhop(:,:,:) + rtrn) + borat(:,:,:) 540 541 END SUBROUTINE anw_infsup 542 543 544 SUBROUTINE p4z_che_solve_hi( p_hini, zhi ) 545 546 ! Universal pH solver that converges from any given initial value, 547 ! determines upper an lower bounds for the solution if required 548 549 ! Argument variables 550 !-------------------- 551 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(IN) :: p_hini 552 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(OUT) :: zhi 553 554 ! Local variables 555 !----------------- 556 INTEGER :: ji, jj, jk, jn 557 REAL(wp) :: zh_ini, zh, zh_prev, zh_lnfactor 558 REAL(wp) :: zdelta, zh_delta 559 REAL(wp) :: zeqn, zdeqndh, zalka 560 REAL(wp) :: aphscale 561 REAL(wp) :: znumer_dic, zdnumer_dic, zdenom_dic, zalk_dic, zdalk_dic 562 REAL(wp) :: znumer_bor, zdnumer_bor, zdenom_bor, zalk_bor, zdalk_bor 563 REAL(wp) :: znumer_po4, zdnumer_po4, zdenom_po4, zalk_po4, zdalk_po4 564 REAL(wp) :: znumer_sil, zdnumer_sil, zdenom_sil, zalk_sil, zdalk_sil 565 REAL(wp) :: znumer_so4, zdnumer_so4, zdenom_so4, zalk_so4, zdalk_so4 566 REAL(wp) :: znumer_flu, zdnumer_flu, zdenom_flu, zalk_flu, zdalk_flu 567 REAL(wp) :: zalk_wat, zdalk_wat 568 REAL(wp) :: zfact, p_alktot, zdic, zbot, zpt, zst, zft, zsit 569 LOGICAL :: l_exitnow 570 REAL(wp), PARAMETER :: pz_exp_threshold = 1.0 571 REAL(wp), POINTER, DIMENSION(:,:,:) :: zalknw_inf, zalknw_sup, rmask, zh_min, zh_max, zeqn_absmin 572 573 IF( nn_timing == 1 ) CALL timing_start('p4z_che_solve_hi') 574 ! Allocate temporary workspace 575 CALL wrk_alloc( jpi, jpj, jpk, zalknw_inf, zalknw_sup, rmask ) 576 CALL wrk_alloc( jpi, jpj, jpk, zh_min, zh_max, zeqn_absmin ) 577 578 CALL anw_infsup( zalknw_inf, zalknw_sup ) 579 580 rmask(:,:,:) = tmask(:,:,:) 581 zhi(:,:,:) = 0. 582 583 ! TOTAL H+ scale: conversion factor for Htot = aphscale * Hfree 584 DO jk = 1, jpk 585 DO jj = 1, jpj 586 DO ji = 1, jpi 587 IF (rmask(ji,jj,jk) == 1.) THEN 588 p_alktot = trb(ji,jj,jk,jptal) * 1000. / (rhop(ji,jj,jk) + rtrn) 589 aphscale = 1. + sulfat(ji,jj,jk)/aks3(ji,jj,jk) 590 zh_ini = p_hini(ji,jj,jk) 591 592 zdelta = (p_alktot-zalknw_inf(ji,jj,jk))**2 + 4.*akw3(ji,jj,jk)/aphscale 593 594 IF(p_alktot >= zalknw_inf(ji,jj,jk)) THEN 595 zh_min(ji,jj,jk) = 2.*akw3(ji,jj,jk) /( p_alktot-zalknw_inf(ji,jj,jk) + SQRT(zdelta) ) 596 ELSE 597 zh_min(ji,jj,jk) = aphscale*(-(p_alktot-zalknw_inf(ji,jj,jk)) + SQRT(zdelta) ) / 2. 598 ENDIF 599 600 zdelta = (p_alktot-zalknw_sup(ji,jj,jk))**2 + 4.*akw3(ji,jj,jk)/aphscale 601 602 IF(p_alktot <= zalknw_sup(ji,jj,jk)) THEN 603 zh_max(ji,jj,jk) = aphscale*(-(p_alktot-zalknw_sup(ji,jj,jk)) + SQRT(zdelta) ) / 2. 604 ELSE 605 zh_max(ji,jj,jk) = 2.*akw3(ji,jj,jk) /( p_alktot-zalknw_sup(ji,jj,jk) + SQRT(zdelta) ) 606 ENDIF 607 608 zhi(ji,jj,jk) = MAX(MIN(zh_max(ji,jj,jk), zh_ini), zh_min(ji,jj,jk)) 609 ENDIF 610 END DO 611 END DO 612 END DO 613 614 zeqn_absmin(:,:,:) = HUGE(1._wp) 615 616 DO jn = 1, jp_maxniter_atgen 617 DO jk = 1, jpk 618 DO jj = 1, jpj 619 DO ji = 1, jpi 620 IF (rmask(ji,jj,jk) == 1.) THEN 621 zfact = rhop(ji,jj,jk) / 1000. + rtrn 622 p_alktot = trb(ji,jj,jk,jptal) / zfact 623 zdic = trb(ji,jj,jk,jpdic) / zfact 624 zbot = borat(ji,jj,jk) 625 zpt = trb(ji,jj,jk,jppo4) / zfact * po4r 626 zsit = trb(ji,jj,jk,jpsil) / zfact 627 zst = sulfat (ji,jj,jk) 628 zft = fluorid(ji,jj,jk) 629 aphscale = 1. + sulfat(ji,jj,jk)/aks3(ji,jj,jk) 630 zh = zhi(ji,jj,jk) 631 zh_prev = zh 632 633 ! H2CO3 - HCO3 - CO3 : n=2, m=0 634 znumer_dic = 2.*ak13(ji,jj,jk)*ak23(ji,jj,jk) + zh*ak13(ji,jj,jk) 635 zdenom_dic = ak13(ji,jj,jk)*ak23(ji,jj,jk) + zh*(ak13(ji,jj,jk) + zh) 636 zalk_dic = zdic * (znumer_dic/zdenom_dic) 637 zdnumer_dic = ak13(ji,jj,jk)*ak13(ji,jj,jk)*ak23(ji,jj,jk) + zh & 638 *(4.*ak13(ji,jj,jk)*ak23(ji,jj,jk) + zh*ak13(ji,jj,jk)) 639 zdalk_dic = -zdic*(zdnumer_dic/zdenom_dic**2) 640 641 642 ! B(OH)3 - B(OH)4 : n=1, m=0 643 znumer_bor = akb3(ji,jj,jk) 644 zdenom_bor = akb3(ji,jj,jk) + zh 645 zalk_bor = zbot * (znumer_bor/zdenom_bor) 646 zdnumer_bor = akb3(ji,jj,jk) 647 zdalk_bor = -zbot*(zdnumer_bor/zdenom_bor**2) 648 649 650 ! H3PO4 - H2PO4 - HPO4 - PO4 : n=3, m=1 651 znumer_po4 = 3.*ak1p3(ji,jj,jk)*ak2p3(ji,jj,jk)*ak3p3(ji,jj,jk) & 652 & + zh*(2.*ak1p3(ji,jj,jk)*ak2p3(ji,jj,jk) + zh* ak1p3(ji,jj,jk)) 653 zdenom_po4 = ak1p3(ji,jj,jk)*ak2p3(ji,jj,jk)*ak3p3(ji,jj,jk) & 654 & + zh*( ak1p3(ji,jj,jk)*ak2p3(ji,jj,jk) + zh*(ak1p3(ji,jj,jk) + zh)) 655 zalk_po4 = zpt * (znumer_po4/zdenom_po4 - 1.) ! Zero level of H3PO4 = 1 656 zdnumer_po4 = ak1p3(ji,jj,jk)*ak2p3(ji,jj,jk)*ak1p3(ji,jj,jk)*ak2p3(ji,jj,jk)*ak3p3(ji,jj,jk) & 657 & + zh*(4.*ak1p3(ji,jj,jk)*ak1p3(ji,jj,jk)*ak2p3(ji,jj,jk)*ak3p3(ji,jj,jk) & 658 & + zh*(9.*ak1p3(ji,jj,jk)*ak2p3(ji,jj,jk)*ak3p3(ji,jj,jk) & 659 & + ak1p3(ji,jj,jk)*ak1p3(ji,jj,jk)*ak2p3(ji,jj,jk) & 660 & + zh*(4.*ak1p3(ji,jj,jk)*ak2p3(ji,jj,jk) + zh * ak1p3(ji,jj,jk) ) ) ) 661 zdalk_po4 = -zpt * (zdnumer_po4/zdenom_po4**2) 662 663 ! H4SiO4 - H3SiO4 : n=1, m=0 664 znumer_sil = aksi3(ji,jj,jk) 665 zdenom_sil = aksi3(ji,jj,jk) + zh 666 zalk_sil = zsit * (znumer_sil/zdenom_sil) 667 zdnumer_sil = aksi3(ji,jj,jk) 668 zdalk_sil = -zsit * (zdnumer_sil/zdenom_sil**2) 669 670 ! HSO4 - SO4 : n=1, m=1 671 aphscale = 1.0 + zst/aks3(ji,jj,jk) 672 znumer_so4 = aks3(ji,jj,jk) * aphscale 673 zdenom_so4 = aks3(ji,jj,jk) * aphscale + zh 674 zalk_so4 = zst * (znumer_so4/zdenom_so4 - 1.) 675 zdnumer_so4 = aks3(ji,jj,jk) 676 zdalk_so4 = -zst * (zdnumer_so4/zdenom_so4**2) 677 678 ! HF - F : n=1, m=1 679 znumer_flu = akf3(ji,jj,jk) 680 zdenom_flu = akf3(ji,jj,jk) + zh 681 zalk_flu = zft * (znumer_flu/zdenom_flu - 1.) 682 zdnumer_flu = akf3(ji,jj,jk) 683 zdalk_flu = -zft * (zdnumer_flu/zdenom_flu**2) 684 685 ! H2O - OH 686 aphscale = 1.0 + zst/aks3(ji,jj,jk) 687 zalk_wat = akw3(ji,jj,jk)/zh - zh/aphscale 688 zdalk_wat = -akw3(ji,jj,jk)/zh**2 - 1./aphscale 689 690 ! CALCULATE [ALK]([CO3--], [HCO3-]) 691 zeqn = zalk_dic + zalk_bor + zalk_po4 + zalk_sil & 692 & + zalk_so4 + zalk_flu & 693 & + zalk_wat - p_alktot 694 695 zalka = p_alktot - (zalk_bor + zalk_po4 + zalk_sil & 696 & + zalk_so4 + zalk_flu + zalk_wat) 697 698 zdeqndh = zdalk_dic + zdalk_bor + zdalk_po4 + zdalk_sil & 699 & + zdalk_so4 + zdalk_flu + zdalk_wat 700 701 ! Adapt bracketing interval 702 IF(zeqn > 0._wp) THEN 703 zh_min(ji,jj,jk) = zh_prev 704 ELSEIF(zeqn < 0._wp) THEN 705 zh_max(ji,jj,jk) = zh_prev 706 ENDIF 707 708 IF(ABS(zeqn) >= 0.5_wp*zeqn_absmin(ji,jj,jk)) THEN 709 ! if the function evaluation at the current point is 710 ! not decreasing faster than with a bisection step (at least linearly) 711 ! in absolute value take one bisection step on [ph_min, ph_max] 712 ! ph_new = (ph_min + ph_max)/2d0 713 ! 714 ! In terms of [H]_new: 715 ! [H]_new = 10**(-ph_new) 716 ! = 10**(-(ph_min + ph_max)/2d0) 717 ! = SQRT(10**(-(ph_min + phmax))) 718 ! = SQRT(zh_max * zh_min) 719 zh = SQRT(zh_max(ji,jj,jk) * zh_min(ji,jj,jk)) 720 zh_lnfactor = (zh - zh_prev)/zh_prev ! Required to test convergence below 721 ELSE 722 ! dzeqn/dpH = dzeqn/d[H] * d[H]/dpH 723 ! = -zdeqndh * LOG(10) * [H] 724 ! \Delta pH = -zeqn/(zdeqndh*d[H]/dpH) = zeqn/(zdeqndh*[H]*LOG(10)) 725 ! 726 ! pH_new = pH_old + \deltapH 727 ! 728 ! [H]_new = 10**(-pH_new) 729 ! = 10**(-pH_old - \Delta pH) 730 ! = [H]_old * 10**(-zeqn/(zdeqndh*[H]_old*LOG(10))) 731 ! = [H]_old * EXP(-LOG(10)*zeqn/(zdeqndh*[H]_old*LOG(10))) 732 ! = [H]_old * EXP(-zeqn/(zdeqndh*[H]_old)) 733 734 zh_lnfactor = -zeqn/(zdeqndh*zh_prev) 735 736 IF(ABS(zh_lnfactor) > pz_exp_threshold) THEN 737 zh = zh_prev*EXP(zh_lnfactor) 738 ELSE 739 zh_delta = zh_lnfactor*zh_prev 740 zh = zh_prev + zh_delta 741 ENDIF 742 743 IF( zh < zh_min(ji,jj,jk) ) THEN 744 ! if [H]_new < [H]_min 745 ! i.e., if ph_new > ph_max then 746 ! take one bisection step on [ph_prev, ph_max] 747 ! ph_new = (ph_prev + ph_max)/2d0 748 ! In terms of [H]_new: 749 ! [H]_new = 10**(-ph_new) 750 ! = 10**(-(ph_prev + ph_max)/2d0) 751 ! = SQRT(10**(-(ph_prev + phmax))) 752 ! = SQRT([H]_old*10**(-ph_max)) 753 ! = SQRT([H]_old * zh_min) 754 zh = SQRT(zh_prev * zh_min(ji,jj,jk)) 755 zh_lnfactor = (zh - zh_prev)/zh_prev ! Required to test convergence below 756 ENDIF 757 758 IF( zh > zh_max(ji,jj,jk) ) THEN 759 ! if [H]_new > [H]_max 760 ! i.e., if ph_new < ph_min, then 761 ! take one bisection step on [ph_min, ph_prev] 762 ! ph_new = (ph_prev + ph_min)/2d0 763 ! In terms of [H]_new: 764 ! [H]_new = 10**(-ph_new) 765 ! = 10**(-(ph_prev + ph_min)/2d0) 766 ! = SQRT(10**(-(ph_prev + ph_min))) 767 ! = SQRT([H]_old*10**(-ph_min)) 768 ! = SQRT([H]_old * zhmax) 769 zh = SQRT(zh_prev * zh_max(ji,jj,jk)) 770 zh_lnfactor = (zh - zh_prev)/zh_prev ! Required to test convergence below 771 ENDIF 772 ENDIF 773 774 zeqn_absmin(ji,jj,jk) = MIN( ABS(zeqn), zeqn_absmin(ji,jj,jk)) 775 776 ! Stop iterations once |\delta{[H]}/[H]| < rdel 777 ! <=> |(zh - zh_prev)/zh_prev| = |EXP(-zeqn/(zdeqndh*zh_prev)) -1| < rdel 778 ! |EXP(-zeqn/(zdeqndh*zh_prev)) -1| ~ |zeqn/(zdeqndh*zh_prev)| 779 780 ! Alternatively: 781 ! |\Delta pH| = |zeqn/(zdeqndh*zh_prev*LOG(10))| 782 ! ~ 1/LOG(10) * |\Delta [H]|/[H] 783 ! < 1/LOG(10) * rdel 784 785 ! Hence |zeqn/(zdeqndh*zh)| < rdel 786 787 ! rdel <-- pp_rdel_ah_target 788 l_exitnow = (ABS(zh_lnfactor) < pp_rdel_ah_target) 789 790 IF(l_exitnow) THEN 791 rmask(ji,jj,jk) = 0. 792 ENDIF 793 794 zhi(ji,jj,jk) = zh 795 796 IF(jn >= jp_maxniter_atgen) THEN 797 zhi(ji,jj,jk) = -1._wp 798 ENDIF 799 800 ENDIF 801 END DO 802 END DO 803 END DO 804 END DO 805 ! 806 CALL wrk_dealloc( jpi, jpj, jpk, zalknw_inf, zalknw_sup, rmask ) 807 CALL wrk_dealloc( jpi, jpj, jpk, zh_min, zh_max, zeqn_absmin ) 808 809 810 IF( nn_timing == 1 ) CALL timing_stop('p4z_che_solve_hi') 811 812 813 END SUBROUTINE p4z_che_solve_hi 333 814 334 815 INTEGER FUNCTION p4z_che_alloc() … … 336 817 !! *** ROUTINE p4z_che_alloc *** 337 818 !!---------------------------------------------------------------------- 338 ALLOCATE( sio3eq(jpi,jpj,jpk), fekeq(jpi,jpj,jpk), chemc(jpi,jpj,3), chemo2(jpi,jpj,jpk), & 339 & tempis(jpi,jpj,jpk), STAT=p4z_che_alloc ) 819 INTEGER :: ierr(3) ! Local variables 820 !!---------------------------------------------------------------------- 821 822 ierr(:) = 0 823 824 ALLOCATE( sio3eq(jpi,jpj,jpk), fekeq(jpi,jpj,jpk), chemc(jpi,jpj,3), chemo2(jpi,jpj,jpk), STAT=ierr(1) ) 825 826 ALLOCATE( akb3(jpi,jpj,jpk) , tempis(jpi, jpj, jpk), & 827 & akw3(jpi,jpj,jpk) , borat (jpi,jpj,jpk) , & 828 & aks3(jpi,jpj,jpk) , akf3(jpi,jpj,jpk) , & 829 & ak1p3(jpi,jpj,jpk) , ak2p3(jpi,jpj,jpk) , & 830 & ak3p3(jpi,jpj,jpk) , aksi3(jpi,jpj,jpk) , & 831 & fluorid(jpi,jpj,jpk) , sulfat(jpi,jpj,jpk) , & 832 & salinprac(jpi,jpj,jpk), STAT=ierr(2) ) 833 834 ALLOCATE( fesol(jpi,jpj,jpk,5), STAT=ierr(3) ) 835 836 !* Variable for chemistry of the CO2 cycle 837 p4z_che_alloc = MAXVAL( ierr ) 340 838 ! 341 839 IF( p4z_che_alloc /= 0 ) CALL ctl_warn('p4z_che_alloc : failed to allocate arrays.') … … 355 853 356 854 !!====================================================================== 357 END MODULE p4zche855 END MODULE p4zche -
branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zflx.F90
r6943 r7607 87 87 REAL(wp) :: zfld, zflu, zfld16, zflu16, zfact 88 88 REAL(wp) :: zvapsw, zsal, zfco2, zxc2, xCO2approx, ztkel, zfugcoeff 89 REAL(wp) :: zph, z ah2, zbot, zdic, zalk, zsch_o2, zalka, zsch_co289 REAL(wp) :: zph, zdic, zsch_o2, zsch_co2 90 90 REAL(wp) :: zyr_dec, zdco2dt 91 91 CHARACTER (len=25) :: charout … … 122 122 #endif 123 123 124 DO jm = 1, 10 125 !CDIR NOVERRCHK 126 DO jj = 1, jpj 127 !CDIR NOVERRCHK 128 DO ji = 1, jpi 129 130 ! DUMMY VARIABLES FOR DIC, H+, AND BORATE 131 zbot = borat(ji,jj,1) 132 zfact = rhop(ji,jj,1) / 1000. + rtrn 133 zdic = trb(ji,jj,1,jpdic) / zfact 134 zph = MAX( hi(ji,jj,1), 1.e-10 ) / zfact 135 zalka = trb(ji,jj,1,jptal) / zfact 136 137 ! CALCULATE [ALK]([CO3--], [HCO3-]) 138 zalk = zalka - ( akw3(ji,jj,1) / zph - zph / aphscale(ji,jj,1) & 139 & + zbot / ( 1.+ zph / akb3(ji,jj,1) ) ) 140 141 ! CALCULATE [H+] AND [H2CO3] 142 zah2 = SQRT( (zdic-zalk)**2 + 4.* ( zalk * ak23(ji,jj,1) & 143 & / ak13(ji,jj,1) ) * ( 2.* zdic - zalk ) ) 144 zah2 = 0.5 * ak13(ji,jj,1) / zalk * ( ( zdic - zalk ) + zah2 ) 145 zh2co3(ji,jj) = ( 2.* zdic - zalk ) / ( 2.+ ak13(ji,jj,1) / zah2 ) * zfact 146 hi(ji,jj,1) = zah2 * zfact 147 END DO 124 DO jj = 1, jpj 125 DO ji = 1, jpi 126 ! DUMMY VARIABLES FOR DIC, H+, AND BORATE 127 zfact = rhop(ji,jj,1) / 1000. + rtrn 128 zdic = trb(ji,jj,1,jpdic) 129 zph = MAX( hi(ji,jj,1), 1.e-10 ) / zfact 130 ! CALCULATE [H2CO3] 131 zh2co3(ji,jj) = zdic/(1. + ak13(ji,jj,1)/zph + ak13(ji,jj,1)*ak23(ji,jj,1)/zph**2) 148 132 END DO 149 133 END DO 150 151 134 152 135 ! -------------- … … 254 237 ENDIF 255 238 ! 239 #if defined key_cpl_carbon_cycle 240 ! change units for carbon cycle coupling 241 oce_co2(:,:) = oce_co2(:,:) / e1e2t(:,:) * rfact2r ! in molC/m2/s 242 #endif 243 ! 256 244 CALL wrk_dealloc( jpi, jpj, zkgco2, zkgo2, zh2co3, zoflx, zpco2atm ) 257 245 ! -
branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zlys.F90
r6943 r7607 23 23 USE sms_pisces ! PISCES Source Minus Sink variables 24 24 USE prtctl_trc ! print control for debugging 25 USE p4zche ! Chemical model 25 26 USE iom ! I/O manager 26 27 … … 60 61 ! 61 62 INTEGER, INTENT(in) :: kt, knt ! ocean time step 62 INTEGER :: ji, jj, jk, jn 63 REAL(wp) :: zalk, zdic, zph, zah2 64 REAL(wp) :: zdispot, zfact, zcalcon, zalka, zaldi 63 INTEGER :: ji, jj, jk 64 REAL(wp) :: zdispot, zfact, zcalcon 65 65 REAL(wp) :: zomegaca, zexcess, zexcess0 66 66 CHARACTER (len=25) :: charout 67 REAL(wp), POINTER, DIMENSION(:,:,:) :: zco3, zco3sat, zcaldiss 67 REAL(wp), POINTER, DIMENSION(:,:,:) :: zco3, zco3sat, zcaldiss, zhinit, zhi 68 68 !!--------------------------------------------------------------------- 69 69 ! 70 70 IF( nn_timing == 1 ) CALL timing_start('p4z_lys') 71 71 ! 72 CALL wrk_alloc( jpi, jpj, jpk, zco3, zco3sat, zcaldiss )72 CALL wrk_alloc( jpi, jpj, jpk, zco3, zco3sat, zcaldiss, zhinit, zhi ) 73 73 ! 74 74 zco3 (:,:,:) = 0. 75 75 zcaldiss(:,:,:) = 0. 76 zhinit(:,:,:) = hi(:,:,:) * 1000. / ( rhop(:,:,:) + rtrn ) 76 77 ! ------------------------------------------- 77 78 ! COMPUTE [CO3--] and [H+] CONCENTRATIONS 78 79 ! ------------------------------------------- 79 80 DO jn = 1, 5 ! BEGIN OF ITERATION 81 ! 82 !CDIR NOVERRCHK 83 DO jk = 1, jpkm1 84 !CDIR NOVERRCHK 85 DO jj = 1, jpj 86 !CDIR NOVERRCHK 87 DO ji = 1, jpi 88 zfact = rhop(ji,jj,jk) / 1000. + rtrn 89 zph = hi(ji,jj,jk) * tmask(ji,jj,jk) / zfact + ( 1.-tmask(ji,jj,jk) ) * 1.e-9 ! [H+] 90 zdic = trb(ji,jj,jk,jpdic) / zfact 91 zalka = trb(ji,jj,jk,jptal) / zfact 92 ! CALCULATE [ALK]([CO3--], [HCO3-]) 93 zalk = zalka - ( akw3(ji,jj,jk) / zph - zph / ( aphscale(ji,jj,jk) + rtrn ) & 94 & + borat(ji,jj,jk) / ( 1. + zph / akb3(ji,jj,jk) ) ) 95 ! CALCULATE [H+] and [CO3--] 96 zaldi = zdic - zalk 97 zah2 = SQRT( zaldi * zaldi + 4.* ( zalk * ak23(ji,jj,jk) / ak13(ji,jj,jk) ) * ( zdic + zaldi ) ) 98 zah2 = 0.5 * ak13(ji,jj,jk) / zalk * ( zaldi + zah2 ) 99 ! 100 zco3(ji,jj,jk) = zalk / ( 2. + zah2 / ak23(ji,jj,jk) ) * zfact 101 hi(ji,jj,jk) = zah2 * zfact 102 END DO 80 81 CALL p4z_che_solve_hi( zhinit, zhi ) 82 83 DO jk = 1, jpkm1 84 DO jj = 1, jpj 85 DO ji = 1, jpi 86 zco3(ji,jj,jk) = trb(ji,jj,jk,jpdic) * ak13(ji,jj,jk) * ak23(ji,jj,jk) / (zhi(ji,jj,jk)**2 & 87 & + ak13(ji,jj,jk) * zhi(ji,jj,jk) + ak13(ji,jj,jk) * ak23(ji,jj,jk) + rtrn ) 88 hi(ji,jj,jk) = zhi(ji,jj,jk) * rhop(ji,jj,jk) / 1000. 103 89 END DO 104 90 END DO 105 ! 106 END DO 91 END DO 107 92 108 93 ! --------------------------------------------------------- … … 138 123 ! AND [SUM(CO2)] DUE TO CACO3 DISSOLUTION/PRECIPITATION 139 124 zcaldiss(ji,jj,jk) = zdispot * rfact2 / rmtss ! calcite dissolution 140 zco3(ji,jj,jk) = zco3(ji,jj,jk) + zcaldiss(ji,jj,jk)141 125 ! 142 126 tra(ji,jj,jk,jptal) = tra(ji,jj,jk,jptal) + 2. * zcaldiss(ji,jj,jk) … … 167 151 ENDIF 168 152 ! 169 CALL wrk_dealloc( jpi, jpj, jpk, zco3, zc o3sat, zcaldiss)153 CALL wrk_dealloc( jpi, jpj, jpk, zco3, zcaldiss, zhinit, zhi, zco3sat ) 170 154 ! 171 155 IF( nn_timing == 1 ) CALL timing_stop('p4z_lys') … … 186 170 !! 187 171 !!---------------------------------------------------------------------- 188 INTEGER :: ji, jj, jk189 172 INTEGER :: ios ! Local integer output status for namelist read 190 REAL(wp) :: zcaralk, zbicarb, zco3191 REAL(wp) :: ztmas, ztmas1192 193 173 NAMELIST/nampiscal/ kdca, nca 194 174 !!---------------------------------------------------------------------- -
branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsms.F90
r6420 r7607 85 85 CALL p4z_che ! initialize the chemical constants 86 86 ! 87 IF( .NOT. ln_rsttr ) THEN ; CALL p4z_ ph_ini ! set PH at kt=nit00087 IF( .NOT. ln_rsttr ) THEN ; CALL p4z_che_ahini( hi ) ! set PH at kt=nit000 88 88 ELSE ; CALL p4z_rst( nittrc000, 'READ' ) !* read or initialize all required fields 89 89 ENDIF … … 310 310 END SUBROUTINE p4z_sms_init 311 311 312 SUBROUTINE p4z_ph_ini313 !!---------------------------------------------------------------------314 !! *** ROUTINE p4z_ini_ph ***315 !!316 !! ** Purpose : Initialization of chemical variables of the carbon cycle317 !!---------------------------------------------------------------------318 INTEGER :: ji, jj, jk319 REAL(wp) :: zcaralk, zbicarb, zco3320 REAL(wp) :: ztmas, ztmas1321 !!---------------------------------------------------------------------322 323 ! Set PH from total alkalinity, borat (???), akb3 (???) and ak23 (???)324 ! --------------------------------------------------------325 DO jk = 1, jpk326 DO jj = 1, jpj327 DO ji = 1, jpi328 ztmas = tmask(ji,jj,jk)329 ztmas1 = 1. - tmask(ji,jj,jk)330 zcaralk = trb(ji,jj,jk,jptal) - borat(ji,jj,jk) / ( 1. + 1.E-8 / ( rtrn + akb3(ji,jj,jk) ) )331 zco3 = ( zcaralk - trb(ji,jj,jk,jpdic) ) * ztmas + 0.5e-3 * ztmas1332 zbicarb = ( 2. * trb(ji,jj,jk,jpdic) - zcaralk )333 hi(ji,jj,jk) = ( ak23(ji,jj,jk) * zbicarb / zco3 ) * ztmas + 1.e-9 * ztmas1334 END DO335 END DO336 END DO337 !338 END SUBROUTINE p4z_ph_ini339 340 312 SUBROUTINE p4z_rst( kt, cdrw ) 341 313 !!--------------------------------------------------------------------- … … 350 322 INTEGER , INTENT(in) :: kt ! ocean time-step 351 323 CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag 352 !353 INTEGER :: ji, jj, jk354 REAL(wp) :: zcaralk, zbicarb, zco3355 REAL(wp) :: ztmas, ztmas1356 324 !!--------------------------------------------------------------------- 357 325 … … 365 333 CALL iom_get( numrtr, jpdom_autoglo, 'PH' , hi(:,:,:) ) 366 334 ELSE 367 ! hi(:,:,:) = 1.e-9 368 CALL p4z_ph_ini 335 CALL p4z_che_ahini( hi ) 369 336 ENDIF 370 337 CALL iom_get( numrtr, jpdom_autoglo, 'Silicalim', xksi(:,:) ) -
branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/TOP_SRC/PISCES/sms_pisces.F90
r6287 r7607 93 93 94 94 !!* Variable for chemistry of the CO2 cycle 95 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: akb3 !: ???96 95 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: ak13 !: ??? 97 96 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: ak23 !: ??? 98 97 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: aksp !: ??? 99 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: akw3 !: ???100 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: borat !: ???101 98 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: hi !: ??? 102 99 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: excess !: ??? … … 153 150 154 151 !* Variable for chemistry of the CO2 cycle 155 ALLOCATE( ak b3(jpi,jpj,jpk) , ak13 (jpi,jpj,jpk) ,&152 ALLOCATE( ak13 (jpi,jpj,jpk) , & 156 153 & ak23(jpi,jpj,jpk) , aksp (jpi,jpj,jpk) , & 157 & akw3(jpi,jpj,jpk) , borat (jpi,jpj,jpk) , &158 154 & hi (jpi,jpj,jpk) , excess(jpi,jpj,jpk) , & 159 155 & aphscale(jpi,jpj,jpk), STAT=ierr(4) ) -
branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/TOP_SRC/trcstp.F90
r7522 r7607 131 131 INTEGER, INTENT(in) :: kt 132 132 INTEGER :: jn 133 REAL(wp) :: zkt 133 REAL(wp) :: zkt, zrec 134 134 CHARACTER(len=1) :: cl1 ! 1 character 135 135 CHARACTER(len=2) :: cl2 ! 2 characters … … 153 153 ! 154 154 ! !* Restart: read in restart file 155 IF( ln_rsttr .AND. iom_varid( numrtr, 'qsr_mean' , ldstop = .FALSE. ) > 0 .AND. & 156 iom_varid( numrtr, 'qsr_arr_1', ldstop = .FALSE. ) > 0 .AND. & 157 iom_varid( numrtr, 'ktdcy' , ldstop = .FALSE. ) > 0 ) THEN 158 CALL iom_get( numrtr, 'ktdcy', zkt ) ! A mean of qsr 155 IF( ln_rsttr .AND. iom_varid( numrtr, 'qsr_mean' , ldstop = .FALSE. ) > 0 & 156 & .AND. iom_varid( numrtr, 'qsr_arr_1', ldstop = .FALSE. ) > 0 & 157 & .AND. iom_varid( numrtr, 'ktdcy' , ldstop = .FALSE. ) > 0 & 158 & .AND. iom_varid( numrtr, 'nrdcy' , ldstop = .FALSE. ) > 0 ) THEN 159 160 CALL iom_get( numrtr, 'ktdcy', zkt ) 159 161 rsecfst = INT( zkt ) * rdttrc(1) 160 162 IF(lwp) WRITE(numout,*) 'trc_qsr_mean: qsr_mean read in the restart file at time-step rsecfst =', rsecfst, ' s ' 161 163 CALL iom_get( numrtr, jpdom_autoglo, 'qsr_mean', qsr_mean ) ! A mean of qsr 162 DO jn = 1, nb_rec_per_day 163 IF( jn <= 9 ) THEN 164 WRITE(cl1,'(i1)') jn 165 CALL iom_get( numrtr, jpdom_autoglo, 'qsr_arr_'//cl1, qsr_arr(:,:,jn) ) ! A mean of qsr 166 ELSE 167 WRITE(cl2,'(i2.2)') jn 168 CALL iom_get( numrtr, jpdom_autoglo, 'qsr_arr_'//cl2, qsr_arr(:,:,jn) ) ! A mean of qsr 169 ENDIF 170 ENDDO 164 CALL iom_get( numrtr, 'nrdcy', zrec ) ! Number of record per days 165 IF( INT( zrec ) == nb_rec_per_day ) THEN 166 DO jn = 1, nb_rec_per_day 167 IF( jn <= 9 ) THEN 168 WRITE(cl1,'(i1)') jn 169 CALL iom_get( numrtr, jpdom_autoglo, 'qsr_arr_'//cl1, qsr_arr(:,:,jn) ) ! A mean of qsr 170 ELSE 171 WRITE(cl2,'(i2.2)') jn 172 CALL iom_get( numrtr, jpdom_autoglo, 'qsr_arr_'//cl2, qsr_arr(:,:,jn) ) ! A mean of qsr 173 ENDIF 174 ENDDO 175 ELSE 176 DO jn = 1, nb_rec_per_day 177 qsr_arr(:,:,jn) = qsr_mean(:,:) 178 ENDDO 179 ENDIF 171 180 ELSE !* no restart: set from nit000 values 172 181 IF(lwp) WRITE(numout,*) 'trc_qsr_mean: qsr_mean set to nit000 values' … … 185 194 llnew = ( rseclast - rsecfst ) .ge. rdt_sampl ! new shortwave to store 186 195 IF( llnew ) THEN 187 IF( lwp ) WRITE(numout,*) ' New shortwave to sample for TOP at time kt = ', kt, &196 IF( lwp .AND. kt < nittrc000 + 100 ) WRITE(numout,*) ' New shortwave to sample for TOP at time kt = ', kt, & 188 197 & ' time = ', rseclast/3600.,'hours ' 189 198 rsecfst = rseclast … … 199 208 IF(lwp) WRITE(numout,*) 'trc_mean_qsr : write qsr_mean in restart file kt =', kt 200 209 IF(lwp) WRITE(numout,*) '~~~~~~~' 201 zkt = REAL( kt, wp ) 202 CALL iom_rstput( kt, nitrst, numrtw, 'ktdcy', zkt ) 210 zkt = REAL( kt, wp ) 211 zrec = REAL( nb_rec_per_day, wp ) 212 CALL iom_rstput( kt, nitrst, numrtw, 'ktdcy', zkt ) 213 CALL iom_rstput( kt, nitrst, numrtw, 'nrdcy', zrec ) 203 214 DO jn = 1, nb_rec_per_day 204 215 IF( jn <= 9 ) THEN
Note: See TracChangeset
for help on using the changeset viewer.
