Changeset 6453 for branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM
- Timestamp:
- 2016-04-08T10:10:11+02:00 (8 years ago)
- Location:
- branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM
- Files:
-
- 28 added
- 24 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zbio.F90
r6204 r6453 24 24 USE p4zmicro ! Sources and sinks of microzooplankton 25 25 USE p4zmeso ! Sources and sinks of mesozooplankton 26 USE p4zrem ! Remineralisation of organic matter 27 USE p4zfechem 26 USE p4zrem ! Remineralisation of organic/inorganic matter 27 USE p4zpoc ! Remineralization of organic particles 28 USE p4zagg ! Aggregation of particles 29 USE p4zlys ! Dissolution of calcite 30 USE p4zfechem ! Iron chemistry 31 USE p4zligand 28 32 USE prtctl_trc ! print control for debugging 29 33 USE iom ! I/O manager … … 76 80 77 81 78 CALL p4z_opt ( kt, knt ) ! Optic: PAR in the water column 79 CALL p4z_sink ( kt, knt ) ! vertical flux of particulate organic matter 80 CALL p4z_fechem(kt, knt ) ! Iron chemistry/scavenging 81 CALL p4z_lim ( kt, knt ) ! co-limitations by the various nutrients 82 CALL p4z_prod ( kt, knt ) ! phytoplankton growth rate over the global ocean. 83 ! ! (for each element : C, Si, Fe, Chl ) 84 CALL p4z_mort ( kt ) ! phytoplankton mortality 85 ! ! zooplankton sources/sinks routines 86 CALL p4z_micro( kt, knt ) ! microzooplankton 87 CALL p4z_meso ( kt, knt ) ! mesozooplankton 88 CALL p4z_rem ( kt, knt ) ! remineralization terms of organic matter+scavenging of Fe 82 CALL p4z_opt ( kt, knt ) ! Optic: PAR in the water column 83 CALL p4z_sink ( kt, knt ) ! vertical flux of particulate organic matter 84 CALL p4z_lys (kt, knt ) ! Dissolution of calcite 85 CALL p4z_fechem(kt, knt ) ! Iron chemistry/scavenging 86 CALL p4z_lim ( kt, knt ) ! co-limitations by the various nutrients 87 CALL p4z_prod ( kt, knt ) ! phytoplankton growth rate over the global ocean. 88 ! ! (for each element : C, Si, Fe, Chl ) 89 CALL p4z_mort ( kt ) ! phytoplankton mortality 90 ! ! zooplankton sources/sinks routines 91 CALL p4z_micro ( kt, knt ) ! microzooplankton 92 CALL p4z_meso ( kt, knt ) ! mesozooplankton 93 CALL p4z_agg ( kt, knt ) ! Aggregation of particles 94 CALL p4z_rem ( kt, knt ) ! remineralization terms of organic matter+scavenging of Fe 95 CALL p4z_poc ( kt, knt ) ! Remineralization of organic particles 96 CALL p4z_ligand(kt, knt ) 89 97 ! ! test if tracers concentrations fall below 0. 90 98 ! ! -
branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zche.F90
r6204 r6453 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 #if defined key_pisces 15 #if defined key_pisces || defined key_pisces_quota 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 … … 26 27 PRIVATE 27 28 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 29 PUBLIC p4z_che ! 30 PUBLIC p4z_che_alloc ! 31 PUBLIC ahini_for_at ! 32 PUBLIC solve_at_general ! 33 34 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: sio3eq ! chemistry of Si 35 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: fekeq ! chemistry of Fe 36 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: chemc ! Solubilities of O2 and CO2 37 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: chemo2 ! Solubilities of O2 and CO2 38 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: fesol ! solubility of Fe 39 40 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: akb3 !: ??? 41 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: akw3 !: ??? 42 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: akf3 !: ??? 43 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: aks3 !: ??? 44 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: ak1p3 !: ??? 45 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: ak2p3 !: ??? 46 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: ak3p3 !: ??? 47 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: aksi3 !: ??? 48 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: borat !: ??? 49 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: fluorid !: ??? 50 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: sulfat !: ??? 51 52 !!* Variable for chemistry of the CO2 cycle 35 53 36 54 REAL(wp), PUBLIC :: atcox = 0.20946 ! units atm 37 55 38 REAL(wp) :: salchl = 1. / 1.80655 ! conversion factor for salinity --> chlorinity (Wooster et al. 1969)39 56 REAL(wp) :: o2atm = 1. / ( 1000. * 0.20946 ) 40 57 41 REAL(wp) :: akcc1 = -171.9065 ! coeff. for apparent solubility equilibrium 42 REAL(wp) :: akcc2 = -0.077993 ! Millero et al. 1995 from Mucci 1983 43 REAL(wp) :: akcc3 = 2839.319 44 REAL(wp) :: akcc4 = 71.595 45 REAL(wp) :: akcc5 = -0.77712 46 REAL(wp) :: akcc6 = 0.00284263 47 REAL(wp) :: akcc7 = 178.34 48 REAL(wp) :: akcc8 = -0.07711 49 REAL(wp) :: akcc9 = 0.0041249 50 51 REAL(wp) :: rgas = 83.143 ! universal gas constants 52 REAL(wp) :: oxyco = 1. / 22.4144 ! converts from liters of an ideal gas to moles 53 54 REAL(wp) :: bor1 = 0.00023 ! borat constants 55 REAL(wp) :: bor2 = 1. / 10.82 56 57 REAL(wp) :: ca0 = -162.8301 ! WEISS & PRICE 1980, units mol/(kg atm) 58 REAL(wp) :: ca1 = 218.2968 59 REAL(wp) :: ca2 = 90.9241 60 REAL(wp) :: ca3 = -1.47696 61 REAL(wp) :: ca4 = 0.025695 62 REAL(wp) :: ca5 = -0.025225 63 REAL(wp) :: ca6 = 0.0049867 64 65 REAL(wp) :: c10 = -3670.7 ! Coeff. for 1. dissoc. of carbonic acid (Edmond and Gieskes, 1970) 66 REAL(wp) :: c11 = 62.008 67 REAL(wp) :: c12 = -9.7944 68 REAL(wp) :: c13 = 0.0118 69 REAL(wp) :: c14 = -0.000116 70 71 REAL(wp) :: c20 = -1394.7 ! coeff. for 2. dissoc. of carbonic acid (Millero, 1995) 72 REAL(wp) :: c21 = -4.777 73 REAL(wp) :: c22 = 0.0184 74 REAL(wp) :: c23 = -0.000118 75 76 REAL(wp) :: st1 = 0.14 ! constants for calculate concentrations for sulfate 77 REAL(wp) :: st2 = 1./96.062 ! (Morris & Riley 1966) 78 REAL(wp) :: ks0 = 141.328 79 REAL(wp) :: ks1 = -4276.1 80 REAL(wp) :: ks2 = -23.093 81 REAL(wp) :: ks3 = -13856. 82 REAL(wp) :: ks4 = 324.57 83 REAL(wp) :: ks5 = -47.986 84 REAL(wp) :: ks6 = 35474. 85 REAL(wp) :: ks7 = -771.54 86 REAL(wp) :: ks8 = 114.723 87 REAL(wp) :: ks9 = -2698. 88 REAL(wp) :: ks10 = 1776. 89 REAL(wp) :: ks11 = 1. 90 REAL(wp) :: ks12 = -0.001005 91 92 REAL(wp) :: ft1 = 0.000067 ! constants for calculate concentrations for fluorides 93 REAL(wp) :: ft2 = 1./18.9984 ! (Dickson & Riley 1979 ) 94 REAL(wp) :: kf0 = -12.641 95 REAL(wp) :: kf1 = 1590.2 96 REAL(wp) :: kf2 = 1.525 97 REAL(wp) :: kf3 = 1.0 98 REAL(wp) :: kf4 = -0.001005 99 100 REAL(wp) :: cb0 = -8966.90 ! Coeff. for 1. dissoc. of boric acid 101 REAL(wp) :: cb1 = -2890.53 ! (Dickson and Goyet, 1994) 102 REAL(wp) :: cb2 = -77.942 103 REAL(wp) :: cb3 = 1.728 104 REAL(wp) :: cb4 = -0.0996 105 REAL(wp) :: cb5 = 148.0248 106 REAL(wp) :: cb6 = 137.1942 107 REAL(wp) :: cb7 = 1.62142 108 REAL(wp) :: cb8 = -24.4344 109 REAL(wp) :: cb9 = -25.085 110 REAL(wp) :: cb10 = -0.2474 111 REAL(wp) :: cb11 = 0.053105 112 113 REAL(wp) :: cw0 = -13847.26 ! Coeff. for dissoc. of water (Dickson and Riley, 1979 ) 114 REAL(wp) :: cw1 = 148.9652 115 REAL(wp) :: cw2 = -23.6521 116 REAL(wp) :: cw3 = 118.67 117 REAL(wp) :: cw4 = -5.977 118 REAL(wp) :: cw5 = 1.0495 119 REAL(wp) :: cw6 = -0.01615 120 121 ! ! volumetric solubility constants for o2 in ml/L 122 REAL(wp) :: ox0 = 2.00856 ! from Table 1 for Eq 8 of Garcia and Gordon, 1992. 123 REAL(wp) :: ox1 = 3.22400 ! corrects for moisture and fugacity, but not total atmospheric pressure 124 REAL(wp) :: ox2 = 3.99063 ! Original PISCES code noted this was a solubility, but 125 REAL(wp) :: ox3 = 4.80299 ! was in fact a bunsen coefficient with units L-O2/(Lsw atm-O2) 126 REAL(wp) :: ox4 = 9.78188e-1 ! Hence, need to divide EXP( zoxy ) by 1000, ml-O2 => L-O2 127 REAL(wp) :: ox5 = 1.71069 ! and atcox = 0.20946 to add the 1/atm dimension. 128 REAL(wp) :: ox6 = -6.24097e-3 129 REAL(wp) :: ox7 = -6.93498e-3 130 REAL(wp) :: ox8 = -6.90358e-3 131 REAL(wp) :: ox9 = -4.29155e-3 132 REAL(wp) :: ox10 = -3.11680e-7 58 REAL(wp) :: rgas = 83.14472 ! universal gas constants 59 REAL(wp) :: oxyco = 1. / 22.4144 ! converts from liters of an ideal gas to moles 133 60 134 61 ! ! coeff. for seawater pressure correction : millero 95 135 62 ! ! AGRIF doesn't like the DATA instruction 136 REAL(wp) :: devk11 = -25.5 137 REAL(wp) :: devk12 = -15.82 138 REAL(wp) :: devk13 = -29.48 139 REAL(wp) :: devk14 = -25.60 140 REAL(wp) :: devk15 = -48.76 63 REAL(wp) :: devk10 = -25.5 64 REAL(wp) :: devk11 = -15.82 65 REAL(wp) :: devk12 = -29.48 66 REAL(wp) :: devk13 = -20.02 67 REAL(wp) :: devk14 = -18.03 68 REAL(wp) :: devk15 = -9.78 69 REAL(wp) :: devk16 = -48.76 70 REAL(wp) :: devk17 = -14.51 71 REAL(wp) :: devk18 = -23.12 72 REAL(wp) :: devk19 = -26.57 73 REAL(wp) :: devk110 = -29.48 141 74 ! 142 REAL(wp) :: devk21 = 0.1271 143 REAL(wp) :: devk22 = -0.0219 144 REAL(wp) :: devk23 = 0.1622 145 REAL(wp) :: devk24 = 0.2324 146 REAL(wp) :: devk25 = 0.5304 75 REAL(wp) :: devk20 = 0.1271 76 REAL(wp) :: devk21 = -0.0219 77 REAL(wp) :: devk22 = 0.1622 78 REAL(wp) :: devk23 = 0.1119 79 REAL(wp) :: devk24 = 0.0466 80 REAL(wp) :: devk25 = -0.0090 81 REAL(wp) :: devk26 = 0.5304 82 REAL(wp) :: devk27 = 0.1211 83 REAL(wp) :: devk28 = 0.1758 84 REAL(wp) :: devk29 = 0.2020 85 REAL(wp) :: devk210 = 0.1622 147 86 ! 87 REAL(wp) :: devk30 = 0. 148 88 REAL(wp) :: devk31 = 0. 149 REAL(wp) :: devk32 = 0. 150 REAL(wp) :: devk33 = 2.608E-3 151 REAL(wp) :: devk34 = -3.6246E-3 152 REAL(wp) :: devk35 = 0. 89 REAL(wp) :: devk32 = 2.608E-3 90 REAL(wp) :: devk33 = -1.409e-3 91 REAL(wp) :: devk34 = 0.316e-3 92 REAL(wp) :: devk35 = -0.942e-3 93 REAL(wp) :: devk36 = 0. 94 REAL(wp) :: devk37 = -0.321e-3 95 REAL(wp) :: devk38 = -2.647e-3 96 REAL(wp) :: devk39 = -3.042e-3 97 REAL(wp) :: devk310 = -2.6080e-3 153 98 ! 154 REAL(wp) :: devk41 = -3.08E-3 155 REAL(wp) :: devk42 = 1.13E-3 156 REAL(wp) :: devk43 = -2.84E-3 157 REAL(wp) :: devk44 = -5.13E-3 158 REAL(wp) :: devk45 = -11.76E-3 99 REAL(wp) :: devk40 = -3.08E-3 100 REAL(wp) :: devk41 = 1.13E-3 101 REAL(wp) :: devk42 = -2.84E-3 102 REAL(wp) :: devk43 = -5.13E-3 103 REAL(wp) :: devk44 = -4.53e-3 104 REAL(wp) :: devk45 = -3.91e-3 105 REAL(wp) :: devk46 = -11.76e-3 106 REAL(wp) :: devk47 = -2.67e-3 107 REAL(wp) :: devk48 = -5.15e-3 108 REAL(wp) :: devk49 = -4.08e-3 109 REAL(wp) :: devk410 = -2.84e-3 159 110 ! 160 REAL(wp) :: devk51 = 0.0877E-3 161 REAL(wp) :: devk52 = -0.1475E-3 162 REAL(wp) :: devk53 = 0. 163 REAL(wp) :: devk54 = 0.0794E-3 164 REAL(wp) :: devk55 = 0.3692E-3 111 REAL(wp) :: devk50 = 0.0877E-3 112 REAL(wp) :: devk51 = -0.1475E-3 113 REAL(wp) :: devk52 = 0. 114 REAL(wp) :: devk53 = 0.0794E-3 115 REAL(wp) :: devk54 = 0.09e-3 116 REAL(wp) :: devk55 = 0.054e-3 117 REAL(wp) :: devk56 = 0.3692E-3 118 REAL(wp) :: devk57 = 0.0427e-3 119 REAL(wp) :: devk58 = 0.09e-3 120 REAL(wp) :: devk59 = 0.0714e-3 121 REAL(wp) :: devk510 = 0.0 122 ! 123 124 125 ! General parameters 126 REAL(wp), PARAMETER :: pp_rdel_ah_target = 1.E-4_wp 127 REAL(wp), PARAMETER :: pp_ln10 = 2.302585092994045684018_wp 128 129 ! Maximum number of iterations for each method 130 INTEGER, PARAMETER :: jp_maxniter_atgen = 20 131 132 ! Bookkeeping variables for each method 133 ! - SOLVE_AT_GENERAL 134 INTEGER :: niter_atgen = jp_maxniter_atgen 135 136 165 137 166 138 !!* Substitution … … 182 154 !!--------------------------------------------------------------------- 183 155 INTEGER :: ji, jj, jk 184 REAL(wp) :: ztkel, zt 156 REAL(wp) :: ztkel, ztkel1, zt , zt2 , zsal , zsal2 , zbuf1 , zbuf2 185 157 REAL(wp) :: ztgg , ztgg2, ztgg3 , ztgg4 , ztgg5 186 158 REAL(wp) :: zpres, ztc , zcl , zcpexp, zoxy , zcpexp2 187 REAL(wp) :: zsqrt, ztr , zlogt , zcek1 159 REAL(wp) :: zsqrt, ztr , zlogt , zcek1, zc1, zplat 188 160 REAL(wp) :: zis , zis2 , zsal15, zisqrt 189 161 REAL(wp) :: zckb , zck1 , zck2 , zckw , zak1 , zak2 , zakb , zaksp0, zakw 162 REAL(wp) :: zck1p, zck2p, zck3p, zcksi, zak1p, zak2p, zak3p, zaksi 190 163 REAL(wp) :: zst , zft , zcks , zckf , zaksp1 164 REAL(wp) :: total2free, free2SWS, total2SWS, SWS2total 165 191 166 !!--------------------------------------------------------------------- 192 167 ! … … 200 175 DO ji = 1, jpi 201 176 ! ! SET ABSOLUTE TEMPERATURE 202 ztkel = tsn(ji,jj,1,jp_tem) + 273.1 6177 ztkel = tsn(ji,jj,1,jp_tem) + 273.15 203 178 zt = ztkel * 0.01 204 zt2 = zt * zt205 179 zsal = tsn(ji,jj,1,jp_sal) + ( 1.- tmask(ji,jj,1) ) * 35. 206 zsal2 = zsal * zsal207 zlogt = LOG( zt )208 180 ! ! LN(K0) OF SOLUBILITY OF CO2 (EQ. 12, WEISS, 1980) 209 181 ! ! AND FOR THE ATMOSPHERE FOR NON IDEAL GAS 210 zcek1 = ca0 + ca1 / zt + ca2 * zlogt + ca3 * zt2 + zsal * ( ca4 + ca5 * zt + ca6 * zt2 ) 211 ! ! LN(K0) OF SOLUBILITY OF O2 and N2 in ml/L (EQ. 8, GARCIA AND GORDON, 1992) 212 ztgg = LOG( ( 298.15 - tsn(ji,jj,1,jp_tem) ) / ztkel ) ! Set the GORDON & GARCIA scaled temperature 213 ztgg2 = ztgg * ztgg 214 ztgg3 = ztgg2 * ztgg 215 ztgg4 = ztgg3 * ztgg 216 ztgg5 = ztgg4 * ztgg 217 zoxy = ox0 + ox1 * ztgg + ox2 * ztgg2 + ox3 * ztgg3 + ox4 * ztgg4 + ox5 * ztgg5 & 218 + zsal * ( ox6 + ox7 * ztgg + ox8 * ztgg2 + ox9 * ztgg3 ) + ox10 * zsal2 219 182 zcek1 = 9345.17/ztkel - 60.2409 + 23.3585 * LOG(zt) + zsal*(0.023517 - 0.00023656*ztkel & 183 & + 0.0047036e-4*ztkel**2) 220 184 ! ! SET SOLUBILITIES OF O2 AND CO2 221 chemc(ji,jj,1) = EXP( zcek1 ) * 1.e-6 * rhop(ji,jj,1) / 1000. ! mol/(L uatm) 222 chemc(ji,jj,2) = ( EXP( zoxy ) * o2atm ) * oxyco ! mol/(L atm) 185 chemc(ji,jj,1) = EXP( zcek1 ) * 1E-6 ! mol/(kg atm) 186 chemc(ji,jj,2) = -1636.75 + 12.0408*ztkel - 0.0327957*ztkel**2 + 0.0000316528*ztkel**3 187 chemc(ji,jj,3) = 57.7 - 0.118*ztkel 223 188 ! 224 189 END DO … … 233 198 !CDIR NOVERRCHK 234 199 DO ji = 1, jpi 235 ztkel = tsn(ji,jj,jk,jp_tem) + 273.1 6200 ztkel = tsn(ji,jj,jk,jp_tem) + 273.15 236 201 zsal = tsn(ji,jj,jk,jp_sal) + ( 1.- tmask(ji,jj,jk) ) * 35. 237 202 zsal2 = zsal * zsal … … 241 206 ztgg4 = ztgg3 * ztgg 242 207 ztgg5 = ztgg4 * ztgg 243 zoxy = ox0 + ox1 * ztgg + ox2 * ztgg2 + ox3 * ztgg3 + ox4 * ztgg4 + ox5 * ztgg5 & 244 + zsal * ( ox6 + ox7 * ztgg + ox8 * ztgg2 + ox9 * ztgg3 ) + ox10 * zsal2 208 209 zoxy = 2.00856 + 3.22400 * ztgg + 3.99063 * ztgg2 + 4.80299 * ztgg3 & 210 & + 9.78188e-1 * ztgg4 + 1.71069 * ztgg5 + zsal * ( -6.24097e-3 & 211 & - 6.93498e-3 * ztgg - 6.90358e-3 * ztgg2 - 4.29155e-3 * ztgg3 ) & 212 & - 3.11680e-7 * zsal2 245 213 chemo2(ji,jj,jk) = ( EXP( zoxy ) * o2atm ) * oxyco * atcox ! mol/(L atm) 246 214 END DO … … 259 227 DO ji = 1, jpi 260 228 261 ! SET PRESSION 262 zpres = 1.025e-1 * fsdept(ji,jj,jk) 229 ! SET PRESSION ACCORDING TO SAUNDER (1980) 230 zplat = SIN ( ABS(gphit(ji,jj)*3.141592654/180.) ) 231 zc1 = 5.92E-3 + zplat**2 * 5.25E-3 232 zpres = ((1-zc1)-SQRT(((1-zc1)**2)-(8.84E-6*fsdept(ji,jj,jk)))) / 4.42E-6 233 zpres = zpres / 10.0 263 234 264 235 ! SET ABSOLUTE TEMPERATURE 265 ztkel = tsn(ji,jj,jk,jp_tem) + 273.1 6236 ztkel = tsn(ji,jj,jk,jp_tem) + 273.15 266 237 zsal = tsn(ji,jj,jk,jp_sal) + ( 1.-tmask(ji,jj,jk) ) * 35. 267 238 zsqrt = SQRT( zsal ) … … 275 246 276 247 ! CHLORINITY (WOOSTER ET AL., 1969) 277 zcl = zsal * salchl248 zcl = zsal / 1.80655 278 249 279 250 ! TOTAL SULFATE CONCENTR. [MOLES/kg soln] 280 zst = st1 * zcl * st2251 zst = 0.14 * zcl /96.062 281 252 282 253 ! TOTAL FLUORIDE CONCENTR. [MOLES/kg soln] 283 zft = ft1 * zcl * ft2254 zft = 0.000067 * zcl /18.9984 284 255 285 256 ! DISSOCIATION CONSTANT FOR SULFATES on free H scale (Dickson 1990) 286 zcks = EXP( ks1 * ztr + ks0 + ks2 * zlogt & 287 & + ( ks3 * ztr + ks4 + ks5 * zlogt ) * zisqrt & 288 & + ( ks6 * ztr + ks7 + ks8 * zlogt ) * zis & 289 & + ks9 * ztr * zis * zisqrt + ks10 * ztr *zis2 + LOG( ks11 + ks12 *zsal ) ) 257 zcks = EXP(-4276.1 * ztr + 141.328 - 23.093 * zlogt & 258 & + (-13856. * ztr + 324.57 - 47.986 * zlogt) * zisqrt & 259 & + (35474. * ztr - 771.54 + 114.723 * zlogt) * zis & 260 & - 2698. * ztr * zis**1.5 + 1776.* ztr * zis2 & 261 & + LOG(1.0 - 0.001005 * zsal)) 262 290 263 291 264 ! DISSOCIATION CONSTANT FOR FLUORIDES on free H scale (Dickson and Riley 79) 292 zckf = EXP( kf1 * ztr + kf0 + kf2 * zisqrt + LOG( kf3 + kf4 * zsal ) ) 293 294 ! DISSOCIATION CONSTANT FOR CARBONATE AND BORATE 295 zckb = ( cb0 + cb1 * zsqrt + cb2 * zsal + cb3 * zsal15 + cb4 * zsal * zsal ) * ztr & 296 & + ( cb5 + cb6 * zsqrt + cb7 * zsal ) & 297 & + ( cb8 + cb9 * zsqrt + cb10 * zsal ) * zlogt + cb11 * zsqrt * ztkel & 298 & + LOG( ( 1.+ zst / zcks + zft / zckf ) / ( 1.+ zst / zcks ) ) 299 300 zck1 = c10 * ztr + c11 + c12 * zlogt + c13 * zsal + c14 * zsal * zsal 301 zck2 = c20 * ztr + c21 + c22 * zsal + c23 * zsal**2 302 303 ! PKW (H2O) (DICKSON AND RILEY, 1979) 304 zckw = cw0 * ztr + cw1 + cw2 * zlogt + ( cw3 * ztr + cw4 + cw5 * zlogt ) * zsqrt + cw6 * zsal 305 265 zckf = EXP( 1590.2*ztr - 12.641 + 1.525*zisqrt & 266 & + LOG(1.0d0 - 0.001005d0*zsal) & 267 & + LOG(1.0d0 + zst/zcks)) 268 269 270 ! DISSOCIATION CONSTANT FOR BORATE 271 zckb= (-8966.90 - 2890.53*zsqrt - 77.942*zsal & 272 & + 1.728*zsal15 - 0.0996*zsal*zsal)*ztr & 273 & + (148.0248 + 137.1942*zsqrt + 1.62142*zsal) & 274 & + (-24.4344 - 25.085*zsqrt - 0.2474*zsal) & 275 & * zlogt + 0.053105*zsqrt*ztkel 276 277 ! DISSOCIATION COEFFICIENT FOR CARBONATE ACCORDING TO 278 ! MEHRBACH (1973) REFIT BY MILLERO (1995), seawater scale 279 zck1 = -1.0*(3633.86*ztr - 61.2172 + 9.6777*zlogt & 280 - 0.011555*zsal + 0.0001152*zsal*zsal) 281 zck2 = -1.0*(471.78*ztr + 25.9290 - 3.16967*zlogt & 282 - 0.01781*zsal + 0.0001122*zsal*zsal) 283 284 285 ! PKW (H2O) (MILLERO, 1995) from composite data 286 zckw = -13847.26 * ztr + 148.9652 - 23.6521 * zlogt + ( 118.67 * ztr & 287 - 5.977 + 1.0495 * zlogt ) * zsqrt - 0.01615 * zsal 288 289 ! CONSTANTS FOR PHOSPHATE (MILLERO, 1995) 290 zck1p = -4576.752*ztr + 115.540 - 18.453*zlogt & 291 & + (-106.736*ztr + 0.69171) * zsqrt & 292 & + (-0.65643*ztr - 0.01844) * zsal 293 294 zck2p = -8814.715*ztr + 172.1033 - 27.927*zlogt & 295 & + (-160.340*ztr + 1.3566)*zsqrt & 296 & + (0.37335*ztr - 0.05778)*zsal 297 298 zck3p = -3070.75*ztr - 18.126 & 299 & + (17.27039*ztr + 2.81197) * zsqrt & 300 & + (-44.99486*ztr - 0.09984) * zsal 301 302 ! CONSTANT FOR SILICATE, MILLERO (1995) 303 zcksi = -8904.2*ztr + 117.400 - 19.334*zlogt & 304 & + (-458.79*ztr + 3.5913) * zisqrt & 305 & + (188.74*ztr - 1.5998) * zis & 306 & + (-12.1652*ztr + 0.07871) * zis2 & 307 & + LOG(1.0 - 0.001005*zsal) 306 308 307 309 ! APPARENT SOLUBILITY PRODUCT K'SP OF CALCITE IN SEAWATER 308 310 ! (S=27-43, T=2-25 DEG C) at pres =0 (atmos. pressure) (MUCCI 1983) 309 zaksp0 = akcc1 + akcc2 * ztkel + akcc3 * ztr + akcc4 * LOG10( ztkel ) & 310 & + ( akcc5 + akcc6 * ztkel + akcc7 * ztr ) * zsqrt + akcc8 * zsal + akcc9 * zsal15 311 zaksp0 = -171.9065 -0.077993*ztkel + 2839.319*ztr + 71.595*LOG10( ztkel ) & 312 & + (-0.77712 + 0.00284263*ztkel + 178.34*ztr) * zsqrt & 313 & - 0.07711*zsal + 0.0041249*zsal15 314 315 ! CONVERT FROM DIFFERENT PH SCALES 316 total2free = 1.0/(1.0 + zst/zcks) 317 free2SWS = 1. + zst/zcks + zft/(zckf*total2free) 318 total2SWS = total2free * free2SWS 319 SWS2total = 1.0 / total2SWS 311 320 312 321 ! K1, K2 OF CARBONIC ACID, KB OF BORIC ACID, KW (H2O) (LIT.?) 313 zak1 = 10**(zck1) 314 zak2 = 10**(zck2) 315 zakb = EXP( zckb )322 zak1 = 10**(zck1) * total2SWS 323 zak2 = 10**(zck2) * total2SWS 324 zakb = EXP( zckb ) * total2SWS 316 325 zakw = EXP( zckw ) 317 326 zaksp1 = 10**(zaksp0) 327 zak1p = exp( zck1p ) 328 zak2p = exp( zck2p ) 329 zak3p = exp( zck3p ) 330 zaksi = exp( zcksi ) 331 zckf = zckf * total2SWS 318 332 319 333 ! FORMULA FOR CPEXP AFTER EDMOND & GIESKES (1970) … … 327 341 ! FORMULA ON P. 1286 IS RIGHT AND CONSISTENT WITH THE 328 342 ! SIGN IN PARTIAL MOLAR VOLUME CHANGE AS SHOWN ON P. 1285)) 329 zcpexp = zpres / (rgas*ztkel)330 zcpexp2 = zpres * z pres/(rgas*ztkel)343 zcpexp = zpres / (rgas*ztkel) 344 zcpexp2 = zpres * zcpexp 331 345 332 346 ! KB OF BORIC ACID, K1,K2 OF CARBONIC ACID PRESSURE … … 334 348 ! (CF. BROECKER ET AL., 1982) 335 349 336 zbuf1 = - ( devk11 + devk21 * ztc + devk31 * ztc * ztc ) 350 zbuf1 = - ( devk10 + devk20 * ztc + devk30 * ztc * ztc ) 351 zbuf2 = 0.5 * ( devk40 + devk50 * ztc ) 352 ak13(ji,jj,jk) = zak1 * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 ) 353 354 zbuf1 = - ( devk11 + devk21 * ztc + devk31 * ztc * ztc ) 337 355 zbuf2 = 0.5 * ( devk41 + devk51 * ztc ) 338 ak 13(ji,jj,jk) = zak1* EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )356 ak23(ji,jj,jk) = zak2 * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 ) 339 357 340 358 zbuf1 = - ( devk12 + devk22 * ztc + devk32 * ztc * ztc ) 341 359 zbuf2 = 0.5 * ( devk42 + devk52 * ztc ) 342 ak 23(ji,jj,jk) = zak2* EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )360 akb3(ji,jj,jk) = zakb * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 ) 343 361 344 362 zbuf1 = - ( devk13 + devk23 * ztc + devk33 * ztc * ztc ) 345 363 zbuf2 = 0.5 * ( devk43 + devk53 * ztc ) 346 ak b3(ji,jj,jk) = zakb* EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )364 akw3(ji,jj,jk) = zakw * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 ) 347 365 348 366 zbuf1 = - ( devk14 + devk24 * ztc + devk34 * ztc * ztc ) 349 367 zbuf2 = 0.5 * ( devk44 + devk54 * ztc ) 350 akw3(ji,jj,jk) = zakw * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 ) 351 368 aks3(ji,jj,jk) = zcks * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 ) 369 370 zbuf1 = - ( devk15 + devk25 * ztc + devk35 * ztc * ztc ) 371 zbuf2 = 0.5 * ( devk45 + devk55 * ztc ) 372 akf3(ji,jj,jk) = zckf * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 ) 373 374 zbuf1 = - ( devk17 + devk27 * ztc + devk37 * ztc * ztc ) 375 zbuf2 = 0.5 * ( devk47 + devk57 * ztc ) 376 ak1p3(ji,jj,jk) = zak1p * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 ) 377 378 zbuf1 = - ( devk18 + devk28 * ztc + devk38 * ztc * ztc ) 379 zbuf2 = 0.5 * ( devk48 + devk58 * ztc ) 380 ak2p3(ji,jj,jk) = zak2p * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 ) 381 382 zbuf1 = - ( devk19 + devk29 * ztc + devk39 * ztc * ztc ) 383 zbuf2 = 0.5 * ( devk49 + devk59 * ztc ) 384 ak3p3(ji,jj,jk) = zak3p * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 ) 385 386 zbuf1 = - ( devk110 + devk210 * ztc + devk310 * ztc * ztc ) 387 zbuf2 = 0.5 * ( devk410 + devk510 * ztc ) 388 aksi3(ji,jj,jk) = zaksi * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 ) 389 390 ! CONVERT FROM DIFFERENT PH SCALES 391 total2free = 1.0/(1.0 + zst/aks3(ji,jj,jk)) 392 free2SWS = 1. + zst/aks3(ji,jj,jk) + zft/akf3(ji,jj,jk) 393 total2SWS = total2free * free2SWS 394 SWS2total = 1.0 / total2SWS 395 396 ! Convert to total scale 397 ak13(ji,jj,jk) = ak13(ji,jj,jk) * SWS2total 398 ak23(ji,jj,jk) = ak23(ji,jj,jk) * SWS2total 399 akb3(ji,jj,jk) = akb3(ji,jj,jk) * SWS2total 400 akw3(ji,jj,jk) = akw3(ji,jj,jk) * SWS2total 401 ak1p3(ji,jj,jk) = ak1p3(ji,jj,jk) * SWS2total 402 ak2p3(ji,jj,jk) = ak2p3(ji,jj,jk) * SWS2total 403 ak3p3(ji,jj,jk) = ak3p3(ji,jj,jk) * SWS2total 404 aksi3(ji,jj,jk) = aksi3(ji,jj,jk) * SWS2total 405 akf3(ji,jj,jk) = akf3(ji,jj,jk) / total2free 352 406 353 407 ! APPARENT SOLUBILITY PRODUCT K'SP OF CALCITE 354 408 ! AS FUNCTION OF PRESSURE FOLLOWING MILLERO 355 409 ! (P. 1285) AND BERNER (1976) 356 zbuf1 = - ( devk1 5 + devk25 * ztc + devk35* ztc * ztc )357 zbuf2 = 0.5 * ( devk4 5 + devk55* ztc )410 zbuf1 = - ( devk16 + devk26 * ztc + devk36 * ztc * ztc ) 411 zbuf2 = 0.5 * ( devk46 + devk56 * ztc ) 358 412 aksp(ji,jj,jk) = zaksp1 * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 ) 359 413 360 ! TOTAL BORATE CONCENTR. [MOLES/L] 361 borat(ji,jj,jk) = bor1 * zcl * bor2 414 ! TOTAL F, S, and BORATE CONCENTR. [MOLES/L] 415 borat(ji,jj,jk) = 0.0002414 * zcl / 10.811 416 sulfat(ji,jj,jk) = zst 417 fluorid(ji,jj,jk) = zft 362 418 363 419 ! Iron and SIO3 saturation concentration from ... 364 420 sio3eq(ji,jj,jk) = EXP( LOG( 10.) * ( 6.44 - 968. / ztkel ) ) * 1.e-6 365 fekeq (ji,jj,jk) = 10**( 17.27 - 1565.7 / ( 273.15 + ztc ) ) 366 421 fekeq (ji,jj,jk) = 10**( 17.27 - 1565.7 / ztkel ) 422 423 ! Liu and Millero (1999) only valid 5 - 50 degC 424 ztkel1 = MAX( 5. , tsn(ji,jj,jk,jp_tem) ) + 273.16 425 fesol(ji,jj,jk,1) = 10**((-13.486) - (0.1856* (zis**0.5)) + (0.3073*zis) + (5254/ztkel1)) 426 fesol(ji,jj,jk,2) = 10**(2.517 - (0.885*(zis**0.5)) + (0.2139 * zis) - (1320/ztkel1) ) 427 fesol(ji,jj,jk,3) = 10**(0.4511 - (0.3305*(ZIS**0.5)) - (1996/ztkel1) ) 428 fesol(ji,jj,jk,4) = 10**(-0.2965 - (0.7881*(zis**0.5)) - (4086/ztkel1) ) 429 fesol(ji,jj,jk,5) = 10**(4.4466 - (0.8505*(zis**0.5)) - (7980/ztkel1) ) 367 430 END DO 368 431 END DO … … 373 436 END SUBROUTINE p4z_che 374 437 438 SUBROUTINE ahini_for_at(p_hini) 439 !!--------------------------------------------------------------------- 440 !! *** ROUTINE ahini_for_at *** 441 !! 442 !! Subroutine returns the root for the 2nd order approximation of the 443 !! DIC -- B_T -- A_CB equation for [H+] (reformulated as a cubic 444 !! polynomial) around the local minimum, if it exists. 445 !! Returns * 1E-03_wp if p_alkcb <= 0 446 !! * 1E-10_wp if p_alkcb >= 2*p_dictot + p_bortot 447 !! * 1E-07_wp if 0 < p_alkcb < 2*p_dictot + p_bortot 448 !! and the 2nd order approximation does not have 449 !! a solution 450 !!--------------------------------------------------------------------- 451 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(OUT) :: p_hini 452 INTEGER :: ji, jj, jk 453 REAL(wp) :: zca1, zba1 454 REAL(wp) :: zd, zsqrtd, zhmin 455 REAL(wp) :: za2, za1, za0 456 REAL(wp) :: p_dictot, p_bortot, p_alkcb 457 458 IF( nn_timing == 1 ) CALL timing_start('ahini_for_at') 459 ! 460 DO jk = 1, jpk 461 DO jj = 1, jpj 462 DO ji = 1, jpj 463 p_alkcb = trb(ji,jj,jk,jptal) * 1000. / (rhop(ji,jj,jk) + rtrn) 464 p_dictot = trb(ji,jj,jk,jpdic) * 1000. / (rhop(ji,jj,jk) + rtrn) 465 p_bortot = borat(ji,jj,jk) 466 IF (p_alkcb <= 0.) THEN 467 p_hini(ji,jj,jk) = 1.e-3 468 ELSEIF (p_alkcb >= (2.*p_dictot + p_bortot)) THEN 469 p_hini(ji,jj,jk) = 1.e-10_wp 470 ELSE 471 zca1 = p_dictot/( p_alkcb + rtrn ) 472 zba1 = p_bortot/ (p_alkcb + rtrn ) 473 ! Coefficients of the cubic polynomial 474 za2 = aKb3(ji,jj,jk)*(1. - zba1) + ak13(ji,jj,jk)*(1.-zca1) 475 za1 = ak13(ji,jj,jk)*akb3(ji,jj,jk)*(1. - zba1 - zca1) & 476 & + ak13(ji,jj,jk)*ak23(ji,jj,jk)*(1. - (zca1+zca1)) 477 za0 = ak13(ji,jj,jk)*ak23(ji,jj,jk)*akb3(ji,jj,jk)*(1. - zba1 - (zca1+zca1)) 478 ! Taylor expansion around the minimum 479 zd = za2*za2 - 3.*za1 ! Discriminant of the quadratic equation 480 ! for the minimum close to the root 481 482 IF(zd > 0.) THEN ! If the discriminant is positive 483 zsqrtd = SQRT(zd) 484 IF(za2 < 0) THEN 485 zhmin = (-za2 + zsqrtd)/3. 486 ELSE 487 zhmin = -za1/(za2 + zsqrtd) 488 ENDIF 489 p_hini(ji,jj,jk) = zhmin + SQRT(-(za0 + zhmin*(za1 + zhmin*(za2 + zhmin)))/zsqrtd) 490 ELSE 491 p_hini(ji,jj,jk) = 1.e-7 492 ENDIF 493 ! 494 ENDIF 495 END DO 496 END DO 497 END DO 498 ! 499 IF( nn_timing == 1 ) CALL timing_stop('ahini_for_at') 500 ! 501 END SUBROUTINE ahini_for_at 502 503 !=============================================================================== 504 SUBROUTINE anw_infsup( p_alknw_inf, p_alknw_sup ) 505 506 ! Subroutine returns the lower and upper bounds of "non-water-selfionization" 507 ! contributions to total alkalinity (the infimum and the supremum), i.e 508 ! inf(TA - [OH-] + [H+]) and sup(TA - [OH-] + [H+]) 509 510 ! Argument variables 511 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(OUT) :: p_alknw_inf 512 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(OUT) :: p_alknw_sup 513 514 p_alknw_inf(:,:,:) = -trb(:,:,:,jppo4) * 1000. / (rhop(:,:,:) + rtrn) - sulfat(:,:,:) & 515 & - fluorid(:,:,:) 516 p_alknw_sup(:,:,:) = (2. * trb(:,:,:,jpdic) + 2. * trb(:,:,:,jppo4) + trb(:,:,:,jpsil) ) & 517 & * 1000. / (rhop(:,:,:) + rtrn) + borat(:,:,:) 518 519 END SUBROUTINE anw_infsup 520 521 522 SUBROUTINE solve_at_general( p_hini, zhi ) 523 524 ! Universal pH solver that converges from any given initial value, 525 ! determines upper an lower bounds for the solution if required 526 527 ! Argument variables 528 !-------------------- 529 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(IN) :: p_hini 530 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(OUT) :: zhi 531 532 ! Local variables 533 !----------------- 534 INTEGER :: ji, jj, jk, jn 535 REAL(wp) :: zh_ini, zh, zh_prev, zh_lnfactor 536 REAL(wp) :: zdelta, zh_delta 537 REAL(wp) :: zeqn, zdeqndh, zalka 538 REAL(wp) :: aphscale 539 REAL(wp) :: znumer_dic, zdnumer_dic, zdenom_dic, zalk_dic, zdalk_dic 540 REAL(wp) :: znumer_bor, zdnumer_bor, zdenom_bor, zalk_bor, zdalk_bor 541 REAL(wp) :: znumer_po4, zdnumer_po4, zdenom_po4, zalk_po4, zdalk_po4 542 REAL(wp) :: znumer_sil, zdnumer_sil, zdenom_sil, zalk_sil, zdalk_sil 543 REAL(wp) :: znumer_nh4, zdnumer_nh4, zdenom_nh4, zalk_nh4, zdalk_nh4 544 REAL(wp) :: znumer_h2s, zdnumer_h2s, zdenom_h2s, zalk_h2s, zdalk_h2s 545 REAL(wp) :: znumer_so4, zdnumer_so4, zdenom_so4, zalk_so4, zdalk_so4 546 REAL(wp) :: znumer_flu, zdnumer_flu, zdenom_flu, zalk_flu, zdalk_flu 547 REAL(wp) :: zalk_wat, zdalk_wat 548 REAL(wp) :: zfact, p_alktot, zdic, zbot, zpt, zst, zft, zsit 549 LOGICAL :: l_exitnow 550 REAL(wp), PARAMETER :: pz_exp_threshold = 1.0 551 REAL(wp), POINTER, DIMENSION(:,:,:) :: zalknw_inf, zalknw_sup, rmask, zh_min, zh_max, zeqn_absmin 552 553 IF( nn_timing == 1 ) CALL timing_start('solve_at_general') 554 ! Allocate temporary workspace 555 CALL wrk_alloc( jpi, jpj, jpk, zalknw_inf, zalknw_sup, rmask ) 556 CALL wrk_alloc( jpi, jpj, jpk, zh_min, zh_max, zeqn_absmin ) 557 558 CALL anw_infsup( zalknw_inf, zalknw_sup ) 559 560 rmask(:,:,:) = tmask(:,:,:) 561 zhi(:,:,:) = 0. 562 563 ! TOTAL H+ scale: conversion factor for Htot = aphscale * Hfree 564 DO jk = 1, jpk 565 DO jj = 1, jpj 566 DO ji = 1, jpj 567 IF (rmask(ji,jj,jk) == 1.) THEN 568 p_alktot = trb(ji,jj,jk,jptal) * 1000. / (rhop(ji,jj,jk) + rtrn) 569 aphscale = 1. + sulfat(ji,jj,jk)/aks3(ji,jj,jk) 570 zh_ini = p_hini(ji,jj,jk) 571 572 zdelta = (p_alktot-zalknw_inf(ji,jj,jk))**2 + 4.*akw3(ji,jj,jk)/aphscale 573 574 IF(p_alktot >= zalknw_inf(ji,jj,jk)) THEN 575 zh_min(ji,jj,jk) = 2.*akw3(ji,jj,jk) /( p_alktot-zalknw_inf(ji,jj,jk) + SQRT(zdelta) ) 576 ELSE 577 zh_min(ji,jj,jk) = aphscale*(-(p_alktot-zalknw_inf(ji,jj,jk)) + SQRT(zdelta) ) / 2. 578 ENDIF 579 580 zdelta = (p_alktot-zalknw_sup(ji,jj,jk))**2 + 4.*akw3(ji,jj,jk)/aphscale 581 582 IF(p_alktot <= zalknw_sup(ji,jj,jk)) THEN 583 zh_max(ji,jj,jk) = aphscale*(-(p_alktot-zalknw_sup(ji,jj,jk)) + SQRT(zdelta) ) / 2. 584 ELSE 585 zh_max(ji,jj,jk) = 2.*akw3(ji,jj,jk) /( p_alktot-zalknw_sup(ji,jj,jk) + SQRT(zdelta) ) 586 ENDIF 587 588 zhi(ji,jj,jk) = MAX(MIN(zh_max(ji,jj,jk), zh_ini), zh_min(ji,jj,jk)) 589 ENDIF 590 END DO 591 END DO 592 END DO 593 594 zeqn_absmin(:,:,:) = HUGE(1._wp) 595 596 DO jn = 1, jp_maxniter_atgen 597 DO jk = 1, jpk 598 DO jj = 1, jpj 599 DO ji = 1, jpj 600 IF (rmask(ji,jj,jk) == 1.) THEN 601 zfact = rhop(ji,jj,jk) / 1000. + rtrn 602 p_alktot = trb(ji,jj,jk,jptal) / zfact 603 zdic = trb(ji,jj,jk,jpdic) / zfact 604 zbot = borat(ji,jj,jk) 605 zpt = trb(ji,jj,jk,jppo4) / zfact * po4r 606 zsit = trb(ji,jj,jk,jpsil) / zfact 607 zst = sulfat (ji,jj,jk) 608 zft = fluorid(ji,jj,jk) 609 aphscale = 1. + sulfat(ji,jj,jk)/aks3(ji,jj,jk) 610 zh = zhi(ji,jj,jk) 611 zh_prev = zh 612 613 ! H2CO3 - HCO3 - CO3 : n=2, m=0 614 znumer_dic = 2.*ak13(ji,jj,jk)*ak23(ji,jj,jk) + zh*ak13(ji,jj,jk) 615 zdenom_dic = ak13(ji,jj,jk)*ak23(ji,jj,jk) + zh*(ak13(ji,jj,jk) + zh) 616 zalk_dic = zdic * (znumer_dic/zdenom_dic) 617 zdnumer_dic = ak13(ji,jj,jk)*ak13(ji,jj,jk)*ak23(ji,jj,jk) + zh & 618 *(4.*ak13(ji,jj,jk)*ak23(ji,jj,jk) + zh*ak13(ji,jj,jk)) 619 zdalk_dic = -zdic*(zdnumer_dic/zdenom_dic**2) 620 621 622 ! B(OH)3 - B(OH)4 : n=1, m=0 623 znumer_bor = akb3(ji,jj,jk) 624 zdenom_bor = akb3(ji,jj,jk) + zh 625 zalk_bor = zbot * (znumer_bor/zdenom_bor) 626 zdnumer_bor = akb3(ji,jj,jk) 627 zdalk_bor = -zbot*(zdnumer_bor/zdenom_bor**2) 628 629 630 ! H3PO4 - H2PO4 - HPO4 - PO4 : n=3, m=1 631 znumer_po4 = 3.*ak1p3(ji,jj,jk)*ak2p3(ji,jj,jk)*ak3p3(ji,jj,jk) & 632 & + zh*(2.*ak1p3(ji,jj,jk)*ak2p3(ji,jj,jk) + zh* ak1p3(ji,jj,jk)) 633 zdenom_po4 = ak1p3(ji,jj,jk)*ak2p3(ji,jj,jk)*ak3p3(ji,jj,jk) & 634 & + zh*( ak1p3(ji,jj,jk)*ak2p3(ji,jj,jk) + zh*(ak1p3(ji,jj,jk) + zh)) 635 zalk_po4 = zpt * (znumer_po4/zdenom_po4 - 1.) ! Zero level of H3PO4 = 1 636 zdnumer_po4 = ak1p3(ji,jj,jk)*ak2p3(ji,jj,jk)*ak1p3(ji,jj,jk)*ak2p3(ji,jj,jk)*ak3p3(ji,jj,jk) & 637 & + zh*(4.*ak1p3(ji,jj,jk)*ak1p3(ji,jj,jk)*ak2p3(ji,jj,jk)*ak3p3(ji,jj,jk) & 638 & + zh*(9.*ak1p3(ji,jj,jk)*ak2p3(ji,jj,jk)*ak3p3(ji,jj,jk) & 639 & + ak1p3(ji,jj,jk)*ak1p3(ji,jj,jk)*ak2p3(ji,jj,jk) & 640 & + zh*(4.*ak1p3(ji,jj,jk)*ak2p3(ji,jj,jk) + zh * ak1p3(ji,jj,jk) ) ) ) 641 zdalk_po4 = -zpt * (zdnumer_po4/zdenom_po4**2) 642 643 ! H4SiO4 - H3SiO4 : n=1, m=0 644 znumer_sil = aksi3(ji,jj,jk) 645 zdenom_sil = aksi3(ji,jj,jk) + zh 646 zalk_sil = zsit * (znumer_sil/zdenom_sil) 647 zdnumer_sil = aksi3(ji,jj,jk) 648 zdalk_sil = -zsit * (zdnumer_sil/zdenom_sil**2) 649 650 ! HSO4 - SO4 : n=1, m=1 651 aphscale = 1.0 + zst/aks3(ji,jj,jk) 652 znumer_so4 = aks3(ji,jj,jk) * aphscale 653 zdenom_so4 = aks3(ji,jj,jk) * aphscale + zh 654 zalk_so4 = zst * (znumer_so4/zdenom_so4 - 1.) 655 zdnumer_so4 = aks3(ji,jj,jk) 656 zdalk_so4 = -zst * (zdnumer_so4/zdenom_so4**2) 657 658 ! HF - F : n=1, m=1 659 znumer_flu = akf3(ji,jj,jk) 660 zdenom_flu = akf3(ji,jj,jk) + zh 661 zalk_flu = zft * (znumer_flu/zdenom_flu - 1.) 662 zdnumer_flu = akf3(ji,jj,jk) 663 zdalk_flu = -zft * (zdnumer_flu/zdenom_flu**2) 664 665 ! H2O - OH 666 aphscale = 1.0 + zst/aks3(ji,jj,jk) 667 zalk_wat = akw3(ji,jj,jk)/zh - zh/aphscale 668 zdalk_wat = -akw3(ji,jj,jk)/zh**2 - 1./aphscale 669 670 ! CALCULATE [ALK]([CO3--], [HCO3-]) 671 zeqn = zalk_dic + zalk_bor + zalk_po4 + zalk_sil & 672 & + zalk_so4 + zalk_flu & 673 & + zalk_wat - p_alktot 674 675 zalka = p_alktot - (zalk_bor + zalk_po4 + zalk_sil & 676 & + zalk_so4 + zalk_flu + zalk_wat) 677 678 zdeqndh = zdalk_dic + zdalk_bor + zdalk_po4 + zdalk_sil & 679 & + zdalk_so4 + zdalk_flu + zdalk_wat 680 681 ! Adapt bracketing interval 682 IF(zeqn > 0._wp) THEN 683 zh_min(ji,jj,jk) = zh_prev 684 ELSEIF(zeqn < 0._wp) THEN 685 zh_max(ji,jj,jk) = zh_prev 686 ENDIF 687 688 IF(ABS(zeqn) >= 0.5_wp*zeqn_absmin(ji,jj,jk)) THEN 689 ! if the function evaluation at the current point is 690 ! not decreasing faster than with a bisection step (at least linearly) 691 ! in absolute value take one bisection step on [ph_min, ph_max] 692 ! ph_new = (ph_min + ph_max)/2d0 693 ! 694 ! In terms of [H]_new: 695 ! [H]_new = 10**(-ph_new) 696 ! = 10**(-(ph_min + ph_max)/2d0) 697 ! = SQRT(10**(-(ph_min + phmax))) 698 ! = SQRT(zh_max * zh_min) 699 zh = SQRT(zh_max(ji,jj,jk) * zh_min(ji,jj,jk)) 700 zh_lnfactor = (zh - zh_prev)/zh_prev ! Required to test convergence below 701 ELSE 702 ! dzeqn/dpH = dzeqn/d[H] * d[H]/dpH 703 ! = -zdeqndh * LOG(10) * [H] 704 ! \Delta pH = -zeqn/(zdeqndh*d[H]/dpH) = zeqn/(zdeqndh*[H]*LOG(10)) 705 ! 706 ! pH_new = pH_old + \deltapH 707 ! 708 ! [H]_new = 10**(-pH_new) 709 ! = 10**(-pH_old - \Delta pH) 710 ! = [H]_old * 10**(-zeqn/(zdeqndh*[H]_old*LOG(10))) 711 ! = [H]_old * EXP(-LOG(10)*zeqn/(zdeqndh*[H]_old*LOG(10))) 712 ! = [H]_old * EXP(-zeqn/(zdeqndh*[H]_old)) 713 714 zh_lnfactor = -zeqn/(zdeqndh*zh_prev) 715 716 IF(ABS(zh_lnfactor) > pz_exp_threshold) THEN 717 zh = zh_prev*EXP(zh_lnfactor) 718 ELSE 719 zh_delta = zh_lnfactor*zh_prev 720 zh = zh_prev + zh_delta 721 ENDIF 722 723 IF( zh < zh_min(ji,jj,jk) ) THEN 724 ! if [H]_new < [H]_min 725 ! i.e., if ph_new > ph_max then 726 ! take one bisection step on [ph_prev, ph_max] 727 ! ph_new = (ph_prev + ph_max)/2d0 728 ! In terms of [H]_new: 729 ! [H]_new = 10**(-ph_new) 730 ! = 10**(-(ph_prev + ph_max)/2d0) 731 ! = SQRT(10**(-(ph_prev + phmax))) 732 ! = SQRT([H]_old*10**(-ph_max)) 733 ! = SQRT([H]_old * zh_min) 734 zh = SQRT(zh_prev * zh_min(ji,jj,jk)) 735 zh_lnfactor = (zh - zh_prev)/zh_prev ! Required to test convergence below 736 ENDIF 737 738 IF( zh > zh_max(ji,jj,jk) ) THEN 739 ! if [H]_new > [H]_max 740 ! i.e., if ph_new < ph_min, then 741 ! take one bisection step on [ph_min, ph_prev] 742 ! ph_new = (ph_prev + ph_min)/2d0 743 ! In terms of [H]_new: 744 ! [H]_new = 10**(-ph_new) 745 ! = 10**(-(ph_prev + ph_min)/2d0) 746 ! = SQRT(10**(-(ph_prev + ph_min))) 747 ! = SQRT([H]_old*10**(-ph_min)) 748 ! = SQRT([H]_old * zhmax) 749 zh = SQRT(zh_prev * zh_max(ji,jj,jk)) 750 zh_lnfactor = (zh - zh_prev)/zh_prev ! Required to test convergence below 751 ENDIF 752 ENDIF 753 754 zeqn_absmin(ji,jj,jk) = MIN( ABS(zeqn), zeqn_absmin(ji,jj,jk)) 755 756 ! Stop iterations once |\delta{[H]}/[H]| < rdel 757 ! <=> |(zh - zh_prev)/zh_prev| = |EXP(-zeqn/(zdeqndh*zh_prev)) -1| < rdel 758 ! |EXP(-zeqn/(zdeqndh*zh_prev)) -1| ~ |zeqn/(zdeqndh*zh_prev)| 759 760 ! Alternatively: 761 ! |\Delta pH| = |zeqn/(zdeqndh*zh_prev*LOG(10))| 762 ! ~ 1/LOG(10) * |\Delta [H]|/[H] 763 ! < 1/LOG(10) * rdel 764 765 ! Hence |zeqn/(zdeqndh*zh)| < rdel 766 767 ! rdel <-- pp_rdel_ah_target 768 l_exitnow = (ABS(zh_lnfactor) < pp_rdel_ah_target) 769 770 IF(l_exitnow) THEN 771 rmask(ji,jj,jk) = 0. 772 ENDIF 773 774 zhi(ji,jj,jk) = zh 775 776 IF(jn >= jp_maxniter_atgen) THEN 777 zhi(ji,jj,jk) = -1._wp 778 ENDIF 779 780 ENDIF 781 END DO 782 END DO 783 END DO 784 END DO 785 ! 786 CALL wrk_dealloc( jpi, jpj, jpk, zalknw_inf, zalknw_sup, rmask ) 787 CALL wrk_dealloc( jpi, jpj, jpk, zh_min, zh_max, zeqn_absmin ) 788 789 790 IF( nn_timing == 1 ) CALL timing_stop('solve_at_general') 791 792 793 END SUBROUTINE solve_at_general 375 794 376 795 INTEGER FUNCTION p4z_che_alloc() … … 378 797 !! *** ROUTINE p4z_che_alloc *** 379 798 !!---------------------------------------------------------------------- 380 ALLOCATE( sio3eq(jpi,jpj,jpk), fekeq(jpi,jpj,jpk), chemc(jpi,jpj,2), chemo2(jpi,jpj,jpk), STAT=p4z_che_alloc ) 799 #if defined key_ligand 800 INTEGER :: ierr(3) ! Local variables 801 #else 802 INTEGER :: ierr(2) ! Local variables 803 #endif 804 !!---------------------------------------------------------------------- 805 806 ierr(:) = 0 807 808 ALLOCATE( sio3eq(jpi,jpj,jpk), fekeq(jpi,jpj,jpk), chemc(jpi,jpj,3), chemo2(jpi,jpj,jpk), STAT=ierr(1) ) 809 810 ALLOCATE( akb3(jpi,jpj,jpk) , & 811 & akw3(jpi,jpj,jpk) , borat (jpi,jpj,jpk) , & 812 & aks3(jpi,jpj,jpk) , akf3(jpi,jpj,jpk) , & 813 & ak1p3(jpi,jpj,jpk) , ak2p3(jpi,jpj,jpk) , & 814 & ak3p3(jpi,jpj,jpk) , aksi3(jpi,jpj,jpk) , & 815 & fluorid(jpi,jpj,jpk) , sulfat(jpi,jpj,jpk) , STAT=ierr(2) ) 816 817 ALLOCATE( fesol(jpi,jpj,jpk,5), STAT=ierr(3) ) 818 819 !* Variable for chemistry of the CO2 cycle 820 p4z_che_alloc = MAXVAL( ierr ) 381 821 ! 382 822 IF( p4z_che_alloc /= 0 ) CALL ctl_warn('p4z_che_alloc : failed to allocate arrays.') … … 396 836 397 837 !!====================================================================== 398 END MODULE p4zche838 END MODULE p4zche -
branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zfechem.F90
r5385 r6453 5 5 !!====================================================================== 6 6 !! History : 3.5 ! 2012-07 (O. Aumont, A. Tagliabue, C. Ethe) Original code 7 !!---------------------------------------------------------------------- 8 #if defined key_pisces 7 !! 3.6 ! 2015-05 (O. Aumont) PISCES quota 8 !!---------------------------------------------------------------------- 9 #if defined key_pisces || defined key_pisces_quota 9 10 !!---------------------------------------------------------------------- 10 11 !! 'key_top' and TOP models 11 !! 'key_pisces ' PISCES bio-model12 !! 'key_pisces*' PISCES bio-model 12 13 !!---------------------------------------------------------------------- 13 14 !! p4z_fechem : Compute remineralization/scavenging of iron … … 33 34 LOGICAL :: ln_fechem !: boolean for complex iron chemistry following Tagliabue and voelker 34 35 LOGICAL :: ln_ligvar !: boolean for variable ligand concentration following Tagliabue and voelker 36 LOGICAL :: ln_fecolloid !: boolean for variable colloidal fraction 35 37 REAL(wp), PUBLIC :: xlam1 !: scavenging rate of Iron 36 38 REAL(wp), PUBLIC :: xlamdust !: scavenging rate of Iron by dust 37 39 REAL(wp), PUBLIC :: ligand !: ligand concentration in the ocean 40 REAL(wp), PUBLIC :: kfep !: rate constant for nanoparticle formation 38 41 39 42 REAL(wp) :: kl1, kl2, kb1, kb2, ks, kpr, spd, con, kth … … 48 51 CONTAINS 49 52 50 SUBROUTINE p4z_fechem( kt, knt )53 SUBROUTINE p4z_fechem( kt, jnt ) 51 54 !!--------------------------------------------------------------------- 52 55 !! *** ROUTINE p4z_fechem *** … … 62 65 !!--------------------------------------------------------------------- 63 66 ! 64 INTEGER, INTENT(in) :: kt, knt ! ocean time step65 ! 66 INTEGER :: ji, jj, jk, jic 67 INTEGER, INTENT(in) :: kt, jnt ! ocean time step 68 ! 69 INTEGER :: ji, jj, jk, jic, jn 67 70 REAL(wp) :: zdep, zlam1a, zlam1b, zlamfac 68 REAL(wp) :: zkeq, zfeequi, zfesatur, zfecoll 71 REAL(wp) :: zkeq, zfeequi, zfesatur, zfecoll, fe3sol 69 72 REAL(wp) :: zdenom1, zscave, zaggdfea, zaggdfeb, zcoag 70 73 REAL(wp) :: ztrc, zdust … … 72 75 REAL(wp) :: zdenom, zdenom2 73 76 #endif 74 REAL(wp), POINTER, DIMENSION(:,:,:) :: zTL1, zFe3, ztotlig 77 REAL(wp), POINTER, DIMENSION(:,:,:) :: zTL1, zFe3, ztotlig, precip 75 78 REAL(wp), POINTER, DIMENSION(:,:,:) :: zFeL1, zFeL2, zTL2, zFe2, zFeP 79 REAL(wp), POINTER, DIMENSION(:,: ) :: zstrn, zstrn2 80 REAL(wp) :: zzFeL1, zzFeL2, zzFe2, zzFeP, zzFe3, zzstrn2 81 REAL(wp) :: zrum, zcodel, zargu, zval, zlight 76 82 REAL(wp) :: zkox, zkph1, zkph2, zph, zionic, ztligand 77 83 REAL(wp) :: za, zb, zc, zkappa1, zkappa2, za0, za1, za2 78 84 REAL(wp) :: zxs, zfunc, zp, zq, zd, zr, zphi, zfff, zp3, zq2 79 REAL(wp) :: ztfe, zoxy 85 REAL(wp) :: ztfe, zoxy, zhplus 80 86 REAL(wp) :: zstep 87 #if defined key_ligand 88 REAL(wp) :: zaggliga, zaggligb 89 REAL(wp) :: dissol, zligco 90 #endif 81 91 CHARACTER (len=25) :: charout 82 92 !!--------------------------------------------------------------------- … … 85 95 ! 86 96 ! Allocate temporary workspace 87 CALL wrk_alloc( jpi, jpj, jpk, zFe3, zFeL1, zTL1, ztotlig )97 CALL wrk_alloc( jpi, jpj, jpk, zFe3, zFeL1, zTL1, ztotlig, precip ) 88 98 zFe3 (:,:,:) = 0. 89 99 zFeL1(:,:,:) = 0. 90 100 zTL1 (:,:,:) = 0. 91 101 IF( ln_fechem ) THEN 102 CALL wrk_alloc( jpi, jpj, zstrn, zstrn2 ) 92 103 CALL wrk_alloc( jpi, jpj, jpk, zFe2, zFeL2, zTL2, zFeP ) 93 104 zFe2 (:,:,:) = 0. … … 104 115 ztotlig(:,:,:) = MIN( ztotlig(:,:,:), 10. ) 105 116 ELSE 117 #if defined key_ligand 118 ztotlig(:,:,:) = trb(:,:,:,jplgw) * 1E9 119 #else 106 120 ztotlig(:,:,:) = ligand * 1E9 121 #endif 107 122 ENDIF 108 123 109 124 IF( ln_fechem ) THEN 125 126 ! compute the day length depending on latitude and the day 127 zrum = REAL( nday_year - 80, wp ) / REAL( nyear_len(1), wp ) 128 zcodel = ASIN( SIN( zrum * rpi * 2._wp ) * SIN( rad * 23.5_wp ) ) 129 130 ! day length in hours 131 zstrn(:,:) = 0. 132 DO jj = 1, jpj 133 DO ji = 1, jpi 134 zargu = TAN( zcodel ) * TAN( gphit(ji,jj) * rad ) 135 zargu = MAX( -1., MIN( 1., zargu ) ) 136 zstrn(ji,jj) = MAX( 0.0, 24. - 2. * ACOS( zargu ) / rad / 15. ) 137 END DO 138 END DO 139 140 ! Maximum light intensity 141 zstrn2(:,:) = zstrn(:,:) / 24. 142 WHERE( zstrn(:,:) < 1.e0 ) zstrn(:,:) = 24. 143 zstrn(:,:) = 24. / zstrn(:,:) 144 110 145 ! ------------------------------------------------------------ 111 146 ! NEW FE CHEMISTRY ROUTINE from Tagliabue and Volker (2009) … … 113 148 ! Chemistry is supposed to be fast enough to be at equilibrium 114 149 ! ------------------------------------------------------------ 150 DO jn = 1, 2 115 151 !CDIR NOVERRCHK 116 152 DO jk = 1, jpkm1 … … 119 155 !CDIR NOVERRCHK 120 156 DO ji = 1, jpi 157 zlight = etot(ji,jj,jk) * zstrn(ji,jj) * float(2-jn) 158 zzstrn2 = zstrn2(ji,jj) * float(2-jn) + (1. - zstrn2(ji,jj) ) * float(jn-1) 121 159 ! Calculate ligand concentrations : assume 2/3rd of excess goes to 122 160 ! strong ligands (L1) and 1/3rd to weak ligands (L2) … … 134 172 zkox = zkox * MAX( 1.e-6, zoxy) / ( chemo2(ji,jj,jk) + rtrn ) 135 173 ! PHOTOREDUCTION of complexed iron : Tagliabue and Arrigo (2006) 136 zkph2 = MAX( 0., 15. * etot(ji,jj,jk) / ( etot(ji,jj,jk)+ 2. ) )174 zkph2 = MAX( 0., 15. * zlight / ( zlight + 2. ) ) 137 175 zkph1 = zkph2 / 5. 138 176 ! pass the dfe concentration from PISCES … … 180 218 ENDIF 181 219 ! solve for the other Fe species 182 z Fe3(ji,jj,jk) = MAX( 0., zxs )183 z Fep(ji,jj,jk) = MAX( 0., ( ks * zFe3(ji,jj,jk)/ kpr ) )220 zzFe3 = MAX( 0., zxs ) 221 zzFep = MAX( 0., ( ks * zzFe3 / kpr ) ) 184 222 zkappa2 = ( kb2 + zkph2 ) / kl2 185 zFeL2(ji,jj,jk) = MAX( 0., ( zFe3(ji,jj,jk) * zTL2(ji,jj,jk) ) / ( zkappa2 + zFe3(ji,jj,jk) ) ) 186 zFeL1(ji,jj,jk) = MAX( 0., ( ztfe / zb - za / zb * zFe3(ji,jj,jk) - zc / zb * zFeL2(ji,jj,jk) ) ) 187 zFe2 (ji,jj,jk) = MAX( 0., ( ( zkph1 * zFeL1(ji,jj,jk) + zkph2 * zFeL2(ji,jj,jk) ) / zkox ) ) 223 zzFeL2 = MAX( 0., ( zzFe3 * zTL2(ji,jj,jk) ) / ( zkappa2 + zzFe3 ) ) 224 zzFeL1 = MAX( 0., ( ztfe / zb - za / zb * zzFe3 - zc / zb * zzFeL2 ) ) 225 zzFe2 = MAX( 0., ( ( zkph1 * zzFeL1 + zkph2 * zzFeL2 ) / zkox ) ) 226 zFe3(ji,jj,jk) = zFe3(ji,jj,jk) + zzFe3 * zzstrn2 227 zFe2(ji,jj,jk) = zFe2(ji,jj,jk) + zzFe2 * zzstrn2 228 zFeL2(ji,jj,jk) = zFeL2(ji,jj,jk) + zzFeL2 * zzstrn2 229 zFeL1(ji,jj,jk) = zFeL1(ji,jj,jk) + zzFeL1 * zzstrn2 230 zFeP(ji,jj,jk) = zFeP(ji,jj,jk) + zzFeP * zzstrn2 188 231 END DO 189 232 END DO 233 END DO 190 234 END DO 191 235 ELSE … … 236 280 zfecoll = ( 0.3 * zFeL1(ji,jj,jk) + 0.5 * zFeL2(ji,jj,jk) ) * 1E-9 237 281 ELSE 238 zfeequi = zFe3(ji,jj,jk) * 1E-9 239 zfecoll = 0.5 * zFeL1(ji,jj,jk) * 1E-9 282 IF (ln_fecolloid) THEN 283 zfeequi = zFe3(ji,jj,jk) * 1E-9 284 zhplus = max( rtrn, hi(ji,jj,jk) ) 285 fe3sol = fesol(ji,jj,jk,1) * ( fesol(ji,jj,jk,2) * zhplus**2 & 286 & + fesol(ji,jj,jk,3) * zhplus + fesol(ji,jj,jk,4) & 287 & + fesol(ji,jj,jk,5) / zhplus ) 288 zfecoll = max( ( 0.1 * zFeL1(ji,jj,jk) * 1E-9 ), ( zFeL1(ji,jj,jk) * 1E-9 -fe3sol ) ) 289 ELSE 290 zfeequi = zFe3(ji,jj,jk) * 1E-9 291 zfecoll = 0.5 * zFeL1(ji,jj,jk) * 1E-9 292 fe3sol = 0. 293 kfep = 0. 294 ENDIF 240 295 ENDIF 241 296 #if defined key_kriest … … 272 327 & + ( 114. * 0.3 * trb(ji,jj,jk,jpdoc) + 5.09E3 * trb(ji,jj,jk,jppoc) ) 273 328 zaggdfea = zlam1a * zstep * zfecoll 329 #if defined key_ligand 330 zligco = max( ( 0.1 * trn(ji,jj,jk,jplgw) ), ( trn(ji,jj,jk,jplgw) - fe3sol ) ) 331 zaggliga = zlam1a * zstep * zligco 332 #endif 333 ! 274 334 #if defined key_kriest 275 335 zaggdfeb = 0. 276 336 ! 337 # if defined key_ligand 338 zaggligb = 0. 339 # endif 340 ! 277 341 tra(ji,jj,jk,jpfer) = tra(ji,jj,jk,jpfer) - zscave - zaggdfea - zaggdfeb - zcoag 278 342 tra(ji,jj,jk,jpsfe) = tra(ji,jj,jk,jpsfe) + zscave * zdenom1 + zaggdfea + zaggdfeb 343 ! 279 344 #else 345 ! 280 346 zlam1b = 3.53E3 * trb(ji,jj,jk,jpgoc) * xdiss(ji,jj,jk) 281 347 zaggdfeb = zlam1b * zstep * zfecoll 282 348 ! 283 tra(ji,jj,jk,jpfer) = tra(ji,jj,jk,jpfer) - zscave - zaggdfea - zaggdfeb - zcoag 349 # if defined key_ligand 350 zaggligb = zlam1b * zstep * zligco 351 # endif 352 ! precipitation of Fe3+, creation of nanoparticles 353 precip(ji,jj,jk) = max( 0., (zfeequi - fe3sol) ) * kfep * zstep 354 ! 355 tra(ji,jj,jk,jpfer) = tra(ji,jj,jk,jpfer) - zscave - zaggdfea - zaggdfeb & 356 & - zcoag - precip(ji,jj,jk) 284 357 tra(ji,jj,jk,jpsfe) = tra(ji,jj,jk,jpsfe) + zscave * zdenom1 + zaggdfea 285 358 tra(ji,jj,jk,jpbfe) = tra(ji,jj,jk,jpbfe) + zscave * zdenom2 + zaggdfeb 359 #endif 360 #if defined key_ligand 361 tra(ji,jj,jk,jpfep) = tra(ji,jj,jk,jpfep) + precip(ji,jj,jk) 362 tra(ji,jj,jk,jplgw) = tra(ji,jj,jk,jplgw) - zaggliga - zaggligb 286 363 #endif 287 364 END DO … … 297 374 ENDIF 298 375 376 #if defined key_ligand 377 IF (.NOT. ln_fechem) THEN 378 plig(:,:,:) = MAX( 0., ( ( zFeL1(:,:,:) * 1E-9 ) / ( trn(:,:,:,jpfer) +rtrn ) ) ) 379 plig(:,:,:) = max(0. , plig(:,:,:) ) 380 ENDIF 381 #endif 382 299 383 ! Output of some diagnostics variables 300 384 ! --------------------------------- 301 IF( lk_iomput .AND. knt == nrdttrc ) THEN 302 IF( iom_use("Fe3") ) CALL iom_put("Fe3" , zFe3 (:,:,:) * tmask(:,:,:) ) ! Fe3+ 303 IF( iom_use("FeL1") ) CALL iom_put("FeL1" , zFeL1 (:,:,:) * tmask(:,:,:) ) ! FeL1 304 IF( iom_use("TL1") ) CALL iom_put("TL1" , zTL1 (:,:,:) * tmask(:,:,:) ) ! TL1 305 IF( iom_use("Totlig") ) CALL iom_put("Totlig" , ztotlig(:,:,:) * tmask(:,:,:) ) ! TL 306 IF( iom_use("Biron") ) CALL iom_put("Biron" , biron (:,:,:) * 1e9 * tmask(:,:,:) ) ! biron 307 IF( ln_fechem ) THEN 308 IF( iom_use("Fe2") ) CALL iom_put("Fe2" , zFe2 (:,:,:) * tmask(:,:,:) ) ! Fe2+ 309 IF( iom_use("FeL2") ) CALL iom_put("FeL2" , zFeL2 (:,:,:) * tmask(:,:,:) ) ! FeL2 310 IF( iom_use("FeP") ) CALL iom_put("FeP" , zFeP (:,:,:) * tmask(:,:,:) ) ! FeP 311 IF( iom_use("TL2") ) CALL iom_put("TL2" , zTL2 (:,:,:) * tmask(:,:,:) ) ! TL2 385 IF( ln_diatrc .AND. lk_iomput ) THEN 386 IF( jnt == nrdttrc ) THEN 387 CALL iom_put("Fe3" , zFe3 (:,:,:) * tmask(:,:,:) ) ! Fe3+ 388 CALL iom_put("FeL1" , zFeL1 (:,:,:) * tmask(:,:,:) ) ! FeL1 389 CALL iom_put("TL1" , zTL1 (:,:,:) * tmask(:,:,:) ) ! TL1 390 CALL iom_put("Totlig" , ztotlig(:,:,:) * tmask(:,:,:) ) ! TL 391 CALL iom_put("Biron" , biron (:,:,:) * 1e9 * tmask(:,:,:) ) ! biron 392 IF( ln_fechem ) THEN 393 CALL iom_put("Fe2" , zFe2 (:,:,:) * tmask(:,:,:) ) ! Fe2+ 394 CALL iom_put("FeL2", zFeL2 (:,:,:) * tmask(:,:,:) ) ! FeL2 395 CALL iom_put("FeP" , zFeP (:,:,:) * tmask(:,:,:) ) ! FeP 396 CALL iom_put("TL2" , zTL2 (:,:,:) * tmask(:,:,:) ) ! TL2 397 ENDIF 312 398 ENDIF 313 399 ENDIF … … 319 405 ENDIF 320 406 ! 321 CALL wrk_dealloc( jpi, jpj, jpk, zFe3, zFeL1, zTL1, ztotlig )407 CALL wrk_dealloc( jpi, jpj, jpk, zFe3, zFeL1, zTL1, ztotlig, precip ) 322 408 IF( ln_fechem ) CALL wrk_dealloc( jpi, jpj, jpk, zFe2, zFeL2, zTL2, zFeP ) 323 409 ! … … 339 425 !! 340 426 !!---------------------------------------------------------------------- 341 NAMELIST/nampisfer/ ln_fechem, ln_ligvar, xlam1, xlamdust, ligand427 NAMELIST/nampisfer/ ln_fechem, ln_ligvar, ln_fecolloid, xlam1, xlamdust, ligand, kfep 342 428 INTEGER :: ios ! Local integer output status for namelist read 343 429 … … 355 441 WRITE(numout,*) ' Namelist parameters for Iron chemistry, nampisfer' 356 442 WRITE(numout,*) ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~' 357 WRITE(numout,*) ' enable complex iron chemistry scheme ln_fechem =', ln_fechem 358 WRITE(numout,*) ' variable concentration of ligand ln_ligvar =', ln_ligvar 359 WRITE(numout,*) ' scavenging rate of Iron xlam1 =', xlam1 360 WRITE(numout,*) ' scavenging rate of Iron by dust xlamdust =', xlamdust 361 WRITE(numout,*) ' ligand concentration in the ocean ligand =', ligand 443 WRITE(numout,*) ' enable complex iron chemistry scheme ln_fechem =', ln_fechem 444 WRITE(numout,*) ' variable concentration of ligand ln_ligvar =', ln_ligvar 445 WRITE(numout,*) ' Variable colloidal fraction of Fe3+ ln_fecolloid =', ln_fecolloid 446 WRITE(numout,*) ' scavenging rate of Iron xlam1 =', xlam1 447 WRITE(numout,*) ' scavenging rate of Iron by dust xlamdust =', xlamdust 448 WRITE(numout,*) ' ligand concentration in the ocean ligand =', ligand 449 WRITE(numout,*) ' rate constant for nanoparticle formation kfep =', kfep 362 450 ENDIF 363 451 ! -
branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zflx.F90
r6204 r6453 10 10 !! 2.0 ! 2007-12 (C. Ethe, G. Madec) F90 11 11 !! ! 2011-02 (J. Simeon, J. Orr) Include total atm P correction 12 !!---------------------------------------------------------------------- 13 #if defined key_pisces 14 !!---------------------------------------------------------------------- 15 !! 'key_pisces' PISCES bio-model 12 !! 3.6 ! 2015-05 (O. Aumont) PISCES quota 13 !!---------------------------------------------------------------------- 14 #if defined key_pisces || defined key_pisces_quota 15 !!---------------------------------------------------------------------- 16 !! 'key_pisces*' PISCES bio-model 16 17 !!---------------------------------------------------------------------- 17 18 !! p4z_flx : CALCULATES GAS EXCHANGE AND CHEMISTRY AT SEA SURFACE … … 52 53 REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:) :: patm ! atmospheric pressure at kt [N/m2] 53 54 TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_patm ! structure of input fields (file informations, fields read) 54 55 LOGICAL, PUBLIC :: ln_presatmco2 !: accounting for spatial atm CO2 in the compuation of carbon flux (T) or not (F) 56 TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_atmco2 ! structure of input fields (file informations, fields read) 55 57 56 58 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: oce_co2 !: ocean carbon flux … … 84 86 ! 85 87 INTEGER :: ji, jj, jm, iind, iindm1 86 REAL(wp) :: ztc, ztc2, ztc3, z ws, zkgwan88 REAL(wp) :: ztc, ztc2, ztc3, ztc4, zws, zkgwan 87 89 REAL(wp) :: zfld, zflu, zfld16, zflu16, zfact 88 REAL(wp) :: zph, zah2, zbot, zdic, zalk, zsch_o2, zalka, zsch_co2 90 REAL(wp) :: zvapsw, zsal, zfco2, zxc2, xCO2approx, ztkel, zfugcoeff 91 REAL(wp) :: zph, zdic, zsch_o2, zsch_co2 89 92 REAL(wp) :: zyr_dec, zdco2dt 90 93 CHARACTER (len=25) :: charout 91 REAL(wp), POINTER, DIMENSION(:,:) :: zkgco2, zkgo2, zh2co3, zoflx, zw2d 94 REAL(wp), POINTER, DIMENSION(:,:) :: zkgco2, zkgo2, zh2co3, zoflx, zw2d 92 95 !!--------------------------------------------------------------------- 93 96 ! … … 103 106 IF( kt /= nit000 .AND. knt == 1 ) CALL p4z_patm( kt ) ! Get sea-level pressure (E&K [1981] climatology) for use in flux calcs 104 107 105 IF( ln_co2int ) THEN108 IF( ln_co2int .AND. .NOT.ln_presatmco2 ) THEN 106 109 ! Linear temporal interpolation of atmospheric pco2. atcco2.txt has annual values. 107 110 ! Caveats: First column of .txt must be in years, decimal years preferably. … … 121 124 #endif 122 125 123 DO jm = 1, 10 124 !CDIR NOVERRCHK 125 DO jj = 1, jpj 126 !CDIR NOVERRCHK 127 DO ji = 1, jpi 128 129 ! DUMMY VARIABLES FOR DIC, H+, AND BORATE 130 zbot = borat(ji,jj,1) 131 zfact = rhop(ji,jj,1) / 1000. + rtrn 132 zdic = trb(ji,jj,1,jpdic) / zfact 133 zph = MAX( hi(ji,jj,1), 1.e-10 ) / zfact 134 zalka = trb(ji,jj,1,jptal) / zfact 135 136 ! CALCULATE [ALK]([CO3--], [HCO3-]) 137 zalk = zalka - ( akw3(ji,jj,1) / zph - zph + zbot / ( 1.+ zph / akb3(ji,jj,1) ) ) 138 139 ! CALCULATE [H+] AND [H2CO3] 140 zah2 = SQRT( (zdic-zalk)**2 + 4.* ( zalk * ak23(ji,jj,1) & 141 & / ak13(ji,jj,1) ) * ( 2.* zdic - zalk ) ) 142 zah2 = 0.5 * ak13(ji,jj,1) / zalk * ( ( zdic - zalk ) + zah2 ) 143 zh2co3(ji,jj) = ( 2.* zdic - zalk ) / ( 2.+ ak13(ji,jj,1) / zah2 ) * zfact 144 hi(ji,jj,1) = zah2 * zfact 145 END DO 126 DO jj = 1, jpj 127 DO ji = 1, jpi 128 ! DUMMY VARIABLES FOR DIC, H+, AND BORATE 129 zfact = rhop(ji,jj,1) / 1000. + rtrn 130 zdic = trb(ji,jj,1,jpdic) / zfact 131 zph = MAX( hi(ji,jj,1), 1.e-10 ) / zfact 132 ! CALCULATE [H2CO3] 133 zh2co3(ji,jj) = zdic/(1. + ak13(ji,jj,1)/zph + ak13(ji,jj,1)*ak23(ji,jj,1)/zph**2)*zfact 146 134 END DO 147 135 END DO 148 149 136 150 137 ! -------------- … … 162 149 ztc2 = ztc * ztc 163 150 ztc3 = ztc * ztc2 151 ztc4 = ztc2 * ztc2 164 152 ! Compute the schmidt Number both O2 and CO2 165 zsch_co2 = 2 073.1 - 125.62 * ztc + 3.6276 * ztc2 - 0.043126 * ztc3166 zsch_o2 = 19 53.4 - 128.0 * ztc + 3.9918 * ztc2 - 0.050091 * ztc3153 zsch_co2 = 2116.8 - 136.25 * ztc + 4.7353 * ztc2 - 0.092307 * ztc3 + 0.0007555 * ztc4 154 zsch_o2 = 1920.4 - 135.6 * ztc + 5.2122 * ztc2 - 0.109390 * ztc3 + 0.0009377 * ztc4 167 155 ! wind speed 168 156 zws = wndm(ji,jj) * wndm(ji,jj) 169 157 ! Compute the piston velocity for O2 and CO2 170 zkgwan = 0. 3 * zws + 2.5 * ( 0.5246 + 0.016256 * ztc + 0.00049946 * ztc2 )158 zkgwan = 0.251 * zws 171 159 zkgwan = zkgwan * xconv * ( 1.- fr_i(ji,jj) ) * tmask(ji,jj,1) 172 160 # if defined key_degrad … … 181 169 DO jj = 1, jpj 182 170 DO ji = 1, jpi 171 ztkel = tsn(ji,jj,1,jp_tem) + 273.15 172 zsal = tsn(ji,jj,1,jp_sal) + ( 1.- tmask(ji,jj,1) ) * 35. 173 zvapsw = exp(24.4543 - 67.4509*(100.0/ztkel) - 4.8489*log(ztkel/100) - 0.000544*zsal) 174 xCO2approx = satmco2(ji,jj) * ( patm(ji,jj) - zvapsw ) * 1.E-6 175 zxc2 = (1.0 - xCO2approx)**2 176 zfugcoeff = exp(patm(ji,jj) * (chemc(ji,jj,2) + 2.0 * zxc2 * chemc(ji,jj,3) ) & 177 & / (82.05736 * ztkel)) 178 zfco2 = satmco2(ji,jj) * ( patm(ji,jj) - zvapsw ) * zfugcoeff 183 179 ! Compute CO2 flux for the sea and air 184 zfld = satmco2(ji,jj) * patm(ji,jj)* tmask(ji,jj,1) * chemc(ji,jj,1) * zkgco2(ji,jj) ! (mol/L) * (m/s)180 zfld = zfco2 * tmask(ji,jj,1) * chemc(ji,jj,1) * zkgco2(ji,jj) ! (mol/L) * (m/s) 185 181 zflu = zh2co3(ji,jj) * tmask(ji,jj,1) * zkgco2(ji,jj) ! (mol/L) (m/s) ? 186 182 oce_co2(ji,jj) = ( zfld - zflu ) * rfact2 * e1e2t(ji,jj) * tmask(ji,jj,1) * 1000. 187 183 ! compute the trend 188 184 tra(ji,jj,1,jpdic) = tra(ji,jj,1,jpdic) + ( zfld - zflu ) * rfact2 / fse3t(ji,jj,1) 189 190 185 ! Compute O2 flux 191 zfld16 = atcox * patm(ji,jj) * chemc(ji,jj,2) * tmask(ji,jj,1) * zkgo2(ji,jj) ! (mol/L) * (m/s)186 zfld16 = patm(ji,jj) * chemo2(ji,jj,1) * tmask(ji,jj,1) * zkgo2(ji,jj) ! (mol/L) * (m/s) 192 187 zflu16 = trb(ji,jj,1,jpoxy) * tmask(ji,jj,1) * zkgo2(ji,jj) 193 188 zoflx(ji,jj) = zfld16 - zflu16 … … 198 193 t_oce_co2_flx = glob_sum( oce_co2(:,:) ) ! Total Flux of Carbon 199 194 t_oce_co2_flx_cum = t_oce_co2_flx_cum + t_oce_co2_flx ! Cumulative Total Flux of Carbon 200 ! t_atm_co2_flx = glob_sum( satmco2(:,:) * e1e2t(:,:) ) ! Total atmospheric pCO2 201 t_atm_co2_flx = atcco2 ! Total atmospheric pCO2 202 195 t_atm_co2_flx = glob_sum( satmco2(:,:) * e1e2t(:,:) ) ! Total atmospheric pCO2 196 203 197 IF(ln_ctl) THEN ! print mean trends (used for debugging) 204 198 WRITE(charout, FMT="('flx ')") … … 208 202 209 203 IF( lk_iomput .AND. knt == nrdttrc ) THEN 210 CALL wrk_alloc( jpi, jpj, zw2d ) 204 CALL wrk_alloc( jpi, jpj, zw2d ) 211 205 IF( iom_use( "Cflx" ) ) THEN 212 206 zw2d(:,:) = oce_co2(:,:) / e1e2t(:,:) * rfact2r 213 CALL iom_put( "Cflx" , zw2d ) 207 CALL iom_put( "Cflx" , zw2d ) 214 208 ENDIF 215 209 IF( iom_use( "Oflx" ) ) THEN … … 222 216 ENDIF 223 217 IF( iom_use( "Dpco2" ) ) THEN 224 zw2d(:,:) = ( satmco2(:,:) * patm(:,:) - zh2co3(:,:) / ( chemc(:,:,1) + rtrn ) ) * tmask(:,:,1)218 zw2d(:,:) = oce_co2(:,:) / e1e2t(:,:) / 1000. * rfact2r * tmask(:,:,1) / ( zkgco2(:,:) * chemc(:,:,1) + rtrn ) 225 219 CALL iom_put( "Dpco2" , zw2d ) 226 220 ENDIF 227 221 IF( iom_use( "Dpo2" ) ) THEN 228 zw2d(:,:) = ( atcox * patm(:,:) - trb(:,:,1,jpoxy) / ( chemc(:,:,2) + rtrn ) )* tmask(:,:,1)222 zw2d(:,:) = ( patm(:,:) - trb(:,:,1,jpoxy) / ( chemo2(:,:,1) + rtrn ) ) * atcox * tmask(:,:,1) 229 223 CALL iom_put( "Dpo2" , zw2d ) 230 224 ENDIF … … 235 229 ELSE 236 230 IF( ln_diatrc ) THEN 237 trc2d(:,:,jp_pcs0_2d ) = oce_co2(:,:) / e1e2t(:,:) * rfact2r 231 trc2d(:,:,jp_pcs0_2d ) = oce_co2(:,:) / e1e2t(:,:) * rfact2r 232 trc2d(:,:,jp_pcs0_2d + 1) = zoflx(:,:) * 1000 * tmask(:,:,1) 233 trc2d(:,:,jp_pcs0_2d + 2) = zkgco2(:,:) * tmask(:,:,1) 234 trc2d(:,:,jp_pcs0_2d + 3) = ( satmco2(:,:) * patm(:,:) - zh2co3(:,:) / ( chemc(:,:,1) + rtrn ) ) * tmask(:,:,1) 235 ENDIF 236 ENDIF 237 238 IF( ln_diatrc ) THEN 239 IF( lk_iomput .AND. knt == nrdttrc ) THEN 240 CALL iom_put( "Cflx" , oce_co2(:,:) / e1e2t(:,:) / rfact ) 241 CALL iom_put( "Oflx" , zoflx(:,:) * 1000 * tmask(:,:,1) ) 242 CALL iom_put( "Kg" , zkgco2(:,:) * tmask(:,:,1) ) 243 CALL iom_put( "Dpco2", ( satmco2(:,:) * patm(:,:) - zh2co3(:,:) / ( chemc(:,:,1) + rtrn ) ) * tmask(:,:,1) ) 244 CALL iom_put( "Dpo2" , ( atcox * patm(:,:) - trb(:,:,1,jpoxy) * atcox / ( chemo2(:,:,1) + rtrn ) ) * tmask(:,:,1) ) 245 ELSE 246 trc2d(:,:,jp_pcs0_2d ) = oce_co2(:,:) / e1e2t(:,:) / rfact 238 247 trc2d(:,:,jp_pcs0_2d + 1) = zoflx(:,:) * 1000 * tmask(:,:,1) 239 248 trc2d(:,:,jp_pcs0_2d + 2) = zkgco2(:,:) * tmask(:,:,1) … … 281 290 WRITE(numout,*) ' ' 282 291 ENDIF 283 IF( .NOT.ln_co2int ) THEN292 IF( .NOT.ln_co2int .AND. .NOT.ln_presatmco2 ) THEN 284 293 IF(lwp) THEN ! control print 285 294 WRITE(numout,*) ' Constant Atmospheric pCO2 value atcco2 =', atcco2 … … 287 296 ENDIF 288 297 satmco2(:,:) = atcco2 ! Initialisation of atmospheric pco2 289 ELSE 298 ELSEIF(ln_co2int .AND. .NOT.ln_presatmco2) THEN 290 299 IF(lwp) THEN 291 300 WRITE(numout,*) ' Atmospheric pCO2 value from file clname =', TRIM( clname ) … … 309 318 END DO 310 319 CLOSE(numco2) 320 ELSEIF(.NOT.ln_co2int .AND. ln_presatmco2) THEN 321 IF(lwp) THEN 322 WRITE(numout,*) ' Spatialized Atmospheric pCO2 from an external file' 323 WRITE(numout,*) ' ' 324 ENDIF 325 ELSE 326 IF(lwp) THEN 327 WRITE(numout,*) ' Spatialized Atmospheric pCO2 from an external file' 328 WRITE(numout,*) ' ' 329 ENDIF 311 330 ENDIF 312 331 ! 313 332 oce_co2(:,:) = 0._wp ! Initialization of Flux of Carbon 333 t_atm_co2_flx = 0._wp 314 334 t_oce_co2_flx = 0._wp 315 t_atm_co2_flx = 0._wp316 335 ! 317 336 CALL p4z_patm( nit000 ) … … 335 354 CHARACTER(len=100) :: cn_dir ! Root directory for location of ssr files 336 355 TYPE(FLD_N) :: sn_patm ! informations about the fields to be read 337 !!338 NAMELIST/nampisatm/ ln_presatm, sn_patm, cn_dir339 356 TYPE(FLD_N) :: sn_atmco2 ! informations about the fields to be read 357 !! 358 NAMELIST/nampisatm/ ln_presatm, ln_presatmco2, sn_patm, sn_atmco2, cn_dir 340 359 ! ! ----------------------- ! 341 360 IF( kt == nit000 ) THEN ! First call kt=nittrc000 ! … … 355 374 WRITE(numout,*) ' Namelist nampisatm : Atmospheric Pressure as external forcing' 356 375 WRITE(numout,*) ' constant atmopsheric pressure (F) or from a file (T) ln_presatm = ', ln_presatm 376 WRITE(numout,*) ' spatial atmopsheric CO2 for flux calcs ln_presatmco2 = ', ln_presatmco2 357 377 WRITE(numout,*) 358 378 ENDIF … … 367 387 ENDIF 368 388 ! 369 IF( .NOT.ln_presatm ) patm(:,:) = 1.e0 ! Initialize patm if no reading from a file 389 IF( ln_presatmco2 ) THEN 390 ALLOCATE( sf_atmco2(1), STAT=ierr ) !* allocate and fill sf_atmco2 (forcing structure) with sn_atmco2 391 IF( ierr > 0 ) CALL ctl_stop( 'STOP', 'p4z_flx: unable to allocate sf_atmco2 structure' ) 392 ! 393 CALL fld_fill( sf_atmco2, (/ sn_atmco2 /), cn_dir, 'p4z_flx', 'Atmospheric co2 partial pressure ', 'nampisatm' ) 394 ALLOCATE( sf_atmco2(1)%fnow(jpi,jpj,1) ) 395 IF( sn_atmco2%ln_tint ) ALLOCATE( sf_atmco2(1)%fdta(jpi,jpj,1,2) ) 396 ENDIF 370 397 ! 371 398 ENDIF … … 374 401 CALL fld_read( kt, 1, sf_patm ) !* input Patm provided at kt + 1/2 375 402 patm(:,:) = sf_patm(1)%fnow(:,:,1) ! atmospheric pressure 403 ELSE 404 patm(:,:) = 1.e0 ! Initialize patm if no reading from a file 405 ENDIF 406 ! 407 IF( ln_presatmco2 ) THEN 408 CALL fld_read( kt, 1, sf_atmco2 ) !* input atmco2 provided at kt + 1/2 409 satmco2(:,:) = sf_atmco2(1)%fnow(:,:,1) ! atmospheric pressure 410 ELSE 411 satmco2(:,:) = atcco2 ! Initialize atmco2 if no reading from a file 412 ENDIF 413 ! 414 IF(lwp) THEN !* control print 415 WRITE(numout,*) ' Min-Max atmospheric CO2 = ', MINVAL(satmco2(:,:)),MAXVAL(satmco2(:,:)) 376 416 ENDIF 377 417 ! … … 400 440 401 441 !!====================================================================== 402 END MODULE p4zflx442 END MODULE p4zflx -
branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zint.F90
r6204 r6453 6 6 !! History : 1.0 ! 2004-03 (O. Aumont) Original code 7 7 !! 2.0 ! 2007-12 (C. Ethe, G. Madec) F90 8 !! 3.6 ! 2015-05 (O. Aumont) PISCES quota 8 9 !!---------------------------------------------------------------------- 9 #if defined key_pisces 10 #if defined key_pisces || defined key_pisces_quota 10 11 !!---------------------------------------------------------------------- 11 !! 'key_pisces ' PISCES bio-model12 !! 'key_pisces*' PISCES bio-model 12 13 !!---------------------------------------------------------------------- 13 14 !! p4z_int : interpolation and computation of various accessory fields … … 81 82 82 83 !!====================================================================== 83 END MODULE p4zint84 END MODULE p4zint -
branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zlim.F90
r6204 r6453 26 26 PUBLIC p4z_lim 27 27 PUBLIC p4z_lim_init 28 PUBLIC p4z_lim_alloc 28 29 29 30 !! * Shared module variables … … 47 48 REAL(wp), PUBLIC :: qdfelim !: optimal Fe quota for diatoms 48 49 REAL(wp), PUBLIC :: caco3r !: mean rainratio 50 51 !!* Phytoplankton limitation terms 52 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: xnanono3 !: ??? 53 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: xdiatno3 !: ??? 54 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: xnanonh4 !: ??? 55 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: xdiatnh4 !: ??? 56 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: xnanopo4 !: ??? 57 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: xdiatpo4 !: ??? 58 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: xlimphy !: ??? 59 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: xlimdia !: ??? 60 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: xlimnfe !: ??? 61 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: xlimdfe !: ??? 62 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: xlimsi !: ??? 63 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: xlimbac !: ?? 64 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: xlimbacl !: ?? 65 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: concdfe !: ??? 66 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: concnfe !: ??? 49 67 50 68 ! Coefficient for iron limitation … … 255 273 END SUBROUTINE p4z_lim_init 256 274 275 INTEGER FUNCTION p4z_lim_alloc() 276 !!---------------------------------------------------------------------- 277 !! *** ROUTINE p5z_lim_alloc *** 278 !!---------------------------------------------------------------------- 279 USE lib_mpp , ONLY: ctl_warn 280 !!---------------------------------------------------------------------- 281 282 !* Biological arrays for phytoplankton growth 283 ALLOCATE( xnanono3(jpi,jpj,jpk), xdiatno3(jpi,jpj,jpk), & 284 & xnanonh4(jpi,jpj,jpk), xdiatnh4(jpi,jpj,jpk), & 285 & xnanopo4(jpi,jpj,jpk), xdiatpo4(jpi,jpj,jpk), & 286 & xlimphy (jpi,jpj,jpk), xlimdia (jpi,jpj,jpk), & 287 & xlimnfe (jpi,jpj,jpk), xlimdfe (jpi,jpj,jpk), & 288 & xlimbac (jpi,jpj,jpk), xlimbacl(jpi,jpj,jpk), & 289 & concnfe (jpi,jpj,jpk), concdfe (jpi,jpj,jpk), & 290 & xlimsi (jpi,jpj,jpk), STAT=p4z_lim_alloc ) 291 ! 292 IF( p4z_lim_alloc /= 0 ) CALL ctl_warn('p4z_lim_alloc : failed to allocate arrays.') 293 ! 294 END FUNCTION p4z_lim_alloc 295 296 257 297 #else 258 298 !!====================================================================== -
branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zlys.F90
r6204 r6453 11 11 !! ! 2011-02 (J. Simeon, J. Orr) Calcon salinity dependence 12 12 !! 3.4 ! 2011-06 (O. Aumont, C. Ethe) Improvment of calcite dissolution 13 !!---------------------------------------------------------------------- 14 #if defined key_pisces 15 !!---------------------------------------------------------------------- 16 !! 'key_pisces' PISCES bio-model 13 !! 3.6 ! 2015-05 (O. Aumont) PISCES quota 14 !!---------------------------------------------------------------------- 15 #if defined key_pisces || defined key_pisces_quota 16 !!---------------------------------------------------------------------- 17 !! 'key_pisces*' PISCES bio-model 17 18 !!---------------------------------------------------------------------- 18 19 !! p4z_lys : Compute the CaCO3 dissolution … … 22 23 USE trc ! passive tracers common variables 23 24 USE sms_pisces ! PISCES Source Minus Sink variables 25 USE p4zche ! Chemical model 24 26 USE prtctl_trc ! print control for debugging 25 27 USE iom ! I/O manager … … 60 62 ! 61 63 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 64 INTEGER :: ji, jj, jk 65 REAL(wp) :: zdispot, zfact, zcalcon 65 66 REAL(wp) :: zomegaca, zexcess, zexcess0 67 REAL(wp) :: zrfact2 66 68 CHARACTER (len=25) :: charout 67 REAL(wp), POINTER, DIMENSION(:,:,:) :: zco3, zcaldiss 69 REAL(wp), POINTER, DIMENSION(:,:,:) :: zco3, zcaldiss, zhinit, zhi 68 70 !!--------------------------------------------------------------------- 69 71 ! 70 72 IF( nn_timing == 1 ) CALL timing_start('p4z_lys') 71 73 ! 72 CALL wrk_alloc( jpi, jpj, jpk, zco3, zcaldiss )74 CALL wrk_alloc( jpi, jpj, jpk, zco3, zcaldiss, zhinit, zhi ) 73 75 ! 74 76 zco3 (:,:,:) = 0. 75 77 zcaldiss(:,:,:) = 0. 78 zhinit(:,:,:) = hi(:,:,:) * 1000. / ( rhop(:,:,:) + rtrn ) 76 79 ! ------------------------------------------- 77 80 ! COMPUTE [CO3--] and [H+] CONCENTRATIONS 78 81 ! ------------------------------------------- 79 82 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 + borat(ji,jj,jk) / ( 1. + zph / akb3(ji,jj,jk) ) ) 94 ! CALCULATE [H+] and [CO3--] 95 zaldi = zdic - zalk 96 zah2 = SQRT( zaldi * zaldi + 4.* ( zalk * ak23(ji,jj,jk) / ak13(ji,jj,jk) ) * ( zdic + zaldi ) ) 97 zah2 = 0.5 * ak13(ji,jj,jk) / zalk * ( zaldi + zah2 ) 98 ! 99 zco3(ji,jj,jk) = zalk / ( 2. + zah2 / ak23(ji,jj,jk) ) * zfact 100 hi(ji,jj,jk) = zah2 * zfact 101 END DO 83 CALL solve_at_general(zhinit, zhi) 84 85 DO jk = 1, jpkm1 86 DO jj = 1, jpj 87 DO ji = 1, jpi 88 zfact = rhop(ji,jj,jk) / 1000. + rtrn 89 zco3(ji,jj,jk) = trb(ji,jj,jk,jpdic) * ak13(ji,jj,jk) * ak23(ji,jj,jk) / (zhi(ji,jj,jk)**2 & 90 & + ak13(ji,jj,jk) * zhi(ji,jj,jk) + ak13(ji,jj,jk) * ak23(ji,jj,jk) + rtrn ) 91 hi(ji,jj,jk) = zhi(ji,jj,jk) * zfact 102 92 END DO 103 93 END DO 104 !105 94 END DO 106 95 … … 125 114 zexcess0 = MAX( 0., excess(ji,jj,jk) ) 126 115 zexcess = zexcess0**nca 127 128 116 ! AMOUNT CACO3 (12C) THAT RE-ENTERS SOLUTION 129 117 ! (ACCORDING TO THIS FORMULATION ALSO SOME PARTICULATE … … 133 121 zdispot = zdispot * facvol(ji,jj,jk) 134 122 # endif 135 ! CHANGE OF [CO3--] , [ALK], PARTICULATE [CACO3], 136 ! AND [SUM(CO2)] DUE TO CACO3 DISSOLUTION/PRECIPITATION 137 zcaldiss(ji,jj,jk) = zdispot * rfact2 / rmtss ! calcite dissolution 138 zco3(ji,jj,jk) = zco3(ji,jj,jk) + zcaldiss(ji,jj,jk) 139 ! 140 tra(ji,jj,jk,jptal) = tra(ji,jj,jk,jptal) + 2. * zcaldiss(ji,jj,jk) 141 tra(ji,jj,jk,jpcal) = tra(ji,jj,jk,jpcal) - zcaldiss(ji,jj,jk) 142 tra(ji,jj,jk,jpdic) = tra(ji,jj,jk,jpdic) + zcaldiss(ji,jj,jk) 123 ! CHANGE OF [CO3--] , [ALK], PARTICULATE [CACO3], 124 ! AND [SUM(CO2)] DUE TO CACO3 DISSOLUTION/PRECIPITATION 125 zcaldiss(ji,jj,jk) = zdispot * rfact2 / rmtss ! calcite dissolution 126 ! 127 tra(ji,jj,jk,jptal) = tra(ji,jj,jk,jptal) + 2. * zcaldiss(ji,jj,jk) 128 tra(ji,jj,jk,jpcal) = tra(ji,jj,jk,jpcal) - zcaldiss(ji,jj,jk) 129 tra(ji,jj,jk,jpdic) = tra(ji,jj,jk,jpdic) + zcaldiss(ji,jj,jk) 143 130 END DO 144 131 END DO 145 132 END DO 146 133 ! 147 134 ! 148 135 IF( lk_iomput .AND. knt == nrdttrc ) THEN 149 136 IF( iom_use( "PH" ) ) CALL iom_put( "PH" , -1. * LOG10( hi(:,:,:) ) * tmask(:,:,:) ) … … 151 138 IF( iom_use( "CO3sat" ) ) CALL iom_put( "CO3sat", aksp(:,:,:) * 1.e+3 / calcon * tmask(:,:,:) ) 152 139 IF( iom_use( "DCAL" ) ) CALL iom_put( "DCAL" , zcaldiss(:,:,:) * 1.e+3 * rfact2r * tmask(:,:,:) ) 140 153 141 ELSE 154 142 IF( ln_diatrc ) THEN … … 158 146 ENDIF 159 147 ENDIF 148 ! 160 149 ! 161 150 IF(ln_ctl) THEN ! print mean trends (used for debugging) … … 165 154 ENDIF 166 155 ! 167 CALL wrk_dealloc( jpi, jpj, jpk, zco3, zcaldiss )156 CALL wrk_dealloc( jpi, jpj, jpk, zco3, zcaldiss, zhinit, zhi ) 168 157 ! 169 158 IF( nn_timing == 1 ) CALL timing_stop('p4z_lys') … … 184 173 !! 185 174 !!---------------------------------------------------------------------- 186 INTEGER :: ji, jj, jk187 175 INTEGER :: ios ! Local integer output status for namelist read 188 REAL(wp) :: zcaralk, zbicarb, zco3189 REAL(wp) :: ztmas, ztmas1190 176 191 177 NAMELIST/nampiscal/ kdca, nca … … 225 211 #endif 226 212 !!====================================================================== 227 END MODULE p4zlys213 END MODULE p4zlys -
branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zmeso.F90
r6204 r6453 73 73 REAL(wp) :: zgraze2 , zdenom, zdenom2 74 74 REAL(wp) :: zfact , zstep, zfood, zfoodlim, zproport 75 REAL(wp) :: zmortzgoc, zfrac, zfracfe, zratio, zratio2 75 REAL(wp) :: zmortzgoc, zfrac, zfracfe, zratio, zratio2, zfracal, zgrazcal 76 76 REAL(wp) :: zepshert, zepsherv, zgrarsig, zgraztot, zgraztotn, zgraztotf 77 77 REAL(wp) :: zgrarem2, zgrafer2, zgrapoc2, zprcaca, zmortz2, zgrasrat, zgrasratn … … 204 204 tra(ji,jj,jk,jpnh4) = tra(ji,jj,jk,jpnh4) + zgrarsig 205 205 tra(ji,jj,jk,jpdoc) = tra(ji,jj,jk,jpdoc) + zgrarem2 - zgrarsig 206 #if defined key_ligand 207 tra(ji,jj,jk,jplgw) = tra(ji,jj,jk,jplgw) + (zgrarem2 - zgrarsig) * ldocz 208 #endif 206 209 tra(ji,jj,jk,jpoxy) = tra(ji,jj,jk,jpoxy) - o2ut * zgrarsig 207 210 tra(ji,jj,jk,jpfer) = tra(ji,jj,jk,jpfer) + zgrafer2 208 211 tra(ji,jj,jk,jpdic) = tra(ji,jj,jk,jpdic) + zgrarsig 209 212 tra(ji,jj,jk,jptal) = tra(ji,jj,jk,jptal) + rno3 * zgrarsig 210 211 213 zmortz2 = ztortz2 + zrespz2 212 214 zmortzgoc = unass2 / ( 1. - epsher2 ) * ztortz2 + zrespz2 … … 222 224 tra(ji,jj,jk,jpdfe) = tra(ji,jj,jk,jpdfe) - zgrazf 223 225 224 ! calcite production225 zprcaca = xfracal(ji,jj,jk) * zgrazn226 prodcal(ji,jj,jk) = prodcal(ji,jj,jk) + zprcaca ! prodcal=prodcal(nanophy)+prodcal(microzoo)+prodcal(mesozoo)227 !228 zprcaca = part2 * zprcaca229 tra(ji,jj,jk,jpdic) = tra(ji,jj,jk,jpdic) - zprcaca230 tra(ji,jj,jk,jptal) = tra(ji,jj,jk,jptal) - 2. * zprcaca231 tra(ji,jj,jk,jpcal) = tra(ji,jj,jk,jpcal) + zprcaca232 226 #if defined key_kriest 233 227 znumpoc = trb(ji,jj,jk,jpnum) / ( trb(ji,jj,jk,jppoc) + rtrn ) … … 237 231 tra(ji,jj,jk,jpsfe) = tra(ji,jj,jk,jpsfe) + ferat3 * zmortzgoc - zgrazfffp - zgrazpof & 238 232 & + zgraztotf * unass2 233 zfracal = trb(ji,jj,jk,jpcal) / (trb(ji,jj,jk,jppoc) + rtrn ) 234 zgrazcal = ( zgrazffep + zgrazpoc ) * (1. - part2) * zfracal 239 235 #else 240 236 tra(ji,jj,jk,jppoc) = tra(ji,jj,jk,jppoc) - zgrazpoc - zgrazffep + zfrac 237 prodpoc(ji,jj,jk) = prodpoc(ji,jj,jk) + zfrac 238 conspoc(ji,jj,jk) = conspoc(ji,jj,jk) - zgrazpoc - zgrazffep 241 239 tra(ji,jj,jk,jpgoc) = tra(ji,jj,jk,jpgoc) + zmortzgoc - zgrazffeg + zgrapoc2 - zfrac 240 prodgoc(ji,jj,jk) = prodgoc(ji,jj,jk) + zmortzgoc + zgrapoc2 241 consgoc(ji,jj,jk) = consgoc(ji,jj,jk) - zgrazffeg - zfrac 242 242 tra(ji,jj,jk,jpsfe) = tra(ji,jj,jk,jpsfe) - zgrazpof - zgrazfffp + zfracfe 243 243 tra(ji,jj,jk,jpbfe) = tra(ji,jj,jk,jpbfe) + ferat3 * zmortzgoc - zgrazfffg & 244 244 & + zgraztotf * unass2 - zfracfe 245 zfracal = trb(ji,jj,jk,jpcal) / (trb(ji,jj,jk,jppoc) + trb(ji,jj,jk,jpgoc) + rtrn ) 246 zgrazcal = zgrazffeg * (1. - part2) * zfracal 245 247 #endif 248 249 ! calcite production 250 zprcaca = xfracal(ji,jj,jk) * zgrazn 251 prodcal(ji,jj,jk) = prodcal(ji,jj,jk) + zprcaca ! prodcal=prodcal(nanophy)+prodcal(microzoo)+prodcal(mesozoo) 252 ! 253 zprcaca = part2 * zprcaca 254 tra(ji,jj,jk,jpdic) = tra(ji,jj,jk,jpdic) + zgrazcal - zprcaca 255 tra(ji,jj,jk,jptal) = tra(ji,jj,jk,jptal) - 2. * ( zgrazcal + zprcaca ) 256 tra(ji,jj,jk,jpcal) = tra(ji,jj,jk,jpcal) - zgrazcal + zprcaca 257 246 258 END DO 247 259 END DO -
branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zmicro.F90
r6204 r6453 20 20 USE p4zlim ! Co-limitations 21 21 USE p4zsink ! vertical flux of particulate matter due to sinking 22 USE p4zint ! interpolation and computation of various fields23 22 USE p4zprod ! production 24 23 USE iom ! I/O manager … … 150 149 tra(ji,jj,jk,jpnh4) = tra(ji,jj,jk,jpnh4) + zgrarsig 151 150 tra(ji,jj,jk,jpdoc) = tra(ji,jj,jk,jpdoc) + zgrarem - zgrarsig 151 #if defined key_ligand 152 tra(ji,jj,jk,jplgw) = tra(ji,jj,jk,jplgw) + (zgrarem - zgrarsig) * ldocz 153 #endif 152 154 tra(ji,jj,jk,jpoxy) = tra(ji,jj,jk,jpoxy) - o2ut * zgrarsig 153 155 tra(ji,jj,jk,jpfer) = tra(ji,jj,jk,jpfer) + zgrafer 154 156 tra(ji,jj,jk,jppoc) = tra(ji,jj,jk,jppoc) + zgrapoc 157 prodpoc(ji,jj,jk) = prodpoc(ji,jj,jk) + zgrapoc 155 158 tra(ji,jj,jk,jpsfe) = tra(ji,jj,jk,jpsfe) + zgraztotf * unass 156 159 tra(ji,jj,jk,jpdic) = tra(ji,jj,jk,jpdic) + zgrarsig … … 172 175 tra(ji,jj,jk,jpdfe) = tra(ji,jj,jk,jpdfe) - zgrazsf 173 176 tra(ji,jj,jk,jppoc) = tra(ji,jj,jk,jppoc) + zmortz - zgrazm 177 prodpoc(ji,jj,jk) = prodpoc(ji,jj,jk) + zmortz 178 conspoc(ji,jj,jk) = conspoc(ji,jj,jk) - zgrazm 174 179 tra(ji,jj,jk,jpsfe) = tra(ji,jj,jk,jpsfe) + ferat3 * zmortz - zgrazmf 175 180 ! -
branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zmort.F90
r6204 r6453 17 17 USE trc ! passive tracers common variables 18 18 USE sms_pisces ! PISCES Source Minus Sink variables 19 USE p4zsink ! vertical flux of particulate matter due to sinking20 19 USE p4zprod ! Primary productivity 20 USE p4zlim ! Phytoplankton limitation terms 21 21 USE prtctl_trc ! print control for debugging 22 22 … … 128 128 tra(ji,jj,jk,jpgoc) = tra(ji,jj,jk,jpgoc) + zfracal * zmortp 129 129 tra(ji,jj,jk,jppoc) = tra(ji,jj,jk,jppoc) + ( 1. - zfracal ) * zmortp 130 prodpoc(ji,jj,jk) = prodpoc(ji,jj,jk) + ( 1. - zfracal ) * zmortp 131 prodgoc(ji,jj,jk) = prodgoc(ji,jj,jk) + zfracal * zmortp 130 132 tra(ji,jj,jk,jpsfe) = tra(ji,jj,jk,jpsfe) + ( 1. - zfracal ) * zmortp * zfactfe 131 133 tra(ji,jj,jk,jpbfe) = tra(ji,jj,jk,jpbfe) + zfracal * zmortp * zfactfe … … 211 213 tra(ji,jj,jk,jpgoc) = tra(ji,jj,jk,jpgoc) + zrespp2 + 0.5 * ztortp2 212 214 tra(ji,jj,jk,jppoc) = tra(ji,jj,jk,jppoc) + 0.5 * ztortp2 215 prodpoc(ji,jj,jk) = prodpoc(ji,jj,jk) + 0.5 * ztortp2 216 prodgoc(ji,jj,jk) = prodgoc(ji,jj,jk) + zrespp2 + 0.5 * ztortp2 213 217 tra(ji,jj,jk,jpsfe) = tra(ji,jj,jk,jpsfe) + 0.5 * ztortp2 * zfactfe 214 218 tra(ji,jj,jk,jpbfe) = tra(ji,jj,jk,jpbfe) + ( zrespp2 + 0.5 * ztortp2 ) * zfactfe -
branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zopt.F90
r6204 r6453 9 9 !! 3.4 ! 2011-06 (O. Aumont, C. Ethe) Improve light availability of nano & diat 10 10 !!---------------------------------------------------------------------- 11 #if defined key_pisces 12 !!---------------------------------------------------------------------- 13 !! 'key_pisces 'PISCES bio-model11 #if defined key_pisces || defined key_pisces_quota 12 !!---------------------------------------------------------------------- 13 !! 'key_pisces*' PISCES bio-model 14 14 !!---------------------------------------------------------------------- 15 15 !! p4z_opt : light availability in the water column … … 43 43 44 44 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: enano, ediat !: PAR for phyto, nano and diat 45 #if defined key_pisces_quota 46 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: epico !: PAR for pico 47 #endif 45 48 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: etot_ndcy !: PAR over 24h in case of diurnal cycle 46 49 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: emoy !: averaged PAR in the mixed layer … … 77 80 REAL(wp) :: zc0 , zc1 , zc2, zc3, z1_dep 78 81 REAL(wp), POINTER, DIMENSION(:,: ) :: zdepmoy, zetmp1, zetmp2, zetmp3, zetmp4, zqsr100 82 #if defined key_pisces_quota 83 REAL(wp), POINTER, DIMENSION(:,: ) :: zetmp5 84 #endif 79 85 REAL(wp), POINTER, DIMENSION(:,:,:) :: zpar, ze0, ze1, ze2, ze3 80 86 !!--------------------------------------------------------------------- … … 85 91 CALL wrk_alloc( jpi, jpj, zqsr100, zdepmoy, zetmp1, zetmp2, zetmp3, zetmp4 ) 86 92 CALL wrk_alloc( jpi, jpj, jpk, zpar, ze0, ze1, ze2, ze3 ) 93 #if defined key_pisces_quota 94 CALL wrk_alloc( jpi, jpj, zetmp5 ) 95 #endif 87 96 88 97 IF( knt == 1 .AND. ln_varpar ) CALL p4z_opt_sbc( kt ) … … 99 108 !CDIR NOVERRCHK 100 109 DO ji = 1, jpi 110 #if defined key_pisces_quota 111 zchl = ( trb(ji,jj,jk,jpnch) + trb(ji,jj,jk,jpdch) + trb(ji,jj,jk,jppch) + rtrn ) * 1.e6 112 #else 101 113 zchl = ( trb(ji,jj,jk,jpnch) + trb(ji,jj,jk,jpdch) + rtrn ) * 1.e6 114 #endif 102 115 zchl = MIN( 10. , MAX( 0.05, zchl ) ) 103 116 irgb = NINT( 41 + 20.* LOG10( zchl ) + rtrn ) … … 121 134 enano (:,:,jk) = 2.1 * ze1(:,:,jk) + 0.42 * ze2(:,:,jk) + 0.4 * ze3(:,:,jk) 122 135 ediat (:,:,jk) = 1.6 * ze1(:,:,jk) + 0.69 * ze2(:,:,jk) + 0.7 * ze3(:,:,jk) 136 #if defined key_pisces_quota 137 epico (:,:,jk) = 2.1 * ze1(:,:,jk) + 0.42 * ze2(:,:,jk) + 0.4 * ze3(:,:,jk) 138 #endif 123 139 END DO 124 140 ! … … 139 155 enano(:,:,jk) = 2.1 * ze1(:,:,jk) + 0.42 * ze2(:,:,jk) + 0.4 * ze3(:,:,jk) 140 156 ediat(:,:,jk) = 1.6 * ze1(:,:,jk) + 0.69 * ze2(:,:,jk) + 0.7 * ze3(:,:,jk) 157 #if defined key_pisces_quota 158 epico(:,:,jk) = 2.1 * ze1(:,:,jk) + 0.42 * ze2(:,:,jk) + 0.4 * ze3(:,:,jk) 159 #endif 141 160 END DO 142 161 etot_ndcy(:,:,:) = etot(:,:,:) … … 177 196 zetmp3 (:,:) = 0.e0 178 197 zetmp4 (:,:) = 0.e0 198 #if defined key_pisces_quota 199 zetmp5 (:,:) = 0.e0 200 #endif 179 201 180 202 DO jk = 1, nksrp … … 188 210 zetmp3 (ji,jj) = zetmp3 (ji,jj) + enano (ji,jj,jk) * fse3t(ji,jj,jk) ! production 189 211 zetmp4 (ji,jj) = zetmp4 (ji,jj) + ediat (ji,jj,jk) * fse3t(ji,jj,jk) ! production 212 #if defined key_pisces_quota 213 zetmp5 (ji,jj) = zetmp5 (ji,jj) + epico (ji,jj,jk) * fse3t(ji,jj,jk) ! production 214 #endif 190 215 zdepmoy(ji,jj) = zdepmoy(ji,jj) + fse3t(ji,jj,jk) 191 216 ENDIF … … 208 233 enano(ji,jj,jk) = zetmp3(ji,jj) * z1_dep 209 234 ediat(ji,jj,jk) = zetmp4(ji,jj) * z1_dep 235 #if defined key_pisces_quota 236 epico(ji,jj,jk) = zetmp5(ji,jj) * z1_dep 237 #endif 210 238 ENDIF 211 239 END DO … … 228 256 CALL wrk_dealloc( jpi, jpj, zqsr100, zdepmoy, zetmp1, zetmp2, zetmp3, zetmp4 ) 229 257 CALL wrk_dealloc( jpi, jpj, jpk, zpar, ze0, ze1, ze2, ze3 ) 258 #if defined key_pisces_quota 259 CALL wrk_dealloc( jpi, jpj, zetmp5 ) 260 #endif 230 261 ! 231 262 IF( nn_timing == 1 ) CALL timing_stop('p4z_opt') … … 410 441 enano (:,:,:) = 0._wp 411 442 ediat (:,:,:) = 0._wp 443 #if defined key_pisces_quota 444 epico (:,:,:) = 0._wp 445 #endif 412 446 IF( ln_qsr_bio ) etot3 (:,:,:) = 0._wp 413 447 ! … … 423 457 ALLOCATE( ekb(jpi,jpj,jpk) , ekr(jpi,jpj,jpk), ekg(jpi,jpj,jpk), & 424 458 & enano(jpi,jpj,jpk) , ediat(jpi,jpj,jpk), & 459 #if defined key_pisces_quota 460 & epico(jpi,jpj,jpk) , & 461 #endif 425 462 & etot_ndcy(jpi,jpj,jpk), emoy (jpi,jpj,jpk), STAT=p4z_opt_alloc ) 426 463 ! -
branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zprod.F90
r6204 r6453 80 80 REAL(wp) :: zratio, zmax, zsilim, ztn, zadap 81 81 REAL(wp) :: zlim, zsilfac2, zsiborn, zprod, zproreg, zproreg2 82 REAL(wp) :: zmxltst, zmxlday, zmaxday 82 REAL(wp) :: zmxltst, zmxlday, zmaxday, zdocprod 83 83 REAL(wp) :: zpislopen , zpislope2n 84 REAL(wp) :: zrum, zcodel, zargu, zval 84 REAL(wp) :: zrum, zcodel, zargu, zval, zfeup 85 85 REAL(wp) :: zfact 86 86 CHARACTER (len=25) :: charout … … 390 390 zproreg = zprorca(ji,jj,jk) - zpronew(ji,jj,jk) 391 391 zproreg2 = zprorcad(ji,jj,jk) - zpronewd(ji,jj,jk) 392 zdocprod = excret2 * zprorcad(ji,jj,jk) + excret * zprorca(ji,jj,jk) 392 393 tra(ji,jj,jk,jppo4) = tra(ji,jj,jk,jppo4) - zprorca(ji,jj,jk) - zprorcad(ji,jj,jk) 393 394 tra(ji,jj,jk,jpno3) = tra(ji,jj,jk,jpno3) - zpronew(ji,jj,jk) - zpronewd(ji,jj,jk) … … 400 401 tra(ji,jj,jk,jpdfe) = tra(ji,jj,jk,jpdfe) + zprofed(ji,jj,jk) * texcret2 401 402 tra(ji,jj,jk,jpdsi) = tra(ji,jj,jk,jpdsi) + zprorcad(ji,jj,jk) * zysopt(ji,jj,jk) * texcret2 402 tra(ji,jj,jk,jpdoc) = tra(ji,jj,jk,jpdoc) + excret2 * zprorcad(ji,jj,jk) + excret * zprorca(ji,jj,jk) 403 tra(ji,jj,jk,jpdoc) = tra(ji,jj,jk,jpdoc) + zdocprod 404 zfeup = texcret * zprofen(ji,jj,jk) + texcret2 * zprofed(ji,jj,jk) 405 #if defined key_ligand 406 tra(ji,jj,jk,jplgw) = tra(ji,jj,jk,jplgw) + zdocprod * ldocp - zfeup * plig(ji,jj,jk) * lthet 407 #endif 403 408 tra(ji,jj,jk,jpoxy) = tra(ji,jj,jk,jpoxy) + o2ut * ( zproreg + zproreg2) & 404 409 & + ( o2ut + o2nit ) * ( zpronew(ji,jj,jk) + zpronewd(ji,jj,jk) ) 405 tra(ji,jj,jk,jpfer) = tra(ji,jj,jk,jpfer) - texcret * zprofen(ji,jj,jk) - texcret2 * zprofed(ji,jj,jk)410 tra(ji,jj,jk,jpfer) = tra(ji,jj,jk,jpfer) - zfeup 406 411 tra(ji,jj,jk,jpsil) = tra(ji,jj,jk,jpsil) - texcret2 * zprorcad(ji,jj,jk) * zysopt(ji,jj,jk) 407 412 tra(ji,jj,jk,jpdic) = tra(ji,jj,jk,jpdic) - zprorca(ji,jj,jk) - zprorcad(ji,jj,jk) -
branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zrem.F90
r5385 r6453 7 7 !! 2.0 ! 2007-12 (C. Ethe, G. Madec) F90 8 8 !! 3.4 ! 2011-06 (O. Aumont, C. Ethe) Quota model for iron 9 !! 3.6 ! 2016-03 (O. Aumont) Quota model and reorganization 9 10 !!---------------------------------------------------------------------- 10 11 #if defined key_pisces … … 22 23 USE p4zopt ! optical model 23 24 USE p4zche ! chemical model 25 USE p4zlim ! Phytoplankton limitation factors 24 26 USE p4zprod ! Growth rate of the 2 phyto groups 25 USE p4zmeso ! Sources and sinks of mesozooplankton26 USE p4zint ! interpolation and computation of various fields27 USE p4zlim28 27 USE prtctl_trc ! print control for debugging 29 28 USE iom ! I/O manager … … 38 37 39 38 !! * Shared module variables 40 REAL(wp), PUBLIC :: xremik !: remineralisation rate of POC 41 REAL(wp), PUBLIC :: xremip !: remineralisation rate of DOC 39 REAL(wp), PUBLIC :: xremik !: remineralisation rate of DOC 42 40 REAL(wp), PUBLIC :: nitrif !: NH4 nitrification rate 43 REAL(wp), PUBLIC :: xsirem !: remineralisation rate of POC44 REAL(wp), PUBLIC :: xsiremlab !: fast remineralisation rate of POC41 REAL(wp), PUBLIC :: xsirem !: remineralisation rate of BSi 42 REAL(wp), PUBLIC :: xsiremlab !: fast remineralisation rate of BSi 45 43 REAL(wp), PUBLIC :: xsilab !: fraction of labile biogenic silica 46 44 REAL(wp), PUBLIC :: oxymin !: halk saturation constant for anoxia 47 45 48 49 46 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: denitr !: denitrification array 50 47 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: denitnh4 !: - - - - - 48 51 49 52 50 !!* Substitution … … 59 57 CONTAINS 60 58 61 SUBROUTINE p4z_rem( kt, knt )59 SUBROUTINE p4z_rem( kt, jnt ) 62 60 !!--------------------------------------------------------------------- 63 61 !! *** ROUTINE p4z_rem *** … … 68 66 !!--------------------------------------------------------------------- 69 67 ! 70 INTEGER, INTENT(in) :: kt, knt ! ocean time step68 INTEGER, INTENT(in) :: kt, jnt ! ocean time step 71 69 ! 72 70 INTEGER :: ji, jj, jk 73 REAL(wp) :: zremi p, zremik, zsiremin71 REAL(wp) :: zremik, zsiremin 74 72 REAL(wp) :: zsatur, zsatur2, znusil, znusil2, zdep, zdepmin, zfactdep 75 REAL(wp) :: zbactfer, zo rem, zorem2, zofer, zolimit73 REAL(wp) :: zbactfer, zolimit 76 74 REAL(wp) :: zosil, ztem 77 #if ! defined key_kriest 78 REAL(wp) :: zofer2 79 #endif 80 REAL(wp) :: zonitr, zstep, zfact 75 REAL(wp) :: zonitr, zstep, zrfact2 81 76 CHARACTER (len=25) :: charout 82 77 REAL(wp), POINTER, DIMENSION(:,: ) :: ztempbac 83 REAL(wp), POINTER, DIMENSION(:,:,:) :: zdepbac, zolimi, zdepprod , zw3d78 REAL(wp), POINTER, DIMENSION(:,:,:) :: zdepbac, zolimi, zdepprod 84 79 !!--------------------------------------------------------------------- 85 80 ! … … 87 82 ! 88 83 ! Allocate temporary workspace 89 CALL wrk_alloc( jpi, jpj, ztempbac 84 CALL wrk_alloc( jpi, jpj, ztempbac ) 90 85 CALL wrk_alloc( jpi, jpj, jpk, zdepbac, zdepprod, zolimi ) 86 87 ! Initialization of local variables 88 ! --------------------------------- 91 89 92 90 ! Initialisation of temprary arrys … … 219 217 zstep = zstep * facvol(ji,jj,jk) 220 218 # endif 221 ! POC disaggregation by turbulence and bacterial activity.222 ! --------------------------------------------------------223 zremip = xremip * zstep * tgfunc(ji,jj,jk) * ( 1.- 0.55 * nitrfac(ji,jj,jk) )224 225 ! POC disaggregation rate is reduced in anoxic zone as shown by226 ! sediment traps data. In oxic area, the exponent of the martin s227 ! law is around -0.87. In anoxic zone, it is around -0.35. This228 ! means a disaggregation constant about 0.5 the value in oxic zones229 ! -----------------------------------------------------------------230 zorem = zremip * trb(ji,jj,jk,jppoc)231 zofer = zremip * trb(ji,jj,jk,jpsfe)232 #if ! defined key_kriest233 zorem2 = zremip * trb(ji,jj,jk,jpgoc)234 zofer2 = zremip * trb(ji,jj,jk,jpbfe)235 #else236 zorem2 = zremip * trb(ji,jj,jk,jpnum)237 #endif238 239 ! Update the appropriate tracers trends240 ! -------------------------------------241 242 tra(ji,jj,jk,jpdoc) = tra(ji,jj,jk,jpdoc) + zorem243 tra(ji,jj,jk,jpfer) = tra(ji,jj,jk,jpfer) + zofer244 #if defined key_kriest245 tra(ji,jj,jk,jppoc) = tra(ji,jj,jk,jppoc) - zorem246 tra(ji,jj,jk,jpnum) = tra(ji,jj,jk,jpnum) - zorem2247 tra(ji,jj,jk,jpsfe) = tra(ji,jj,jk,jpsfe) - zofer248 #else249 tra(ji,jj,jk,jppoc) = tra(ji,jj,jk,jppoc) + zorem2 - zorem250 tra(ji,jj,jk,jpgoc) = tra(ji,jj,jk,jpgoc) - zorem2251 tra(ji,jj,jk,jpsfe) = tra(ji,jj,jk,jpsfe) + zofer2 - zofer252 tra(ji,jj,jk,jpbfe) = tra(ji,jj,jk,jpbfe) - zofer2253 #endif254 255 END DO256 END DO257 END DO258 259 IF(ln_ctl) THEN ! print mean trends (used for debugging)260 WRITE(charout, FMT="('rem3')")261 CALL prt_ctl_trc_info(charout)262 CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm)263 ENDIF264 265 DO jk = 1, jpkm1266 DO jj = 1, jpj267 DO ji = 1, jpi268 zstep = xstep269 # if defined key_degrad270 zstep = zstep * facvol(ji,jj,jk)271 # endif272 219 ! Remineralization rate of BSi depedant on T and saturation 273 220 ! --------------------------------------------------------- … … 297 244 298 245 IF(ln_ctl) THEN ! print mean trends (used for debugging) 299 WRITE(charout, FMT="('rem 4')")246 WRITE(charout, FMT="('rem3')") 300 247 CALL prt_ctl_trc_info(charout) 301 248 CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm) … … 315 262 END DO 316 263 317 IF( knt == nrdttrc ) THEN 318 CALL wrk_alloc( jpi, jpj, jpk, zw3d ) 319 zfact = 1.e+3 * rfact2r ! conversion from mol/l/kt to mol/m3/s 320 ! 321 IF( iom_use( "REMIN" ) ) THEN 322 zw3d(:,:,:) = zolimi(:,:,:) * tmask(:,:,:) * zfact ! Remineralisation rate 323 CALL iom_put( "REMIN" , zw3d ) 324 ENDIF 325 IF( iom_use( "DENIT" ) ) THEN 326 zw3d(:,:,:) = denitr(:,:,:) * rdenit * rno3 * tmask(:,:,:) * zfact ! Denitrification 327 CALL iom_put( "DENIT" , zw3d ) 328 ENDIF 329 ! 330 CALL wrk_dealloc( jpi, jpj, jpk, zw3d ) 331 ENDIF 264 IF( ln_diatrc .AND. lk_iomput .AND. jnt == nrdttrc ) THEN 265 zrfact2 = 1.e3 * rfact2r 266 CALL iom_put( "REMIN" , zolimi(:,:,:) * tmask(:,:,:) * zrfact2 ) ! Remineralisation rate 267 CALL iom_put( "DENIT" , denitr(:,:,:) * rdenit * rno3 * tmask(:,:,:) * zrfact2 ) ! Denitrification 268 ENDIF 332 269 333 270 IF(ln_ctl) THEN ! print mean trends (used for debugging) 334 WRITE(charout, FMT="('rem 6')")271 WRITE(charout, FMT="('rem4')") 335 272 CALL prt_ctl_trc_info(charout) 336 273 CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm) 337 274 ENDIF 338 275 ! 339 CALL wrk_dealloc( jpi, jpj, ztempbac 276 CALL wrk_dealloc( jpi, jpj, ztempbac ) 340 277 CALL wrk_dealloc( jpi, jpj, jpk, zdepbac, zdepprod, zolimi ) 341 278 ! … … 357 294 !! 358 295 !!---------------------------------------------------------------------- 359 NAMELIST/nampisrem/ xremik, xremip,nitrif, xsirem, xsiremlab, xsilab, &296 NAMELIST/nampisrem/ xremik, nitrif, xsirem, xsiremlab, xsilab, & 360 297 & oxymin 361 298 INTEGER :: ios ! Local integer output status for namelist read … … 374 311 WRITE(numout,*) ' Namelist parameters for remineralization, nampisrem' 375 312 WRITE(numout,*) ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~' 376 WRITE(numout,*) ' remineralisation rate of POC xremip =', xremip377 313 WRITE(numout,*) ' remineralization rate of DOC xremik =', xremik 378 314 WRITE(numout,*) ' remineralization rate of Si xsirem =', xsirem … … 386 322 denitr (:,:,:) = 0._wp 387 323 denitnh4(:,:,:) = 0._wp 388 ! 324 389 325 END SUBROUTINE p4z_rem_init 390 326 -
branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsbc.F90
r6204 r6453 6 6 !! History : 3.5 ! 2012-07 (O. Aumont, C. Ethe) Original code 7 7 !!---------------------------------------------------------------------- 8 #if defined key_pisces 8 #if defined key_pisces || defined key_pisces_quota 9 9 !!---------------------------------------------------------------------- 10 10 !! 'key_pisces' PISCES bio-model … … 42 42 REAL(wp), PUBLIC :: concfediaz !: Fe half-saturation Cste for diazotrophs 43 43 REAL(wp) :: hratio !: Fe:3He ratio assumed for vent iron supply 44 #if defined key_ligand 45 REAL(wp), PUBLIC :: fep_rats !: Fep/Fer ratio from sed sources 46 REAL(wp), PUBLIC :: fep_rath !: Fep/Fer ratio from hydro sources 47 #endif 44 48 45 49 LOGICAL , PUBLIC :: ll_sbc … … 72 76 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: rivdic, rivalk !: river input fields 73 77 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: rivdin, rivdip !: river input fields 78 #if defined key_pisces_quota 79 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: rivdon, rivdop !: river input fields 80 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: rivdoc !: river input fields 81 #endif 74 82 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: rivdsi !: river input fields 75 83 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: nitdep !: atmospheric N deposition … … 80 88 REAL(wp), PUBLIC :: rivdininput, rivdipinput, rivdsiinput 81 89 82 83 90 !!* Substitution 84 91 # include "top_substitute.h90" 85 92 !!---------------------------------------------------------------------- 86 93 !! NEMO/TOP 3.3 , NEMO Consortium (2010) 87 !! $ Id$94 !! $Header:$ 88 95 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 89 96 !!---------------------------------------------------------------------- … … 140 147 DO jj = 1, jpj 141 148 DO ji = 1, jpi 142 zcoef = ryyss * e1e2t(ji,jj) * h_rnf(ji,jj)149 zcoef = ryyss * cvol(ji,jj,1) 143 150 rivalk(ji,jj) = sf_river(jr_dic)%fnow(ji,jj,1) & 144 & * 1.E3 / ( 12. * zcoef + rtrn ) 151 & * 1.E3 / ( 12. * zcoef + rtrn ) * tmask(ji,jj,1) 152 #if defined key_pisces_quota 153 rivdic(ji,jj) = ( sf_river(jr_dic)%fnow(ji,jj,1) ) & 154 & * 1.E3 / ( 12. * zcoef + rtrn ) * tmask(ji,jj,1) 155 rivdin(ji,jj) = ( sf_river(jr_din)%fnow(ji,jj,1) ) & 156 & * 1.E3 / rno3 / ( 14. * zcoef + rtrn ) * tmask(ji,jj,1) 157 rivdip(ji,jj) = ( sf_river(jr_dip)%fnow(ji,jj,1) ) & 158 & * 1.E3 / po4r / ( 31. * zcoef + rtrn ) * tmask(ji,jj,1) 159 rivdoc(ji,jj) = ( sf_river(jr_doc)%fnow(ji,jj,1) ) & 160 & * 1.E3 / ( 12. * zcoef + rtrn ) * tmask(ji,jj,1) 161 rivdon(ji,jj) = ( sf_river(jr_don)%fnow(ji,jj,1) ) & 162 & * 1.E3 / rno3 / ( 14. * zcoef + rtrn ) * tmask(ji,jj,1) 163 rivdop(ji,jj) = ( sf_river(jr_dop)%fnow(ji,jj,1) ) & 164 & * 1.E3 / po4r / ( 31. * zcoef + rtrn ) * tmask(ji,jj,1) 165 #else 145 166 rivdic(ji,jj) = ( sf_river(jr_dic)%fnow(ji,jj,1) + sf_river(jr_doc)%fnow(ji,jj,1) ) & 146 & * 1.E3 / ( 12. * zcoef + rtrn)167 & * 1.E3 / ( 12. * zcoef + rtrn ) * tmask(ji,jj,1) 147 168 rivdin(ji,jj) = ( sf_river(jr_din)%fnow(ji,jj,1) + sf_river(jr_don)%fnow(ji,jj,1) ) & 148 & * 1.E3 / rno3 / ( 14. * zcoef + rtrn ) 169 & * 1.E3 / rno3 / ( 14. * zcoef + rtrn ) * tmask(ji,jj,1) 149 170 rivdip(ji,jj) = ( sf_river(jr_dip)%fnow(ji,jj,1) + sf_river(jr_dop)%fnow(ji,jj,1) ) & 150 & * 1.E3 / po4r / ( 31. * zcoef + rtrn ) 171 & * 1.E3 / po4r / ( 31. * zcoef + rtrn ) * tmask(ji,jj,1) 172 #endif 151 173 rivdsi(ji,jj) = sf_river(jr_dsi)%fnow(ji,jj,1) & 152 & * 1.E3 / ( 28.1 * zcoef + rtrn)174 & * 1.E3 / ( 28.1 * zcoef + rtrn ) * tmask(ji,jj,1) 153 175 END DO 154 176 END DO … … 192 214 INTEGER :: ios ! Local integer output status for namelist read 193 215 INTEGER :: ik50 ! last level where depth less than 50 m 194 INTEGER :: isrow ! index for ORCA1 starting row195 216 REAL(wp) :: zexpide, zdenitide, zmaskt 196 217 REAL(wp) :: ztimes_dust, ztimes_riv, ztimes_ndep … … 208 229 & sn_riverdip, sn_riverdop, sn_riverdsi, sn_ndepo, sn_ironsed, sn_hydrofe, & 209 230 & ln_dust, ln_solub, ln_river, ln_ndepo, ln_ironsed, ln_ironice, ln_hydrofe, & 210 & sedfeinput, dustsolub, icefeinput, wdust, mfrac, nitrfix, diazolight, concfediaz, hratio 231 & sedfeinput, dustsolub, icefeinput, wdust, mfrac, nitrfix, diazolight, concfediaz, & 232 #if defined key_ligand 233 & fep_rats, fep_rath, & 234 #endif 235 & hratio 211 236 !!---------------------------------------------------------------------- 212 237 ! … … 222 247 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nampissbc in configuration namelist', lwp ) 223 248 IF(lwm) WRITE ( numonp, nampissbc ) 224 225 IF ( ( nn_ice_tr >= 0 ) .AND. ln_ironice ) THEN226 IF(lwp) THEN227 WRITE(numout,*) ' ln_ironice incompatible with nn_ice_tr = ', nn_ice_tr228 WRITE(numout,*) ' Specify your sea ice iron concentration in nampisice instead '229 WRITE(numout,*) ' ln_ironice is forced to .FALSE. '230 ln_ironice = .FALSE.231 ENDIF232 ENDIF233 249 234 250 IF(lwp) THEN … … 252 268 WRITE(numout,*) ' fe half-saturation cste for diazotrophs concfediaz = ', concfediaz 253 269 WRITE(numout,*) ' Fe to 3He ratio assumed for vent iron supply hratio = ', hratio 270 #if defined key_ligand 271 WRITE(numout,*) ' Fep/Fer ratio from sed sources = ', fep_rats 272 WRITE(numout,*) ' Fep/Fer ratio from sed hydro sources = ', fep_rath 273 #endif 254 274 END IF 255 275 … … 337 357 ! 338 358 ALLOCATE( rivdic(jpi,jpj), rivalk(jpi,jpj), rivdin(jpi,jpj), rivdip(jpi,jpj), rivdsi(jpi,jpj) ) 359 #if defined key_pisces_quota 360 ALLOCATE( rivdon(jpi,jpj), rivdop(jpi,jpj), rivdoc(jpi,jpj) ) 361 #endif 339 362 ! 340 363 ALLOCATE( sf_river(jpriv), rivinput(jpriv), STAT=ierr1 ) !* allocate and fill sf_river (forcing structure) with sn_river_ … … 380 403 rivalkinput = 0._wp 381 404 END IF 405 382 406 ! nutrient input from dust 383 407 ! ------------------------ … … 449 473 END DO 450 474 END DO 451 ! 475 IF( cp_cfg == 'orca' .AND. jp_cfg == 2 ) THEN 476 ii0 = 176 ; ii1 = 176 ! Southern Island : Kerguelen 477 ij0 = 37 ; ij1 = 37 ; zcmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1), 1:jpk ) = 0.3_wp 478 ! 479 ii0 = 119 ; ii1 = 119 ! South Georgia 480 ij0 = 29 ; ij1 = 29 ; zcmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1), 1:jpk ) = 0.3_wp 481 ! 482 ii0 = 111 ; ii1 = 111 ! Falklands 483 ij0 = 35 ; ij1 = 35 ; zcmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1), 1:jpk ) = 0.3_wp 484 ! 485 ii0 = 168 ; ii1 = 168 ! Crozet 486 ij0 = 40 ; ij1 = 40 ; zcmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1), 1:jpk ) = 0.3_wp 487 ! 488 ii0 = 119 ; ii1 = 119 ! South Orkney 489 ij0 = 28 ; ij1 = 28 ; zcmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1), 1:jpk ) = 0.3_wp 490 ! 491 ii0 = 140 ; ii1 = 140 ! Bouvet Island 492 ij0 = 33 ; ij1 = 33 ; zcmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1), 1:jpk ) = 0.3_wp 493 ! 494 ii0 = 178 ; ii1 = 178 ! Prince edwards 495 ij0 = 34 ; ij1 = 34 ; zcmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1), 1:jpk ) = 0.3_wp 496 ! 497 ii0 = 43 ; ii1 = 43 ! Balleny islands 498 ij0 = 21 ; ij1 = 21 ; zcmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1), 1:jpk ) = 0.3_wp 499 ENDIF 452 500 CALL lbc_lnk( zcmask , 'T', 1. ) ! lateral boundary conditions on cmask (sign unchanged) 453 !454 501 DO jk = 1, jpk 455 502 DO jj = 1, jpj … … 457 504 zexpide = MIN( 8.,( fsdept(ji,jj,jk) / 500. )**(-1.5) ) 458 505 zdenitide = -0.9543 + 0.7662 * LOG( zexpide ) - 0.235 * LOG( zexpide )**2 459 zcmask(ji,jj,jk) = zcmask(ji,jj,jk) * MIN( 1., EXP( zdenitide ) / 0.5 ) 506 zcmask(ji,jj,jk) = zcmask(ji,jj,jk) * MIN( 1., EXP( zdenitide ) / 0.5 ) * tmask(ji,jj,jk) 460 507 END DO 461 508 END DO … … 465 512 ironsed(:,:,jpk) = 0._wp 466 513 DO jk = 1, jpkm1 467 ironsed(:,:,jk) = sedfeinput * zcmask(:,:,jk) / ( fse3t(:,:,jk) * rday ) 514 ironsed(:,:,jk) = sedfeinput * zcmask(:,:,jk) / ( fse3t(:,:,jk) * rday ) 468 515 END DO 469 516 DEALLOCATE( zcmask) … … 483 530 CALL iom_close( numhydro ) 484 531 ! 485 hydrofe(:,:,:) = ( hydrofe(:,:,:) * hratio ) / ( cvol(:,:,:) * ryyss + rtrn ) / 1000._wp 532 hydrofe(:,:,:) = ( hydrofe(:,:,:) * hratio ) / ( cvol(:,:,:) * ryyss + rtrn ) / 1000._wp * tmask(:,:,:) 486 533 ! 487 534 ENDIF … … 519 566 520 567 !!====================================================================== 521 END MODULE p4zsbc568 END MODULE p4zsbc -
branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsed.F90
r6204 r6453 75 75 REAL(wp), POINTER, DIMENSION(:,: ) :: zwsbio3, zwsbio4, zwscal 76 76 REAL(wp), POINTER, DIMENSION(:,:,:) :: zirondep, zsoufer 77 #if defined key_ligand 78 REAL(wp) :: zwssfep 79 REAL(wp), POINTER, DIMENSION(:,: ) :: zwsfep 80 #endif 77 81 !!--------------------------------------------------------------------- 78 82 ! … … 85 89 CALL wrk_alloc( jpi, jpj, zwsbio3, zwsbio4, zwscal ) 86 90 CALL wrk_alloc( jpi, jpj, jpk, zsoufer ) 91 #if defined key_ligand 92 CALL wrk_alloc( jpi, jpj, zwsfep ) 93 #endif 87 94 88 95 zdenit2d(:,:) = 0.e0 … … 186 193 IF( ln_ironsed ) THEN 187 194 tra(:,:,:,jpfer) = tra(:,:,:,jpfer) + ironsed(:,:,:) * rfact2 195 #if defined key_ligand 196 trn(:,:,:,jpfep) = trn(:,:,:,jpfep) + (ironsed(:,:,:) * fep_rats ) * rfact2 197 #endif 188 198 ! 189 199 IF( lk_iomput .AND. knt == nrdttrc .AND. iom_use( "Ironsed" ) ) & … … 195 205 IF( ln_hydrofe ) THEN 196 206 tra(:,:,:,jpfer) = tra(:,:,:,jpfer) + hydrofe(:,:,:) * rfact2 207 #if defined key_ligand 208 trn(:,:,:,jpfep) = trn(:,:,:,jpfep) + ( hydrofe(:,:,:) * fep_rath ) * rfact2 209 trn(:,:,:,jplgw) = trn(:,:,:,jplgw) + ( hydrofe(:,:,:) * 0.5 ) * rfact2 210 #endif 197 211 ! 198 212 IF( lk_iomput .AND. knt == nrdttrc .AND. iom_use( "HYDR" ) ) & … … 210 224 zwscal (ji,jj) = MIN( 0.99 * zdep, wscal (ji,jj,ikt) ) 211 225 zwsbio3(ji,jj) = MIN( 0.99 * zdep, wsbio3(ji,jj,ikt) ) 226 #if defined key_ligand 227 zwsfep(ji,jj) = MIN( 0.99 * zdep, wsfep(ji,jj,ikt) ) 228 #endif 212 229 END DO 213 230 END DO … … 308 325 zws4 = zwsbio4(ji,jj) * zdep 309 326 zws3 = zwsbio3(ji,jj) * zdep 327 #if defined key_ligand 328 zwssfep = zwsfep(ji,jj) * zdep 329 #endif 310 330 zrivno3 = 1. - zbureff(ji,jj) 311 331 # if ! defined key_kriest … … 315 335 tra(ji,jj,ikt,jpsfe) = tra(ji,jj,ikt,jpsfe) - trb(ji,jj,ikt,jpsfe) * zws3 316 336 zwstpoc = trb(ji,jj,ikt,jpgoc) * zws4 + trb(ji,jj,ikt,jppoc) * zws3 337 # if defined key_ligand 338 trn(ji,jj,ikt,jpfep) = trn(ji,jj,ikt,jpfep) - trn(ji,jj,ikt,jpfep) * zwssfep 339 # endif 317 340 # else 318 341 tra(ji,jj,ikt,jpnum) = tra(ji,jj,ikt,jpnum) - trb(ji,jj,ikt,jpnum) * zws4 … … 320 343 tra(ji,jj,ikt,jpsfe) = tra(ji,jj,ikt,jpsfe) - trb(ji,jj,ikt,jpsfe) * zws3 321 344 zwstpoc = trb(ji,jj,ikt,jppoc) * zws3 345 # if defined key_ligand 346 trn(ji,jj,ikt,jpfep) = trn(ji,jj,ikt,jpfep) - trn(ji,jj,ikt,jpfep) * zwssfep 347 # endif 322 348 # endif 323 349 … … 407 433 CALL wrk_dealloc( jpi, jpj, zwsbio3, zwsbio4, zwscal ) 408 434 CALL wrk_dealloc( jpi, jpj, jpk, zsoufer ) 435 #if defined key_ligand 436 CALL wrk_dealloc( jpi, jpj, zwsfep ) 437 #endif 409 438 ! 410 439 IF( nn_timing == 1 ) CALL timing_stop('p4z_sed') -
branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsink.F90
r6204 r6453 40 40 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: sinkfer2 !: Big iron sinking fluxes 41 41 #endif 42 #if defined key_ligand 43 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: sinkfep !: Fep sinking fluxes 44 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: wsfep 45 #endif 42 46 43 47 INTEGER :: ik100 … … 91 95 INTEGER :: ji, jj, jk, jit 92 96 INTEGER :: iiter1, iiter2 93 REAL(wp) :: zagg1, zagg2, zagg3, zagg494 REAL(wp) :: zagg , zaggfe, zaggdoc, zaggdoc2, zaggdoc395 97 REAL(wp) :: zfact, zwsmax, zmax, zstep 96 98 CHARACTER (len=25) :: charout … … 100 102 ! 101 103 IF( nn_timing == 1 ) CALL timing_start('p4z_sink') 104 102 105 ! 103 106 ! Sinking speeds of detritus is increased with depth as shown … … 108 111 DO ji = 1,jpi 109 112 zmax = MAX( heup(ji,jj), hmld(ji,jj) ) 110 zfact = MAX( 0., fsdepw(ji,jj,jk+1) - zmax ) / 5000._wp111 wsbio4(ji,jj,jk) = wsbio2 + ( 200.- wsbio2) * zfact113 zfact = MAX( 0., fsdepw(ji,jj,jk+1) - zmax ) / wsbio2scale 114 wsbio4(ji,jj,jk) = wsbio2 + MAX(0., ( wsbio2max - wsbio2 )) * zfact 112 115 END DO 113 116 END DO … … 117 120 wsbio3(:,:,:) = wsbio 118 121 wscal (:,:,:) = wsbio4(:,:,:) 122 #if defined key_ligand 123 wsfep (:,:,:) = wfep 124 #endif 119 125 ! 120 126 ! OA This is (I hope) a temporary solution for the problem that may … … 159 165 wsbio3(ji,jj,jk) = MIN( wsbio3(ji,jj,jk), zwsmax * FLOAT( iiter1 ) ) 160 166 wsbio4(ji,jj,jk) = MIN( wsbio4(ji,jj,jk), zwsmax * FLOAT( iiter2 ) ) 167 #if defined key_ligand 168 wsfep(ji,jj,jk) = MIN( wsfep(ji,jj,jk), zwsmax * FLOAT( iiter1 ) ) 169 #endif 161 170 ENDIF 162 171 END DO … … 172 181 sinksil (:,:,:) = 0.e0 173 182 sinkfer2(:,:,:) = 0.e0 183 #if defined key_ligand 184 sinkfep(:,:,:) = 0.e0 185 #endif 174 186 175 187 ! Compute the sedimentation term using p4zsink2 for all the sinking particles … … 178 190 CALL p4z_sink2( wsbio3, sinking , jppoc, iiter1 ) 179 191 CALL p4z_sink2( wsbio3, sinkfer , jpsfe, iiter1 ) 192 #if defined key_ligand 193 CALL p4z_sink2( wsfep , sinkfep , jpfep, iiter1 ) 194 #endif 180 195 END DO 181 196 … … 186 201 CALL p4z_sink2( wscal , sinkcal , jpcal, iiter2 ) 187 202 END DO 188 189 ! Exchange between organic matter compartments due to coagulation/disaggregation190 ! ---------------------------------------------------191 DO jk = 1, jpkm1192 DO jj = 1, jpj193 DO ji = 1, jpi194 !195 zstep = xstep196 # if defined key_degrad197 zstep = zstep * facvol(ji,jj,jk)198 # endif199 zfact = zstep * xdiss(ji,jj,jk)200 ! Part I : Coagulation dependent on turbulence201 zagg1 = 25.9 * zfact * trb(ji,jj,jk,jppoc) * trb(ji,jj,jk,jppoc)202 zagg2 = 4452. * zfact * trb(ji,jj,jk,jppoc) * trb(ji,jj,jk,jpgoc)203 204 ! Part II : Differential settling205 206 ! Aggregation of small into large particles207 zagg3 = 47.1 * zstep * trb(ji,jj,jk,jppoc) * trb(ji,jj,jk,jpgoc)208 zagg4 = 3.3 * zstep * trb(ji,jj,jk,jppoc) * trb(ji,jj,jk,jppoc)209 210 zagg = zagg1 + zagg2 + zagg3 + zagg4211 zaggfe = zagg * trb(ji,jj,jk,jpsfe) / ( trb(ji,jj,jk,jppoc) + rtrn )212 213 ! Aggregation of DOC to POC :214 ! 1st term is shear aggregation of DOC-DOC215 ! 2nd term is shear aggregation of DOC-POC216 ! 3rd term is differential settling of DOC-POC217 zaggdoc = ( ( 0.369 * 0.3 * trb(ji,jj,jk,jpdoc) + 102.4 * trb(ji,jj,jk,jppoc) ) * zfact &218 & + 2.4 * zstep * trb(ji,jj,jk,jppoc) ) * 0.3 * trb(ji,jj,jk,jpdoc)219 ! transfer of DOC to GOC :220 ! 1st term is shear aggregation221 ! 2nd term is differential settling222 zaggdoc2 = ( 3.53E3 * zfact + 0.1 * zstep ) * trb(ji,jj,jk,jpgoc) * 0.3 * trb(ji,jj,jk,jpdoc)223 ! tranfer of DOC to POC due to brownian motion224 zaggdoc3 = ( 5095. * trb(ji,jj,jk,jppoc) + 114. * 0.3 * trb(ji,jj,jk,jpdoc) ) *zstep * 0.3 * trb(ji,jj,jk,jpdoc)225 226 ! Update the trends227 tra(ji,jj,jk,jppoc) = tra(ji,jj,jk,jppoc) - zagg + zaggdoc + zaggdoc3228 tra(ji,jj,jk,jpgoc) = tra(ji,jj,jk,jpgoc) + zagg + zaggdoc2229 tra(ji,jj,jk,jpsfe) = tra(ji,jj,jk,jpsfe) - zaggfe230 tra(ji,jj,jk,jpbfe) = tra(ji,jj,jk,jpbfe) + zaggfe231 tra(ji,jj,jk,jpdoc) = tra(ji,jj,jk,jpdoc) - zaggdoc - zaggdoc2 - zaggdoc3232 !233 END DO234 END DO235 END DO236 237 203 238 204 ! Total carbon export per year … … 341 307 ! 342 308 INTEGER :: ji, jj, jk, jit, niter1, niter2 343 REAL(wp) :: zagg1, zagg2, zagg3, zagg4, zagg5, zfract, zaggsi, zaggsh344 REAL(wp) :: zagg , zaggdoc, zaggdoc1, znumdoc345 309 REAL(wp) :: znum , zeps, zfm, zgm, zsm 346 310 REAL(wp) :: zdiv , zdiv1, zdiv2, zdiv3, zdiv4, zdiv5 347 REAL(wp) :: zval1, zval2, zval3 , zval4311 REAL(wp) :: zval1, zval2, zval3 348 312 REAL(wp) :: zfact 349 313 INTEGER :: ik1 … … 395 359 396 360 wscal(:,:,:) = MAX( wsbio3(:,:,:), 30._wp ) 361 #if defined key_ligand 362 wsfep (:,:,:) = wfep 363 #endif 397 364 398 365 ! INITIALIZE TO ZERO ALL THE SINKING ARRAYS … … 404 371 sinkfer (:,:,:) = 0.e0 405 372 sinksil (:,:,:) = 0.e0 373 #if defined key_ligand 374 sinkfep(:,:,:) = 0.e0 375 #endif 406 376 407 377 ! Compute the sedimentation term using p4zsink2 for all the sinking particles … … 416 386 CALL p4z_sink2( wscal , sinksil , jpgsi, niter1 ) 417 387 CALL p4z_sink2( wscal , sinkcal , jpcal, niter1 ) 388 #if defined key_ligand 389 CALL p4z_sink2( wsfep, sinkfep , jpfep, iiter1 ) 390 #endif 418 391 END DO 419 392 … … 422 395 END DO 423 396 424 ! Exchange between organic matter compartments due to coagulation/disaggregation 425 ! --------------------------------------------------- 426 427 zval1 = 1. + xkr_zeta 428 zval2 = 1. + xkr_eta 429 zval3 = 3. + xkr_eta 430 zval4 = 4. + xkr_eta 431 432 DO jk = 1,jpkm1 433 DO jj = 1,jpj 434 DO ji = 1,jpi 435 IF( tmask(ji,jj,jk) /= 0.e0 ) THEN 436 437 znum = trb(ji,jj,jk,jppoc)/(trb(ji,jj,jk,jpnum)+rtrn) / xkr_massp 438 !-------------- To avoid sinking speed over 50 m/day ------- 439 znum = min(xnumm(jk),znum) 440 znum = MAX( 1.1,znum) 441 !------------------------------------------------------------ 442 zeps = ( zval1 * znum - 1.) / ( znum - 1.) 443 zdiv = MAX( 1.e-4, ABS( zeps - zval3) ) * SIGN( 1., zeps - zval3 ) 444 zdiv1 = MAX( 1.e-4, ABS( zeps - 4. ) ) * SIGN( 1., zeps - 4. ) 445 zdiv2 = zeps - 2. 446 zdiv3 = zeps - 3. 447 zdiv4 = zeps - zval2 448 zdiv5 = 2.* zeps - zval4 449 zfm = xkr_frac**( 1.- zeps ) 450 zsm = xkr_frac**xkr_eta 451 452 ! Part I : Coagulation dependant on turbulence 453 ! ---------------------------------------------- 454 455 zagg1 = 0.163 * trb(ji,jj,jk,jpnum)**2 & 456 & * 2.*( (zfm-1.)*(zfm*xkr_mass_max**3-xkr_mass_min**3) & 457 & * (zeps-1)/zdiv1 + 3.*(zfm*xkr_mass_max-xkr_mass_min) & 458 & * (zfm*xkr_mass_max**2-xkr_mass_min**2) & 459 & * (zeps-1.)**2/(zdiv2*zdiv3)) 460 zagg2 = 2*0.163*trb(ji,jj,jk,jpnum)**2*zfm* & 461 & ((xkr_mass_max**3+3.*(xkr_mass_max**2 & 462 & *xkr_mass_min*(zeps-1.)/zdiv2 & 463 & +xkr_mass_max*xkr_mass_min**2*(zeps-1.)/zdiv3) & 464 & +xkr_mass_min**3*(zeps-1)/zdiv1) & 465 & -zfm*xkr_mass_max**3*(1.+3.*((zeps-1.)/ & 466 & (zeps-2.)+(zeps-1.)/zdiv3)+(zeps-1.)/zdiv1)) 467 468 zagg3 = 0.163*trb(ji,jj,jk,jpnum)**2*zfm**2*8. * xkr_mass_max**3 469 470 ! Aggregation of small into large particles 471 ! Part II : Differential settling 472 ! ---------------------------------------------- 473 474 zagg4 = 2.*3.141*0.125*trb(ji,jj,jk,jpnum)**2* & 475 & xkr_wsbio_min*(zeps-1.)**2 & 476 & *(xkr_mass_min**2*((1.-zsm*zfm)/(zdiv3*zdiv4) & 477 & -(1.-zfm)/(zdiv*(zeps-1.)))- & 478 & ((zfm*zfm*xkr_mass_max**2*zsm-xkr_mass_min**2) & 479 & *xkr_eta)/(zdiv*zdiv3*zdiv5) ) 480 481 zagg5 = 2.*3.141*0.125*trb(ji,jj,jk,jpnum)**2 & 482 & *(zeps-1.)*zfm*xkr_wsbio_min & 483 & *(zsm*(xkr_mass_min**2-zfm*xkr_mass_max**2) & 484 & /zdiv3-(xkr_mass_min**2-zfm*zsm*xkr_mass_max**2) & 485 & /zdiv) 486 487 ! 488 ! Fractionnation by swimming organisms 489 ! ------------------------------------ 490 491 zfract = 2.*3.141*0.125*trb(ji,jj,jk,jpmes)*12./0.12/0.06**3*trb(ji,jj,jk,jpnum) & 492 & * (0.01/xkr_mass_min)**(1.-zeps)*0.1**2 & 493 & * 10000.*xstep 494 495 ! Aggregation of DOC to small particles 496 ! -------------------------------------- 497 498 zaggdoc = 0.83 * trb(ji,jj,jk,jpdoc) * xstep * xdiss(ji,jj,jk) * trb(ji,jj,jk,jpdoc) & 499 & + 0.005 * 231. * trb(ji,jj,jk,jpdoc) * xstep * trb(ji,jj,jk,jpdoc) 500 zaggdoc1 = 271. * trb(ji,jj,jk,jppoc) * xstep * xdiss(ji,jj,jk) * trb(ji,jj,jk,jpdoc) & 501 & + 0.02 * 16706. * trb(ji,jj,jk,jppoc) * xstep * trb(ji,jj,jk,jpdoc) 502 503 # if defined key_degrad 504 zagg1 = zagg1 * facvol(ji,jj,jk) 505 zagg2 = zagg2 * facvol(ji,jj,jk) 506 zagg3 = zagg3 * facvol(ji,jj,jk) 507 zagg4 = zagg4 * facvol(ji,jj,jk) 508 zagg5 = zagg5 * facvol(ji,jj,jk) 509 zaggdoc = zaggdoc * facvol(ji,jj,jk) 510 zaggdoc1 = zaggdoc1 * facvol(ji,jj,jk) 511 # endif 512 zaggsh = ( zagg1 + zagg2 + zagg3 ) * rfact2 * xdiss(ji,jj,jk) / 1000. 513 zaggsi = ( zagg4 + zagg5 ) * xstep / 10. 514 zagg = 0.5 * xkr_stick * ( zaggsh + zaggsi ) 515 ! 516 znumdoc = trb(ji,jj,jk,jpnum) / ( trb(ji,jj,jk,jppoc) + rtrn ) 517 tra(ji,jj,jk,jppoc) = tra(ji,jj,jk,jppoc) + zaggdoc + zaggdoc1 518 tra(ji,jj,jk,jpnum) = tra(ji,jj,jk,jpnum) + zfract + zaggdoc / xkr_massp - zagg 519 tra(ji,jj,jk,jpdoc) = tra(ji,jj,jk,jpdoc) - zaggdoc - zaggdoc1 520 521 ENDIF 522 END DO 523 END DO 524 END DO 525 526 ! Total primary production per year 527 t_oce_co2_exp = t_oce_co2_exp + glob_sum( ( sinking(:,:,ik100) * e1e2t(:,:) * tmask(:,:,1) ) 397 IF( iom_use( "tcexp" ) .OR. ( ln_check_mass .AND. kt == nitend .AND. knt == nrdttrc ) ) & 398 & t_oce_co2_exp = glob_sum( sinking(:,:,ik100) * e1e2t(:,:) * tmask(:,:,1) ) 528 399 ! 529 400 IF( lk_iomput ) THEN … … 755 626 ik100 = 10 ! last level where depth less than 100 m 756 627 DO jk = jpkm1, 1, -1 757 IF( gdept_1d(jk) > 100. ) ik sed= jk - 1628 IF( gdept_1d(jk) > 100. ) ik100 = jk - 1 758 629 END DO 759 630 IF (lwp) WRITE(numout,*) … … 897 768 & sinkfer2(jpi,jpj,jpk) , & 898 769 #endif 770 #if defined key_ligand 771 & wsfep(jpi,jpj,jpk) , sinkfep(jpi,jpj,jpk) , & 772 #endif 899 773 & sinkfer(jpi,jpj,jpk) , STAT=p4z_sink_alloc ) 900 774 ! -
branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsms.F90
r5546 r6453 19 19 USE p4zbio ! Biological model 20 20 USE p4zche ! Chemical model 21 USE p4zlys ! Calcite saturation22 21 USE p4zflx ! Gas exchange 23 22 USE p4zsbc ! External source of nutrients … … 85 84 CALL p4z_che ! initialize the chemical constants 86 85 ! 87 IF( .NOT. ln_rsttr ) THEN ; CALL p4z_ph_ini ! set PH at kt=nit00086 IF( .NOT. ln_rsttr ) THEN ; CALL ahini_for_at(hi) ! set PH at kt=nit000 88 87 ELSE ; CALL p4z_rst( nittrc000, 'READ' ) !* read or initialize all required fields 89 88 ENDIF … … 133 132 ! 134 133 CALL p4z_bio( kt, jnt ) ! Biology 134 CALL p4z_flx( kt, jnt ) ! Compute surface fluxes 135 135 CALL p4z_sed( kt, jnt ) ! Sedimentation 136 CALL p4z_lys( kt, jnt ) ! Compute CaCO3 saturation137 CALL p4z_flx( kt, jnt ) ! Compute surface fluxes138 136 ! 139 137 xnegtr(:,:,:) = 1.e0 … … 216 214 !! natkriest ("key_kriest") 217 215 !!---------------------------------------------------------------------- 218 NAMELIST/nampisbio/ nrdttrc, wsbio, xkmort, ferat3, wsbio2, niter1max, niter2max 216 #if defined key_ligand 217 NAMELIST/nampisbio/ nrdttrc, wsbio, xkmort, ferat3, wsbio2, wsbio2max, wsbio2scale, & 218 & niter1max, niter2max, wfep, ldocp, ldocz, lthet 219 #else 220 NAMELIST/nampisbio/ nrdttrc, wsbio, xkmort, ferat3, wsbio2, wsbio2max, wsbio2scale, niter1max, niter2max 221 #endif 219 222 #if defined key_kriest 220 223 NAMELIST/nampiskrp/ xkr_eta, xkr_zeta, xkr_ncontent, xkr_mass_min, xkr_mass_max … … 241 244 WRITE(numout,*) ' Fe/C in zooplankton ferat3 =', ferat3 242 245 WRITE(numout,*) ' Big particles sinking speed wsbio2 =', wsbio2 246 WRITE(numout,*) ' Big particles maximum sinking speed wsbio2max =', wsbio2max 247 WRITE(numout,*) ' Big particles sinking speed length scale wsbio2scale =', wsbio2scale 243 248 WRITE(numout,*) ' Maximum number of iterations for POC niter1max =', niter1max 244 249 WRITE(numout,*) ' Maximum number of iterations for GOC niter2max =', niter2max 250 #if defined key_ligand 251 WRITE(numout,*) ' FeP sinking speed wfep =', wfep 252 WRITE(numout,*) ' Phyto ligand production per unit doc ldocp =', ldocp 253 WRITE(numout,*) ' Zoo ligand production per unit doc ldocz =', ldocz 254 WRITE(numout,*) ' Proportional loss of ligands due to Fe uptake lthet =', lthet 255 #endif 245 256 ENDIF 246 257 … … 310 321 END SUBROUTINE p4z_sms_init 311 322 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 323 SUBROUTINE p4z_rst( kt, cdrw ) 341 324 !!--------------------------------------------------------------------- … … 350 333 INTEGER , INTENT(in) :: kt ! ocean time-step 351 334 CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag 352 !353 INTEGER :: ji, jj, jk354 REAL(wp) :: zcaralk, zbicarb, zco3355 REAL(wp) :: ztmas, ztmas1356 335 !!--------------------------------------------------------------------- 357 336 … … 365 344 CALL iom_get( numrtr, jpdom_autoglo, 'PH' , hi(:,:,:) ) 366 345 ELSE 367 ! hi(:,:,:) = 1.e-9 368 CALL p4z_ph_ini 346 CALL ahini_for_at(hi) 369 347 ENDIF 370 348 CALL iom_get( numrtr, jpdom_autoglo, 'Silicalim', xksi(:,:) ) … … 550 528 & + trn(:,:,:,jpbfe) & 551 529 #endif 530 #if defined key_ligand 531 & + trn(:,:,:,jpfep) & 532 #endif 552 533 & + trn(:,:,:,jpsfe) & 553 534 & + trn(:,:,:,jpzoo) * ferat3 & -
branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/par_pisces.F90
r5385 r6453 18 18 !!--------------------------------------------------------------------- 19 19 LOGICAL, PUBLIC, PARAMETER :: lk_pisces = .TRUE. !: PISCES flag 20 LOGICAL, PUBLIC, PARAMETER :: lk_p4z = .FALSE.!: p4z flag20 INTEGER, PUBLIC, PARAMETER :: nn_p4z = 1 !: p4z flag 21 21 INTEGER, PUBLIC, PARAMETER :: jp_pisces = 6 !: number of passive tracers 22 22 INTEGER, PUBLIC, PARAMETER :: jp_pisces_2d = 19 !: additional 2d output … … 36 36 INTEGER, PUBLIC, PARAMETER :: jpkbm1 = jpkb - 1 !: first vertical layers where biology is active 37 37 38 #elif defined key_pisces && defined key_kriest38 #elif defined key_pisces 39 39 !!--------------------------------------------------------------------- 40 40 !! 'key_pisces' & 'key_kriest' PISCES bio-model + ??? 41 41 !!--------------------------------------------------------------------- 42 LOGICAL, PUBLIC, PARAMETER :: lk_pisces = .TRUE. !: PISCES flag 43 LOGICAL, PUBLIC, PARAMETER :: lk_p4z = .TRUE. !: p4z flag 44 LOGICAL, PUBLIC, PARAMETER :: lk_kriest = .TRUE. !: Kriest flag 45 INTEGER, PUBLIC, PARAMETER :: jp_pisces = 23 !: number of passive tracers 46 INTEGER, PUBLIC, PARAMETER :: jp_pisces_2d = 13 !: additional 2d output 47 INTEGER, PUBLIC, PARAMETER :: jp_pisces_3d = 18 !: additional 3d output 48 INTEGER, PUBLIC, PARAMETER :: jp_pisces_trd = 1 !: number of sms trends for PISCES 42 LOGICAL, PUBLIC, PARAMETER :: lk_pisces = .TRUE. !: PISCES flag 43 INTEGER, PUBLIC, PARAMETER :: nn_p4z = 2 !: p4z flag 44 # if defined key_kriest 45 LOGICAL, PUBLIC, PARAMETER :: lk_kriest = .TRUE. !: Kriest flag 46 INTEGER, PUBLIC, PARAMETER :: jp_kriest = 1 47 INTEGER, PUBLIC, PARAMETER :: jp_kriest_diag = 7 48 # else 49 LOGICAL, PUBLIC, PARAMETER :: lk_kriest = .FALSE. !: Kriest flag 50 INTEGER, PUBLIC, PARAMETER :: jp_kriest = 2 51 INTEGER, PUBLIC, PARAMETER :: jp_kriest_diag = 0 52 # endif 53 # if defined key_ligand 54 LOGICAL, PUBLIC, PARAMETER :: lk_ligand = .TRUE. !: Kriest flag 55 INTEGER, PUBLIC, PARAMETER :: jp_ligand = 2 56 INTEGER, PUBLIC, PARAMETER :: jp_ligand_diag = 0 57 # else 58 LOGICAL, PUBLIC, PARAMETER :: lk_ligand = .FALSE. !: Kriest flag 59 INTEGER, PUBLIC, PARAMETER :: jp_ligand = 0 60 INTEGER, PUBLIC, PARAMETER :: jp_ligand_diag = 0 61 # endif 62 INTEGER, PUBLIC, PARAMETER :: jp_pisces = 22 + jp_kriest + jp_ligand !: number of passive tracers 63 INTEGER, PUBLIC, PARAMETER :: jp_pisces_2d = 13 !: additional 2d output 64 INTEGER, PUBLIC, PARAMETER :: jp_pisces_3d = 13 + jp_kriest_diag + jp_ligand_diag !: additional 3d output 65 INTEGER, PUBLIC, PARAMETER :: jp_pisces_trd = 1 !: number of sms trends for PISCES 49 66 50 67 ! assign an index in trc arrays for each LOBSTER prognostic variables … … 63 80 INTEGER, PUBLIC, PARAMETER :: jpdia = 11 !: Diatoms Concentration 64 81 INTEGER, PUBLIC, PARAMETER :: jpmes = 12 !: Mesozooplankton Concentration 65 INTEGER, PUBLIC, PARAMETER :: jpdsi = 13 !: DiatomsSilicate Concentration82 INTEGER, PUBLIC, PARAMETER :: jpdsi = 13 !: (big) Silicate Concentration 66 83 INTEGER, PUBLIC, PARAMETER :: jpfer = 14 !: Iron Concentration 67 INTEGER, PUBLIC, PARAMETER :: jpnum = 15 !: Big iron particles Concentration 68 INTEGER, PUBLIC, PARAMETER :: jpsfe = 16 !: number of particulate organic phosphate concentration 69 INTEGER, PUBLIC, PARAMETER :: jpdfe = 17 !: Diatoms iron Concentration 70 INTEGER, PUBLIC, PARAMETER :: jpgsi = 18 !: (big) Silicate Concentration 71 INTEGER, PUBLIC, PARAMETER :: jpnfe = 19 !: Nano iron Concentration 72 INTEGER, PUBLIC, PARAMETER :: jpnch = 20 !: Nano Chlorophyll Concentration 73 INTEGER, PUBLIC, PARAMETER :: jpdch = 21 !: Diatoms Chlorophyll Concentration 74 INTEGER, PUBLIC, PARAMETER :: jpno3 = 22 !: Nitrates Concentration 75 INTEGER, PUBLIC, PARAMETER :: jpnh4 = 23 !: Ammonium Concentration 76 77 #elif defined key_pisces 78 !!--------------------------------------------------------------------- 79 !! 'key_pisces' : standard PISCES bio-model 84 INTEGER, PUBLIC, PARAMETER :: jpsfe = 15 !: number of particulate organic phosphate concentration 85 INTEGER, PUBLIC, PARAMETER :: jpdfe = 16 !: Diatoms iron Concentration 86 INTEGER, PUBLIC, PARAMETER :: jpgsi = 17 !: Diatoms Silicate Concentration 87 INTEGER, PUBLIC, PARAMETER :: jpnfe = 18 !: Nano iron Concentration 88 INTEGER, PUBLIC, PARAMETER :: jpnch = 19 !: Nano Chlorophyll Concentration 89 INTEGER, PUBLIC, PARAMETER :: jpdch = 20 !: Diatoms Chlorophyll Concentration 90 INTEGER, PUBLIC, PARAMETER :: jpno3 = 21 !: Nitrates Concentration 91 INTEGER, PUBLIC, PARAMETER :: jpnh4 = 22 !: Ammonium Concentration 92 # if defined key_kriest 93 INTEGER, PUBLIC, PARAMETER :: jpnum = 23 !: Number of particles in the aggregates 94 # else 95 INTEGER, PUBLIC, PARAMETER :: jpgoc = 23 !: Big OM particles Concentration 96 INTEGER, PUBLIC, PARAMETER :: jpbfe = 24 !: Big iron particles Concentration 97 # endif 98 # if defined key_ligand 99 INTEGER, PUBLIC, PARAMETER :: jplgw = 22 + jp_kriest + 1 !: Weak Ligands 100 INTEGER, PUBLIC, PARAMETER :: jpfep = 22 + jp_kriest + 2 !: Fe nanoparticle 101 # endif 102 103 #elif defined key_pisces_quota 104 !!--------------------------------------------------------------------- 105 !! 'key_pisces' & 'key_kriest' PISCES bio-model + ??? 80 106 !!--------------------------------------------------------------------- 81 107 LOGICAL, PUBLIC, PARAMETER :: lk_pisces = .TRUE. !: PISCES flag 82 LOGICAL, PUBLIC, PARAMETER :: lk_p4z = .TRUE. !: p4z flag 108 INTEGER, PUBLIC, PARAMETER :: nn_p4z = 3 !: p4z flag 109 # if defined key_kriest 110 LOGICAL, PUBLIC, PARAMETER :: lk_kriest = .TRUE. !: Kriest flag 111 INTEGER, PUBLIC, PARAMETER :: jp_kriest = 1 112 INTEGER, PUBLIC, PARAMETER :: jp_kriest_diag = 7 113 # else 83 114 LOGICAL, PUBLIC, PARAMETER :: lk_kriest = .FALSE. !: Kriest flag 84 INTEGER, PUBLIC, PARAMETER :: jp_pisces = 24 !: number of PISCES passive tracers 85 INTEGER, PUBLIC, PARAMETER :: jp_pisces_2d = 13 !: additional 2d output 86 INTEGER, PUBLIC, PARAMETER :: jp_pisces_3d = 11 !: additional 3d output 87 INTEGER, PUBLIC, PARAMETER :: jp_pisces_trd = 1 !: number of sms trends for PISCES 115 INTEGER, PUBLIC, PARAMETER :: jp_kriest = 4 116 INTEGER, PUBLIC, PARAMETER :: jp_kriest_diag = 7 117 # endif 118 # if defined key_ligand 119 LOGICAL, PUBLIC, PARAMETER :: lk_ligand = .TRUE. !: Kriest flag 120 INTEGER, PUBLIC, PARAMETER :: jp_ligand = 2 121 INTEGER, PUBLIC, PARAMETER :: jp_ligand_diag = 0 122 # else 123 LOGICAL, PUBLIC, PARAMETER :: lk_ligand = .FALSE. !: Kriest flag 124 INTEGER, PUBLIC, PARAMETER :: jp_ligand = 0 125 INTEGER, PUBLIC, PARAMETER :: jp_ligand_diag = 0 126 # endif 127 INTEGER, PUBLIC, PARAMETER :: jp_pisces = 35 + jp_kriest + jp_ligand !: number of passive tracers 128 INTEGER, PUBLIC, PARAMETER :: jp_pisces_2d = 13 !: additional 2d output 129 INTEGER, PUBLIC, PARAMETER :: jp_pisces_3d = 13 + jp_kriest_diag + jp_ligand_diag !: additional 3d output 130 INTEGER, PUBLIC, PARAMETER :: jp_pisces_trd = 1 !: number of sms trends for PISCES 88 131 89 132 ! assign an index in trc arrays for each LOBSTER prognostic variables … … 102 145 INTEGER, PUBLIC, PARAMETER :: jpdia = 11 !: Diatoms Concentration 103 146 INTEGER, PUBLIC, PARAMETER :: jpmes = 12 !: Mesozooplankton Concentration 104 INTEGER, PUBLIC, PARAMETER :: jpdsi = 13 !: DiatomsSilicate Concentration147 INTEGER, PUBLIC, PARAMETER :: jpdsi = 13 !: (big) Silicate Concentration 105 148 INTEGER, PUBLIC, PARAMETER :: jpfer = 14 !: Iron Concentration 106 INTEGER, PUBLIC, PARAMETER :: jpbfe = 15 !: Big iron particles Concentration 107 INTEGER, PUBLIC, PARAMETER :: jpgoc = 16 !: big particulate organic phosphate concentration 108 INTEGER, PUBLIC, PARAMETER :: jpsfe = 17 !: Small iron particles Concentration 109 INTEGER, PUBLIC, PARAMETER :: jpdfe = 18 !: Diatoms iron Concentration 110 INTEGER, PUBLIC, PARAMETER :: jpgsi = 19 !: (big) Silicate Concentration 111 INTEGER, PUBLIC, PARAMETER :: jpnfe = 20 !: Nano iron Concentration 112 INTEGER, PUBLIC, PARAMETER :: jpnch = 21 !: Nano Chlorophyll Concentration 113 INTEGER, PUBLIC, PARAMETER :: jpdch = 22 !: Diatoms Chlorophyll Concentration 114 INTEGER, PUBLIC, PARAMETER :: jpno3 = 23 !: Nitrates Concentration 115 INTEGER, PUBLIC, PARAMETER :: jpnh4 = 24 !: Ammonium Concentration 149 INTEGER, PUBLIC, PARAMETER :: jpsfe = 15 !: number of particulate organic phosphate concentration 150 INTEGER, PUBLIC, PARAMETER :: jpdfe = 16 !: Diatoms iron Concentration 151 INTEGER, PUBLIC, PARAMETER :: jpgsi = 17 !: Diatoms Silicate Concentration 152 INTEGER, PUBLIC, PARAMETER :: jpnfe = 18 !: Nano iron Concentration 153 INTEGER, PUBLIC, PARAMETER :: jpnch = 19 !: Nano Chlorophyll Concentration 154 INTEGER, PUBLIC, PARAMETER :: jpdch = 20 !: Diatoms Chlorophyll Concentration 155 INTEGER, PUBLIC, PARAMETER :: jpno3 = 21 !: Nitrates Concentration 156 INTEGER, PUBLIC, PARAMETER :: jpnh4 = 22 !: Ammonium Concentration 157 INTEGER, PUBLIC, PARAMETER :: jpdon = 23 !: dissolved organic nitrogen concentration 158 INTEGER, PUBLIC, PARAMETER :: jpdop = 24 !: dissolved organic phosphorus concentration 159 INTEGER, PUBLIC, PARAMETER :: jppon = 25 !: small particulate organic nitrogen concentration 160 INTEGER, PUBLIC, PARAMETER :: jppop = 26 !: small particulate organic phosphorus concentration 161 INTEGER, PUBLIC, PARAMETER :: jpnph = 27 !: small particulate organic phosphorus concentration 162 INTEGER, PUBLIC, PARAMETER :: jppph = 28 !: small particulate organic phosphorus concentration 163 INTEGER, PUBLIC, PARAMETER :: jpndi = 29 !: small particulate organic phosphorus concentration 164 INTEGER, PUBLIC, PARAMETER :: jppdi = 30 !: small particulate organic phosphorus concentration 165 INTEGER, PUBLIC, PARAMETER :: jppic = 31 !: small particulate organic phosphorus concentration 166 INTEGER, PUBLIC, PARAMETER :: jpnpi = 32 !: small particulate organic phosphorus concentration 167 INTEGER, PUBLIC, PARAMETER :: jpppi = 33 !: small particulate organic phosphorus concentration 168 INTEGER, PUBLIC, PARAMETER :: jppfe = 34 !: small particulate organic phosphorus concentration 169 INTEGER, PUBLIC, PARAMETER :: jppch = 35 !: small particulate organic phosphorus concentration 170 # if defined key_kriest 171 INTEGER, PUBLIC, PARAMETER :: jpnum = 36 !: Number of particles in the aggregates 172 # else 173 INTEGER, PUBLIC, PARAMETER :: jpgoc = 36 !: Big carbon particles Concentration 174 INTEGER, PUBLIC, PARAMETER :: jpgon = 37 !: Big nitrogen particles Concentration 175 INTEGER, PUBLIC, PARAMETER :: jpgop = 38 !: Big phosphorus particles Concentration 176 INTEGER, PUBLIC, PARAMETER :: jpbfe = 39 !: Big iron particles Concentration 177 # endif 178 # if defined key_ligand 179 INTEGER, PUBLIC, PARAMETER :: jplgw = 35 + jp_kriest + 1 !: Weak Ligands 180 INTEGER, PUBLIC, PARAMETER :: jpfep = 35 + jp_kriest + 2 !: Fe nanoparticle 181 # endif 116 182 117 183 #else … … 120 186 !!--------------------------------------------------------------------- 121 187 LOGICAL, PUBLIC, PARAMETER :: lk_pisces = .FALSE. !: PISCES flag 122 LOGICAL, PUBLIC, PARAMETER :: lk_p4z = .FALSE.!: p4z flag188 INTEGER, PUBLIC, PARAMETER :: nn_p4z = 0 !: p4z flag 123 189 INTEGER, PUBLIC, PARAMETER :: jp_pisces = 0 !: No CFC tracers 124 190 INTEGER, PUBLIC, PARAMETER :: jp_pisces_2d = 0 !: No CFC additional 2d output arrays -
branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/sms_pisces.F90
r5385 r6453 6 6 !! History : 1.0 ! 2000-02 (O. Aumont) original code 7 7 !! 3.2 ! 2009-04 (C. Ethe & NEMO team) style 8 !! 3.6 ! 2015-05 (O. Aumont) PISCES quota 8 9 !!---------------------------------------------------------------------- 9 #if defined key_pisces || defined key_pisces_reduced 10 #if defined key_pisces || defined key_pisces_reduced || defined key_pisces_quota 10 11 !!---------------------------------------------------------------------- 11 !! 'key_pisces ' PISCES model12 !! 'key_pisces*' PISCES model 12 13 !!---------------------------------------------------------------------- 13 14 USE par_oce … … 22 23 23 24 !!* Biological fluxes for light : variables shared by pisces & lobster 24 INTEGER , ALLOCATABLE, SAVE, DIMENSION(:,:) :: neln !: number of T-levels + 1 in the euphotic layer25 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: heup !: euphotic layer depth26 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,: ,:) :: etot !: par (photosynthetic available radiation)27 !28 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: xksi !: LOBSTER : zooplakton closure29 ! !: PISCES : silicon dependant half saturation25 INTEGER , ALLOCATABLE, SAVE, DIMENSION(:,:) :: neln !: number of T-levels + 1 in the euphotic layer 26 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: heup !: euphotic layer depth 27 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: heup_01 !: Absolute euphotic layer depth 28 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: etot !: par (photosynthetic available radiation) 29 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: xksi !: LOBSTER : zooplakton closure 30 ! !: PISCES : silicon dependant half saturation 30 31 31 #if defined key_pisces 32 #if defined key_pisces || defined key_pisces_quota 32 33 !!* Time variables 33 34 INTEGER :: nrdttrc !: ??? … … 38 39 REAL(wp) :: ryyss !: number of seconds per year 39 40 REAL(wp) :: r1_ryyss !: inverse number of seconds per year 40 41 41 42 42 !!* Biological parameters … … 49 49 REAL(wp) :: o2nit !: ??? 50 50 REAL(wp) :: wsbio, wsbio2 !: ??? 51 REAL(wp) :: wsbio2max !: ??? 52 REAL(wp) :: wsbio2scale !: ??? 51 53 REAL(wp) :: xkmort !: ??? 52 54 REAL(wp) :: ferat3 !: ??? 55 #if defined key_ligand 56 REAL(wp) :: wfep !: ??? 57 REAL(wp) :: ldocp !: ??? 58 REAL(wp) :: ldocz !: ??? 59 REAL(wp) :: lthet !: ??? 60 #endif 61 53 62 54 63 !!* diagnostic parameters … … 68 77 !!* Biological fluxes for primary production 69 78 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: xksimax !: ??? 70 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: xnanono3 !: ???71 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: xdiatno3 !: ???72 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: xnanonh4 !: ???73 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: xdiatnh4 !: ???74 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: xnanopo4 !: ???75 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: xdiatpo4 !: ???76 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: xlimphy !: ???77 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: xlimdia !: ???78 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: concdfe !: ???79 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: concnfe !: ???80 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: xlimnfe !: ???81 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: xlimdfe !: ???82 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: xlimsi !: ???83 79 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: biron !: bioavailable fraction of iron 84 80 #if defined key_ligand 81 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: plig !: proportion of iron organically complexed 82 #endif 85 83 86 84 !!* SMS for the organic matter 87 85 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: xfracal !: ?? 88 86 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: nitrfac !: ?? 89 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: xlimbac !: ?? 90 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: xlimbacl !: ?? 87 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: orem !: ?? 91 88 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: xdiss !: ?? 92 89 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: prodcal !: Calcite production 90 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: prodpoc !: Calcite production 91 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: conspoc !: Calcite production 92 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: prodgoc !: Calcite production 93 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: consgoc !: Calcite production 94 93 95 94 96 !!* Variable for chemistry of the CO2 cycle 95 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: akb3 !: ???96 97 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: ak13 !: ??? 97 98 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: ak23 !: ??? 98 99 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: aksp !: ??? 99 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: akw3 !: ???100 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: borat !: ???101 100 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: hi !: ??? 102 101 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: excess !: ??? … … 116 115 117 116 #endif 117 118 # if defined key_pisces_quota 119 !!* Biological parameters 120 REAL(wp) :: no3rat3 !: ??? 121 REAL(wp) :: po4rat3 !: ??? 122 123 !!* SMS for the organic matter 124 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: sizen !: size of diatoms 125 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: sizep !: size of diatoms 126 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: sized !: size of diatoms 127 128 #endif 118 129 !!---------------------------------------------------------------------- 119 130 !! NEMO/TOP 3.3 , NEMO Consortium (2010) … … 128 139 !!---------------------------------------------------------------------- 129 140 USE lib_mpp , ONLY: ctl_warn 141 #if defined key_pisces_quota 142 INTEGER :: ierr(6) ! Local variables 143 #else 130 144 INTEGER :: ierr(5) ! Local variables 145 #endif 131 146 !!---------------------------------------------------------------------- 132 147 ierr(:) = 0 133 148 !* Biological fluxes for light : shared variables for pisces & lobster 134 ALLOCATE( etot(jpi,jpj,jpk), neln(jpi,jpj), heup(jpi,jpj), xksi(jpi,jpj), STAT=ierr(1) ) 149 ALLOCATE( etot(jpi,jpj,jpk), neln(jpi,jpj), heup(jpi,jpj), & 150 & heup_01(jpi,jpj), xksi(jpi,jpj), STAT=ierr(1) ) 135 151 ! 136 #if defined key_pisces 152 #if defined key_pisces || defined key_pisces_quota 137 153 !* Biological fluxes for primary production 138 154 ALLOCATE( xksimax(jpi,jpj) , biron (jpi,jpj,jpk), & 139 & xnanono3(jpi,jpj,jpk), xdiatno3(jpi,jpj,jpk), & 140 & xnanonh4(jpi,jpj,jpk), xdiatnh4(jpi,jpj,jpk), & 141 & xnanopo4(jpi,jpj,jpk), xdiatpo4(jpi,jpj,jpk), & 142 & xlimphy (jpi,jpj,jpk), xlimdia (jpi,jpj,jpk), & 143 & xlimnfe (jpi,jpj,jpk), xlimdfe (jpi,jpj,jpk), & 144 & xlimsi (jpi,jpj,jpk), concdfe (jpi,jpj,jpk), & 145 & concnfe (jpi,jpj,jpk), STAT=ierr(2) ) 155 #if defined key_ligand 156 & plig(jpi,jpj,jpk) , & 157 #endif 158 & STAT=ierr(2) ) 146 159 ! 147 160 !* SMS for the organic matter 148 ALLOCATE( xfracal (jpi,jpj,jpk), nitrfac(jpi,jpj,jpk), & 149 & xlimbac (jpi,jpj,jpk), xdiss (jpi,jpj,jpk), & 150 & xlimbacl(jpi,jpj,jpk), prodcal(jpi,jpj,jpk), STAT=ierr(3) ) 161 ALLOCATE( xfracal (jpi,jpj,jpk), nitrfac (jpi,jpj,jpk), & 162 & orem (jpi,jpj,jpk), & 163 & prodcal(jpi,jpj,jpk), xdiss (jpi,jpj,jpk), & 164 & prodpoc(jpi,jpj,jpk) , conspoc(jpi,jpj,jpk), & 165 & prodgoc(jpi,jpj,jpk) , consgoc(jpi,jpj,jpk), STAT=ierr(3) ) 151 166 152 167 !* Variable for chemistry of the CO2 cycle 153 ALLOCATE( ak b3(jpi,jpj,jpk) , ak13 (jpi,jpj,jpk) ,&168 ALLOCATE( ak13 (jpi,jpj,jpk) , & 154 169 & ak23(jpi,jpj,jpk) , aksp (jpi,jpj,jpk) , & 155 & akw3(jpi,jpj,jpk) , borat (jpi,jpj,jpk) , &156 170 & hi (jpi,jpj,jpk) , excess(jpi,jpj,jpk) , STAT=ierr(4) ) 157 171 ! … … 160 174 ! 161 175 #endif 176 #if defined key_pisces_quota 177 !* Size of phytoplankton cells 178 ALLOCATE( sizen (jpi,jpj,jpk), sizep (jpi,jpj,jpk), & 179 & sized (jpi,jpj,jpk), STAT=ierr(6) ) 180 ! 181 #endif 182 162 183 ! 163 184 sms_pisces_alloc = MAXVAL( ierr ) -
branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/trcice_pisces.F90
r5724 r6453 37 37 !!---------------------------------------------------------------------- 38 38 39 IF( lk_p4z ) THEN ; CALL p4z_ice_ini ! PISCES40 ELSE ; CALL p2z_ice_ini ! LOBSTER39 IF( nn_p4z == 1 ) THEN ; CALL p2z_ice_ini ! PISCES 40 ELSE ; CALL p4z_ice_ini ! LOBSTER 41 41 ENDIF 42 42 … … 128 128 zpisc(jpno3,1) = 5.79e-6_wp / rno3 129 129 zpisc(jpnh4,1) = 3.22e-7_wp / rno3 130 #if defined key_pisces_quota 131 zpisc(jppic,1) = 9.57e-8_wp 132 zpisc(jpnpi,1) = 9.57e-8_wp 133 zpisc(jpppi,1) = 9.57e-8_wp 134 zpisc(jppfe,1) = 1.76e-11_wp 135 zpisc(jppch,1) = 1.67e-7_wp 136 zpisc(jpnph,1) = 9.57e-8_wp 137 zpisc(jppph,1) = 9.57e-8_wp 138 zpisc(jpndi,1) = 4.24e-7_wp 139 zpisc(jppdi,1) = 4.24e-7_wp 140 zpisc(jppon,1) = 9.57e-8_wp 141 zpisc(jppop,1) = 9.57e-8_wp 142 zpisc(jpdon,1) = 2.04e-5_wp 143 zpisc(jpdop,1) = 2.04e-5_wp 144 zpisc(jpgon,1) = 5.23e-8_wp 145 zpisc(jpgop,1) = 5.23e-8_wp 146 #endif 130 147 131 148 !--- Arctic specificities (dissolved inorganic & DOM) … … 158 175 zpisc(jpno3,2) = 3.51e-06_wp / rno3 159 176 zpisc(jpnh4,2) = 6.15e-08_wp / rno3 177 #if defined key_pisces_quota 178 zpisc(jppic,1) = 5.25e-7_wp 179 zpisc(jpnpi,1) = 5.25e-7_wp 180 zpisc(jpppi,1) = 5.25e-7_wp 181 zpisc(jppfe,1) = 1.75e-11_wp 182 zpisc(jppch,1) = 1.46e-07_wp 183 zpisc(jpnph,1) = 5.25e-7_wp 184 zpisc(jppph,1) = 5.25e-7_wp 185 zpisc(jpndi,1) = 7.75e-7_wp 186 zpisc(jppdi,1) = 7.75e-7_wp 187 zpisc(jppon,1) = 4.05e-7_wp 188 zpisc(jppop,1) = 4.05e-7_wp 189 zpisc(jpdon,1) = 6.00e-6_wp 190 zpisc(jpdop,1) = 6.00e-6_wp 191 zpisc(jpgon,1) = 2.84e-8_wp 192 zpisc(jpgop,1) = 2.84e-8_wp 193 #endif 194 160 195 161 196 !--- Antarctic specificities (dissolved inorganic & DOM) … … 188 223 zpisc(jpno3,3) = 2.64e-5_wp / rno3 189 224 zpisc(jpnh4,3) = 3.39e-7_wp / rno3 225 #if defined key_pisces_quota 226 zpisc(jppic,1) = 8.10e-7_wp 227 zpisc(jpnpi,1) = 8.10e-7_wp 228 zpisc(jpppi,1) = 8.10e-7_wp 229 zpisc(jppfe,1) = 1.48e-11_wp 230 zpisc(jppch,1) = 2.02e-7_wp 231 zpisc(jpnph,1) = 9.57e-8_wp 232 zpisc(jppph,1) = 9.57e-8_wp 233 zpisc(jpndi,1) = 5.77e-7_wp 234 zpisc(jppdi,1) = 5.77e-7_wp 235 zpisc(jppon,1) = 1.13e-6_wp 236 zpisc(jppop,1) = 1.13e-6_wp 237 zpisc(jpdon,1) = 7.02e-6_wp 238 zpisc(jpdop,1) = 7.02e-6_wp 239 zpisc(jpgon,1) = 2.89e-8_wp 240 zpisc(jpgop,1) = 2.89e-8_wp 241 #endif 242 190 243 191 244 !--- Baltic Sea particular case for ORCA configurations … … 218 271 zpisc(jpno3,4) = 5.36e-5_wp / rno3 219 272 zpisc(jpnh4,4) = 7.18e-7_wp / rno3 273 #if defined key_pisces_quota 274 zpisc(jppic,1) = 6.64e-7_wp 275 zpisc(jpnpi,1) = 6.64e-7_wp 276 zpisc(jpppi,1) = 6.64e-7_wp 277 zpisc(jppfe,1) = 3.89e-11_wp 278 zpisc(jppch,1) = 1.17e-7_wp 279 zpisc(jpnph,1) = 6.64e-7_wp 280 zpisc(jppph,1) = 6.64e-7_wp 281 zpisc(jpndi,1) = 3.41e-7_wp 282 zpisc(jppdi,1) = 3.41e-7_wp 283 zpisc(jppon,1) = 4.84e-7_wp 284 zpisc(jppop,1) = 4.84e-7_wp 285 zpisc(jpdon,1) = 1.06e-5_wp 286 zpisc(jpdop,1) = 1.06e-5_wp 287 zpisc(jpgon,1) = 1.05e-8_wp 288 zpisc(jpgop,1) = 1.05e-8_wp 289 #endif 290 220 291 221 292 DO jn = jp_pcs0, jp_pcs1 -
branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/trcini_pisces.F90
r5385 r6453 10 10 !! 2.0 ! 2007-12 (C. Ethe, G. Madec) from trcini.pisces.h90 11 11 !! 3.5 ! 2012-05 (C. Ethe) Merge PISCES-LOBSTER 12 !!---------------------------------------------------------------------- 13 #if defined key_pisces || defined key_pisces_reduced 14 !!---------------------------------------------------------------------- 15 !! 'key_pisces' PISCES bio-model 12 !! 3.6 ! 2015-05 (O. Aumont) PISCES quota 13 !!---------------------------------------------------------------------- 14 #if defined key_pisces || defined key_pisces_reduced || defined key_pisces_quota 15 !!---------------------------------------------------------------------- 16 !! 'key_pisces*' PISCES bio-model 16 17 !!---------------------------------------------------------------------- 17 18 !! trc_ini_pisces : PISCES biochemical model initialisation … … 43 44 !!---------------------------------------------------------------------- 44 45 45 IF( lk_p4z ) THEN ; CALL p4z_ini ! PISCES 46 ELSE ; CALL p2z_ini ! LOBSTER 47 ENDIF 46 SELECT CASE ( nn_p4z ) 47 ! 48 CASE(1) ; CALL p2z_ini ! LOBSTER 49 CASE(2) ; CALL p4z_ini ! PISCES 50 CASE(3) ; CALL p5z_ini ! PISCES QUOTA 51 52 END SELECT 48 53 49 54 END SUBROUTINE trc_ini_pisces 55 56 SUBROUTINE p5z_ini 57 !!---------------------------------------------------------------------- 58 !! *** ROUTINE p5z_ini *** 59 !! 60 !! ** Purpose : Initialisation of the PISCES biochemical model 61 !! with variable stoichiometry 62 !!---------------------------------------------------------------------- 63 #if defined key_pisces_quota 64 ! 65 USE p5zsms ! Main P4Z routine 66 USE p4zche ! Chemical model 67 USE p5zsink ! vertical flux of particulate matter due to sinking 68 USE p4zopt ! optical model 69 USE p4zsbc ! Boundary conditions 70 USE p4zfechem ! Iron chemistry 71 USE p5zrem ! Remineralisation of organic matter 72 USE p5zpoc ! Remineralization of organic particles 73 USE p4zligand ! Remineralization of organic ligands 74 USE p4zflx ! Gas exchange 75 USE p5zlim ! Co-limitations of differents nutrients 76 USE p5zprod ! Growth rate of the 2 phyto groups 77 USE p5zmicro ! Sources and sinks of microzooplankton 78 USE p5zmeso ! Sources and sinks of mesozooplankton 79 USE p5zmort ! Mortality terms for phytoplankton 80 USE p4zlys ! Calcite saturation 81 USE p5zsed ! Sedimentation & burial 82 ! 83 REAL(wp), SAVE :: sco2 = 2.312e-3_wp 84 REAL(wp), SAVE :: alka0 = 2.423e-3_wp 85 REAL(wp), SAVE :: oxyg0 = 177.6e-6_wp 86 REAL(wp), SAVE :: po4 = 2.174e-6_wp 87 REAL(wp), SAVE :: bioma0 = 1.000e-8_wp 88 REAL(wp), SAVE :: silic1 = 91.65e-6_wp 89 REAL(wp), SAVE :: no3 = 31.04e-6_wp * 7.625_wp 90 ! 91 INTEGER :: ierr 92 !!---------------------------------------------------------------------- 93 94 IF(lwp) WRITE(numout,*) 95 IF(lwp) WRITE(numout,*) ' p5z_ini : PISCES biochemical model initialisation' 96 IF(lwp) WRITE(numout,*) ' With variable stoichiometry' 97 IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~~~~' 98 99 ! Allocate PISCES arrays 100 ierr = sms_pisces_alloc() 101 ierr = ierr + p4z_che_alloc() 102 ierr = ierr + p5z_sink_alloc() 103 ierr = ierr + p4z_opt_alloc() 104 ierr = ierr + p5z_lim_alloc() 105 ierr = ierr + p5z_prod_alloc() 106 ierr = ierr + p5z_rem_alloc() 107 ierr = ierr + p4z_flx_alloc() 108 ierr = ierr + p5z_sed_alloc() 109 ! 110 IF( lk_mpp ) CALL mpp_sum( ierr ) 111 IF( ierr /= 0 ) CALL ctl_stop( 'STOP', 'pisces_alloc: unable to allocate PISCES arrays' ) 112 ! 113 ryyss = nyear_len(1) * rday ! number of seconds per year 114 r1_ryyss = 1. / ryyss 115 ! 116 CALL p5z_sms_init ! Maint routine 117 ! ! Time-step 118 ! Set biological ratios 119 ! --------------------- 120 rno3 = 16._wp / 122._wp 121 po4r = 1._wp / 122._wp 122 o2nit = 32._wp / 122._wp 123 rdenit = 105._wp / 16._wp 124 rdenita = 3._wp / 5._wp 125 o2ut = 133._wp / 122._wp 126 no3rat3 = no3rat3 / rno3 127 po4rat3 = po4rat3 / po4r 128 129 ! Initialization of tracer concentration in case of no restart 130 !-------------------------------------------------------------- 131 IF( .NOT. ln_rsttr ) THEN 132 133 trn(:,:,:,jpdic) = sco2 134 trn(:,:,:,jpdoc) = bioma0 135 trn(:,:,:,jpdon) = bioma0 136 trn(:,:,:,jpdop) = bioma0 137 trn(:,:,:,jptal) = alka0 138 trn(:,:,:,jpoxy) = oxyg0 139 trn(:,:,:,jpcal) = bioma0 140 trn(:,:,:,jppo4) = po4 / po4r 141 trn(:,:,:,jppoc) = bioma0 142 trn(:,:,:,jppon) = bioma0 143 trn(:,:,:,jppop) = bioma0 144 # if ! defined key_kriest 145 trn(:,:,:,jpgoc) = bioma0 146 trn(:,:,:,jpgon) = bioma0 147 trn(:,:,:,jpgop) = bioma0 148 trn(:,:,:,jpbfe) = bioma0 * 5.e-6 149 # else 150 trn(:,:,:,jpnum) = bioma0 / ( 6. * xkr_massp ) 151 # endif 152 trn(:,:,:,jpsil) = silic1 153 trn(:,:,:,jpdsi) = bioma0 * 0.15 154 trn(:,:,:,jpgsi) = bioma0 * 5.e-6 155 trn(:,:,:,jpphy) = bioma0 156 trn(:,:,:,jpnph) = bioma0 157 trn(:,:,:,jppph) = bioma0 158 trn(:,:,:,jppic) = bioma0 159 trn(:,:,:,jpnpi) = bioma0 160 trn(:,:,:,jpppi) = bioma0 161 trn(:,:,:,jpdia) = bioma0 162 trn(:,:,:,jpndi) = bioma0 163 trn(:,:,:,jppdi) = bioma0 164 trn(:,:,:,jpzoo) = bioma0 165 trn(:,:,:,jpmes) = bioma0 166 trn(:,:,:,jpfer) = 0.6E-9 167 trn(:,:,:,jpsfe) = bioma0 * 5.e-6 168 trn(:,:,:,jppfe) = bioma0 * 5.e-6 169 trn(:,:,:,jpdfe) = bioma0 * 5.e-6 170 trn(:,:,:,jpnfe) = bioma0 * 5.e-6 171 trn(:,:,:,jpnch) = bioma0 * 12. / 55. 172 trn(:,:,:,jppch) = bioma0 * 12. / 55. 173 trn(:,:,:,jpdch) = bioma0 * 12. / 55. 174 trn(:,:,:,jpno3) = no3 175 trn(:,:,:,jpnh4) = bioma0 176 #if defined key_ligand 177 trn(:,:,:,jplgw) = 0.6E-9 178 trn(:,:,:,jpfep) = 0. * 5.e-6 179 #endif 180 181 ! initialize the half saturation constant for silicate 182 ! ---------------------------------------------------- 183 xksi(:,:) = 2.e-6 184 xksimax(:,:) = xksi(:,:) 185 END IF 186 187 CALL p5z_sink_init ! vertical flux of particulate organic matter 188 CALL p4z_opt_init ! Optic: PAR in the water column 189 CALL p5z_lim_init ! co-limitations by the various nutrients 190 CALL p5z_prod_init ! phytoplankton growth rate over the global ocean. 191 CALL p4z_sbc_init ! boundary conditions 192 CALL p4z_fechem_init ! Iron chemistry 193 CALL p5z_rem_init ! remineralisation 194 CALL p5z_poc_init ! remineralisation of organic particles 195 #if defined key_ligand 196 CALL p4z_ligand_init ! remineralisation of organic ligands 197 #endif 198 CALL p5z_mort_init ! phytoplankton mortality 199 CALL p5z_micro_init ! microzooplankton 200 CALL p5z_meso_init ! mesozooplankton 201 CALL p4z_lys_init ! calcite saturation 202 CALL p4z_flx_init ! gas exchange 203 204 ndayflxtr = 0 205 206 IF(lwp) WRITE(numout,*) 207 IF(lwp) WRITE(numout,*) 'Initialization of PISCES tracers done' 208 IF(lwp) WRITE(numout,*) 209 #endif 210 ! 211 END SUBROUTINE p5z_ini 50 212 51 213 SUBROUTINE p4z_ini … … 64 226 USE p4zfechem ! Iron chemistry 65 227 USE p4zrem ! Remineralisation of organic matter 228 USE p4zpoc ! Remineralisation of organic particles 229 USE p4zligand ! Remineralization of organic ligands 66 230 USE p4zflx ! Gas exchange 67 231 USE p4zlim ! Co-limitations of differents nutrients … … 74 238 ! 75 239 REAL(wp), SAVE :: sco2 = 2.312e-3_wp 76 REAL(wp), SAVE :: alka0 = 2.426e-3_wp 77 REAL(wp), SAVE :: oxyg0 = 177.6e-6_wp 78 REAL(wp), SAVE :: po4 = 2.165e-6_wp 79 REAL(wp), SAVE :: bioma0 = 1.000e-8_wp 80 REAL(wp), SAVE :: silic1 = 91.51e-6_wp 81 REAL(wp), SAVE :: no3 = 30.9e-6_wp * 7.625_wp 82 ! 83 INTEGER :: ji, jj, jk, ierr 84 REAL(wp) :: zcaralk, zbicarb, zco3 85 REAL(wp) :: ztmas, ztmas1 240 REAL(wp), SAVE :: alka0 = 2.423e-3_wp 241 REAL(wp), SAVE :: oxyg0 = 177.6e-6_wp 242 REAL(wp), SAVE :: po4 = 2.174e-6_wp 243 REAL(wp), SAVE :: bioma0 = 1.000e-8_wp 244 REAL(wp), SAVE :: silic1 = 91.65e-6_wp 245 REAL(wp), SAVE :: no3 = 31.04e-6_wp * 7.625_wp 246 ! 247 INTEGER :: ierr 86 248 !!---------------------------------------------------------------------- 87 249 … … 91 253 92 254 ! Allocate PISCES arrays 93 ierr = sms_pisces_alloc() 255 ierr = sms_pisces_alloc() 94 256 ierr = ierr + p4z_che_alloc() 95 257 ierr = ierr + p4z_sink_alloc() 96 258 ierr = ierr + p4z_opt_alloc() 259 ierr = ierr + p4z_lim_alloc() 97 260 ierr = ierr + p4z_prod_alloc() 98 261 ierr = ierr + p4z_rem_alloc() … … 109 272 CALL p4z_sms_init ! Maint routine 110 273 ! ! Time-step 111 112 274 ! Set biological ratios 113 275 ! --------------------- … … 121 283 ! Initialization of tracer concentration in case of no restart 122 284 !-------------------------------------------------------------- 123 IF( .NOT. ln_rsttr ) THEN 124 285 IF( .NOT. ln_rsttr ) THEN 286 125 287 trn(:,:,:,jpdic) = sco2 126 288 trn(:,:,:,jpdoc) = bioma0 … … 151 313 trn(:,:,:,jpno3) = no3 152 314 trn(:,:,:,jpnh4) = bioma0 315 #if defined key_ligand 316 trn(:,:,:,jplgw) = 0.6E-9 317 trn(:,:,:,jpfep) = 0. * 5.e-6 318 #endif 153 319 154 320 ! initialize the half saturation constant for silicate … … 157 323 xksimax(:,:) = xksi(:,:) 158 324 END IF 159 160 325 161 326 CALL p4z_sink_init ! vertical flux of particulate organic matter … … 166 331 CALL p4z_fechem_init ! Iron chemistry 167 332 CALL p4z_rem_init ! remineralisation 333 CALL p4z_poc_init ! remineralization of organic particles 334 #if defined key_ligand 335 CALL p4z_ligand_init ! remineralisation of organic ligands 336 #endif 168 337 CALL p4z_mort_init ! phytoplankton mortality 169 338 CALL p4z_micro_init ! microzooplankton … … 174 343 ndayflxtr = 0 175 344 176 IF(lwp) WRITE(numout,*) 345 IF(lwp) WRITE(numout,*) 177 346 IF(lwp) WRITE(numout,*) 'Initialization of PISCES tracers done' 178 IF(lwp) WRITE(numout,*) 347 IF(lwp) WRITE(numout,*) 179 348 #endif 180 349 ! -
branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/trcnam_pisces.F90
r4990 r6453 8 8 !! 1.0 ! 2003-08 (C. Ethe) module F90 9 9 !! 2.0 ! 2007-12 (C. Ethe, G. Madec) from trcnam.pisces.h90 10 !! 3.6 ! 2015-05 (O. Aumont) PISCES quota 10 11 !!---------------------------------------------------------------------- 11 #if defined key_pisces || defined key_pisces_reduced 12 #if defined key_pisces || defined key_pisces_reduced || defined key_pisces_quota 12 13 !!---------------------------------------------------------------------- 13 !! 'key_pisces ': PISCES bio-model14 !! 'key_pisces*' : PISCES bio-model 14 15 !!---------------------------------------------------------------------- 15 16 !! trc_nam_pisces : PISCES model namelist read … … 66 67 #if defined key_pisces 67 68 IF(lwp) WRITE(numout,*) ' trc_nam_pisces : read PISCES namelist' 69 #elif defined key_pisces_quota 70 IF(lwp) WRITE(numout,*) ' trc_nam_pisces : read PISCES-QUOTA namelist' 68 71 #else 69 72 IF(lwp) WRITE(numout,*) ' trc_nam_pisces : read LOBSTER namelist' … … 123 126 #if defined key_pisces_reduced 124 127 125 IF( ( .NOT.lk_iomput .AND. ln_diabio ) .OR. lk_trdm xl_trc ) THEN128 IF( ( .NOT.lk_iomput .AND. ln_diabio ) .OR. lk_trdmld_trc ) THEN 126 129 ! 127 130 ! Namelist nampisdbi -
branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/trcsms_pisces.F90
r4147 r6453 6 6 !! History : 1.0 ! 2004-03 (O. Aumont) Original code 7 7 !! 2.0 ! 2007-12 (C. Ethe, G. Madec) F90 8 !! 3.6 ! 2015-05 (O. Aumont) PISCES quota 8 9 !!---------------------------------------------------------------------- 9 #if defined key_pisces || defined key_pisces_reduced 10 #if defined key_pisces || defined key_pisces_reduced || defined key_pisces_quota 10 11 !!---------------------------------------------------------------------- 11 !! 'key_pisces ' PISCES bio-model12 !! 'key_pisces*' PISCES bio-model 12 13 !!---------------------------------------------------------------------- 13 14 !! trcsms_pisces : Time loop of passive tracers sms 14 15 !!---------------------------------------------------------------------- 15 16 USE par_pisces 17 USE p5zsms 16 18 USE p4zsms 17 19 USE p2zsms … … 48 50 !!--------------------------------------------------------------------- 49 51 ! 50 IF( lk_p4z ) THEN ; CALL p4z_sms( kt ) ! PISCES 51 ELSE ; CALL p2z_sms( kt ) ! LOBSTER 52 ENDIF 52 SELECT CASE ( nn_p4z ) 53 54 CASE(1) ; CALL p2z_sms( kt ) ! LOBSTER 55 CASE(2) ; CALL p4z_sms( kt ) ! PISCES 56 CASE(3) ; CALL p5z_sms( kt ) ! PISCES with variable stoichiometry 57 58 END SELECT 53 59 54 60 ! -
branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/trcwri_pisces.F90
r4996 r6453 5 5 !!====================================================================== 6 6 !! History : 1.0 ! 2009-05 (C. Ethe) Original code 7 !! 3.6 ! 2015-05 (O. Aumont) PISCES quota 7 8 !!---------------------------------------------------------------------- 8 #if defined key_top && defined key_iomput && ( defined key_pisces || defined key_pisces_reduced )9 #if defined key_top && defined key_iomput && ( defined key_pisces || defined key_pisces_reduced || defined key_pisces_quota) 9 10 !!---------------------------------------------------------------------- 10 !! 'key_pisces or key_pisces_reduced' PISCES model11 !! 'key_pisces*' PISCES model 11 12 !!---------------------------------------------------------------------- 12 13 !! trc_wri_pisces : outputs of concentration fields … … 30 31 !! ** Purpose : output passive tracers fields 31 32 !!--------------------------------------------------------------------- 32 CHARACTER (len=20) :: cltra 33 REAL(wp) :: zfact 34 INTEGER :: ji, jj, jk, jn 35 REAL(wp), DIMENSION(jpi,jpj) :: zdic, zo2min, zdepo2min 33 CHARACTER (len=20) :: cltra 34 REAL(wp) :: zrfact 35 INTEGER :: jn 36 36 !!--------------------------------------------------------------------- 37 37 … … 41 41 DO jn = jp_pcs0, jp_pcs1 42 42 cltra = TRIM( ctrcnm(jn) ) ! short title for tracer 43 CALL iom_put( cltra, trn(:,:,:,jn) ) 43 IF( lk_vvl ) THEN 44 CALL iom_put( cltra, trn(:,:,:,jn) * fse3t_n(:,:,:) ) 45 ELSE 46 CALL iom_put( cltra, trn(:,:,:,jn) ) 47 ENDIF 48 CALL iom_put( cltra, trn(:,:,:,jn) * zrfact ) 44 49 END DO 45 50 #else 46 51 DO jn = jp_pcs0, jp_pcs1 47 z fact = 1.0e+648 IF( jn == jpno3 .OR. jn == jpnh4 ) z fact = rno3 * 1.0e+649 IF( jn == jppo4 ) z fact = po4r * 1.0e+652 zrfact = 1.0e+6 53 IF( jn == jpno3 .OR. jn == jpnh4 ) zrfact = rno3 * 1.0e+6 54 IF( jn == jppo4 ) zrfact = po4r * 1.0e+6 50 55 cltra = TRIM( ctrcnm(jn) ) ! short title for tracer 51 IF( iom_use( cltra ) ) CALL iom_put( cltra, trn(:,:,:,jn) * zfact ) 56 IF( lk_vvl ) THEN 57 CALL iom_put( cltra, trn(:,:,:,jn) * fse3t_n(:,:,:) * zrfact ) 58 ELSE 59 CALL iom_put( cltra, trn(:,:,:,jn) * zrfact ) 60 ENDIF 52 61 END DO 53 54 IF( iom_use( "INTDIC" ) ) THEN ! DIC content in kg/m255 zdic(:,:) = 0.56 DO jk = 1, jpkm157 zdic(:,:) = zdic(:,:) + trn(:,:,jk,jpdic) * fse3t(:,:,jk) * tmask(:,:,jk) * 12.58 ENDDO59 CALL iom_put( 'INTDIC', zdic )60 ENDIF61 !62 IF( iom_use( "O2MIN" ) .OR. iom_use ( "ZO2MIN" ) ) THEN ! Oxygen minimum concentration and depth63 zo2min (:,:) = trn(:,:,1,jpoxy) * tmask(:,:,1)64 zdepo2min(:,:) = fsdepw(:,:,1) * tmask(:,:,1)65 DO jk = 2, jpkm166 DO jj = 1, jpj67 DO ji = 1, jpi68 IF( tmask(ji,jj,jk) == 1 ) then69 IF( trn(ji,jj,jk,jpoxy) < zo2min(ji,jj) ) then70 zo2min (ji,jj) = trn(ji,jj,jk,jpoxy)71 zdepo2min(ji,jj) = fsdepw(ji,jj,jk)72 ENDIF73 ENDIF74 END DO75 END DO76 END DO77 !78 CALL iom_put('O2MIN' , zo2min ) ! oxygen minimum concentration79 CALL iom_put('ZO2MIN', zdepo2min ) ! depth of oxygen minimum concentration80 !81 ENDIF82 62 #endif 83 63 !
Note: See TracChangeset
for help on using the changeset viewer.