Changeset 7698 for trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk.F90
- Timestamp:
- 2017-02-18T10:02:03+01:00 (7 years ago)
- File:
-
- 1 edited
-
trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk.F90 (modified) (25 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk.F90
r7646 r7698 316 316 #if defined key_cice 317 317 IF( MOD( kt - 1, nn_fsbc ) == 0 ) THEN 318 qlw_ice(:,:,1) = sf(jp_qlw )%fnow(:,:,1) 319 IF( ln_dm2dc ) THEN ; qsr_ice(:,:,1) = sbc_dcy( sf(jp_qsr)%fnow(:,:,1) ) 320 ELSE ; qsr_ice(:,:,1) = sf(jp_qsr)%fnow(:,:,1) 318 !$OMP PARALLEL DO schedule(static) private(jj, ji) 319 DO jj = 1, jpj 320 DO ji = 1, jpi 321 qlw_ice(ji,jj,1) = sf(jp_qlw)%fnow(ji,jj,1) 322 END DO 323 END DO 324 IF( ln_dm2dc ) THEN ; qsr_ice(:,:,1) = sbc_dcy( sf(jp_qsr)%fnow(:,:,1) ) 325 ELSE 326 !$OMP PARALLEL DO schedule(static) private(jj, ji) 327 DO jj = 1, jpj 328 DO ji = 1, jpi 329 qsr_ice(ji,jj,1) = sf(jp_qsr)%fnow(ji,jj,1) 330 END DO 331 END DO 321 332 ENDIF 322 tatm_ice(:,:) = sf(jp_tair)%fnow(:,:,1) 323 qatm_ice(:,:) = sf(jp_humi)%fnow(:,:,1) 324 tprecip(:,:) = sf(jp_prec)%fnow(:,:,1) * rn_pfac 325 sprecip(:,:) = sf(jp_snow)%fnow(:,:,1) * rn_pfac 326 wndi_ice(:,:) = sf(jp_wndi)%fnow(:,:,1) 327 wndj_ice(:,:) = sf(jp_wndj)%fnow(:,:,1) 333 !$OMP PARALLEL DO schedule(static) private(jj, ji) 334 DO jj = 1, jpj 335 DO ji = 1, jpi 336 tatm_ice(ji,jj) = sf(jp_tair)%fnow(ji,jj,1) 337 qatm_ice(ji,jj) = sf(jp_humi)%fnow(ji,jj,1) 338 tprecip(ji,jj) = sf(jp_prec)%fnow(ji,jj,1) * rn_pfac 339 sprecip(ji,jj) = sf(jp_snow)%fnow(ji,jj,1) * rn_pfac 340 wndi_ice(ji,jj) = sf(jp_wndi)%fnow(ji,jj,1) 341 wndj_ice(ji,jj) = sf(jp_wndj)%fnow(ji,jj,1) 342 END DO 343 END DO 328 344 ENDIF 329 345 #endif … … 382 398 ! 383 399 384 ! local scalars ( place there for vector optimisation purposes) 385 zst(:,:) = pst(:,:) + rt0 ! convert SST from Celcius to Kelvin (and set minimum value far above 0 K) 386 400 !$OMP PARALLEL DO schedule(static) private(jj, ji) 401 DO jj = 1, jpj 402 DO ji = 1, jpi 403 ! local scalars ( place there for vector optimisation purposes) 404 zst(ji,jj) = pst(ji,jj) + rt0 ! convert SST from Celcius to Kelvin (and set minimum value far above 0 K) 405 406 ! ... components ( U10m - U_oce ) at T-point (unmasked) 407 !!gm move zwnd_i (_j) set to zero inside the key_cyclone ??? 408 zwnd_i(ji,jj) = 0._wp 409 zwnd_j(ji,jj) = 0._wp 410 END DO 411 END DO 387 412 ! ----------------------------------------------------------------------------- ! 388 413 ! 0 Wind components and module at T-point relative to the moving ocean ! 389 414 ! ----------------------------------------------------------------------------- ! 390 415 391 ! ... components ( U10m - U_oce ) at T-point (unmasked)392 !!gm move zwnd_i (_j) set to zero inside the key_cyclone ???393 zwnd_i(:,:) = 0._wp394 zwnd_j(:,:) = 0._wp395 416 #if defined key_cyclone 396 417 CALL wnd_cyc( kt, zwnd_i, zwnd_j ) ! add analytical tropical cyclone (Vincent et al. JGR 2012) 418 !$OMP PARALLEL DO schedule(static) private(jj, ji) 397 419 DO jj = 2, jpjm1 398 420 DO ji = fs_2, fs_jpim1 ! vect. opt. … … 402 424 END DO 403 425 #endif 426 !$OMP PARALLEL DO schedule(static) private(jj, ji) 404 427 DO jj = 2, jpjm1 405 428 DO ji = fs_2, fs_jpim1 ! vect. opt. … … 411 434 CALL lbc_lnk( zwnd_j(:,:) , 'T', -1. ) 412 435 ! ... scalar wind ( = | U10m - U_oce | ) at T-point (masked) 413 wndm(:,:) = SQRT( zwnd_i(:,:) * zwnd_i(:,:) & 414 & + zwnd_j(:,:) * zwnd_j(:,:) ) * tmask(:,:,1) 415 436 !$OMP PARALLEL DO schedule(static) private(jj, ji) 437 DO jj = 1, jpj 438 DO ji = 1, jpi 439 wndm(ji,jj) = SQRT( zwnd_i(ji,jj) * zwnd_i(ji,jj) & 440 & + zwnd_j(ji,jj) * zwnd_j(ji,jj) ) * tmask(ji,jj,1) 441 442 END DO 443 END DO 416 444 ! ----------------------------------------------------------------------------- ! 417 445 ! I Radiative FLUXES ! … … 421 449 zztmp = 1. - albo 422 450 IF( ln_dm2dc ) THEN ; qsr(:,:) = zztmp * sbc_dcy( sf(jp_qsr)%fnow(:,:,1) ) * tmask(:,:,1) 423 ELSE ; qsr(:,:) = zztmp * sf(jp_qsr)%fnow(:,:,1) * tmask(:,:,1) 451 ELSE 452 !$OMP PARALLEL DO schedule(static) private(jj, ji) 453 DO jj = 1, jpj 454 DO ji = 1, jpi 455 qsr(ji,jj) = zztmp * sf(jp_qsr)%fnow(ji,jj,1) * tmask(ji,jj,1) 456 END DO 457 END DO 424 458 ENDIF 425 459 426 zqlw(:,:) = ( sf(jp_qlw)%fnow(:,:,1) - Stef * zst(:,:)*zst(:,:)*zst(:,:)*zst(:,:) ) * tmask(:,:,1) ! Long Wave 460 !$OMP PARALLEL DO schedule(static) private(jj, ji) 461 DO jj = 1, jpj 462 DO ji = 1, jpi 463 zqlw(ji,jj) = ( sf(jp_qlw)%fnow(ji,jj,1) - Stef * zst(ji,jj)*zst(ji,jj)*zst(ji,jj)*zst(ji,jj) ) * tmask(ji,jj,1) ! Long Wave 464 END DO 465 END DO 427 466 428 467 … … 461 500 END IF 462 501 463 Cd_oce(:,:) = Cd(:,:) ! record value of pure ocean-atm. drag (clem) 464 502 !$OMP PARALLEL 503 !$OMP DO schedule(static) private(jj, ji) 504 DO jj = 1, jpj 505 DO ji = 1, jpi 506 Cd_oce(ji,jj) = Cd(ji,jj) ! record value of pure ocean-atm. drag (clem) 507 END DO 508 END DO 509 510 !$OMP DO schedule(static) private(jj, ji) 465 511 DO jj = 1, jpj ! tau module, i and j component 466 512 DO ji = 1, jpi … … 471 517 END DO 472 518 END DO 519 !$OMP END PARALLEL 473 520 474 521 ! ! add the HF tau contribution to the wind stress module 475 IF( lhftau ) taum(:,:) = taum(:,:) + sf(jp_tdif)%fnow(:,:,1) 522 IF( lhftau ) THEN 523 !$OMP PARALLEL DO schedule(static) private(jj, ji) 524 DO jj = 1, jpj 525 DO ji = 1, jpi 526 taum(ji,jj) = taum(ji,jj) + sf(jp_tdif)%fnow(ji,jj,1) 527 END DO 528 END DO 529 END IF 476 530 477 531 CALL iom_put( "taum_oce", taum ) ! output wind stress module … … 480 534 ! Note the use of 0.5*(2-umask) in order to unmask the stress along coastlines 481 535 ! Note the use of MAX(tmask(i,j),tmask(i+1,j) is to mask tau over ice shelves 536 !$OMP PARALLEL DO schedule(static) private(jj, ji) 482 537 DO jj = 1, jpjm1 483 538 DO ji = 1, fs_jpim1 … … 496 551 497 552 ! zqla used as temporary array, for rho*U (common term of bulk formulae): 498 zqla(:,:) = zrhoa(:,:) * zU_zu(:,:) 553 !$OMP PARALLEL DO schedule(static) private(jj, ji) 554 DO jj = 1, jpj 555 DO ji = 1, jpi 556 zqla(ji,jj) = zrhoa(ji,jj) * zU_zu(ji,jj) 557 END DO 558 END DO 499 559 500 560 IF( ABS( rn_zu - rn_zqt) < 0.01_wp ) THEN 501 561 !! q_air and t_air are given at 10m (wind reference height) 502 zevap(:,:) = rn_efac*MAX( 0._wp, zqla(:,:)*Ce(:,:)*(zsq(:,:) - sf(jp_humi)%fnow(:,:,1)) ) ! Evaporation, using bulk wind speed 503 zqsb (:,:) = cp_air(sf(jp_humi)%fnow(:,:,1))*zqla(:,:)*Ch(:,:)*(zst(:,:) - ztpot(:,:) ) ! Sensible Heat, using bulk wind speed 562 !$OMP PARALLEL DO schedule(static) private(jj, ji) 563 DO jj = 1, jpj 564 DO ji = 1, jpi 565 zevap(ji,jj) = rn_efac*MAX( 0._wp, zqla(ji,jj)*Ce(ji,jj)*(zsq(ji,jj) - sf(jp_humi)%fnow(ji,jj,1)) ) ! Evaporation, using bulk wind speed 566 END DO 567 END DO 568 zqsb (:,:) = cp_air(sf(jp_humi)%fnow(:,:,1))*zqla(:,:)*Ch(:,:)*(zst(:,:) - ztpot(:,:) ) ! Sensible Heat, using bulk wind speed 504 569 ELSE 505 570 !! q_air and t_air are not given at 10m (wind reference height) 506 571 ! Values of temp. and hum. adjusted to height of wind during bulk algorithm iteration must be used!!! 507 zevap(:,:) = rn_efac*MAX( 0._wp, zqla(:,:)*Ce(:,:)*(zsq(:,:) - zq_zu(:,:) ) ) ! Evaporation ! using bulk wind speed 572 !$OMP PARALLEL DO schedule(static) private(jj, ji) 573 DO jj = 1, jpj 574 DO ji = 1, jpi 575 zevap(ji,jj) = rn_efac*MAX( 0._wp, zqla(ji,jj)*Ce(ji,jj)*(zsq(ji,jj) - zq_zu(ji,jj) ) ) ! Evaporation ! using bulk wind speed 576 END DO 577 END DO 508 578 zqsb (:,:) = cp_air(sf(jp_humi)%fnow(:,:,1))*zqla(:,:)*Ch(:,:)*(zst(:,:) - zt_zu(:,:) ) ! Sensible Heat ! using bulk wind speed 509 579 ENDIF … … 527 597 ! ----------------------------------------------------------------------------- ! 528 598 ! 529 emp (:,:) = ( zevap(:,:) & ! mass flux (evap. - precip.) 530 & - sf(jp_prec)%fnow(:,:,1) * rn_pfac ) * tmask(:,:,1) 531 ! 532 qns(:,:) = zqlw(:,:) - zqsb(:,:) - zqla(:,:) & ! Downward Non Solar 533 & - sf(jp_snow)%fnow(:,:,1) * rn_pfac * lfus & ! remove latent melting heat for solid precip 534 & - zevap(:,:) * pst(:,:) * rcp & ! remove evap heat content at SST 535 & + ( sf(jp_prec)%fnow(:,:,1) - sf(jp_snow)%fnow(:,:,1) ) * rn_pfac & ! add liquid precip heat content at Tair 536 & * ( sf(jp_tair)%fnow(:,:,1) - rt0 ) * rcp & 537 & + sf(jp_snow)%fnow(:,:,1) * rn_pfac & ! add solid precip heat content at min(Tair,Tsnow) 538 & * ( MIN( sf(jp_tair)%fnow(:,:,1), rt0_snow ) - rt0 ) * cpic * tmask(:,:,1) 539 ! 599 !$OMP PARALLEL DO schedule(static) private(jj, ji) 600 DO jj = 1, jpj 601 DO ji = 1, jpi 602 emp (ji,jj) = ( zevap(ji,jj) & ! mass flux (evap. - precip.) 603 & - sf(jp_prec)%fnow(ji,jj,1) * rn_pfac ) * tmask(ji,jj,1) 604 ! 605 qns(ji,jj) = zqlw(ji,jj) - zqsb(ji,jj) - zqla(ji,jj) & ! Downward Non Solar 606 & - sf(jp_snow)%fnow(ji,jj,1) * rn_pfac * lfus & ! remove latent melting heat for solid precip 607 & - zevap(ji,jj) * pst(ji,jj) * rcp & ! remove evap heat content at SST 608 & + ( sf(jp_prec)%fnow(ji,jj,1) - sf(jp_snow)%fnow(ji,jj,1) ) * rn_pfac & ! add liquid precip heat content at Tair 609 & * ( sf(jp_tair)%fnow(ji,jj,1) - rt0 ) * rcp & 610 & + sf(jp_snow)%fnow(ji,jj,1) * rn_pfac & ! add solid precip heat content at min(Tair,Tsnow) 611 & * ( MIN( sf(jp_tair)%fnow(ji,jj,1), rt0_snow ) - rt0 ) * cpic * tmask(ji,jj,1) 612 ! 540 613 #if defined key_lim3 541 qns_oce(:,:) = zqlw(:,:) - zqsb(:,:) - zqla(:,:) ! non solar without emp (only needed by LIM3)542 qsr_oce(:,:) = qsr(:,:)614 qns_oce(ji,jj) = zqlw(ji,jj) - zqsb(ji,jj) - zqla(ji,jj) ! non solar without emp (only needed by LIM3) 615 qsr_oce(ji,jj) = qsr(ji,jj) 543 616 #endif 617 END DO 618 END DO 544 619 ! 545 620 IF ( nn_ice == 0 ) THEN … … 551 626 CALL iom_put( "qsr_oce" , qsr ) ! output downward solar heat over the ocean 552 627 CALL iom_put( "qt_oce" , qns+qsr ) ! output total downward heat over the ocean 553 tprecip(:,:) = sf(jp_prec)%fnow(:,:,1) * rn_pfac ! output total precipitation [kg/m2/s] 554 sprecip(:,:) = sf(jp_snow)%fnow(:,:,1) * rn_pfac ! output solid precipitation [kg/m2/s] 628 !$OMP PARALLEL DO schedule(static) private(jj, ji) 629 DO jj = 1, jpj 630 DO ji = 1, jpi 631 tprecip(ji,jj) = sf(jp_prec)%fnow(ji,jj,1) * rn_pfac ! output total precipitation [kg/m2/s] 632 sprecip(ji,jj) = sf(jp_snow)%fnow(ji,jj,1) * rn_pfac ! output solid precipitation [kg/m2/s] 633 END DO 634 END DO 555 635 CALL iom_put( 'snowpre', sprecip * 86400. ) ! Snow 556 636 CALL iom_put( 'precip' , tprecip * 86400. ) ! Total precipitation … … 599 679 CALL wrk_alloc( jpi,jpj, Cd ) 600 680 601 Cd(:,:) = Cd_ice 681 !$OMP PARALLEL DO schedule(static) private(jj, ji) 682 DO jj = 1, jpj 683 DO ji = 1, jpi 684 Cd(ji,jj) = Cd_ice 685 END DO 686 END DO 602 687 603 688 ! Make ice-atm. drag dependent on ice concentration (see Lupkes et al. 2012) (clem) … … 613 698 zrhoa (:,:) = rho_air(sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1)) 614 699 615 !!gm brutal.... 616 utau_ice (:,:) = 0._wp 617 vtau_ice (:,:) = 0._wp 618 wndm_ice (:,:) = 0._wp 619 !!gm end 700 !$OMP PARALLEL DO schedule(static) private(jj, ji) 701 DO jj = 1, jpj 702 DO ji = 1, jpi 703 !!gm brutal.... 704 utau_ice (ji,jj) = 0._wp 705 vtau_ice (ji,jj) = 0._wp 706 wndm_ice (ji,jj) = 0._wp 707 !!gm end 708 END DO 709 END DO 620 710 621 711 ! ----------------------------------------------------------------------------- ! … … 625 715 CASE( 'I' ) ! B-grid ice dynamics : I-point (i.e. F-point with sea-ice indexation) 626 716 ! and scalar wind at T-point ( = | U10m - U_ice | ) (masked) 717 !$OMP PARALLEL DO schedule(static) private(jj,ji,zwndi_f,zwndj_f,zwnorm_f,zwndi_t,zwndj_t) 627 718 DO jj = 2, jpjm1 628 719 DO ji = 2, jpim1 ! B grid : NO vector opt … … 649 740 ! 650 741 CASE( 'C' ) ! C-grid ice dynamics : U & V-points (same as ocean) 742 !$OMP PARALLEL DO schedule(static) private(jj,ji,zwndi_t,zwndj_t) 651 743 DO jj = 2, jpj 652 744 DO ji = fs_2, jpi ! vect. opt. … … 656 748 END DO 657 749 END DO 750 !$OMP PARALLEL DO schedule(static) private(jj,ji) 658 751 DO jj = 2, jpjm1 659 752 DO ji = fs_2, fs_jpim1 ! vect. opt. … … 700 793 REAL(wp) :: zztmp, z1_lsub ! - - 701 794 REAL(wp), DIMENSION(:,:,:), POINTER :: z_qlw ! long wave heat flux over ice 795 REAL(wp), DIMENSION(:,:,:), POINTER :: zevap_ice3d, zqns_ice3d, zqsr_ice3d 702 796 REAL(wp), DIMENSION(:,:,:), POINTER :: z_qsb ! sensible heat flux over ice 703 797 REAL(wp), DIMENSION(:,:,:), POINTER :: z_dqlw ! long wave heat sensitivity over ice 704 798 REAL(wp), DIMENSION(:,:,:), POINTER :: z_dqsb ! sensible heat sensitivity over ice 705 799 REAL(wp), DIMENSION(:,:) , POINTER :: zevap, zsnw ! evaporation and snw distribution after wind blowing (LIM3) 800 REAL(wp), DIMENSION(:,:) , POINTER :: zevap_ice2d, zqns_ice2d, zqsr_ice2d 706 801 REAL(wp), DIMENSION(:,:) , POINTER :: zrhoa 707 802 REAL(wp), DIMENSION(:,:) , POINTER :: Cd ! transfer coefficient for momentum (tau) … … 710 805 IF( nn_timing == 1 ) CALL timing_start('blk_ice_flx') 711 806 ! 712 CALL wrk_alloc( jpi,jpj,jpl, z_qlw, z_qsb, z_dqlw, z_dqsb )713 CALL wrk_alloc( jpi,jpj, zrhoa )807 CALL wrk_alloc( jpi,jpj,jpl, z_qlw, z_qsb, z_dqlw, z_dqsb, zevap_ice3d, zqns_ice3d, zqsr_ice3d ) 808 CALL wrk_alloc( jpi,jpj, zrhoa, zevap_ice2d, zqns_ice2d, zqsr_ice2d) 714 809 CALL wrk_alloc( jpi,jpj, Cd ) 715 810 716 Cd(:,:) = Cd_ice 811 !$OMP PARALLEL DO schedule(static) private(jj, ji) 812 DO jj = 1, jpj 813 DO ji = 1, jpi 814 Cd(ji,jj) = Cd_ice 815 END DO 816 END DO 717 817 718 818 ! Make ice-atm. drag dependent on ice concentration (see Lupkes et al. 2012) (clem) … … 731 831 ! 732 832 zztmp = 1. / ( 1. - albo ) 733 ! ! ========================== ! 734 DO jl = 1, jpl ! Loop over ice categories ! 735 ! ! ========================== ! 833 !$OMP PARALLEL 834 !$OMP DO schedule(static) private(jl,jj,ji,zst2,zst3) ! ========================== ! 835 DO jl = 1, jpl ! Loop over ice categories ! 836 ! ! ========================== ! 736 837 DO jj = 1 , jpj 737 838 DO ji = 1, jpi … … 781 882 END DO 782 883 ! 783 tprecip(:,:) = sf(jp_prec)%fnow(:,:,1) * rn_pfac ! total precipitation [kg/m2/s] 784 sprecip(:,:) = sf(jp_snow)%fnow(:,:,1) * rn_pfac ! solid precipitation [kg/m2/s] 884 !$OMP DO schedule(static) private(jj, ji) 885 DO jj = 1, jpj 886 DO ji = 1, jpi 887 tprecip(ji,jj) = sf(jp_prec)%fnow(ji,jj,1) * rn_pfac ! total precipitation [kg/m2/s] 888 sprecip(ji,jj) = sf(jp_snow)%fnow(ji,jj,1) * rn_pfac ! solid precipitation [kg/m2/s] 889 END DO 890 END DO 891 !$OMP END PARALLEL 785 892 CALL iom_put( 'snowpre', sprecip * 86400. ) ! Snow precipitation 786 893 CALL iom_put( 'precip' , tprecip * 86400. ) ! Total precipitation … … 791 898 ! --- evaporation --- ! 792 899 z1_lsub = 1._wp / Lsub 793 evap_ice (:,:,:) = rn_efac * qla_ice (:,:,:) * z1_lsub ! sublimation 794 devap_ice(:,:,:) = rn_efac * dqla_ice(:,:,:) * z1_lsub ! d(sublimation)/dT 795 zevap (:,:) = rn_efac * ( emp(:,:) + tprecip(:,:) ) ! evaporation over ocean 796 797 ! --- evaporation minus precipitation --- ! 798 zsnw(:,:) = 0._wp 900 !$OMP PARALLEL 901 !$OMP DO schedule(static) private(jl,jj,ji) 902 DO jl = 1, jpl 903 DO jj = 1 , jpj 904 DO ji = 1, jpi 905 evap_ice (ji,jj,jl) = rn_efac * qla_ice (ji,jj,jl) * z1_lsub ! sublimation 906 devap_ice(ji,jj,jl) = rn_efac * dqla_ice(ji,jj,jl) * z1_lsub ! d(sublimation)/dT 907 END DO 908 END DO 909 END DO 910 ! 911 !$OMP DO schedule(static) private(jj, ji) 912 DO jj = 1, jpj 913 DO ji = 1, jpi 914 zevap (ji,jj) = rn_efac * ( emp(ji,jj) + tprecip(ji,jj) ) ! evaporation over ocean 915 916 ! --- evaporation minus precipitation --- ! 917 zsnw(ji,jj) = 0._wp 918 END DO 919 END DO 920 !$OMP END PARALLEL 799 921 CALL lim_thd_snwblow( pfrld, zsnw ) ! snow distribution over ice after wind blowing 800 emp_oce(:,:) = pfrld(:,:) * zevap(:,:) - ( tprecip(:,:) - sprecip(:,:) ) - sprecip(:,:) * (1._wp - zsnw ) 801 emp_ice(:,:) = SUM( a_i_b(:,:,:) * evap_ice(:,:,:), dim=3 ) - sprecip(:,:) * zsnw 802 emp_tot(:,:) = emp_oce(:,:) + emp_ice(:,:) 803 804 ! --- heat flux associated with emp --- ! 805 qemp_oce(:,:) = - pfrld(:,:) * zevap(:,:) * sst_m(:,:) * rcp & ! evap at sst 806 & + ( tprecip(:,:) - sprecip(:,:) ) * ( sf(jp_tair)%fnow(:,:,1) - rt0 ) * rcp & ! liquid precip at Tair 807 & + sprecip(:,:) * ( 1._wp - zsnw ) * & ! solid precip at min(Tair,Tsnow) 808 & ( ( MIN( sf(jp_tair)%fnow(:,:,1), rt0_snow ) - rt0 ) * cpic * tmask(:,:,1) - lfus ) 809 qemp_ice(:,:) = sprecip(:,:) * zsnw * & ! solid precip (only) 810 & ( ( MIN( sf(jp_tair)%fnow(:,:,1), rt0_snow ) - rt0 ) * cpic * tmask(:,:,1) - lfus ) 811 812 ! --- total solar and non solar fluxes --- ! 813 qns_tot(:,:) = pfrld(:,:) * qns_oce(:,:) + SUM( a_i_b(:,:,:) * qns_ice(:,:,:), dim=3 ) + qemp_ice(:,:) + qemp_oce(:,:) 814 qsr_tot(:,:) = pfrld(:,:) * qsr_oce(:,:) + SUM( a_i_b(:,:,:) * qsr_ice(:,:,:), dim=3 ) 815 816 ! --- heat content of precip over ice in J/m3 (to be used in 1D-thermo) --- ! 817 qprec_ice(:,:) = rhosn * ( ( MIN( sf(jp_tair)%fnow(:,:,1), rt0_snow ) - rt0 ) * cpic * tmask(:,:,1) - lfus ) 922 !$OMP PARALLEL 923 !$OMP DO schedule(static) private(jj,ji) 924 DO jj = 1, jpj 925 DO ji = 1, jpi 926 emp_oce(ji,jj) = pfrld(ji,jj) * zevap(ji,jj) - ( tprecip(ji,jj) - sprecip(ji,jj) ) - sprecip(ji,jj) * (1._wp - zsnw(ji,jj)) 927 END DO 928 END DO 929 !$OMP END DO NOWAIT 930 !$OMP DO schedule(static) private(jl,jj,ji) 931 DO jl = 1, jpl 932 DO jj = 1 , jpj 933 DO ji = 1, jpi 934 zevap_ice3d(ji,jj,jl) = a_i_b(ji,jj,jl) * evap_ice(ji,jj,jl) 935 zqns_ice3d(ji,jj,jl) = a_i_b(ji,jj,jl) * qns_ice(ji,jj,jl) 936 zqsr_ice3d(ji,jj,jl) = a_i_b(ji,jj,jl) * qsr_ice(ji,jj,jl) 937 END DO 938 END DO 939 END DO 940 !$OMP END DO NOWAIT 941 !$OMP DO schedule(static) private(jj,ji) 942 DO jj = 1, jpj 943 DO ji = 1, jpi 944 zevap_ice2d(ji,jj) = 0._wp 945 zqns_ice2d(ji,jj) = 0._wp 946 zqsr_ice2d(ji,jj) = 0._wp 947 END DO 948 END DO 949 DO jl = 1, jpl 950 !$OMP DO schedule(static) private(jj,ji) 951 DO jj = 1 , jpj 952 DO ji = 1, jpi 953 zevap_ice2d(ji,jj) = zevap_ice2d(ji,jj) + zevap_ice3d(ji,jj,jl) 954 zqns_ice2d(ji,jj) = zqns_ice2d(ji,jj) + zqns_ice3d(ji,jj,jl) 955 zqsr_ice2d(ji,jj) = zqsr_ice2d(ji,jj) + zqsr_ice3d(ji,jj,jl) 956 END DO 957 END DO 958 END DO 959 !$OMP DO schedule(static) private(jj,ji) 960 DO jj = 1 , jpj 961 DO ji = 1, jpi 962 emp_ice(ji,jj) = zevap_ice2d(ji,jj) - sprecip(ji,jj) * zsnw(ji,jj) 963 emp_tot(ji,jj) = emp_oce(ji,jj) + emp_ice(ji,jj) 964 965 ! --- heat flux associated with emp --- ! 966 qemp_oce(ji,jj) = - pfrld(ji,jj) * zevap(ji,jj) * sst_m(ji,jj) * rcp & ! evap at sst 967 & + ( tprecip(ji,jj) - sprecip(ji,jj) ) * ( sf(jp_tair)%fnow(ji,jj,1) - rt0 ) * rcp & ! liquid precip at Tair 968 & + sprecip(ji,jj) * ( 1._wp - zsnw(ji,jj) ) * & ! solid precip at min(Tair,Tsnow) 969 & ( ( MIN( sf(jp_tair)%fnow(ji,jj,1), rt0_snow ) - rt0 ) * cpic * tmask(ji,jj,1) - lfus ) 970 qemp_ice(ji,jj) = sprecip(ji,jj) * zsnw(ji,jj) * & ! solid precip (only) 971 & ( ( MIN( sf(jp_tair)%fnow(ji,jj,1), rt0_snow ) - rt0 ) * cpic * tmask(ji,jj,1) - lfus ) 972 973 ! --- total solar and non solar fluxes --- ! 974 qns_tot(ji,jj) = pfrld(ji,jj) * qns_oce(ji,jj) + zqns_ice2d(ji,jj) + qemp_ice(ji,jj) + qemp_oce(ji,jj) 975 qsr_tot(ji,jj) = pfrld(ji,jj) * qsr_oce(ji,jj) + zqsr_ice2d(ji,jj) 976 977 ! --- heat content of precip over ice in J/m3 (to be used in 1D-thermo) --- ! 978 qprec_ice(ji,jj) = rhosn * ( ( MIN( sf(jp_tair)%fnow(ji,jj,1), rt0_snow ) - rt0 ) * cpic * tmask(ji,jj,1) - lfus ) 979 END DO 980 END DO 981 !$OMP END DO NOWAIT 818 982 819 983 ! --- heat content of evap over ice in W/m2 (to be used in 1D-thermo) --- 984 !$OMP DO schedule(static) private(jl,jj,ji) 820 985 DO jl = 1, jpl 821 qevap_ice(:,:,jl) = 0._wp ! should be -evap_ice(:,:,jl)*( ( Tice - rt0 ) * cpic * tmask(:,:,1) ) 822 ! But we do not have Tice => consider it at 0degC => evap=0 823 END DO 986 DO jj = 1, jpj 987 DO ji = 1, jpi 988 qevap_ice(ji,jj,jl) = 0._wp ! should be -evap_ice(:,:,jl)*( ( Tice - rt0 ) * cpic * tmask(:,:,1) ) 989 ! But we do not have Tice => consider it at 0degC => evap=0 990 END DO 991 END DO 992 END DO 993 !$OMP END PARALLEL 824 994 825 995 CALL wrk_dealloc( jpi,jpj, zevap, zsnw ) … … 831 1001 ! ( Maykut and Untersteiner, 1971 ; Ebert and Curry, 1993 ) 832 1002 ! 833 fr1_i0(:,:) = ( 0.18 * ( 1.0 - cldf_ice ) + 0.35 * cldf_ice ) 834 fr2_i0(:,:) = ( 0.82 * ( 1.0 - cldf_ice ) + 0.65 * cldf_ice ) 1003 !$OMP PARALLEL DO schedule(static) private(jj,ji) 1004 DO jj = 1, jpj 1005 DO ji = 1, jpi 1006 fr1_i0(ji,jj) = ( 0.18 * ( 1.0 - cldf_ice ) + 0.35 * cldf_ice ) 1007 fr2_i0(ji,jj) = ( 0.82 * ( 1.0 - cldf_ice ) + 0.65 * cldf_ice ) 1008 END DO 1009 END DO 835 1010 ! 836 1011 ! … … 844 1019 ENDIF 845 1020 846 CALL wrk_dealloc( jpi,jpj,jpl, z_qlw, z_qsb, z_dqlw, z_dqsb )1021 CALL wrk_dealloc( jpi,jpj,jpl, z_qlw, z_qsb, z_dqlw, z_dqsb, zevap_ice3d, zqns_ice3d, zqsr_ice3d ) 847 1022 CALL wrk_dealloc( jpi,jpj, zrhoa ) 848 CALL wrk_dealloc( jpi,jpj, Cd 1023 CALL wrk_dealloc( jpi,jpj, Cd, zevap_ice2d, zqns_ice2d, zqsr_ice2d) 849 1024 ! 850 1025 IF( nn_timing == 1 ) CALL timing_stop('blk_ice_flx') … … 908 1083 !!---------------------------------------------------------------------------------- 909 1084 ! 1085 !$OMP PARALLEL DO schedule(static) private(jj,ji,ztmp,ze_sat) 910 1086 DO jj = 1, jpj 911 1087 DO ji = 1, jpi … … 944 1120 !!---------------------------------------------------------------------------------- 945 1121 ! 1122 !$OMP PARALLEL DO schedule(static) private(jj,ji,zrv,ziRT) 946 1123 DO jj = 1, jpj 947 1124 DO ji = 1, jpi
Note: See TracChangeset
for help on using the changeset viewer.
