- Timestamp:
- 2020-09-14T17:40:34+02:00 (4 years ago)
- Location:
- NEMO/branches/2019/dev_r11351_fldread_with_XIOS
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r11351_fldread_with_XIOS
- Property svn:externals
-
old new 3 3 ^/utils/build/mk@HEAD mk 4 4 ^/utils/tools@HEAD tools 5 ^/vendors/AGRIF/dev @HEADext/AGRIF5 ^/vendors/AGRIF/dev_r12970_AGRIF_CMEMS ext/AGRIF 6 6 ^/vendors/FCM@HEAD ext/FCM 7 7 ^/vendors/IOIPSL@HEAD ext/IOIPSL 8 9 # SETTE 10 ^/utils/CI/sette@13382 sette
-
- Property svn:externals
-
NEMO/branches/2019/dev_r11351_fldread_with_XIOS/src/TOP/PISCES/P4Z/p4zflx.F90
r10425 r13463 19 19 USE sms_pisces ! PISCES Source Minus Sink variables 20 20 USE p4zche ! Chemical model 21 USE prtctl _trc! print control for debugging21 USE prtctl ! print control for debugging 22 22 USE iom ! I/O manager 23 23 USE fldread ! read input fields … … 52 52 REAL(wp) :: xconv = 0.01_wp / 3600._wp !: coefficients for conversion 53 53 54 !! * Substitutions 55 # include "do_loop_substitute.h90" 56 # include "domzgr_substitute.h90" 54 57 !!---------------------------------------------------------------------- 55 58 !! NEMO/TOP 4.0 , NEMO Consortium (2018) … … 59 62 CONTAINS 60 63 61 SUBROUTINE p4z_flx ( kt, knt )64 SUBROUTINE p4z_flx ( kt, knt, Kbb, Kmm, Krhs ) 62 65 !!--------------------------------------------------------------------- 63 66 !! *** ROUTINE p4z_flx *** … … 71 74 !!--------------------------------------------------------------------- 72 75 INTEGER, INTENT(in) :: kt, knt ! 76 INTEGER, INTENT(in) :: Kbb, Kmm, Krhs ! time level indices 73 77 ! 74 78 INTEGER :: ji, jj, jm, iind, iindm1 … … 80 84 CHARACTER (len=25) :: charout 81 85 REAL(wp), DIMENSION(jpi,jpj) :: zkgco2, zkgo2, zh2co3, zoflx, zpco2atm 82 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zw2d83 86 !!--------------------------------------------------------------------- 84 87 ! … … 107 110 IF( l_co2cpl ) satmco2(:,:) = atm_co2(:,:) 108 111 109 DO jj = 1, jpj 110 DO ji = 1, jpi 111 ! DUMMY VARIABLES FOR DIC, H+, AND BORATE 112 zfact = rhop(ji,jj,1) / 1000. + rtrn 113 zdic = trb(ji,jj,1,jpdic) 114 zph = MAX( hi(ji,jj,1), 1.e-10 ) / zfact 115 ! CALCULATE [H2CO3] 116 zh2co3(ji,jj) = zdic/(1. + ak13(ji,jj,1)/zph + ak13(ji,jj,1)*ak23(ji,jj,1)/zph**2) 117 END DO 118 END DO 112 DO_2D( 1, 1, 1, 1 ) 113 ! DUMMY VARIABLES FOR DIC, H+, AND BORATE 114 zfact = rhop(ji,jj,1) / 1000. + rtrn 115 zdic = tr(ji,jj,1,jpdic,Kbb) 116 zph = MAX( hi(ji,jj,1), 1.e-10 ) / zfact 117 ! CALCULATE [H2CO3] 118 zh2co3(ji,jj) = zdic/(1. + ak13(ji,jj,1)/zph + ak13(ji,jj,1)*ak23(ji,jj,1)/zph**2) 119 END_2D 119 120 120 121 ! -------------- … … 125 126 ! ------------------------------------------- 126 127 127 DO jj = 1, jpj 128 DO ji = 1, jpi 129 ztc = MIN( 35., tsn(ji,jj,1,jp_tem) ) 130 ztc2 = ztc * ztc 131 ztc3 = ztc * ztc2 132 ztc4 = ztc2 * ztc2 133 ! Compute the schmidt Number both O2 and CO2 134 zsch_co2 = 2116.8 - 136.25 * ztc + 4.7353 * ztc2 - 0.092307 * ztc3 + 0.0007555 * ztc4 135 zsch_o2 = 1920.4 - 135.6 * ztc + 5.2122 * ztc2 - 0.109390 * ztc3 + 0.0009377 * ztc4 136 ! wind speed 137 zws = wndm(ji,jj) * wndm(ji,jj) 138 ! Compute the piston velocity for O2 and CO2 139 zkgwan = 0.251 * zws 140 zkgwan = zkgwan * xconv * ( 1.- fr_i(ji,jj) ) * tmask(ji,jj,1) 141 ! compute gas exchange for CO2 and O2 142 zkgco2(ji,jj) = zkgwan * SQRT( 660./ zsch_co2 ) 143 zkgo2 (ji,jj) = zkgwan * SQRT( 660./ zsch_o2 ) 144 END DO 145 END DO 146 147 148 DO jj = 1, jpj 149 DO ji = 1, jpi 150 ztkel = tempis(ji,jj,1) + 273.15 151 zsal = salinprac(ji,jj,1) + ( 1.- tmask(ji,jj,1) ) * 35. 152 zvapsw = EXP(24.4543 - 67.4509*(100.0/ztkel) - 4.8489*LOG(ztkel/100) - 0.000544*zsal) 153 zpco2atm(ji,jj) = satmco2(ji,jj) * ( patm(ji,jj) - zvapsw ) 154 zxc2 = ( 1.0 - zpco2atm(ji,jj) * 1E-6 )**2 155 zfugcoeff = EXP( patm(ji,jj) * (chemc(ji,jj,2) + 2.0 * zxc2 * chemc(ji,jj,3) ) & 156 & / ( 82.05736 * ztkel )) 157 zfco2 = zpco2atm(ji,jj) * zfugcoeff 158 159 ! Compute CO2 flux for the sea and air 160 zfld = zfco2 * chemc(ji,jj,1) * zkgco2(ji,jj) ! (mol/L) * (m/s) 161 zflu = zh2co3(ji,jj) * zkgco2(ji,jj) ! (mol/L) (m/s) ? 162 oce_co2(ji,jj) = ( zfld - zflu ) * rfact2 * e1e2t(ji,jj) * tmask(ji,jj,1) * 1000. 163 ! compute the trend 164 tra(ji,jj,1,jpdic) = tra(ji,jj,1,jpdic) + ( zfld - zflu ) * rfact2 / e3t_n(ji,jj,1) * tmask(ji,jj,1) 165 166 ! Compute O2 flux 167 zfld16 = patm(ji,jj) * chemo2(ji,jj,1) * zkgo2(ji,jj) ! (mol/L) * (m/s) 168 zflu16 = trb(ji,jj,1,jpoxy) * zkgo2(ji,jj) 169 zoflx(ji,jj) = ( zfld16 - zflu16 ) * tmask(ji,jj,1) 170 tra(ji,jj,1,jpoxy) = tra(ji,jj,1,jpoxy) + zoflx(ji,jj) * rfact2 / e3t_n(ji,jj,1) 171 END DO 172 END DO 128 DO_2D( 1, 1, 1, 1 ) 129 ztc = MIN( 35., ts(ji,jj,1,jp_tem,Kmm) ) 130 ztc2 = ztc * ztc 131 ztc3 = ztc * ztc2 132 ztc4 = ztc2 * ztc2 133 ! Compute the schmidt Number both O2 and CO2 134 zsch_co2 = 2116.8 - 136.25 * ztc + 4.7353 * ztc2 - 0.092307 * ztc3 + 0.0007555 * ztc4 135 zsch_o2 = 1920.4 - 135.6 * ztc + 5.2122 * ztc2 - 0.109390 * ztc3 + 0.0009377 * ztc4 136 ! wind speed 137 zws = wndm(ji,jj) * wndm(ji,jj) 138 ! Compute the piston velocity for O2 and CO2 139 zkgwan = 0.251 * zws 140 zkgwan = zkgwan * xconv * ( 1.- fr_i(ji,jj) ) * tmask(ji,jj,1) 141 ! compute gas exchange for CO2 and O2 142 zkgco2(ji,jj) = zkgwan * SQRT( 660./ zsch_co2 ) 143 zkgo2 (ji,jj) = zkgwan * SQRT( 660./ zsch_o2 ) 144 END_2D 145 146 147 DO_2D( 1, 1, 1, 1 ) 148 ztkel = tempis(ji,jj,1) + 273.15 149 zsal = salinprac(ji,jj,1) + ( 1.- tmask(ji,jj,1) ) * 35. 150 zvapsw = EXP(24.4543 - 67.4509*(100.0/ztkel) - 4.8489*LOG(ztkel/100) - 0.000544*zsal) 151 zpco2atm(ji,jj) = satmco2(ji,jj) * ( patm(ji,jj) - zvapsw ) 152 zxc2 = ( 1.0 - zpco2atm(ji,jj) * 1E-6 )**2 153 zfugcoeff = EXP( patm(ji,jj) * (chemc(ji,jj,2) + 2.0 * zxc2 * chemc(ji,jj,3) ) & 154 & / ( 82.05736 * ztkel )) 155 zfco2 = zpco2atm(ji,jj) * zfugcoeff 156 157 ! Compute CO2 flux for the sea and air 158 zfld = zfco2 * chemc(ji,jj,1) * zkgco2(ji,jj) ! (mol/L) * (m/s) 159 zflu = zh2co3(ji,jj) * zkgco2(ji,jj) ! (mol/L) (m/s) ? 160 oce_co2(ji,jj) = ( zfld - zflu ) * tmask(ji,jj,1) 161 ! compute the trend 162 tr(ji,jj,1,jpdic,Krhs) = tr(ji,jj,1,jpdic,Krhs) + oce_co2(ji,jj) * rfact2 / e3t(ji,jj,1,Kmm) 163 164 ! Compute O2 flux 165 zfld16 = patm(ji,jj) * chemo2(ji,jj,1) * zkgo2(ji,jj) ! (mol/L) * (m/s) 166 zflu16 = tr(ji,jj,1,jpoxy,Kbb) * zkgo2(ji,jj) 167 zoflx(ji,jj) = ( zfld16 - zflu16 ) * tmask(ji,jj,1) 168 tr(ji,jj,1,jpoxy,Krhs) = tr(ji,jj,1,jpoxy,Krhs) + zoflx(ji,jj) * rfact2 / e3t(ji,jj,1,Kmm) 169 END_2D 173 170 174 171 IF( iom_use("tcflx") .OR. iom_use("tcflxcum") .OR. kt == nitrst & 175 172 & .OR. (ln_check_mass .AND. kt == nitend) ) & 176 t_oce_co2_flx = glob_sum( 'p4zflx', oce_co2(:,:) ) ! Total Flux of Carbon173 t_oce_co2_flx = glob_sum( 'p4zflx', oce_co2(:,:) * e1e2t(:,:) * 1000. ) ! Total Flux of Carbon 177 174 t_oce_co2_flx_cum = t_oce_co2_flx_cum + t_oce_co2_flx ! Cumulative Total Flux of Carbon 178 175 ! t_atm_co2_flx = glob_sum( 'p4zflx', satmco2(:,:) * e1e2t(:,:) ) ! Total atmospheric pCO2 179 176 t_atm_co2_flx = atcco2 ! Total atmospheric pCO2 180 177 181 IF( ln_ctl) THEN ! print mean trends (used for debugging)178 IF(sn_cfctl%l_prttrc) THEN ! print mean trends (used for debugging) 182 179 WRITE(charout, FMT="('flx ')") 183 CALL prt_ctl_ trc_info(charout)184 CALL prt_ctl _trc(tab4d=tra, mask=tmask, clinfo=ctrcnm)180 CALL prt_ctl_info( charout, cdcomp = 'top' ) 181 CALL prt_ctl(tab4d_1=tr(:,:,:,:,Krhs), mask1=tmask, clinfo=ctrcnm) 185 182 ENDIF 186 183 187 184 IF( lk_iomput .AND. knt == nrdttrc ) THEN 188 ALLOCATE( zw2d(jpi,jpj) ) 189 IF( iom_use( "Cflx" ) ) THEN 190 zw2d(:,:) = oce_co2(:,:) / e1e2t(:,:) * rfact2r 191 CALL iom_put( "Cflx" , zw2d ) 192 ENDIF 193 IF( iom_use( "Oflx" ) ) THEN 194 zw2d(:,:) = zoflx(:,:) * 1000 * tmask(:,:,1) 195 CALL iom_put( "Oflx" , zw2d ) 196 ENDIF 197 IF( iom_use( "Kg" ) ) THEN 198 zw2d(:,:) = zkgco2(:,:) * tmask(:,:,1) 199 CALL iom_put( "Kg" , zw2d ) 200 ENDIF 201 IF( iom_use( "Dpco2" ) ) THEN 202 zw2d(:,:) = ( zpco2atm(:,:) - zh2co3(:,:) / ( chemc(:,:,1) + rtrn ) ) * tmask(:,:,1) 203 CALL iom_put( "Dpco2" , zw2d ) 204 ENDIF 205 IF( iom_use( "Dpo2" ) ) THEN 206 zw2d(:,:) = ( atcox * patm(:,:) - atcox * trb(:,:,1,jpoxy) / ( chemo2(:,:,1) + rtrn ) ) * tmask(:,:,1) 207 CALL iom_put( "Dpo2" , zw2d ) 208 ENDIF 209 CALL iom_put( "tcflx" , t_oce_co2_flx * rfact2r ) ! molC/s 210 CALL iom_put( "tcflxcum" , t_oce_co2_flx_cum ) ! molC 211 ! 212 DEALLOCATE( zw2d ) 185 CALL iom_put( "AtmCo2" , satmco2(:,:) * tmask(:,:,1) ) ! Atmospheric CO2 concentration 186 CALL iom_put( "Cflx" , oce_co2(:,:) * 1000. ) 187 CALL iom_put( "Oflx" , zoflx(:,:) * 1000. ) 188 CALL iom_put( "Kg" , zkgco2(:,:) * tmask(:,:,1) ) 189 CALL iom_put( "Dpco2" , ( zpco2atm(:,:) - zh2co3(:,:) / ( chemc(:,:,1) + rtrn ) ) * tmask(:,:,1) ) 190 CALL iom_put( "pCO2sea" , ( zh2co3(:,:) / ( chemc(:,:,1) + rtrn ) ) * tmask(:,:,1) ) 191 CALL iom_put( "Dpo2" , ( atcox * patm(:,:) - atcox * tr(:,:,1,jpoxy,Kbb) / ( chemo2(:,:,1) + rtrn ) ) * tmask(:,:,1) ) 192 CALL iom_put( "tcflx" , t_oce_co2_flx ) ! molC/s 193 CALL iom_put( "tcflxcum", t_oce_co2_flx_cum ) ! molC 213 194 ENDIF 214 195 ! … … 239 220 ENDIF 240 221 ! 241 REWIND( numnatp_ref ) ! Namelist nampisext in reference namelist : Pisces atm. conditions242 222 READ ( numnatp_ref, nampisext, IOSTAT = ios, ERR = 901) 243 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nampisext in reference namelist', lwp ) 244 REWIND( numnatp_cfg ) ! Namelist nampisext in configuration namelist : Pisces atm. conditions 223 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nampisext in reference namelist' ) 245 224 READ ( numnatp_cfg, nampisext, IOSTAT = ios, ERR = 902 ) 246 902 IF( ios > 0 ) CALL ctl_nam ( ios , 'nampisext in configuration namelist' , lwp)225 902 IF( ios > 0 ) CALL ctl_nam ( ios , 'nampisext in configuration namelist' ) 247 226 IF(lwm) WRITE ( numonp, nampisext ) 248 227 ! … … 320 299 ENDIF 321 300 ! 322 REWIND( numnatp_ref ) ! Namelist nampisatm in reference namelist : Pisces atm. sea level pressure file323 301 READ ( numnatp_ref, nampisatm, IOSTAT = ios, ERR = 901) 324 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nampisatm in reference namelist', lwp ) 325 REWIND( numnatp_cfg ) ! Namelist nampisatm in configuration namelist : Pisces atm. sea level pressure file 302 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nampisatm in reference namelist' ) 326 303 READ ( numnatp_cfg, nampisatm, IOSTAT = ios, ERR = 902 ) 327 902 IF( ios > 0 ) CALL ctl_nam ( ios , 'nampisatm in configuration namelist' , lwp)304 902 IF( ios > 0 ) CALL ctl_nam ( ios , 'nampisatm in configuration namelist' ) 328 305 IF(lwm) WRITE ( numonp, nampisatm ) 329 306 !
Note: See TracChangeset
for help on using the changeset viewer.