[7041] | 1 | MODULE trcatm_c14 |
---|
| 2 | !!====================================================================== |
---|
| 3 | !! *** MODULE trcatm_c14 *** |
---|
| 4 | !! TOP : read and manages atmospheric values for radiocarbon model |
---|
| 5 | !!===================================================================== |
---|
| 6 | !! History: Based on trcini_c14b & trcsms_c14b : |
---|
| 7 | !! Anne Mouchet |
---|
| 8 | !!---------------------------------------------------------------------- |
---|
| 9 | !! trc_atm_c14_ini : initialize c14atm & pco2atm |
---|
| 10 | !! trc_atm_c14 : read and time interpolate c14atm & pco2atm |
---|
| 11 | !!---------------------------------------------------------------------- |
---|
| 12 | USE par_trc ! passive tracers parameters |
---|
| 13 | USE oce_trc ! shared variables between ocean and passive tracers |
---|
| 14 | USE trc ! passive tracers common variables |
---|
| 15 | USE sms_c14 ! c14 simulation type, atm default values... |
---|
| 16 | |
---|
| 17 | IMPLICIT NONE |
---|
| 18 | PRIVATE |
---|
| 19 | |
---|
| 20 | PUBLIC trc_atm_c14 ! called in trcsms_c14.F90 |
---|
| 21 | PUBLIC trc_atm_c14_ini ! called in trcini_c14.F90 |
---|
| 22 | ! |
---|
| 23 | !!---------------------------------------------------------------------- |
---|
[10068] | 24 | !! NEMO/TOP 4.0 , NEMO Consortium (2018) |
---|
[10888] | 25 | !! $Id$ |
---|
[10068] | 26 | !! Software governed by the CeCILL license (see ./LICENSE) |
---|
[7041] | 27 | !!---------------------------------------------------------------------- |
---|
| 28 | CONTAINS |
---|
| 29 | |
---|
| 30 | SUBROUTINE trc_atm_c14_ini |
---|
| 31 | !!---------------------------------------------------------------------- |
---|
| 32 | !! *** ROUTINE trc_c14_ini *** |
---|
| 33 | !! |
---|
| 34 | !! ** Purpose : initialisation of sbc for radiocarbon |
---|
| 35 | !! |
---|
| 36 | !! ** Method : |
---|
| 37 | !!---------------------------------------------------------------------- |
---|
| 38 | ! |
---|
| 39 | CHARACTER (len=20) :: clfile ! forcing file name |
---|
| 40 | INTEGER :: ji,jj,jn ! dummy loop indice |
---|
| 41 | INTEGER :: ierr1,ierr2,ierr3,izco2 ! temporary integers |
---|
| 42 | INTEGER :: inum1,inum2,incom,iyear ! temporary integers |
---|
| 43 | REAL(wp) :: ys40 = -40. ! 40 degrees south |
---|
| 44 | REAL(wp) :: ys20 = -20. ! 20 degrees south |
---|
| 45 | REAL(wp) :: yn20 = 20. ! 20 degrees north |
---|
| 46 | REAL(wp) :: yn40 = 40. ! 40 degrees north |
---|
[9125] | 47 | REAL(wp), ALLOCATABLE, DIMENSION(:) :: zco2, zyrco2 ! temporary arrays for swap |
---|
[7041] | 48 | ! |
---|
| 49 | !!---------------------------------------------------------------------- |
---|
| 50 | ! |
---|
[9124] | 51 | IF( ln_timing ) CALL timing_start('trc_atm_c14_ini') |
---|
[7041] | 52 | ! |
---|
| 53 | |
---|
| 54 | IF( lwp ) WRITE(numout,*) ' ' |
---|
| 55 | IF( lwp ) WRITE(numout,*) ' trc_atm_c14_ini : initialize atm CO2 & C14-ratio ' |
---|
| 56 | IF( lwp ) WRITE(numout,*) ' ' |
---|
| 57 | ! |
---|
| 58 | tyrc14_now = 0._wp ! initialize |
---|
| 59 | ! |
---|
| 60 | IF(kc14typ >= 1) THEN ! Transient atmospheric forcing: CO2 |
---|
| 61 | ! |
---|
| 62 | clfile = TRIM( cfileco2 ) |
---|
| 63 | IF(lwp) WRITE(numout,*) 'Read CO2 atmospheric concentrations file ',clfile |
---|
| 64 | CALL ctl_opn( inum1, clfile, 'OLD', 'FORMATTED', 'SEQUENTIAL', -1, numout, .FALSE. ) |
---|
| 65 | REWIND(inum1) |
---|
| 66 | |
---|
| 67 | READ(inum1,*) nrecco2,incom |
---|
| 68 | DO jn = 1, incom ! Skip over descriptor lines |
---|
| 69 | READ(inum1,'(1x)') |
---|
| 70 | END DO |
---|
| 71 | ALLOCATE( spco2(nrecco2), tyrco2(nrecco2) , STAT=ierr1 ) |
---|
| 72 | IF( ierr1 /= 0 ) CALL ctl_stop( 'STOP', 'trc_atm_c14_ini: unable to allocate co2 arrays' ) |
---|
| 73 | ! get CO2 data |
---|
| 74 | DO jn = 1, nrecco2 |
---|
| 75 | READ(inum1, *) tyrco2(jn), spco2(jn) |
---|
| 76 | END DO |
---|
| 77 | CLOSE(inum1) |
---|
| 78 | ! |
---|
| 79 | IF(kc14typ==2) THEN |
---|
[9125] | 80 | ALLOCATE( zco2(nrecco2), zyrco2(nrecco2) ) |
---|
[7041] | 81 | zco2(:)=spco2(:) |
---|
| 82 | zyrco2(:)=tyrco2(:) |
---|
| 83 | ! Set CO2 times on AD time scale & swap records : In CO2 file : youngest first |
---|
| 84 | DO jn = 1, nrecco2 |
---|
| 85 | izco2=nrecco2-jn+1 |
---|
| 86 | spco2(izco2)=zco2(jn) |
---|
| 87 | tyrco2(izco2)=1950._wp-zyrco2(jn) ! BP to AD dates |
---|
| 88 | END DO |
---|
[9125] | 89 | DEALLOCATE( zco2,zyrco2 ) |
---|
[7041] | 90 | ENDIF |
---|
| 91 | ! |
---|
| 92 | ! ! Transient atmospheric forcing: Bomb C14 & Paleo C14 : open file |
---|
| 93 | ! |
---|
| 94 | clfile = TRIM( cfilec14 ) |
---|
| 95 | IF (lwp) WRITE(numout,*) 'Read C-14 atmospheric concentrations file ',clfile |
---|
| 96 | CALL ctl_opn( inum2, clfile, 'OLD', 'FORMATTED', 'SEQUENTIAL', -1, numout, .FALSE. ) |
---|
| 97 | REWIND(inum2) |
---|
| 98 | ! |
---|
| 99 | ! Bomb C14: 3 zones for atm C14 ! |
---|
| 100 | IF(kc14typ == 1) THEN ! Transient atmospheric forcing: Bomb C14 |
---|
| 101 | ! |
---|
| 102 | READ(inum2,*) nrecc14,incom |
---|
| 103 | DO jn = 1, incom ! Skip over descriptor lines |
---|
| 104 | READ(inum2,'(1x)') |
---|
| 105 | END DO |
---|
| 106 | ALLOCATE( bomb(nrecc14,nc14zon), tyrc14(nrecc14) , STAT=ierr2 ) |
---|
| 107 | IF( ierr2 /= 0 ) CALL ctl_stop( 'STOP', 'trc_atm_c14_ini: unable to allocate c14 arrays' ) |
---|
| 108 | ! get bomb c14 data |
---|
| 109 | DO jn = 1, nrecc14 |
---|
| 110 | READ(inum2,*) tyrc14(jn), bomb(jn,1), bomb(jn,2), bomb(jn,3) |
---|
| 111 | END DO |
---|
| 112 | CLOSE(inum2) |
---|
| 113 | |
---|
| 114 | ! Linear interpolation of the C-14 source fonction |
---|
| 115 | ! in linear latitude bands (20N,40N) and (20S,40S) |
---|
| 116 | !------------------------------------------------------ |
---|
| 117 | ALLOCATE( fareaz (jpi,jpj ,nc14zon) , STAT=ierr3 ) |
---|
| 118 | IF( ierr3 /= 0 ) CALL ctl_stop( 'STOP', 'trc_atm_c14_ini: unable to allocate fareaz' ) |
---|
| 119 | ! |
---|
| 120 | DO jj = 1 , jpj ! from C14b package |
---|
| 121 | DO ji = 1 , jpi |
---|
| 122 | IF( gphit(ji,jj) >= yn40 ) THEN |
---|
| 123 | fareaz(ji,jj,1) = 0. |
---|
| 124 | fareaz(ji,jj,2) = 0. |
---|
| 125 | fareaz(ji,jj,3) = 1. |
---|
| 126 | ELSE IF( gphit(ji,jj ) <= ys40) THEN |
---|
| 127 | fareaz(ji,jj,1) = 1. |
---|
| 128 | fareaz(ji,jj,2) = 0. |
---|
| 129 | fareaz(ji,jj,3) = 0. |
---|
| 130 | ELSE IF( gphit(ji,jj) >= yn20 ) THEN |
---|
| 131 | fareaz(ji,jj,1) = 0. |
---|
| 132 | fareaz(ji,jj,2) = 2. * ( 1. - gphit(ji,jj) / yn40 ) |
---|
| 133 | fareaz(ji,jj,3) = 2. * gphit(ji,jj) / yn40 - 1. |
---|
| 134 | ELSE IF( gphit(ji,jj) <= ys20 ) THEN |
---|
| 135 | fareaz(ji,jj,1) = 2. * gphit(ji,jj) / ys40 - 1. |
---|
| 136 | fareaz(ji,jj,2) = 2. * ( 1. - gphit(ji,jj) / ys40 ) |
---|
| 137 | fareaz(ji,jj,3) = 0. |
---|
| 138 | ELSE |
---|
| 139 | fareaz(ji,jj,1) = 0. |
---|
| 140 | fareaz(ji,jj,2) = 1. |
---|
| 141 | fareaz(ji,jj,3) = 0. |
---|
| 142 | ENDIF |
---|
| 143 | END DO |
---|
| 144 | END DO |
---|
| 145 | ! |
---|
| 146 | ENDIF |
---|
| 147 | ! |
---|
| 148 | ! Paleo C14: 1 zone for atm C14 ! |
---|
| 149 | IF(kc14typ == 2) THEN ! Transient atmospheric forcing: Paleo C14 |
---|
| 150 | ! |
---|
| 151 | READ(inum2,*) nrecc14,incom |
---|
| 152 | DO jn = 1, incom ! Skip over descriptor lines |
---|
| 153 | READ(inum2,'(1x)') |
---|
| 154 | END DO |
---|
| 155 | ALLOCATE( atmc14(nrecc14), tyrc14(nrecc14) , STAT=ierr2 ) |
---|
| 156 | IF( ierr2 /= 0 ) CALL ctl_stop( 'STOP', 'trc_atm_c14_ini: unable to allocate c14 arrays' ) |
---|
| 157 | ! get past c14 data |
---|
| 158 | DO jn = 1, nrecc14 |
---|
| 159 | READ(inum2,*) iyear,incom,incom,atmc14(jn) |
---|
[7192] | 160 | tyrc14(jn)=1950._wp-REAL(iyear,wp) ! BP to AD dates |
---|
[7041] | 161 | END DO |
---|
| 162 | CLOSE(inum2) |
---|
| 163 | ! |
---|
| 164 | ENDIF |
---|
| 165 | ! |
---|
| 166 | ! Note on dates: |
---|
| 167 | ! In files dates have dimension yr; either AD or BP; if BP date is changed into AD here |
---|
| 168 | ! When dealing with dates previous to 0. AD one needs to set tyrc14_beg to the actual starting year |
---|
| 169 | ! Do not forget to appropriately set nn_date0 and nn_rstctl in namelist |
---|
| 170 | ! AND nn_rsttr in namelist_top if offline run |
---|
| 171 | ! All details are given in NEMO-C14.pdf report |
---|
| 172 | ! |
---|
| 173 | tyrc14_now=nyear ! actual initial yr - Bomb |
---|
| 174 | if(kc14typ == 2) tyrc14_now=nyear+tyrc14_beg-1 ! actual initial yr - Paleo |
---|
| 175 | ! ! we suppose we start on tyrc14_now/01/01 @ 0h |
---|
| 176 | m1_c14= 1 |
---|
| 177 | m1_co2= 1 |
---|
| 178 | DO jn = 1,nrecco2 |
---|
| 179 | IF ( tyrc14_now >= tyrco2(jn) ) m1_co2 = jn ! index of first co2 record to consider |
---|
| 180 | END DO |
---|
| 181 | DO jn = 1,nrecc14 |
---|
| 182 | IF ( tyrc14_now >= tyrc14(jn) ) m1_c14 = jn ! index of first c14 record to consider |
---|
| 183 | END DO |
---|
| 184 | IF (lwp) WRITE(numout,*) 'Initial yr for experiment', tyrc14_now |
---|
| 185 | IF (lwp) WRITE(numout,*) ' CO2 & C14 start years:', tyrco2(m1_co2),tyrc14(m1_c14) |
---|
| 186 | ! |
---|
| 187 | m2_c14= m1_c14 |
---|
| 188 | m2_co2= m1_co2 |
---|
| 189 | ! |
---|
| 190 | ENDIF |
---|
| 191 | ! |
---|
[9124] | 192 | IF( ln_timing ) CALL timing_stop('trc_atm_c14_ini') |
---|
[7041] | 193 | ! |
---|
| 194 | END SUBROUTINE trc_atm_c14_ini |
---|
| 195 | |
---|
| 196 | |
---|
| 197 | SUBROUTINE trc_atm_c14( kt, co2sbc, c14sbc ) |
---|
| 198 | !!---------------------------------------------------------------------- |
---|
| 199 | !! *** ROUTINE trc_flx *** |
---|
| 200 | !! |
---|
| 201 | !! ** Purpose : provides sbc for co2 & c14 at kt |
---|
| 202 | !! |
---|
| 203 | !! ** Method : read files |
---|
| 204 | !! |
---|
| 205 | !! ** Action : atmospheric values interpolated at time-step kt |
---|
| 206 | !!---------------------------------------------------------------------- |
---|
| 207 | INTEGER , INTENT(in ) :: kt ! ocean time-step |
---|
| 208 | REAL(wp), DIMENSION(:,:), INTENT( out) :: c14sbc ! atm c14 ratio |
---|
| 209 | REAL(wp), INTENT( out) :: co2sbc ! atm co2 p |
---|
| 210 | INTEGER :: jz ! dummy loop indice |
---|
| 211 | REAL(wp) :: zdint,zint ! work |
---|
| 212 | REAL(wp), DIMENSION(nc14zon) :: zonbc14 ! work |
---|
| 213 | ! |
---|
| 214 | !!---------------------------------------------------------------------- |
---|
| 215 | ! |
---|
[9124] | 216 | IF( ln_timing ) CALL timing_start('trc_atm_c14') |
---|
[7041] | 217 | ! |
---|
| 218 | IF( kc14typ == 0) THEN |
---|
| 219 | co2sbc=pco2at |
---|
| 220 | c14sbc(:,:)=rc14at |
---|
| 221 | ENDIF |
---|
| 222 | ! |
---|
| 223 | IF(kc14typ >= 1) THEN ! Transient C14 & CO2 |
---|
| 224 | ! |
---|
| 225 | tyrc14_now = tyrc14_now + ( rdt / ( rday * nyear_len(1)) ) ! current time step in yr relative to tyrc14_beg |
---|
| 226 | ! |
---|
| 227 | ! CO2 -------------------------------------------------------- |
---|
| 228 | ! |
---|
| 229 | ! time interpolation of CO2 concentrations ! if out of record keeps first/last value |
---|
| 230 | IF( tyrc14_now > tyrco2(m2_co2) ) THEN ! next interval |
---|
| 231 | m1_co2 = m2_co2 |
---|
| 232 | m2_co2 = MIN ( m2_co2 + 1 , nrecco2 ) |
---|
| 233 | ENDIF |
---|
| 234 | ! |
---|
| 235 | zdint = tyrco2(m2_co2) - tyrco2(m1_co2) |
---|
| 236 | co2sbc = spco2(m2_co2) ! if out of record keeps first/last value |
---|
| 237 | zint = 0._wp |
---|
| 238 | IF ( zdint > 0._wp ) THEN ! if within record interpolate: |
---|
| 239 | zint = ( tyrco2(m2_co2) - tyrc14_now ) / zdint |
---|
| 240 | co2sbc = spco2(m2_co2) + zint * ( spco2(m1_co2) - spco2(m2_co2) ) |
---|
| 241 | ENDIF |
---|
| 242 | ! |
---|
| 243 | IF( lwp .AND. kt == nitend ) THEN |
---|
| 244 | WRITE(numout, '(3(A,F12.4))') 't1/tn/t2:',tyrco2(m1_co2),'/', tyrc14_now,'/',tyrco2(m2_co2) |
---|
| 245 | WRITE(numout, *) 'CO2:',spco2(m1_co2),co2sbc ,spco2(m2_co2) |
---|
| 246 | ENDIF |
---|
| 247 | ! |
---|
| 248 | ! C14 -------------------------------------------------------- |
---|
| 249 | ! |
---|
| 250 | ! time interpolation of C14 concentrations |
---|
| 251 | IF ( tyrc14_now > tyrc14(m2_c14) ) THEN ! next interval |
---|
| 252 | m1_c14 = m2_c14 |
---|
| 253 | m2_c14 = MIN ( m2_c14 + 1 , nrecc14 ) |
---|
| 254 | ENDIF |
---|
| 255 | zdint = tyrc14(m2_c14) - tyrc14(m1_c14) |
---|
| 256 | zint=0._wp |
---|
| 257 | IF ( zdint > 0._wp ) zint = ( tyrc14(m2_c14) - tyrc14_now ) / zdint ! if within record |
---|
| 258 | IF( lwp .AND. kt == nitend ) & |
---|
| 259 | & WRITE(numout,'(3(A,F12.4))') 't1/tn/t2:',tyrc14(m1_c14),'/', tyrc14_now,'/',tyrc14(m2_c14) |
---|
| 260 | ! |
---|
| 261 | ! ------- Bomb C14 ------------------------------------------ |
---|
| 262 | ! |
---|
| 263 | IF( kc14typ == 1) THEN |
---|
| 264 | ! ! time interpolation |
---|
| 265 | zonbc14(:) = bomb(m2_c14,:) ! if out of record keeps first/last value |
---|
| 266 | ! ! if within record interpolate: |
---|
| 267 | IF ( zdint > 0._wp ) zonbc14(:) = bomb(m2_c14,:) + zint * ( bomb(m1_c14,:) - bomb(m2_c14,:) ) |
---|
| 268 | ! |
---|
| 269 | IF(lwp .AND. kt == nitend ) & |
---|
| 270 | & WRITE(numout, *) 'C14:',bomb(m1_c14,1),zonbc14(1),bomb(m2_c14,1) |
---|
| 271 | ! Transform DeltaC14 --> C14 ratio |
---|
| 272 | zonbc14(:) = 1._wp + zonbc14(:)/1.d03 |
---|
| 273 | ! |
---|
| 274 | ! For each (i,j)-box, with information from the fractional area |
---|
| 275 | ! (zonmean), computes area-weighted mean to give the atmospheric C-14 |
---|
| 276 | ! ---------------------------------------------------------------- |
---|
| 277 | c14sbc(:,:) = zonbc14(1) * fareaz(:,:,1) & |
---|
| 278 | & + zonbc14(2) * fareaz(:,:,2) & |
---|
| 279 | & + zonbc14(3) * fareaz(:,:,3) |
---|
| 280 | ENDIF |
---|
| 281 | ! |
---|
| 282 | ! ------- Paleo C14 ----------------------------------------- |
---|
| 283 | ! |
---|
| 284 | IF( kc14typ == 2 ) THEN |
---|
| 285 | ! ! time interpolation |
---|
| 286 | zonbc14(1) = atmc14(m2_c14) ! if out of record keeps first/last value |
---|
| 287 | ! ! if within record interpolate: |
---|
| 288 | IF ( zdint > 0._wp ) zonbc14(1) = atmc14(m2_c14) + zint * ( atmc14(m1_c14) - atmc14(m2_c14) ) |
---|
| 289 | IF(lwp .AND. kt == nitend ) & |
---|
| 290 | & WRITE(numout, *) 'C14: ',atmc14(m1_c14),zonbc14(1),atmc14(m2_c14) |
---|
| 291 | ! Transform DeltaC14 --> C14 ratio |
---|
| 292 | c14sbc(:,:) = 1._wp + zonbc14(1)/1.d03 |
---|
| 293 | ENDIF |
---|
| 294 | ! |
---|
| 295 | ENDIF |
---|
| 296 | ! |
---|
[9124] | 297 | IF( ln_timing ) CALL timing_stop('trc_atm_c14') |
---|
[7041] | 298 | ! |
---|
| 299 | END SUBROUTINE trc_atm_c14 |
---|
| 300 | |
---|
| 301 | !!====================================================================== |
---|
| 302 | END MODULE trcatm_c14 |
---|