Changeset 9190 for branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/ICB
- Timestamp:
- 2018-01-06T15:18:23+01:00 (6 years ago)
- Location:
- branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/ICB
- Files:
-
- 11 edited
-
icb_oce.F90 (modified) (1 diff)
-
icbclv.F90 (modified) (7 diffs)
-
icbdia.F90 (modified) (23 diffs)
-
icbdyn.F90 (modified) (13 diffs)
-
icbini.F90 (modified) (4 diffs)
-
icblbc.F90 (modified) (4 diffs)
-
icbrst.F90 (modified) (5 diffs)
-
icbstp.F90 (modified) (1 diff)
-
icbthm.F90 (modified) (10 diffs)
-
icbtrj.F90 (modified) (6 diffs)
-
icbutl.F90 (modified) (10 diffs)
Legend:
- Unmodified
- Added
- Removed
-
branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/ICB/icb_oce.F90
r9019 r9190 1 1 MODULE icb_oce 2 3 2 !!====================================================================== 4 3 !! *** MODULE icb_oce *** -
branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/ICB/icbclv.F90
r5215 r9190 1 1 MODULE icbclv 2 3 2 !!====================================================================== 4 3 !! *** MODULE icbclv *** … … 11 10 !! - ! 2011-05 (Alderson) budgets into separate module 12 11 !!---------------------------------------------------------------------- 12 13 13 !!---------------------------------------------------------------------- 14 14 !! icb_clv_flx : transfer input flux of ice into iceberg classes … … 45 45 !! 46 46 !!---------------------------------------------------------------------- 47 INTEGER, INTENT(in) ::kt47 INTEGER, INTENT(in) :: kt 48 48 ! 49 REAL(wp) ::zcalving_used, zdist, zfact50 INTEGER ::jn, ji, jj ! loop counters51 INTEGER ::imx ! temporary integer for max berg class52 LOGICAL, SAVE ::ll_first_call = .TRUE.49 REAL(wp) :: zcalving_used, zdist, zfact 50 INTEGER :: jn, ji, jj ! loop counters 51 INTEGER :: imx ! temporary integer for max berg class 52 LOGICAL, SAVE :: ll_first_call = .TRUE. 53 53 !!---------------------------------------------------------------------- 54 54 ! … … 70 70 DO jj = 2, jpjm1 71 71 DO ji = 2, jpim1 72 IF( berg_grid%calving(ji,jj) /= 0._wp ) & ! Need units of J73 berg_grid%stored_heat(ji,jj) = SUM( berg_grid%stored_ice(ji,jj,:) ) * & ! initial stored ice in kg74 berg_grid%calving_hflx(ji,jj) * e1e2t(ji,jj) / & ! J/s/m2 x m^2 = J/s75 berg_grid%calving(ji,jj) !/calving in kg/s72 IF( berg_grid%calving(ji,jj) /= 0._wp ) & ! Need units of J 73 berg_grid%stored_heat(ji,jj) = SUM( berg_grid%stored_ice(ji,jj,:) ) * & ! initial stored ice in kg 74 & berg_grid%calving_hflx(ji,jj) * e1e2t(ji,jj) / berg_grid%calving(ji,jj) ! J/s/m2 x m^2 75 ! ! = J/s/calving in kg/s 76 76 END DO 77 77 END DO … … 80 80 ! assume that all calving flux must be distributed even if distribution array does not sum 81 81 ! to one - this may not be what is intended, but it's what you've got 82 DO jj = 1, jpj83 DO ji = 1, jpi82 DO jj = 1, jpj 83 DO ji = 1, jpi 84 84 imx = berg_grid%maxclass(ji,jj) 85 85 zdist = SUM( rn_distribution(1:nclasses) ) / SUM( rn_distribution(1:imx) ) 86 86 DO jn = 1, imx 87 berg_grid%stored_ice(ji,jj,jn) = berg_grid%stored_ice(ji,jj,jn) +&88 87 berg_grid%stored_ice(ji,jj,jn) = berg_grid%stored_ice(ji,jj,jn) & 88 & + berg_dt * berg_grid%calving(ji,jj) * rn_distribution(jn) * zdist 89 89 END DO 90 90 END DO … … 98 98 ! 99 99 END SUBROUTINE icb_clv_flx 100 100 101 101 102 SUBROUTINE icb_clv() … … 171 172 END DO 172 173 ! 173 DO jn = 1, nclasses174 DO jn = 1, nclasses 174 175 CALL lbc_lnk( berg_grid%stored_ice(:,:,jn), 'T', 1._wp ) 175 176 END DO -
branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/ICB/icbdia.F90
r5836 r9190 1 1 MODULE icbdia 2 3 2 !!====================================================================== 4 3 !! *** MODULE icbdia *** … … 13 12 !! - ! from the right points in the code 14 13 !!---------------------------------------------------------------------- 14 15 15 !!---------------------------------------------------------------------- 16 !! icb_dia_init : initialise iceberg budgeting 16 !! icb_dia_init : initialise iceberg budgeting 17 !! icb_dia : global iceberg diagnostics 18 !! icb_dia_step : reset at the beginning of each timestep 19 !! icb_dia_put : output (via iom_put) iceberg fields 20 !! icb_dia_calve : 21 !! icb_dia_income: 22 !! icb_dia_size : 23 !! icb_dia_speed : 24 !! icb_dia_melt : 25 !! report_state : 26 !! report_consistant : 27 !! report_budget : 28 !! report_istate : 29 !! report_ibudget: 17 30 !!---------------------------------------------------------------------- 18 31 USE par_oce ! ocean parameters … … 85 98 !!---------------------------------------------------------------------- 86 99 ! 87 IF( .NOT. ln_bergdia )RETURN100 IF( .NOT.ln_bergdia ) RETURN 88 101 89 102 ALLOCATE( berg_melt (jpi,jpj) ) ; berg_melt (:,:) = 0._wp … … 136 149 137 150 floating_mass_start = icb_utl_mass( first_berg ) 138 bergs_mass_start = icb_utl_mass( first_berg, justbergs=. true. )139 bits_mass_start = icb_utl_mass( first_berg, justbits =.true. )151 bergs_mass_start = icb_utl_mass( first_berg, justbergs=.TRUE. ) 152 bits_mass_start = icb_utl_mass( first_berg, justbits =.TRUE. ) 140 153 IF( lk_mpp ) THEN 141 154 ALLOCATE( rsumbuf(23) ) ; rsumbuf(:) = 0._wp … … 159 172 !! for this we pack variables into buffer so we only need one mpp_sum 160 173 !!---------------------------------------------------------------------- 161 LOGICAL, INTENT(in) :: ld_budge 162 ! 163 INTEGER 164 REAL(wp) 165 !!---------------------------------------------------------------------- 166 ! 167 IF( .NOT. 174 LOGICAL, INTENT(in) :: ld_budge ! 175 ! 176 INTEGER :: ik 177 REAL(wp):: zunused_calving, ztmpsum, zgrdd_berg_mass, zgrdd_bits_mass 178 !!---------------------------------------------------------------------- 179 ! 180 IF( .NOT.ln_bergdia ) RETURN 168 181 169 182 zunused_calving = SUM( berg_grid%calving(:,:) ) … … 181 194 ztmpsum = SUM( berg_grid%calving_hflx(:,:) * e1e2t(:,:) * tmask_i(:,:) ) 182 195 calving_out_heat_net = calving_out_heat_net + ztmpsum * berg_dt ! Units of J 183 196 ! 184 197 IF( ld_budge ) THEN 185 198 stored_end = SUM( berg_grid%stored_ice(:,:,:) ) 186 199 stored_heat_end = SUM( berg_grid%stored_heat(:,:) ) 187 200 floating_mass_end = icb_utl_mass( first_berg ) 188 bergs_mass_end = icb_utl_mass( first_berg,justbergs=. true. )189 bits_mass_end = icb_utl_mass( first_berg,justbits =.true. )201 bergs_mass_end = icb_utl_mass( first_berg,justbergs=.TRUE. ) 202 bits_mass_end = icb_utl_mass( first_berg,justbits =.TRUE. ) 190 203 floating_heat_end = icb_utl_heat( first_berg ) 191 204 ! 192 205 nbergs_end = icb_utl_count() 193 206 zgrdd_berg_mass = SUM( berg_mass(:,:)*e1e2t(:,:)*tmask_i(:,:) ) 194 207 zgrdd_bits_mass = SUM( bits_mass(:,:)*e1e2t(:,:)*tmask_i(:,:) ) 195 208 ! 196 209 IF( lk_mpp ) THEN 197 210 rsumbuf( 1) = stored_end … … 218 231 rsumbuf(22) = zgrdd_berg_mass 219 232 rsumbuf(23) = zgrdd_bits_mass 220 233 ! 221 234 CALL mpp_sum( rsumbuf(1:23), 23) 222 235 ! 223 236 stored_end = rsumbuf( 1) 224 237 stored_heat_end = rsumbuf( 2) … … 244 257 zgrdd_berg_mass = rsumbuf(22) 245 258 zgrdd_bits_mass = rsumbuf(23) 246 259 ! 247 260 nsumbuf(1) = nbergs_end 248 261 nsumbuf(2) = nbergs_calved … … 252 265 nsumbuf(4+ik) = nbergs_calved_by_class(ik) 253 266 END DO 254 255 267 CALL mpp_sum( nsumbuf(1:nclasses+4), nclasses+4 ) 256 268 ! 257 269 nbergs_end = nsumbuf(1) 258 270 nbergs_calved = nsumbuf(2) … … 261 273 DO ik = 1,nclasses 262 274 nbergs_calved_by_class(ik)= nsumbuf(4+ik) 263 END DO264 275 END DO 276 ! 265 277 ENDIF 266 267 CALL report_state ( 'stored ice','kg','',stored_start,'',stored_end,'')268 CALL report_state ( 'floating','kg','',floating_mass_start,'',floating_mass_end,'',nbergs_end)269 CALL report_state ( 'icebergs','kg','',bergs_mass_start,'',bergs_mass_end,'')270 CALL report_state ( 'bits','kg','',bits_mass_start,'',bits_mass_end,'')271 CALL report_istate ( 'berg #','',nbergs_start,'',nbergs_end,'')278 ! 279 CALL report_state ( 'stored ice','kg','',stored_start,'',stored_end,'') 280 CALL report_state ( 'floating','kg','',floating_mass_start,'',floating_mass_end,'',nbergs_end ) 281 CALL report_state ( 'icebergs','kg','',bergs_mass_start,'',bergs_mass_end,'') 282 CALL report_state ( 'bits','kg','',bits_mass_start,'',bits_mass_end,'') 283 CALL report_istate ( 'berg #','',nbergs_start,'',nbergs_end,'') 272 284 CALL report_ibudget( 'berg #','calved',nbergs_calved, & 273 274 285 & 'melted',nbergs_melted, & 286 & '#',nbergs_start,nbergs_end) 275 287 CALL report_budget( 'stored mass','kg','calving used',calving_used_net, & 276 277 288 & 'bergs',calving_to_bergs_net, & 289 & 'stored mass',stored_start,stored_end) 278 290 CALL report_budget( 'floating mass','kg','calving used',calving_to_bergs_net, & 279 280 291 & 'bergs',melt_net, & 292 & 'stored mass',floating_mass_start,floating_mass_end) 281 293 CALL report_budget( 'berg mass','kg','calving',calving_to_bergs_net, & 282 283 294 & 'melt+eros',berg_melt_net, & 295 & 'berg mass',bergs_mass_start,bergs_mass_end) 284 296 CALL report_budget( 'bits mass','kg','eros used',bits_src_net, & 285 286 297 & 'bergs',bits_melt_net, & 298 & 'stored mass',bits_mass_start,bits_mass_end) 287 299 CALL report_budget( 'net mass','kg','recvd',calving_rcv_net, & 288 289 290 300 & 'rtrnd',calving_ret_net, & 301 & 'net mass',stored_start+floating_mass_start, & 302 & stored_end+floating_mass_end) 291 303 CALL report_consistant( 'iceberg mass','kg','gridded',zgrdd_berg_mass,'bergs',bergs_mass_end) 292 304 CALL report_consistant( 'bits mass','kg','gridded',zgrdd_bits_mass,'bits',bits_mass_end) 293 305 CALL report_state( 'net heat','J','',stored_heat_start+floating_heat_start,'', & 294 306 & stored_heat_end+floating_heat_end,'') 295 307 CALL report_state( 'stored heat','J','',stored_heat_start,'',stored_heat_end,'') 296 308 CALL report_state( 'floating heat','J','',floating_heat_start,'',floating_heat_end,'') 297 309 CALL report_budget( 'net heat','J','net heat',calving_src_heat_net, & 298 299 300 310 & 'net heat',calving_out_heat_net, & 311 & 'net heat',stored_heat_start+floating_heat_start, & 312 & stored_heat_end+floating_heat_end) 301 313 CALL report_budget( 'stored heat','J','calving used',calving_src_heat_used_net, & 302 303 314 & 'bergs',heat_to_bergs_net, & 315 & 'net heat',stored_heat_start,stored_heat_end) 304 316 CALL report_budget( 'flting heat','J','calved',heat_to_bergs_net, & 305 306 317 & 'melt',heat_to_ocean_net, & 318 & 'net heat',floating_heat_start,floating_heat_end) 307 319 IF (nn_verbose_level >= 1) THEN 308 320 CALL report_consistant( 'top interface','kg','from SIS',calving_src_net, & 309 321 & 'received',calving_rcv_net) 310 322 CALL report_consistant( 'bot interface','kg','sent',calving_out_net, & 311 323 & 'returned',calving_ret_net) 312 324 ENDIF 313 325 WRITE( numicb, '("calved by class = ",i6,20(",",i6))') (nbergs_calved_by_class(ik),ik=1,nclasses) 314 IF ( nspeeding_tickets > 0 )WRITE( numicb, '("speeding tickets issued = ",i6)') nspeeding_tickets315 326 IF( nspeeding_tickets > 0 ) WRITE( numicb, '("speeding tickets issued = ",i6)') nspeeding_tickets 327 ! 316 328 nbergs_start = nbergs_end 317 329 stored_start = stored_end … … 350 362 !!---------------------------------------------------------------------- 351 363 ! 352 IF( .NOT. ln_bergdia )RETURN353 berg_melt 354 buoy_melt 355 eros_melt 356 conv_melt 357 bits_src 358 bits_melt 359 bits_mass 360 berg_mass 361 virtual_area 362 real_calving 364 IF( .NOT.ln_bergdia ) RETURN 365 berg_melt (:,:) = 0._wp 366 buoy_melt (:,:) = 0._wp 367 eros_melt (:,:) = 0._wp 368 conv_melt (:,:) = 0._wp 369 bits_src (:,:) = 0._wp 370 bits_melt (:,:) = 0._wp 371 bits_mass (:,:) = 0._wp 372 berg_mass (:,:) = 0._wp 373 virtual_area(:,:) = 0._wp 374 real_calving(:,:,:) = 0._wp 363 375 ! 364 376 END SUBROUTINE icb_dia_step … … 369 381 !!---------------------------------------------------------------------- 370 382 ! 371 IF( .NOT. 383 IF( .NOT.ln_bergdia ) RETURN !!gm useless iom will control whether it is output or not 372 384 ! 373 385 CALL iom_put( "berg_melt" , berg_melt (:,:) ) ! Melt rate of icebergs [kg/m2/s] … … 388 400 !!---------------------------------------------------------------------- 389 401 !!---------------------------------------------------------------------- 390 INTEGER ,INTENT(in) :: ki, kj, kn402 INTEGER , INTENT(in) :: ki, kj, kn 391 403 REAL(wp), INTENT(in) :: pcalved 392 404 REAL(wp), INTENT(in) :: pheated … … 411 423 !!---------------------------------------------------------------------- 412 424 ! 413 IF( .NOT. ln_bergdia )RETURN425 IF( .NOT.ln_bergdia ) RETURN 414 426 ! 415 427 IF( kt == nit000 ) THEN … … 437 449 !!---------------------------------------------------------------------- 438 450 !!---------------------------------------------------------------------- 439 INTEGER , INTENT(in) ::ki, kj440 REAL(wp), INTENT(in) :: pWn, pLn, pAbits, pmass_scale, pMnew, pnMbits, pz1_e1e2441 !!---------------------------------------------------------------------- 442 ! 443 IF( .NOT. ln_bergdia )RETURN451 INTEGER , INTENT(in) :: ki, kj 452 REAL(wp), INTENT(in) :: pWn, pLn, pAbits, pmass_scale, pMnew, pnMbits, pz1_e1e2 453 !!---------------------------------------------------------------------- 454 ! 455 IF( .NOT.ln_bergdia ) RETURN 444 456 virtual_area(ki,kj) = virtual_area(ki,kj) + ( pWn * pLn + pAbits ) * pmass_scale ! m^2 445 457 berg_mass(ki,kj) = berg_mass(ki,kj) + pMnew * pz1_e1e2 ! kg/m2 … … 453 465 !!---------------------------------------------------------------------- 454 466 ! 455 IF( .NOT. ln_bergdia )RETURN467 IF( .NOT.ln_bergdia ) RETURN 456 468 nspeeding_tickets = nspeeding_tickets + 1 457 469 ! … … 459 471 460 472 461 SUBROUTINE icb_dia_melt(ki, kj, pmnew, pheat, pmass_scale, &462 & pdM, pdMbitsE, pdMbitsM, pdMb, pdMe, &463 & pdMv, pz1_dt_e1e2 )473 SUBROUTINE icb_dia_melt(ki, kj, pmnew, pheat, pmass_scale, & 474 & pdM, pdMbitsE, pdMbitsM, pdMb, pdMe, & 475 & pdMv, pz1_dt_e1e2 ) 464 476 !!---------------------------------------------------------------------- 465 477 !!---------------------------------------------------------------------- … … 469 481 !!---------------------------------------------------------------------- 470 482 ! 471 IF( .NOT. ln_bergdia )RETURN472 483 IF( .NOT.ln_bergdia ) RETURN 484 ! 473 485 berg_melt (ki,kj) = berg_melt (ki,kj) + pdM * pz1_dt_e1e2 ! kg/m2/s 474 486 bits_src (ki,kj) = bits_src (ki,kj) + pdMbitsE * pz1_dt_e1e2 ! mass flux into bergy bitskg/m2/s … … 492 504 !!---------------------------------------------------------------------- 493 505 ! 494 IF 506 IF( PRESENT(kbergs) ) THEN 495 507 WRITE(numicb,100) cd_budgetstr // ' state:', & 496 497 498 499 508 & cd_startstr // ' start', pstartval, cd_budgetunits, & 509 & cd_endstr // ' end', pendval, cd_budgetunits, & 510 & 'Delta ' // cd_delstr, pendval-pstartval, cd_budgetunits, & 511 & '# of bergs', kbergs 500 512 ELSE 501 513 WRITE(numicb,100) cd_budgetstr // ' state:', & 502 503 504 514 & cd_startstr // ' start', pstartval, cd_budgetunits, & 515 & cd_endstr // ' end', pendval, cd_budgetunits, & 516 & cd_delstr // 'Delta', pendval-pstartval, cd_budgetunits 505 517 ENDIF 506 100 FORMAT(a19,3(a18,"=",es14.7,x,a2,:,","),a12,i8) 518 100 FORMAT(a19,3(a18,"=",es14.7,x,a2,:,","),a12,i8) 519 ! 507 520 END SUBROUTINE report_state 508 521 … … 516 529 ! 517 530 WRITE(numicb,200) cd_budgetstr // ' check:', & 518 cd_startstr, pstartval, cd_budgetunits, & 519 cd_endstr, pendval, cd_budgetunits, & 520 'error', (pendval-pstartval)/((pendval+pstartval)+1e-30), 'nd' 521 200 FORMAT(a19,10(a18,"=",es14.7,x,a2,:,",")) 531 & cd_startstr, pstartval, cd_budgetunits, & 532 & cd_endstr, pendval, cd_budgetunits, & 533 & 'error', (pendval-pstartval)/((pendval+pstartval)+1e-30), 'nd' 534 200 FORMAT(a19,10(a18,"=",es14.7,x,a2,:,",")) 535 ! 522 536 END SUBROUTINE report_consistant 523 537 … … 530 544 REAL(wp), INTENT(in) :: pinval, poutval, pstartval, pendval 531 545 ! 532 REAL(wp) ::zval546 REAL(wp) :: zval 533 547 !!---------------------------------------------------------------------- 534 548 ! 535 549 zval = ( ( pendval - pstartval ) - ( pinval - poutval ) ) / & 536 & MAX( 1.e-30, MAX( abs( pendval - pstartval ) , ABS( pinval - poutval ) ) )537 550 & MAX( 1.e-30, MAX( ABS( pendval - pstartval ) , ABS( pinval - poutval ) ) ) 551 ! 538 552 WRITE(numicb,200) cd_budgetstr // ' budget:', & 539 553 & cd_instr // ' in', pinval, cd_budgetunits, & … … 549 563 !!---------------------------------------------------------------------- 550 564 !!---------------------------------------------------------------------- 551 CHARACTER*(*), INTENT(in) :: cd_budgetstr, cd_startstr, cd_endstr, cd_delstr 552 INTEGER, INTENT(in) :: pstartval, pendval 565 CHARACTER*(*), INTENT(in) :: cd_budgetstr, cd_startstr, cd_endstr, cd_delstr 566 INTEGER , INTENT(in) :: pstartval, pendval 567 !!---------------------------------------------------------------------- 553 568 ! 554 569 WRITE(numicb,100) cd_budgetstr // ' state:', & … … 570 585 ! 571 586 WRITE(numicb,200) cd_budgetstr // ' budget:', & 572 cd_instr // ' in', pinval, & 573 cd_outstr // ' out', poutval, & 574 'Delta ' // cd_delstr, pinval-poutval, & 575 'error', ( ( pendval - pstartval ) - ( pinval - poutval ) ) 576 200 FORMAT(a19,10(a18,"=",i14,x,:,",")) 587 & cd_instr // ' in', pinval, & 588 & cd_outstr // ' out', poutval, & 589 & 'Delta ' // cd_delstr, pinval-poutval, & 590 & 'error', ( ( pendval - pstartval ) - ( pinval - poutval ) ) 591 200 FORMAT(a19,10(a18,"=",i14,x,:,",")) 592 ! 577 593 END SUBROUTINE report_ibudget 578 594 -
branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/ICB/icbdyn.F90
r5215 r9190 1 1 MODULE icbdyn 2 3 2 !!====================================================================== 4 3 !! *** MODULE icbdyn *** 5 4 !! Iceberg: time stepping routine for iceberg tracking 6 5 !!====================================================================== 7 !! History : 3.3.1 ! 2010-01 (Martin&Adcroft) Original code 8 !! - ! 2011-03 (Madec) Part conversion to NEMO form 9 !! - ! Removal of mapping from another grid 10 !! - ! 2011-04 (Alderson) Split into separate modules 11 !! - ! 2011-05 (Alderson) Replace broken grounding routine 12 !! - ! with one of Gurvan's suggestions (just like 13 !! - ! the broken one) 6 !! History : 3.3 ! 2010-01 (Martin&Adcroft) Original code 7 !! - ! 2011-03 (Madec) Part conversion to NEMO form 8 !! - ! Removal of mapping from another grid 9 !! - ! 2011-04 (Alderson) Split into separate modules 10 !! - ! 2011-05 (Alderson) Replace broken grounding routine with one of 11 !! - ! Gurvan's suggestions (just like the broken one) 14 12 !!---------------------------------------------------------------------- 15 13 USE par_oce ! NEMO parameters … … 41 39 !! ** Method : - See Martin & Adcroft, Ocean Modelling 34, 2010 42 40 !!---------------------------------------------------------------------- 43 REAL(wp) :: zuvel1 , zvvel1 , zu1, zv1, zax1, zay1, zxi1 , zyj1 44 REAL(wp) :: zuvel2 , zvvel2 , zu2, zv2, zax2, zay2, zxi2 , zyj2 45 REAL(wp) :: zuvel3 , zvvel3 , zu3, zv3, zax3, zay3, zxi3 , zyj3 46 REAL(wp) :: zuvel4 , zvvel4 , zu4, zv4, zax4, zay4, zxi4 , zyj4 47 REAL(wp) :: zuvel_n, zvvel_n, zxi_n , zyj_n 48 REAL(wp) :: zdt, zdt_2, zdt_6, ze1, ze2 49 LOGICAL :: ll_bounced 50 TYPE(iceberg), POINTER :: berg 51 TYPE(point) , POINTER :: pt 52 INTEGER :: kt 53 !!---------------------------------------------------------------------- 54 41 INTEGER, INTENT(in) :: kt ! 42 ! 43 LOGICAL :: ll_bounced 44 REAL(wp) :: zuvel1 , zvvel1 , zu1, zv1, zax1, zay1, zxi1 , zyj1 45 REAL(wp) :: zuvel2 , zvvel2 , zu2, zv2, zax2, zay2, zxi2 , zyj2 46 REAL(wp) :: zuvel3 , zvvel3 , zu3, zv3, zax3, zay3, zxi3 , zyj3 47 REAL(wp) :: zuvel4 , zvvel4 , zu4, zv4, zax4, zay4, zxi4 , zyj4 48 REAL(wp) :: zuvel_n, zvvel_n, zxi_n , zyj_n 49 REAL(wp) :: zdt, zdt_2, zdt_6, ze1, ze2 50 TYPE(iceberg), POINTER :: berg 51 TYPE(point) , POINTER :: pt 52 !!---------------------------------------------------------------------- 53 ! 55 54 ! 4th order Runge-Kutta to solve: d/dt X = V, d/dt V = A 56 55 ! with I.C.'s: X=X1 and V=V1 … … 75 74 pt => berg%current_point 76 75 77 ll_bounced = . false.76 ll_bounced = .FALSE. 78 77 79 78 … … 99 98 ! 100 99 CALL icb_ground( zxi2, zxi1, zu1, & 101 &zyj2, zyj1, zv1, ll_bounced )100 & zyj2, zyj1, zv1, ll_bounced ) 102 101 103 102 ! !** A2 = A(X2,V2) … … 115 114 ! 116 115 CALL icb_ground( zxi3, zxi1, zu3, & 117 & zyj3, zyj1, zv3, ll_bounced )116 & zyj3, zyj1, zv3, ll_bounced ) 118 117 119 118 ! !** A3 = A(X3,V3) … … 131 130 132 131 CALL icb_ground( zxi4, zxi1, zu4, & 133 &zyj4, zyj1, zv4, ll_bounced )132 & zyj4, zyj1, zv4, ll_bounced ) 134 133 135 134 ! !** A4 = A(X4,V4) … … 150 149 151 150 CALL icb_ground( zxi_n, zxi1, zuvel_n, & 152 &zyj_n, zyj1, zvvel_n, ll_bounced )151 & zyj_n, zyj1, zvvel_n, ll_bounced ) 153 152 154 153 pt%uvel = zuvel_n !** save in berg structure … … 169 168 170 169 SUBROUTINE icb_ground( pi, pi0, pu, & 171 & 170 & pj, pj0, pv, ld_bounced ) 172 171 !!---------------------------------------------------------------------- 173 172 !! *** ROUTINE icb_ground *** … … 216 215 ibounce_method = 2 217 216 SELECT CASE ( ibounce_method ) 218 219 220 221 222 223 224 225 226 227 228 229 230 231 217 CASE ( 1 ) 218 pi = pi0 219 pj = pj0 220 pu = 0._wp 221 pv = 0._wp 222 CASE ( 2 ) 223 IF( ii0 /= ii ) THEN 224 pi = pi0 ! return back to the initial position 225 pu = 0._wp ! zeroing of velocity in the direction of the grounding 226 ENDIF 227 IF( ij0 /= ij ) THEN 228 pj = pj0 ! return back to the initial position 229 pv = 0._wp ! zeroing of velocity in the direction of the grounding 230 ENDIF 232 231 END SELECT 233 232 ! … … 259 258 ! 260 259 INTEGER :: itloop 261 REAL(wp) :: zuo, zvo, zui, zvi, zua, zva, zuwave, zvwave, zssh_x, zssh_y, zsst, zcn, zhi 260 REAL(wp) :: zuo, zui, zua, zuwave, zssh_x, zsst, zcn, zhi 261 REAL(wp) :: zvo, zvi, zva, zvwave, zssh_y 262 262 REAL(wp) :: zff, zT, zD, zW, zL, zM, zF 263 263 REAL(wp) :: zdrag_ocn, zdrag_atm, zdrag_ice, zwave_rad … … 339 339 zaxe = zaxe - zdrag_ocn*(puvel -zuo) - zdrag_atm*(puvel -zua) -zdrag_ice*(puvel -zui) 340 340 zaye = zaye - zdrag_ocn*(pvvel -zvo) - zdrag_atm*(pvvel -zva) -zdrag_ice*(pvvel -zvi) 341 endif341 ENDIF 342 342 343 343 ! Solve for implicit accelerations … … 349 349 pax = zdetA * ( zA11*zaxe + zA12*zaye ) 350 350 pay = zdetA * ( zA11*zaye - zA12*zaxe ) 351 else351 ELSE 352 352 pax = zaxe ; pay = zaye 353 endif353 ENDIF 354 354 355 355 zuveln = puvel0 + pdt*pax … … 362 362 IF( zspeed > 0._wp ) THEN 363 363 zloc_dx = MIN( pe1, pe2 ) ! minimum grid spacing 364 zspeed_new = zloc_dx / pdt * rn_speed_limit ! Speed limit as a factor of dx / dt364 zspeed_new = zloc_dx / pdt * rn_speed_limit ! Speed limit as a factor of dx / dt 365 365 IF( zspeed_new < zspeed ) THEN 366 366 zuveln = zuveln * ( zspeed_new / zspeed ) ! Scale velocity to reduce speed -
branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/ICB/icbini.F90
r9168 r9190 361 361 !!---------------------------------------------------------------------- 362 362 363 #if !defined key_agrif 363 #if defined key_agrif 364 IF(lwp) THEN 365 WRITE(numout,*) 366 WRITE(numout,*) 'icb_nam : AGRIF is not compatible with namelist namberg : ' 367 WRITE(numout,*) '~~~~~~~ definition of rn_initial_mass(nclasses) with nclasses as PARAMETER ' 368 WRITE(numout,*) 369 WRITE(numout,*) ' ==>>> force NO icebergs used. The namelist namberg is not read' 370 ENDIF 371 ln_icebergs = .false. 372 RETURN 373 #else 374 IF(lwp) THEN 375 WRITE(numout,*) 376 WRITE(numout,*) 'icb_nam : iceberg initialization through namberg namelist read' 377 WRITE(numout,*) '~~~~~~~~ ' 378 ENDIF 379 #endif 380 ! !== read namelist ==! 364 381 REWIND( numnam_ref ) ! Namelist namberg in reference namelist : Iceberg parameters 365 382 READ ( numnam_ref, namberg, IOSTAT = ios, ERR = 901) … … 369 386 902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namberg in configuration namelist', lwp ) 370 387 IF(lwm) WRITE ( numond, namberg ) 371 #else 372 IF(lwp) THEN 373 WRITE(numout,*) 374 WRITE(numout,*) 'icbini : AGRIF is not compatible with namelist namberg : ' 375 WRITE(numout,*) '~~~~~~ definition of rn_initial_mass(nclasses) with nclasses as PARAMETER ' 376 WRITE(numout,*) ' ==>>> force NO icebergs used. The namelist namberg is not read' 377 ENDIF 378 ln_icebergs = .false. 379 #endif 380 IF( .NOT. ln_icebergs ) THEN ! no icebergs 381 IF(lwp) THEN 382 WRITE(numout,*) 383 WRITE(numout,*) 'icbini : Namelist namberg ln_icebergs = F , NO icebergs used' 384 WRITE(numout,*) '~~~~~~ ' 385 ENDIF 388 ! 389 IF(lwp) WRITE(numout,*) 390 IF( ln_icebergs ) THEN 391 IF(lwp) WRITE(numout,*) ' ==>>> icebergs are used' 392 ELSE 393 IF(lwp) WRITE(numout,*) ' ==>>> No icebergs used' 386 394 RETURN 387 395 ENDIF 388 389 390 ! IF( lk_lim3 .AND. ln_icebergs ) THEN 391 ! CALL ctl_stop( 'icb_nam: the use of ICB with LIM3 not allowed. ice thickness missing in ICB' ) 392 ! ENDIF 393 396 ! 397 IF( nn_test_icebergs > nclasses ) THEN 398 IF(lwp) WRITE(numout,*) 399 IF(lwp) WRITE(numout,*) ' ==>>> Resetting of nn_test_icebergs to ', nclasses 400 nn_test_icebergs = nclasses 401 ENDIF 402 ! 394 403 IF(lwp) THEN ! control print 395 404 WRITE(numout,*) … … 399 408 WRITE(numout,*) ' Period between sampling of position for trajectory storage nn_sample_rate = ', nn_sample_rate 400 409 WRITE(numout,*) ' Mass thresholds between iceberg classes (kg) rn_initial_mass =' 401 DO jn =1,nclasses402 WRITE(numout,'(a,f15.2)') ' ', rn_initial_mass(jn)410 DO jn = 1, nclasses 411 WRITE(numout,'(a,f15.2)') ' ', rn_initial_mass(jn) 403 412 ENDDO 404 413 WRITE(numout,*) ' Fraction of calving to apply to this class (non-dim) rn_distribution =' 405 414 DO jn = 1, nclasses 406 WRITE(numout,'(a,f10.4)') ' ', rn_distribution(jn)415 WRITE(numout,'(a,f10.4)') ' ', rn_distribution(jn) 407 416 END DO 408 417 WRITE(numout,*) ' Ratio between effective and real iceberg mass (non-dim) rn_mass_scaling = ' 409 418 DO jn = 1, nclasses 410 WRITE(numout,'(a,f10.2)') ' ', rn_mass_scaling(jn)419 WRITE(numout,'(a,f10.2)') ' ', rn_mass_scaling(jn) 411 420 END DO 412 421 WRITE(numout,*) ' Total thickness of newly calved bergs (m) rn_initial_thickness = ' 413 422 DO jn = 1, nclasses 414 WRITE(numout,'(a,f10.2)') ' ', rn_initial_thickness(jn)423 WRITE(numout,'(a,f10.2)') ' ', rn_initial_thickness(jn) 415 424 END DO 416 425 WRITE(numout,*) ' Timesteps between verbose messages nn_verbose_write = ', nn_verbose_write … … 435 444 ENDIF 436 445 ! 437 IF( nn_test_icebergs > nclasses ) THEN438 IF(lwp) WRITE(numout,*) ' ==>>> Resetting of nn_test_icebergs to ', nclasses439 nn_test_icebergs = nclasses440 ENDIF441 442 446 ! ensure that the sum of berg input distribution is equal to one 443 447 zfact = SUM( rn_distribution ) -
branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/ICB/icblbc.F90
r5215 r9190 1 1 MODULE icblbc 2 3 2 !!====================================================================== 4 3 !! *** MODULE icblbc *** 5 4 !! Ocean physics: routines to handle boundary exchanges for icebergs 6 5 !!====================================================================== 7 !! History : 3.3.1 ! 2010-01 (Martin&Adcroft) Original code 8 !! - ! 2011-03 (Madec) Part conversion to NEMO form 9 !! - ! Removal of mapping from another grid 10 !! - ! 2011-04 (Alderson) Split into separate modules 11 !! - ! 2011-05 (Alderson) MPP exchanges written based on lib_mpp 12 !! - ! 2011-05 (Alderson) MPP and single processor boundary 13 !! - ! conditions added 6 !! History : 3.3 ! 2010-01 (Martin&Adcroft) Original code 7 !! - ! 2011-03 (Madec) Part conversion to NEMO form 8 !! - ! Removal of mapping from another grid 9 !! - ! 2011-04 (Alderson) Split into separate modules 10 !! - ! 2011-05 (Alderson) MPP exchanges written based on lib_mpp 11 !! - ! 2011-05 (Alderson) MPP and single processor boundary conditions added 14 12 !!---------------------------------------------------------------------- 13 15 14 !!---------------------------------------------------------------------- 16 15 !! icb_lbc : - Pass icebergs across cyclic boundaries … … 27 26 !! nicbfldpts - packed i,j point in exchanging processor 28 27 !!---------------------------------------------------------------------- 29 30 28 USE par_oce ! ocean parameters 31 29 USE dom_oce ! ocean domain … … 45 43 46 44 TYPE, PUBLIC :: buffer 47 INTEGER :: size =048 REAL(wp), DIMENSION(:,:), POINTER :: data45 INTEGER :: size = 0 46 REAL(wp), DIMENSION(:,:), POINTER :: data 49 47 END TYPE buffer 50 48 51 TYPE(buffer), POINTER 52 TYPE(buffer), POINTER 53 TYPE(buffer), POINTER 54 TYPE(buffer), POINTER 49 TYPE(buffer), POINTER :: obuffer_n=>NULL() , ibuffer_n=>NULL() 50 TYPE(buffer), POINTER :: obuffer_s=>NULL() , ibuffer_s=>NULL() 51 TYPE(buffer), POINTER :: obuffer_e=>NULL() , ibuffer_e=>NULL() 52 TYPE(buffer), POINTER :: obuffer_w=>NULL() , ibuffer_w=>NULL() 55 53 56 54 ! north fold exchange buffers 57 TYPE(buffer), POINTER 58 59 INTEGER, PARAMETER, PRIVATE 60 INTEGER, PARAMETER, PRIVATE 55 TYPE(buffer), POINTER :: obuffer_f=>NULL() , ibuffer_f=>NULL() 56 57 INTEGER, PARAMETER, PRIVATE :: jp_delta_buf = 25 ! Size by which to increment buffers 58 INTEGER, PARAMETER, PRIVATE :: jp_buffer_width = 15+nkounts ! items to store for each berg 61 59 62 60 #endif … … 926 924 WRITE(numout,*) 'icb_lbc_mpp: You should not have seen this message!!' 927 925 END SUBROUTINE icb_lbc_mpp 928 929 926 #endif 930 927 -
branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/ICB/icbrst.F90
r6623 r9190 1 1 MODULE icbrst 2 3 2 !!====================================================================== 4 3 !! *** MODULE icbrst *** … … 15 14 !! read single restart files 16 15 !!---------------------------------------------------------------------- 16 17 17 !!---------------------------------------------------------------------- 18 18 !! icb_rst_read : read restart file … … 110 110 CALL iom_get( ncid, 'mass_of_bits' , localpt%mass_of_bits , ktime=jn ) 111 111 CALL iom_get( ncid, 'heat_density' , localpt%heat_density , ktime=jn ) 112 113 112 ! 114 113 CALL icb_utl_add( localberg, localpt ) 115 114 ! 116 115 ENDIF 117 116 ! 118 117 END DO 119 118 ! 120 119 ENDIF 121 120 … … 144 143 CALL iom_close( ncid ) 145 144 ! 146 IF( lwp . and. nn_verbose_level >= 0) WRITE(numout,'(a)') 'icebergs, read_restart_bergs: completed'145 IF( lwp .AND. nn_verbose_level >= 0) WRITE(numout,'(a)') 'icebergs, read_restart_bergs: completed' 147 146 ! 148 147 END SUBROUTINE icb_rst_read … … 361 360 END SUBROUTINE icb_rst_write 362 361 ! 362 !!====================================================================== 363 363 END MODULE icbrst -
branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/ICB/icbstp.F90
r9124 r9190 171 171 END SUBROUTINE icb_end 172 172 173 !! -------------------------------------------------------------------------173 !!====================================================================== 174 174 END MODULE icbstp -
branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/ICB/icbthm.F90
r5836 r9190 1 1 MODULE icbthm 2 3 2 !!====================================================================== 4 3 !! *** MODULE icbthm *** … … 31 30 PUBLIC icb_thm ! routine called in icbstp.F90 module 32 31 32 !!---------------------------------------------------------------------- 33 !! NEMO/OPA 3.3 , NEMO Consortium (2011) 33 34 !! $Id$ 35 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 36 !!---------------------------------------------------------------------- 34 37 CONTAINS 35 38 … … 55 58 ! 56 59 z1_rday = 1._wp / rday 57 60 ! 58 61 ! we're either going to ignore berg fresh water melt flux and associated heat 59 62 ! or we pass it into the ocean, so at this point we set them both to zero, … … 63 66 berg_grid%floating_melt(:,:) = 0._wp 64 67 berg_grid%calving_hflx(:,:) = 0._wp 65 68 ! 66 69 this => first_berg 67 DO WHILE( associated(this) )70 DO WHILE( ASSOCIATED(this) ) 68 71 ! 69 72 pt => this%current_point 70 73 nknberg = this%number(1) 71 CALL icb_utl_interp( pt%xi, pt%e1, pt%uo, pt%ui, pt%ua, pt%ssh_x, &72 & pt%yj, pt%e2, pt%vo, pt%vi, pt%va, pt%ssh_y,&73 & pt%sst, pt%cn, pt%hi, zff )74 CALL icb_utl_interp( pt%xi, pt%e1, pt%uo, pt%ui, pt%ua, pt%ssh_x, & 75 & pt%yj, pt%e2, pt%vo, pt%vi, pt%va, pt%ssh_y, & 76 & pt%sst, pt%cn, pt%hi, zff ) 74 77 ! 75 78 zSST = pt%sst … … 98 101 zMv = MAX( 7.62e-3*zSST+1.29e-3*(zSST**2) , 0._wp ) * z1_rday ! Buoyant convection at sides (eqn M.A10) 99 102 zMb = MAX( 0.58*(zdvo**0.8)*(zSST+4.0)/(zL**0.2) , 0._wp ) * z1_rday ! Basal turbulent melting (eqn M.A7 ) 100 zMe = MAX( 1./12.*(zSST+2.)*zSs*(1+ cos(rpi*(zIC**3))) , 0._wp ) * z1_rday ! Wave erosion (eqn M.A8 )103 zMe = MAX( 1./12.*(zSST+2.)*zSs*(1+COS(rpi*(zIC**3))) , 0._wp ) * z1_rday ! Wave erosion (eqn M.A8 ) 101 104 102 105 IF( ln_operator_splitting ) THEN ! Operator split update of volume/mass 103 106 zTn = MAX( zT - zMb*zdt , 0._wp ) ! new total thickness (m) 104 znVol = zTn * zW * zL 105 zMnew1 = (znVol/zVol) * zM 107 znVol = zTn * zW * zL ! new volume (m^3) 108 zMnew1 = (znVol/zVol) * zM ! new mass (kg) 106 109 zdMb = zM - zMnew1 ! mass lost to basal melting (>0) (kg) 107 110 ! 108 111 zLn = MAX( zL - zMv*zdt , 0._wp ) ! new length (m) 109 112 zWn = MAX( zW - zMv*zdt , 0._wp ) ! new width (m) 110 znVol = zTn * zWn * zLn 111 zMnew2 = (znVol/zVol) * zM 113 znVol = zTn * zWn * zLn ! new volume (m^3) 114 zMnew2 = (znVol/zVol) * zM ! new mass (kg) 112 115 zdMv = zMnew1 - zMnew2 ! mass lost to buoyant convection (>0) (kg) 113 116 ! 114 117 zLn = MAX( zLn - zMe*zdt , 0._wp ) ! new length (m) 115 118 zWn = MAX( zWn - zMe*zdt , 0._wp ) ! new width (m) 116 znVol = zTn * zWn * zLn 117 zMnew = ( znVol / zVol ) * zM 119 znVol = zTn * zWn * zLn ! new volume (m^3) 120 zMnew = ( znVol / zVol ) * zM ! new mass (kg) 118 121 zdMe = zMnew2 - zMnew ! mass lost to erosion (>0) (kg) 119 122 zdM = zM - zMnew ! mass lost to all erosion and melting (>0) (kg) 120 123 ! 121 124 ELSE ! Update dimensions of berg 122 zLn = MAX( zL -(zMv+zMe)*zdt ,0._wp ) 123 zWn = MAX( zW -(zMv+zMe)*zdt ,0._wp ) 125 zLn = MAX( zL -(zMv+zMe)*zdt ,0._wp ) ! (m) 126 zWn = MAX( zW -(zMv+zMe)*zdt ,0._wp ) ! (m) 124 127 zTn = MAX( zT - zMb *zdt ,0._wp ) ! (m) 125 128 ! Update volume and mass of berg 126 znVol = zTn*zWn*zLn 127 zMnew = (znVol/zVol)*zM 129 znVol = zTn*zWn*zLn ! (m^3) 130 zMnew = (znVol/zVol)*zM ! (kg) 128 131 zdM = zM - zMnew ! (kg) 129 zdMb = (zM/zVol) * (zW* zL ) *zMb*zdt 130 zdMe = (zM/zVol) * (zT*(zW+zL)) *zMe*zdt 131 zdMv = (zM/zVol) * (zT*(zW+zL)) *zMv*zdt 132 ENDIF 133 134 IF( rn_bits_erosion_fraction > 0._wp ) THEN 132 zdMb = (zM/zVol) * (zW* zL ) *zMb*zdt ! approx. mass loss to basal melting (kg) 133 zdMe = (zM/zVol) * (zT*(zW+zL)) *zMe*zdt ! approx. mass lost to erosion (kg) 134 zdMv = (zM/zVol) * (zT*(zW+zL)) *zMv*zdt ! approx. mass loss to buoyant convection (kg) 135 ENDIF 136 137 IF( rn_bits_erosion_fraction > 0._wp ) THEN ! Bergy bits 135 138 ! 136 139 zMbits = pt%mass_of_bits ! mass of bergy bits (kg) 137 zdMbitsE = rn_bits_erosion_fraction * zdMe ! change in mass of bits (kg)138 znMbits = zMbits + zdMbitsE 139 zLbits = MIN( zL, zW, zT, 40._wp ) 140 zAbits = ( zMbits / rn_rho_bergs ) / zLbits ! Effective bottom area (assuming T=Lbits)141 zMbb = MAX( 0.58*(zdvo**0.8)*(zSST+2.0)/(zLbits**0.2), 0.) * z1_rday 142 zMbb = rn_rho_bergs * zAbits * zMbb ! in kg/s143 zdMbitsM = MIN( zMbb*zdt , znMbits ) 144 znMbits = znMbits-zdMbitsM 140 zdMbitsE = rn_bits_erosion_fraction * zdMe ! change in mass of bits (kg) 141 znMbits = zMbits + zdMbitsE ! add new bergy bits to mass (kg) 142 zLbits = MIN( zL, zW, zT, 40._wp ) ! assume bergy bits are smallest dimension or 40 meters 143 zAbits = ( zMbits / rn_rho_bergs ) / zLbits ! Effective bottom area (assuming T=Lbits) 144 zMbb = MAX( 0.58*(zdvo**0.8)*(zSST+2.0)/(zLbits**0.2), 0.) * z1_rday ! Basal turbulent melting (for bits) 145 zMbb = rn_rho_bergs * zAbits * zMbb ! in kg/s 146 zdMbitsM = MIN( zMbb*zdt , znMbits ) ! bergy bits mass lost to melting (kg) 147 znMbits = znMbits-zdMbitsM ! remove mass lost to bergy bits melt 145 148 IF( zMnew == 0._wp ) THEN ! if parent berg has completely melted then 146 zdMbitsM = zdMbitsM + znMbits 149 zdMbitsM = zdMbitsM + znMbits ! instantly melt all the bergy bits 147 150 znMbits = 0._wp 148 151 ENDIF … … 163 166 berg_grid%calving_hflx (ii,ij) = berg_grid%calving_hflx (ii,ij) + zheat * z1_e1e2 ! W/m2 164 167 CALL icb_dia_melt( ii, ij, zMnew, zheat, this%mass_scaling, & 165 &zdM, zdMbitsE, zdMbitsM, zdMb, zdMe, &166 &zdMv, z1_dt_e1e2 )168 & zdM, zdMbitsE, zdMbitsM, zdMb, zdMe, & 169 & zdMv, z1_dt_e1e2 ) 167 170 ELSE 168 171 WRITE(numout,*) 'icb_thm: berg ',this%number(:),' appears to have grounded at ',narea,ii,ij … … 178 181 zTn = zWn 179 182 zWn = zT 180 endif183 ENDIF 181 184 182 185 ! Store the new state of iceberg (with L>W) … … 184 187 pt%mass_of_bits = znMbits 185 188 pt%thickness = zTn 186 pt%width = min(zWn,zLn)187 pt%length = max(zWn,zLn)189 pt%width = MIN( zWn , zLn ) 190 pt%length = MAX( zWn , zLn ) 188 191 189 192 next=>this%next … … 197 200 z1_e1e2 = r1_e1e2t(ii,ij) * this%mass_scaling 198 201 CALL icb_dia_size( ii, ij, zWn, zLn, zAbits, & 199 & this%mass_scaling, zMnew, znMbits, z1_e1e2)202 & this%mass_scaling, zMnew, znMbits, z1_e1e2 ) 200 203 ENDIF 201 204 ! … … 203 206 ! 204 207 END DO 205 208 206 209 ! now use melt and associated heat flux in ocean (or not) 207 210 ! -
branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/ICB/icbtrj.F90
r9019 r9190 1 1 MODULE icbtrj 2 3 2 !!====================================================================== 4 3 !! *** MODULE icbtrj *** 5 4 !! Ocean physics: trajectory I/O routines 6 5 !!====================================================================== 7 !! History : 3.3.1 ! 2010-01 (Martin&Adcroft) Original code 8 !! - ! 2011-03 (Madec) Part conversion to NEMO form 9 !! - ! Removal of mapping from another grid 10 !! - ! 2011-05 (Alderson) New module to handle trajectory output 11 !!---------------------------------------------------------------------- 12 !!---------------------------------------------------------------------- 13 !! icb_trj_init : 6 !! History : 3.3 ! 2010-01 (Martin&Adcroft) Original code 7 !! - ! 2011-03 (Madec) Part conversion to NEMO form 8 !! - ! Removal of mapping from another grid 9 !! - ! 2011-05 (Alderson) New module to handle trajectory output 10 !!---------------------------------------------------------------------- 11 12 !!---------------------------------------------------------------------- 13 !! icb_trj_init : 14 !! icb_trj_write : 15 !! icb_trj_sync : 16 !! icb_trj_end : 14 17 !!---------------------------------------------------------------------- 15 18 USE par_oce ! NEMO parameters … … 49 52 !!---------------------------------------------------------------------- 50 53 CONTAINS 51 52 !!-------------------------------------------------------------------------53 54 54 55 SUBROUTINE icb_trj_init( ktend ) … … 252 253 END SUBROUTINE icb_trj_write 253 254 254 !!-------------------------------------------------------------------------255 255 256 256 SUBROUTINE icb_trj_sync() … … 260 260 !! ** Purpose : 261 261 !!---------------------------------------------------------------------- 262 INTEGER ::iret262 INTEGER :: iret 263 263 !!---------------------------------------------------------------------- 264 264 ! flush to file … … 270 270 271 271 SUBROUTINE icb_trj_end() 272 ! Local variables273 INTEGER ::iret272 !!---------------------------------------------------------------------- 273 INTEGER :: iret 274 274 !!---------------------------------------------------------------------- 275 275 ! Finish up … … 279 279 END SUBROUTINE icb_trj_end 280 280 281 !!------------------------------------------------------------------------- 282 281 !!====================================================================== 283 282 END MODULE icbtrj -
branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/ICB/icbutl.F90
r9019 r9190 9 9 !! - ! 2011-04 (Alderson) Split into separate modules 10 10 !!---------------------------------------------------------------------- 11 11 12 !!---------------------------------------------------------------------- 12 13 !! icb_utl_interp : … … 48 49 !! $Id$ 49 50 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 50 !!------------------------------------------------------------------------- 51 51 !!---------------------------------------------------------------------- 52 52 CONTAINS 53 53 … … 65 65 ! and ssh which is used to calculate gradients 66 66 67 uo_e(:,:) = 0._wp ; uo_e(1:jpi,1:jpj) = ssu_m(:,:) * umask(:,:,1)68 vo_e(:,:) = 0._wp ; vo_e(1:jpi,1:jpj) = ssv_m(:,:) * vmask(:,:,1)69 ff_e(:,:) = 0._wp ; ff_e(1:jpi,1:jpj) = ff_f (:,:)70 tt_e(:,:) = 0._wp ; tt_e(1:jpi,1:jpj) = sst_m(:,:)71 fr_e(:,:) = 0._wp ; fr_e(1:jpi,1:jpj) = fr_i (:,:)72 ua_e(:,:) = 0._wp ; ua_e(1:jpi,1:jpj) = utau (:,:) * umask(:,:,1) ! maybe mask useless because mask applied in sbcblk73 va_e(:,:) = 0._wp ; va_e(1:jpi,1:jpj) = vtau (:,:) * vmask(:,:,1) ! maybe mask useless because mask applied in sbcblk74 67 uo_e(:,:) = 0._wp ; uo_e(1:jpi,1:jpj) = ssu_m(:,:) * umask(:,:,1) 68 vo_e(:,:) = 0._wp ; vo_e(1:jpi,1:jpj) = ssv_m(:,:) * vmask(:,:,1) 69 ff_e(:,:) = 0._wp ; ff_e(1:jpi,1:jpj) = ff_f (:,:) 70 tt_e(:,:) = 0._wp ; tt_e(1:jpi,1:jpj) = sst_m(:,:) 71 fr_e(:,:) = 0._wp ; fr_e(1:jpi,1:jpj) = fr_i (:,:) 72 ua_e(:,:) = 0._wp ; ua_e(1:jpi,1:jpj) = utau (:,:) * umask(:,:,1) ! maybe mask useless because mask applied in sbcblk 73 va_e(:,:) = 0._wp ; va_e(1:jpi,1:jpj) = vtau (:,:) * vmask(:,:,1) ! maybe mask useless because mask applied in sbcblk 74 ! 75 75 CALL lbc_lnk_icb( uo_e, 'U', -1._wp, 1, 1 ) 76 76 CALL lbc_lnk_icb( vo_e, 'V', -1._wp, 1, 1 ) … … 84 84 ui_e(:,:) = 0._wp ; ui_e(1:jpi, 1:jpj) = u_ice(:,:) 85 85 vi_e(:,:) = 0._wp ; vi_e(1:jpi, 1:jpj) = v_ice(:,:) 86 CALL lbc_lnk_icb(hicth, 'T', +1._wp, 1, 1 ) 87 CALL lbc_lnk_icb( ui_e, 'U', -1._wp, 1, 1 ) 88 CALL lbc_lnk_icb( vi_e, 'V', -1._wp, 1, 1 ) 86 ! 87 CALL lbc_lnk_icb( hicth, 'T', +1._wp, 1, 1 ) 88 CALL lbc_lnk_icb( ui_e , 'U', -1._wp, 1, 1 ) 89 CALL lbc_lnk_icb( vi_e , 'V', -1._wp, 1, 1 ) 89 90 #endif 90 91 … … 149 150 150 151 #if defined key_lim3 151 pui = icb_utl_bilin_h( ui_e , pi, pj, 'U' ) ! sea-ice velocities152 pvi = icb_utl_bilin_h( vi_e , pi, pj, 'V' )153 phi = icb_utl_bilin_h( hicth, pi, pj, 'T' ) ! ice thickness152 pui = icb_utl_bilin_h( ui_e , pi, pj, 'U' ) ! sea-ice velocities 153 pvi = icb_utl_bilin_h( vi_e , pi, pj, 'V' ) 154 phi = icb_utl_bilin_h( hicth, pi, pj, 'T' ) ! ice thickness 154 155 #else 155 156 pui = 0._wp … … 160 161 ! Estimate SSH gradient in i- and j-direction (centred evaluation) 161 162 pssh_i = ( icb_utl_bilin_h( ssh_e, pi+0.1_wp, pj, 'T' ) - & 162 &icb_utl_bilin_h( ssh_e, pi-0.1_wp, pj, 'T' ) ) / ( 0.2_wp * pe1 )163 & icb_utl_bilin_h( ssh_e, pi-0.1_wp, pj, 'T' ) ) / ( 0.2_wp * pe1 ) 163 164 pssh_j = ( icb_utl_bilin_h( ssh_e, pi, pj+0.1_wp, 'T' ) - & 164 &icb_utl_bilin_h( ssh_e, pi, pj-0.1_wp, 'T' ) ) / ( 0.2_wp * pe2 )165 & icb_utl_bilin_h( ssh_e, pi, pj-0.1_wp, 'T' ) ) / ( 0.2_wp * pe2 ) 165 166 ! 166 167 END SUBROUTINE icb_utl_interp … … 181 182 REAL(wp) , INTENT(in) :: pi, pj ! targeted coordinates in (i,j) referential 182 183 CHARACTER(len=1) , INTENT(in) :: cd_type ! type of pfld array grid-points: = T , U , V or F points 184 ! 185 INTEGER :: ii, ij ! local integer 186 REAL(wp) :: zi, zj ! local real 187 !!---------------------------------------------------------------------- 188 ! 189 SELECT CASE ( cd_type ) 190 CASE ( 'T' ) 191 ! note that here there is no +0.5 added 192 ! since we're looking for four T points containing quadrant we're in of 193 ! current T cell 194 ii = MAX(1, INT( pi )) 195 ij = MAX(1, INT( pj )) ! T-point 196 zi = pi - REAL(ii,wp) 197 zj = pj - REAL(ij,wp) 198 CASE ( 'U' ) 199 ii = MAX(1, INT( pi-0.5 )) 200 ij = MAX(1, INT( pj )) ! U-point 201 zi = pi - 0.5 - REAL(ii,wp) 202 zj = pj - REAL(ij,wp) 203 CASE ( 'V' ) 204 ii = MAX(1, INT( pi )) 205 ij = MAX(1, INT( pj-0.5 )) ! V-point 206 zi = pi - REAL(ii,wp) 207 zj = pj - 0.5 - REAL(ij,wp) 208 CASE ( 'F' ) 209 ii = MAX(1, INT( pi-0.5 )) 210 ij = MAX(1, INT( pj-0.5 )) ! F-point 211 zi = pi - 0.5 - REAL(ii,wp) 212 zj = pj - 0.5 - REAL(ij,wp) 213 END SELECT 214 ! 215 ! find position in this processor. Prevent near edge problems (see #1389) 216 ! 217 IF ( ii < mig( 1 ) ) THEN ; ii = 1 218 ELSEIF( ii > mig(jpi) ) THEN ; ii = jpi 219 ELSE ; ii = mi1(ii) 220 ENDIF 221 IF ( ij < mjg( 1 ) ) THEN ; ij = 1 222 ELSEIF( ij > mjg(jpj) ) THEN ; ij = jpj 223 ELSE ; ij = mj1(ij) 224 ENDIF 225 ! 226 IF( ii == jpi ) ii = ii-1 227 IF( ij == jpj ) ij = ij-1 228 ! 229 icb_utl_bilin_h = ( pfld(ii,ij ) * (1.-zi) + pfld(ii+1,ij ) * zi ) * (1.-zj) & 230 & + ( pfld(ii,ij+1) * (1.-zi) + pfld(ii+1,ij+1) * zi ) * zj 231 ! 232 END FUNCTION icb_utl_bilin_h 233 234 235 REAL(wp) FUNCTION icb_utl_bilin( pfld, pi, pj, cd_type ) 236 !!---------------------------------------------------------------------- 237 !! *** FUNCTION icb_utl_bilin *** 238 !! 239 !! ** Purpose : bilinear interpolation at berg location depending on the grid-point type 240 !! 241 !! !!gm CAUTION an optional argument should be added to handle 242 !! the slip/no-slip conditions ==>>> to be done later 243 !! 244 !!---------------------------------------------------------------------- 245 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pfld ! field to be interpolated 246 REAL(wp) , INTENT(in) :: pi, pj ! targeted coordinates in (i,j) referential 247 CHARACTER(len=1) , INTENT(in) :: cd_type ! type of pfld array grid-points: = T , U , V or F points 183 248 ! 184 249 INTEGER :: ii, ij ! local integer … … 213 278 ! 214 279 ! find position in this processor. Prevent near edge problems (see #1389) 215 216 if (ii.lt.mig(1)) then 217 ii = 1 218 else if (ii.gt.mig(jpi)) then 219 ii = jpi 220 else 221 ii = mi1( ii ) 222 end if 223 224 if (ij.lt.mjg(1)) then 225 ij = 1 226 else if (ij.gt.mjg(jpj)) then 227 ij = jpj 228 else 229 ij = mj1( ij ) 230 end if 231 232 if (ij.eq.jpj) ij=ij-1 233 if (ii.eq.jpi) ii=ii-1 234 235 ! 236 icb_utl_bilin_h = ( pfld(ii,ij ) * (1.-zi) + pfld(ii+1,ij ) * zi ) * (1.-zj) & 237 & + ( pfld(ii,ij+1) * (1.-zi) + pfld(ii+1,ij+1) * zi ) * zj 238 ! 239 END FUNCTION icb_utl_bilin_h 240 241 242 REAL(wp) FUNCTION icb_utl_bilin( pfld, pi, pj, cd_type ) 243 !!---------------------------------------------------------------------- 244 !! *** FUNCTION icb_utl_bilin *** 245 !! 246 !! ** Purpose : bilinear interpolation at berg location depending on the grid-point type 247 !! 248 !! !!gm CAUTION an optional argument should be added to handle 249 !! the slip/no-slip conditions ==>>> to be done later 250 !! 251 !!---------------------------------------------------------------------- 252 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pfld ! field to be interpolated 253 REAL(wp) , INTENT(in) :: pi, pj ! targeted coordinates in (i,j) referential 254 CHARACTER(len=1) , INTENT(in) :: cd_type ! type of pfld array grid-points: = T , U , V or F points 255 ! 256 INTEGER :: ii, ij ! local integer 257 REAL(wp) :: zi, zj ! local real 258 !!---------------------------------------------------------------------- 259 ! 260 SELECT CASE ( cd_type ) 261 CASE ( 'T' ) 262 ! note that here there is no +0.5 added 263 ! since we're looking for four T points containing quadrant we're in of 264 ! current T cell 265 ii = MAX(1, INT( pi )) 266 ij = MAX(1, INT( pj )) ! T-point 267 zi = pi - REAL(ii,wp) 268 zj = pj - REAL(ij,wp) 269 CASE ( 'U' ) 270 ii = MAX(1, INT( pi-0.5 )) 271 ij = MAX(1, INT( pj )) ! U-point 272 zi = pi - 0.5 - REAL(ii,wp) 273 zj = pj - REAL(ij,wp) 274 CASE ( 'V' ) 275 ii = MAX(1, INT( pi )) 276 ij = MAX(1, INT( pj-0.5 )) ! V-point 277 zi = pi - REAL(ii,wp) 278 zj = pj - 0.5 - REAL(ij,wp) 279 CASE ( 'F' ) 280 ii = MAX(1, INT( pi-0.5 )) 281 ij = MAX(1, INT( pj-0.5 )) ! F-point 282 zi = pi - 0.5 - REAL(ii,wp) 283 zj = pj - 0.5 - REAL(ij,wp) 284 END SELECT 285 ! 286 ! find position in this processor. Prevent near edge problems (see #1389) 287 288 if (ii.lt.mig(1)) then 289 ii = 1 290 else if (ii.gt.mig(jpi)) then 291 ii = jpi 292 else 293 ii = mi1( ii ) 294 end if 295 296 if (ij.lt.mjg(1)) then 297 ij = 1 298 else if (ij.gt.mjg(jpj)) then 299 ij = jpj 300 else 301 ij = mj1( ij ) 302 end if 303 304 if (ij.eq.jpj) ij=ij-1 305 if (ii.eq.jpi) ii=ii-1 306 280 IF ( ii < mig( 1 ) ) THEN ; ii = 1 281 ELSEIF( ii > mig(jpi) ) THEN ; ii = jpi 282 ELSE ; ii = mi1(ii) 283 ENDIF 284 IF ( ij < mjg( 1 ) ) THEN ; ij = 1 285 ELSEIF( ij > mjg(jpj) ) THEN ; ij = jpj 286 ELSE ; ij = mj1(ij) 287 ENDIF 288 ! 289 IF( ii == jpi ) ii = ii-1 290 IF( ij == jpj ) ij = ij-1 291 ! 307 292 icb_utl_bilin = ( pfld(ii,ij ) * (1.-zi) + pfld(ii+1,ij ) * zi ) * (1.-zj) & 308 293 & + ( pfld(ii,ij+1) * (1.-zi) + pfld(ii+1,ij+1) * zi ) * zj … … 340 325 ! 341 326 ! find position in this processor. Prevent near edge problems (see #1389) 342 343 if (ii.lt.mig(1)) then 344 ii = 1 345 else if (ii.gt.mig(jpi)) then 346 ii = jpi 347 else 348 ii = mi1( ii ) 349 end if 350 351 if (ij.lt.mjg(1)) then 352 ij = 1 353 else if (ij.gt.mjg(jpj)) then 354 ij = jpj 355 else 356 ij = mj1( ij ) 357 end if 358 359 if (ij.eq.jpj) ij=ij-1 360 if (ii.eq.jpi) ii=ii-1 361 327 IF ( ii < mig( 1 ) ) THEN ; ii = 1 328 ELSEIF( ii > mig(jpi) ) THEN ; ii = jpi 329 ELSE ; ii = mi1(ii) 330 ENDIF 331 IF ( ij < mjg( 1 ) ) THEN ; ij = 1 332 ELSEIF( ij > mjg(jpj) ) THEN ; ij = jpj 333 ELSE ; ij = mj1(ij) 334 ENDIF 335 ! 336 IF( ii == jpi ) ii = ii-1 337 IF( ij == jpj ) ij = ij-1 338 ! 362 339 z4(1) = pfld(ii ,ij ) 363 340 z4(2) = pfld(ii+1,ij ) … … 408 385 409 386 ! find position in this processor. Prevent near edge problems (see #1389) 410 411 if (ii.lt.mig(1)) then 412 ii = 1 413 else if (ii.gt.mig(jpi)) then 414 ii = jpi 415 else 416 ii = mi1( ii ) 417 end if 418 419 if (ij.lt.mjg(1)) then 420 ij = 1 421 else if (ij.gt.mjg(jpj)) then 422 ij = jpj 423 else 424 ij = mj1( ij ) 425 end if 426 427 if (ij.eq.jpj) ij=ij-1 428 if (ii.eq.jpi) ii=ii-1 429 387 IF ( ii < mig( 1 ) ) THEN ; ii = 1 388 ELSEIF( ii > mig(jpi) ) THEN ; ii = jpi 389 ELSE ; ii = mi1(ii) 390 ENDIF 391 IF ( ij < mjg( 1 ) ) THEN ; ij = 1 392 ELSEIF( ij > mjg(jpj) ) THEN ; ij = jpj 393 ELSE ; ij = mj1(ij) 394 ENDIF 395 ! 396 IF( ii == jpi ) ii = ii-1 397 IF( ij == jpj ) ij = ij-1 398 ! 430 399 IF( 0.0_wp <= zi .AND. zi < 0.5_wp ) THEN 431 400 IF( 0.0_wp <= zj .AND. zj < 0.5_wp ) THEN ! NE quadrant
Note: See TracChangeset
for help on using the changeset viewer.
