New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 7753 for trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk.F90 – NEMO

Ignore:
Timestamp:
2017-03-03T12:46:59+01:00 (6 years ago)
Author:
mocavero
Message:

Reverting trunk to remove OpenMP

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk.F90

    r7698 r7753  
    316316#if defined key_cice 
    317317      IF( MOD( kt - 1, nn_fsbc ) == 0 )   THEN 
    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 
     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)  
    332321         ENDIF  
    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 
     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) 
    344328      ENDIF 
    345329#endif 
     
    398382      ! 
    399383 
    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 
     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 
    412387      ! ----------------------------------------------------------------------------- ! 
    413388      !      0   Wind components and module at T-point relative to the moving ocean   ! 
    414389      ! ----------------------------------------------------------------------------- ! 
    415390 
     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._wp 
     394      zwnd_j(:,:) = 0._wp 
    416395#if defined key_cyclone 
    417396      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) 
    419397      DO jj = 2, jpjm1 
    420398         DO ji = fs_2, fs_jpim1   ! vect. opt. 
     
    424402      END DO 
    425403#endif 
    426 !$OMP PARALLEL DO schedule(static) private(jj, ji) 
    427404      DO jj = 2, jpjm1 
    428405         DO ji = fs_2, fs_jpim1   ! vect. opt. 
     
    434411      CALL lbc_lnk( zwnd_j(:,:) , 'T', -1. ) 
    435412      ! ... scalar wind ( = | U10m - U_oce | ) at T-point (masked) 
    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 
     413      wndm(:,:) = SQRT(  zwnd_i(:,:) * zwnd_i(:,:)   & 
     414         &             + zwnd_j(:,:) * zwnd_j(:,:)  ) * tmask(:,:,1) 
     415 
    444416      ! ----------------------------------------------------------------------------- ! 
    445417      !      I   Radiative FLUXES                                                     ! 
     
    449421      zztmp = 1. - albo 
    450422      IF( ln_dm2dc ) THEN   ;   qsr(:,:) = zztmp * sbc_dcy( 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 
    458       ENDIF 
    459  
    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 
     423      ELSE                  ;   qsr(:,:) = zztmp *          sf(jp_qsr)%fnow(:,:,1)   * tmask(:,:,1) 
     424      ENDIF 
     425 
     426      zqlw(:,:) = (  sf(jp_qlw)%fnow(:,:,1) - Stef * zst(:,:)*zst(:,:)*zst(:,:)*zst(:,:)  ) * tmask(:,:,1)   ! Long  Wave 
    466427 
    467428 
     
    500461      END IF 
    501462 
    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) 
     463      Cd_oce(:,:) = Cd(:,:)  ! record value of pure ocean-atm. drag (clem) 
     464 
    511465      DO jj = 1, jpj             ! tau module, i and j component 
    512466         DO ji = 1, jpi 
     
    517471         END DO 
    518472      END DO 
    519 !$OMP END PARALLEL 
    520473 
    521474      !                          ! add the HF tau contribution to the wind stress module 
    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 
     475      IF( lhftau )   taum(:,:) = taum(:,:) + sf(jp_tdif)%fnow(:,:,1) 
    530476 
    531477      CALL iom_put( "taum_oce", taum )   ! output wind stress module 
     
    534480      !     Note the use of 0.5*(2-umask) in order to unmask the stress along coastlines 
    535481      !     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) 
    537482      DO jj = 1, jpjm1 
    538483         DO ji = 1, fs_jpim1 
     
    551496 
    552497      ! zqla used as temporary array, for rho*U (common term of bulk formulae): 
    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 
     498      zqla(:,:) = zrhoa(:,:) * zU_zu(:,:) 
    559499 
    560500      IF( ABS( rn_zu - rn_zqt) < 0.01_wp ) THEN 
    561501         !! q_air and t_air are given at 10m (wind reference height) 
    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 
     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 
    569504      ELSE 
    570505         !! q_air and t_air are not given at 10m (wind reference height) 
    571506         ! Values of temp. and hum. adjusted to height of wind during bulk algorithm iteration must be used!!! 
    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 
     507         zevap(:,:) = rn_efac*MAX( 0._wp,             zqla(:,:)*Ce(:,:)*(zsq(:,:) - zq_zu(:,:) ) ) ! Evaporation ! using bulk wind speed 
    578508         zqsb (:,:) = cp_air(sf(jp_humi)%fnow(:,:,1))*zqla(:,:)*Ch(:,:)*(zst(:,:) - zt_zu(:,:) )   ! Sensible Heat ! using bulk wind speed 
    579509      ENDIF 
     
    597527      ! ----------------------------------------------------------------------------- ! 
    598528      ! 
    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             ! 
     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      ! 
    613540#if defined key_lim3 
    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) 
     541      qns_oce(:,:) = zqlw(:,:) - zqsb(:,:) - zqla(:,:)                                ! non solar without emp (only needed by LIM3) 
     542      qsr_oce(:,:) = qsr(:,:) 
    616543#endif 
    617          END DO 
    618       END DO 
    619544      ! 
    620545      IF ( nn_ice == 0 ) THEN 
     
    626551         CALL iom_put( "qsr_oce" ,   qsr  )                 ! output downward solar heat over the ocean 
    627552         CALL iom_put( "qt_oce"  ,   qns+qsr )              ! output total downward heat over the ocean 
    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 
     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] 
    635555         CALL iom_put( 'snowpre', sprecip * 86400. )        ! Snow 
    636556         CALL iom_put( 'precip' , tprecip * 86400. )        ! Total precipitation 
     
    679599      CALL wrk_alloc( jpi,jpj, Cd ) 
    680600 
    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 
     601      Cd(:,:) = Cd_ice 
    687602 
    688603      ! Make ice-atm. drag dependent on ice concentration (see Lupkes et al. 2012) (clem) 
     
    698613      zrhoa (:,:) =  rho_air(sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1)) 
    699614 
    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 
     615      !!gm brutal.... 
     616      utau_ice  (:,:) = 0._wp 
     617      vtau_ice  (:,:) = 0._wp 
     618      wndm_ice  (:,:) = 0._wp 
     619      !!gm end 
    710620 
    711621      ! ----------------------------------------------------------------------------- ! 
     
    715625      CASE( 'I' )                  ! B-grid ice dynamics :   I-point (i.e. F-point with sea-ice indexation) 
    716626         !                           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) 
    718627         DO jj = 2, jpjm1 
    719628            DO ji = 2, jpim1   ! B grid : NO vector opt 
     
    740649         ! 
    741650      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) 
    743651         DO jj = 2, jpj 
    744652            DO ji = fs_2, jpi   ! vect. opt. 
     
    748656            END DO 
    749657         END DO 
    750 !$OMP PARALLEL DO schedule(static) private(jj,ji) 
    751658         DO jj = 2, jpjm1 
    752659            DO ji = fs_2, fs_jpim1   ! vect. opt. 
     
    793700      REAL(wp) ::   zztmp, z1_lsub           !   -      - 
    794701      REAL(wp), DIMENSION(:,:,:), POINTER ::   z_qlw         ! long wave heat flux over ice 
    795       REAL(wp), DIMENSION(:,:,:), POINTER ::   zevap_ice3d, zqns_ice3d, zqsr_ice3d  
    796702      REAL(wp), DIMENSION(:,:,:), POINTER ::   z_qsb         ! sensible  heat flux over ice 
    797703      REAL(wp), DIMENSION(:,:,:), POINTER ::   z_dqlw        ! long wave heat sensitivity over ice 
    798704      REAL(wp), DIMENSION(:,:,:), POINTER ::   z_dqsb        ! sensible  heat sensitivity over ice 
    799705      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 
    801706      REAL(wp), DIMENSION(:,:)  , POINTER ::   zrhoa 
    802707      REAL(wp), DIMENSION(:,:)  , POINTER ::   Cd            ! transfer coefficient for momentum      (tau) 
     
    805710      IF( nn_timing == 1 )  CALL timing_start('blk_ice_flx') 
    806711      ! 
    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) 
     712      CALL wrk_alloc( jpi,jpj,jpl,   z_qlw, z_qsb, z_dqlw, z_dqsb ) 
     713      CALL wrk_alloc( jpi,jpj,       zrhoa) 
    809714      CALL wrk_alloc( jpi,jpj, Cd ) 
    810715 
    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 
     716      Cd(:,:) = Cd_ice 
    817717 
    818718      ! Make ice-atm. drag dependent on ice concentration (see Lupkes et al.  2012) (clem) 
     
    831731      ! 
    832732      zztmp = 1. / ( 1. - albo ) 
    833 !$OMP PARALLEL 
    834 !$OMP DO schedule(static) private(jl,jj,ji,zst2,zst3)            ! ========================== ! 
    835       DO jl = 1, jpl                                             !  Loop over ice categories  ! 
    836          !                                                       ! ========================== ! 
     733      !                                     ! ========================== ! 
     734      DO jl = 1, jpl                        !  Loop over ice categories  ! 
     735         !                                  ! ========================== ! 
    837736         DO jj = 1 , jpj 
    838737            DO ji = 1, jpi 
     
    882781      END DO 
    883782      ! 
    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 
     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] 
    892785      CALL iom_put( 'snowpre', sprecip * 86400. )                  ! Snow precipitation 
    893786      CALL iom_put( 'precip' , tprecip * 86400. )                  ! Total precipitation 
     
    898791      ! --- evaporation --- ! 
    899792      z1_lsub = 1._wp / Lsub 
    900 !$OMP PARALLEL 
    901 !$OMP DO schedule(static) private(jl,jj,ji) 
     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 
     799      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 ) 
     818 
     819      ! --- heat content of evap over ice in W/m2 (to be used in 1D-thermo) --- 
    902820      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 
     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  
    909823      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 
    921       CALL lim_thd_snwblow( pfrld, zsnw )  ! snow distribution over ice after wind blowing 
    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 
    982  
    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) 
    985       DO jl = 1, jpl 
    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 
    994824 
    995825      CALL wrk_dealloc( jpi,jpj,   zevap, zsnw ) 
     
    1001831      ! ( Maykut and Untersteiner, 1971 ; Ebert and Curry, 1993 ) 
    1002832      ! 
    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 
     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 ) 
    1010835      ! 
    1011836      ! 
     
    1019844      ENDIF 
    1020845 
    1021       CALL wrk_dealloc( jpi,jpj,jpl,   z_qlw, z_qsb, z_dqlw, z_dqsb, zevap_ice3d, zqns_ice3d, zqsr_ice3d ) 
     846      CALL wrk_dealloc( jpi,jpj,jpl,   z_qlw, z_qsb, z_dqlw, z_dqsb ) 
    1022847      CALL wrk_dealloc( jpi,jpj,       zrhoa ) 
    1023       CALL wrk_dealloc( jpi,jpj, Cd, zevap_ice2d, zqns_ice2d, zqsr_ice2d) 
     848      CALL wrk_dealloc( jpi,jpj, Cd ) 
    1024849      ! 
    1025850      IF( nn_timing == 1 )  CALL timing_stop('blk_ice_flx') 
     
    1083908      !!---------------------------------------------------------------------------------- 
    1084909      ! 
    1085 !$OMP PARALLEL DO schedule(static) private(jj,ji,ztmp,ze_sat) 
    1086910      DO jj = 1, jpj 
    1087911         DO ji = 1, jpi 
     
    1120944      !!---------------------------------------------------------------------------------- 
    1121945      ! 
    1122 !$OMP PARALLEL DO schedule(static) private(jj,ji,zrv,ziRT) 
    1123946      DO jj = 1, jpj 
    1124947         DO ji = 1, jpi 
Note: See TracChangeset for help on using the changeset viewer.