| 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 | !!---------------------------------------------------------------------- |
|---|
| 24 | !! NEMO/OPA 3.3 , NEMO Consortium (2010) |
|---|
| 25 | !! $Id$ |
|---|
| 26 | !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) |
|---|
| 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 |
|---|
| 47 | REAL(wp), POINTER, DIMENSION(:) :: zco2, zyrco2 ! temporary arrays for swap |
|---|
| 48 | ! |
|---|
| 49 | !!---------------------------------------------------------------------- |
|---|
| 50 | ! |
|---|
| 51 | IF( nn_timing == 1 ) CALL timing_start('trc_atm_c14_ini') |
|---|
| 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 |
|---|
| 80 | CALL wrk_alloc( nrecco2,zco2) |
|---|
| 81 | CALL wrk_alloc( nrecco2,zyrco2) |
|---|
| 82 | zco2(:)=spco2(:) |
|---|
| 83 | zyrco2(:)=tyrco2(:) |
|---|
| 84 | ! Set CO2 times on AD time scale & swap records : In CO2 file : youngest first |
|---|
| 85 | DO jn = 1, nrecco2 |
|---|
| 86 | izco2=nrecco2-jn+1 |
|---|
| 87 | spco2(izco2)=zco2(jn) |
|---|
| 88 | tyrco2(izco2)=1950._wp-zyrco2(jn) ! BP to AD dates |
|---|
| 89 | END DO |
|---|
| 90 | CALL wrk_dealloc(nrecco2,zco2) |
|---|
| 91 | CALL wrk_dealloc(nrecco2,zyrco2) |
|---|
| 92 | ENDIF |
|---|
| 93 | ! |
|---|
| 94 | ! ! Transient atmospheric forcing: Bomb C14 & Paleo C14 : open file |
|---|
| 95 | ! |
|---|
| 96 | clfile = TRIM( cfilec14 ) |
|---|
| 97 | IF (lwp) WRITE(numout,*) 'Read C-14 atmospheric concentrations file ',clfile |
|---|
| 98 | CALL ctl_opn( inum2, clfile, 'OLD', 'FORMATTED', 'SEQUENTIAL', -1, numout, .FALSE. ) |
|---|
| 99 | REWIND(inum2) |
|---|
| 100 | ! |
|---|
| 101 | ! Bomb C14: 3 zones for atm C14 ! |
|---|
| 102 | IF(kc14typ == 1) THEN ! Transient atmospheric forcing: Bomb C14 |
|---|
| 103 | ! |
|---|
| 104 | READ(inum2,*) nrecc14,incom |
|---|
| 105 | DO jn = 1, incom ! Skip over descriptor lines |
|---|
| 106 | READ(inum2,'(1x)') |
|---|
| 107 | END DO |
|---|
| 108 | ALLOCATE( bomb(nrecc14,nc14zon), tyrc14(nrecc14) , STAT=ierr2 ) |
|---|
| 109 | IF( ierr2 /= 0 ) CALL ctl_stop( 'STOP', 'trc_atm_c14_ini: unable to allocate c14 arrays' ) |
|---|
| 110 | ! get bomb c14 data |
|---|
| 111 | DO jn = 1, nrecc14 |
|---|
| 112 | READ(inum2,*) tyrc14(jn), bomb(jn,1), bomb(jn,2), bomb(jn,3) |
|---|
| 113 | END DO |
|---|
| 114 | CLOSE(inum2) |
|---|
| 115 | |
|---|
| 116 | ! Linear interpolation of the C-14 source fonction |
|---|
| 117 | ! in linear latitude bands (20N,40N) and (20S,40S) |
|---|
| 118 | !------------------------------------------------------ |
|---|
| 119 | ALLOCATE( fareaz (jpi,jpj ,nc14zon) , STAT=ierr3 ) |
|---|
| 120 | IF( ierr3 /= 0 ) CALL ctl_stop( 'STOP', 'trc_atm_c14_ini: unable to allocate fareaz' ) |
|---|
| 121 | ! |
|---|
| 122 | DO jj = 1 , jpj ! from C14b package |
|---|
| 123 | DO ji = 1 , jpi |
|---|
| 124 | IF( gphit(ji,jj) >= yn40 ) THEN |
|---|
| 125 | fareaz(ji,jj,1) = 0. |
|---|
| 126 | fareaz(ji,jj,2) = 0. |
|---|
| 127 | fareaz(ji,jj,3) = 1. |
|---|
| 128 | ELSE IF( gphit(ji,jj ) <= ys40) THEN |
|---|
| 129 | fareaz(ji,jj,1) = 1. |
|---|
| 130 | fareaz(ji,jj,2) = 0. |
|---|
| 131 | fareaz(ji,jj,3) = 0. |
|---|
| 132 | ELSE IF( gphit(ji,jj) >= yn20 ) THEN |
|---|
| 133 | fareaz(ji,jj,1) = 0. |
|---|
| 134 | fareaz(ji,jj,2) = 2. * ( 1. - gphit(ji,jj) / yn40 ) |
|---|
| 135 | fareaz(ji,jj,3) = 2. * gphit(ji,jj) / yn40 - 1. |
|---|
| 136 | ELSE IF( gphit(ji,jj) <= ys20 ) THEN |
|---|
| 137 | fareaz(ji,jj,1) = 2. * gphit(ji,jj) / ys40 - 1. |
|---|
| 138 | fareaz(ji,jj,2) = 2. * ( 1. - gphit(ji,jj) / ys40 ) |
|---|
| 139 | fareaz(ji,jj,3) = 0. |
|---|
| 140 | ELSE |
|---|
| 141 | fareaz(ji,jj,1) = 0. |
|---|
| 142 | fareaz(ji,jj,2) = 1. |
|---|
| 143 | fareaz(ji,jj,3) = 0. |
|---|
| 144 | ENDIF |
|---|
| 145 | END DO |
|---|
| 146 | END DO |
|---|
| 147 | ! |
|---|
| 148 | ENDIF |
|---|
| 149 | ! |
|---|
| 150 | ! Paleo C14: 1 zone for atm C14 ! |
|---|
| 151 | IF(kc14typ == 2) THEN ! Transient atmospheric forcing: Paleo C14 |
|---|
| 152 | ! |
|---|
| 153 | READ(inum2,*) nrecc14,incom |
|---|
| 154 | DO jn = 1, incom ! Skip over descriptor lines |
|---|
| 155 | READ(inum2,'(1x)') |
|---|
| 156 | END DO |
|---|
| 157 | ALLOCATE( atmc14(nrecc14), tyrc14(nrecc14) , STAT=ierr2 ) |
|---|
| 158 | IF( ierr2 /= 0 ) CALL ctl_stop( 'STOP', 'trc_atm_c14_ini: unable to allocate c14 arrays' ) |
|---|
| 159 | ! get past c14 data |
|---|
| 160 | DO jn = 1, nrecc14 |
|---|
| 161 | READ(inum2,*) iyear,incom,incom,atmc14(jn) |
|---|
| 162 | tyrc14(jn)=1950._wp-REAL(iyear,wp) ! BP to AD dates |
|---|
| 163 | END DO |
|---|
| 164 | CLOSE(inum2) |
|---|
| 165 | ! |
|---|
| 166 | ENDIF |
|---|
| 167 | ! |
|---|
| 168 | ! Note on dates: |
|---|
| 169 | ! In files dates have dimension yr; either AD or BP; if BP date is changed into AD here |
|---|
| 170 | ! When dealing with dates previous to 0. AD one needs to set tyrc14_beg to the actual starting year |
|---|
| 171 | ! Do not forget to appropriately set nn_date0 and nn_rstctl in namelist |
|---|
| 172 | ! AND nn_rsttr in namelist_top if offline run |
|---|
| 173 | ! All details are given in NEMO-C14.pdf report |
|---|
| 174 | ! |
|---|
| 175 | tyrc14_now=nyear ! actual initial yr - Bomb |
|---|
| 176 | if(kc14typ == 2) tyrc14_now=nyear+tyrc14_beg-1 ! actual initial yr - Paleo |
|---|
| 177 | ! ! we suppose we start on tyrc14_now/01/01 @ 0h |
|---|
| 178 | m1_c14= 1 |
|---|
| 179 | m1_co2= 1 |
|---|
| 180 | DO jn = 1,nrecco2 |
|---|
| 181 | IF ( tyrc14_now >= tyrco2(jn) ) m1_co2 = jn ! index of first co2 record to consider |
|---|
| 182 | END DO |
|---|
| 183 | DO jn = 1,nrecc14 |
|---|
| 184 | IF ( tyrc14_now >= tyrc14(jn) ) m1_c14 = jn ! index of first c14 record to consider |
|---|
| 185 | END DO |
|---|
| 186 | IF (lwp) WRITE(numout,*) 'Initial yr for experiment', tyrc14_now |
|---|
| 187 | IF (lwp) WRITE(numout,*) ' CO2 & C14 start years:', tyrco2(m1_co2),tyrc14(m1_c14) |
|---|
| 188 | ! |
|---|
| 189 | m2_c14= m1_c14 |
|---|
| 190 | m2_co2= m1_co2 |
|---|
| 191 | ! |
|---|
| 192 | ENDIF |
|---|
| 193 | ! |
|---|
| 194 | IF( nn_timing == 1 ) CALL timing_stop('trc_atm_c14_ini') |
|---|
| 195 | ! |
|---|
| 196 | END SUBROUTINE trc_atm_c14_ini |
|---|
| 197 | |
|---|
| 198 | |
|---|
| 199 | SUBROUTINE trc_atm_c14( kt, co2sbc, c14sbc ) |
|---|
| 200 | !!---------------------------------------------------------------------- |
|---|
| 201 | !! *** ROUTINE trc_flx *** |
|---|
| 202 | !! |
|---|
| 203 | !! ** Purpose : provides sbc for co2 & c14 at kt |
|---|
| 204 | !! |
|---|
| 205 | !! ** Method : read files |
|---|
| 206 | !! |
|---|
| 207 | !! ** Action : atmospheric values interpolated at time-step kt |
|---|
| 208 | !!---------------------------------------------------------------------- |
|---|
| 209 | INTEGER , INTENT(in ) :: kt ! ocean time-step |
|---|
| 210 | REAL(wp), DIMENSION(:,:), INTENT( out) :: c14sbc ! atm c14 ratio |
|---|
| 211 | REAL(wp), INTENT( out) :: co2sbc ! atm co2 p |
|---|
| 212 | INTEGER :: jz ! dummy loop indice |
|---|
| 213 | REAL(wp) :: zdint,zint ! work |
|---|
| 214 | REAL(wp), DIMENSION(nc14zon) :: zonbc14 ! work |
|---|
| 215 | ! |
|---|
| 216 | !!---------------------------------------------------------------------- |
|---|
| 217 | ! |
|---|
| 218 | IF( nn_timing == 1 ) CALL timing_start('trc_atm_c14') |
|---|
| 219 | ! |
|---|
| 220 | IF( kc14typ == 0) THEN |
|---|
| 221 | co2sbc=pco2at |
|---|
| 222 | c14sbc(:,:)=rc14at |
|---|
| 223 | ENDIF |
|---|
| 224 | ! |
|---|
| 225 | IF(kc14typ >= 1) THEN ! Transient C14 & CO2 |
|---|
| 226 | ! |
|---|
| 227 | tyrc14_now = tyrc14_now + ( rdt / ( rday * nyear_len(1)) ) ! current time step in yr relative to tyrc14_beg |
|---|
| 228 | ! |
|---|
| 229 | ! CO2 -------------------------------------------------------- |
|---|
| 230 | ! |
|---|
| 231 | ! time interpolation of CO2 concentrations ! if out of record keeps first/last value |
|---|
| 232 | IF( tyrc14_now > tyrco2(m2_co2) ) THEN ! next interval |
|---|
| 233 | m1_co2 = m2_co2 |
|---|
| 234 | m2_co2 = MIN ( m2_co2 + 1 , nrecco2 ) |
|---|
| 235 | ENDIF |
|---|
| 236 | ! |
|---|
| 237 | zdint = tyrco2(m2_co2) - tyrco2(m1_co2) |
|---|
| 238 | co2sbc = spco2(m2_co2) ! if out of record keeps first/last value |
|---|
| 239 | zint = 0._wp |
|---|
| 240 | IF ( zdint > 0._wp ) THEN ! if within record interpolate: |
|---|
| 241 | zint = ( tyrco2(m2_co2) - tyrc14_now ) / zdint |
|---|
| 242 | co2sbc = spco2(m2_co2) + zint * ( spco2(m1_co2) - spco2(m2_co2) ) |
|---|
| 243 | ENDIF |
|---|
| 244 | ! |
|---|
| 245 | IF( lwp .AND. kt == nitend ) THEN |
|---|
| 246 | WRITE(numout, '(3(A,F12.4))') 't1/tn/t2:',tyrco2(m1_co2),'/', tyrc14_now,'/',tyrco2(m2_co2) |
|---|
| 247 | WRITE(numout, *) 'CO2:',spco2(m1_co2),co2sbc ,spco2(m2_co2) |
|---|
| 248 | ENDIF |
|---|
| 249 | ! |
|---|
| 250 | ! C14 -------------------------------------------------------- |
|---|
| 251 | ! |
|---|
| 252 | ! time interpolation of C14 concentrations |
|---|
| 253 | IF ( tyrc14_now > tyrc14(m2_c14) ) THEN ! next interval |
|---|
| 254 | m1_c14 = m2_c14 |
|---|
| 255 | m2_c14 = MIN ( m2_c14 + 1 , nrecc14 ) |
|---|
| 256 | ENDIF |
|---|
| 257 | zdint = tyrc14(m2_c14) - tyrc14(m1_c14) |
|---|
| 258 | zint=0._wp |
|---|
| 259 | IF ( zdint > 0._wp ) zint = ( tyrc14(m2_c14) - tyrc14_now ) / zdint ! if within record |
|---|
| 260 | IF( lwp .AND. kt == nitend ) & |
|---|
| 261 | & WRITE(numout,'(3(A,F12.4))') 't1/tn/t2:',tyrc14(m1_c14),'/', tyrc14_now,'/',tyrc14(m2_c14) |
|---|
| 262 | ! |
|---|
| 263 | ! ------- Bomb C14 ------------------------------------------ |
|---|
| 264 | ! |
|---|
| 265 | IF( kc14typ == 1) THEN |
|---|
| 266 | ! ! time interpolation |
|---|
| 267 | zonbc14(:) = bomb(m2_c14,:) ! if out of record keeps first/last value |
|---|
| 268 | ! ! if within record interpolate: |
|---|
| 269 | IF ( zdint > 0._wp ) zonbc14(:) = bomb(m2_c14,:) + zint * ( bomb(m1_c14,:) - bomb(m2_c14,:) ) |
|---|
| 270 | ! |
|---|
| 271 | IF(lwp .AND. kt == nitend ) & |
|---|
| 272 | & WRITE(numout, *) 'C14:',bomb(m1_c14,1),zonbc14(1),bomb(m2_c14,1) |
|---|
| 273 | ! Transform DeltaC14 --> C14 ratio |
|---|
| 274 | zonbc14(:) = 1._wp + zonbc14(:)/1.d03 |
|---|
| 275 | ! |
|---|
| 276 | ! For each (i,j)-box, with information from the fractional area |
|---|
| 277 | ! (zonmean), computes area-weighted mean to give the atmospheric C-14 |
|---|
| 278 | ! ---------------------------------------------------------------- |
|---|
| 279 | c14sbc(:,:) = zonbc14(1) * fareaz(:,:,1) & |
|---|
| 280 | & + zonbc14(2) * fareaz(:,:,2) & |
|---|
| 281 | & + zonbc14(3) * fareaz(:,:,3) |
|---|
| 282 | ENDIF |
|---|
| 283 | ! |
|---|
| 284 | ! ------- Paleo C14 ----------------------------------------- |
|---|
| 285 | ! |
|---|
| 286 | IF( kc14typ == 2 ) THEN |
|---|
| 287 | ! ! time interpolation |
|---|
| 288 | zonbc14(1) = atmc14(m2_c14) ! if out of record keeps first/last value |
|---|
| 289 | ! ! if within record interpolate: |
|---|
| 290 | IF ( zdint > 0._wp ) zonbc14(1) = atmc14(m2_c14) + zint * ( atmc14(m1_c14) - atmc14(m2_c14) ) |
|---|
| 291 | IF(lwp .AND. kt == nitend ) & |
|---|
| 292 | & WRITE(numout, *) 'C14: ',atmc14(m1_c14),zonbc14(1),atmc14(m2_c14) |
|---|
| 293 | ! Transform DeltaC14 --> C14 ratio |
|---|
| 294 | c14sbc(:,:) = 1._wp + zonbc14(1)/1.d03 |
|---|
| 295 | ENDIF |
|---|
| 296 | ! |
|---|
| 297 | ENDIF |
|---|
| 298 | ! |
|---|
| 299 | IF( nn_timing == 1 ) CALL timing_stop('trc_atm_c14') |
|---|
| 300 | ! |
|---|
| 301 | END SUBROUTINE trc_atm_c14 |
|---|
| 302 | |
|---|
| 303 | !!====================================================================== |
|---|
| 304 | END MODULE trcatm_c14 |
|---|
