- Timestamp:
- 2016-04-08T10:10:11+02:00 (8 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
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
Note: See TracChangeset
for help on using the changeset viewer.
