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 12706 for NEMO/branches/UKMO/NEMO_4.0.2_ENHANCE-02_ISF_nemo/src/OCE/SBC – NEMO

Ignore:
Timestamp:
2020-04-07T18:34:56+02:00 (4 years ago)
Author:
mathiot
Message:

NEMO_4.0.2_ENHANCE-02_ISF_nemo: in sync with trunk right before release_4.0-HEAD was created (svn merge -r 12072:12367 /NEMO/trunk)

Location:
NEMO/branches/UKMO/NEMO_4.0.2_ENHANCE-02_ISF_nemo/src/OCE/SBC
Files:
8 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/UKMO/NEMO_4.0.2_ENHANCE-02_ISF_nemo/src/OCE/SBC/cpl_oasis3.F90

    r10582 r12706  
    306306      ! End of definition phase 
    307307      !------------------------------------------------------------------ 
    308        
     308      !      
     309#if defined key_agrif 
     310      IF( agrif_fixed() == Agrif_Nb_Fine_Grids() ) THEN 
     311#endif 
    309312      CALL oasis_enddef(nerror) 
    310313      IF( nerror /= OASIS_Ok )   CALL oasis_abort ( ncomp_id, 'cpl_define', 'Failure in oasis_enddef') 
     314#if defined key_agrif 
     315      ENDIF 
     316#endif 
    311317      ! 
    312318      IF ( ltmp_wapatch ) THEN 
     
    357363                     WRITE(numout,*) 'oasis_put:  kstep ', kstep 
    358364                     WRITE(numout,*) 'oasis_put:   info ', kinfo 
    359                      WRITE(numout,*) '     - Minimum value is ', MINVAL(pdata(:,:,jc)) 
    360                      WRITE(numout,*) '     - Maximum value is ', MAXVAL(pdata(:,:,jc)) 
    361                      WRITE(numout,*) '     -     Sum value is ', SUM(pdata(:,:,jc)) 
     365                     WRITE(numout,*) '     - Minimum value is ', MINVAL(pdata(nldi:nlei,nldj:nlej,jc)) 
     366                     WRITE(numout,*) '     - Maximum value is ', MAXVAL(pdata(nldi:nlei,nldj:nlej,jc)) 
     367                     WRITE(numout,*) '     -     Sum value is ', SUM(pdata(nldi:nlei,nldj:nlej,jc)) 
    362368                     WRITE(numout,*) '****************' 
    363369                  ENDIF 
  • NEMO/branches/UKMO/NEMO_4.0.2_ENHANCE-02_ISF_nemo/src/OCE/SBC/fldread.F90

    r12143 r12706  
    833833 
    834834      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   pdta_read       ! data read in bdy file 
    835       REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pdta_read_z     ! depth of the data read in bdy file 
    836       REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pdta_read_dz    ! thickness of the levels in bdy file 
     835      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   pdta_read_z     ! depth of the data read in bdy file 
     836      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   pdta_read_dz    ! thickness of the levels in bdy file 
    837837      REAL(wp), DIMENSION(:,:,:), INTENT(  out) ::   pdta            ! output field on model grid (2 dimensional) 
    838838      REAL(wp)                  , INTENT(in   ) ::   pfv             ! fillvalue of the data read in bdy file 
     
    841841      INTEGER                   , INTENT(in   ) ::   kbdy            ! bdy number 
    842842      !! 
    843       INTEGER                                   ::   ipi             ! length of boundary data on local process 
    844       INTEGER                                   ::   ipkb            ! number of vertical levels in boundary data file 
    845       INTEGER                                   ::   jb, ji, jj, jk, jkb   ! loop counters 
    846       REAL(wp)                                  ::   zcoef 
    847       REAL(wp)                                  ::   zl, zi, zh      ! tmp variable for current depth and interpolation factor 
    848       REAL(wp)                                  ::   zfv_alt, ztrans, ztrans_new ! fillvalue and alternative -ABS(pfv) 
    849       REAL(wp), DIMENSION(jpk)                  ::   zdepth, zdhalf  ! level and half-level depth 
     843      INTEGER                  ::   ipi                 ! length of boundary data on local process 
     844      INTEGER                  ::   ipkb                ! number of vertical levels in boundary data file 
     845      INTEGER                  ::   ipkmax              ! number of vertical levels in boundary data file where no mask 
     846      INTEGER                  ::   jb, ji, jj, jk, jkb ! loop counters 
     847      REAL(wp)                 ::   zcoef, zi           !  
     848      REAL(wp)                 ::   ztrans, ztrans_new  ! transports 
     849      REAL(wp), DIMENSION(jpk) ::   zdepth, zdhalf      ! level and half-level depth 
    850850      !!--------------------------------------------------------------------- 
    851851      
     
    853853      ipkb = SIZE( pdta_read, 3 ) 
    854854       
    855       zfv_alt = -ABS(pfv)  ! set _FillValue < 0 as we make use of MAXVAL and MAXLOC later 
    856       ! 
    857       WHERE( pdta_read == pfv ) 
    858          pdta_read_z  = zfv_alt ! safety: put fillvalue into external depth field so consistent with data 
    859          pdta_read_dz = 0._wp   ! safety: put 0._wp into external thickness factors to ensure transport is correct 
    860       ENDWHERE 
    861        
    862855      DO jb = 1, ipi 
    863856         ji = idx_bdy(kbdy)%nbi(jb,kgrd) 
    864857         jj = idx_bdy(kbdy)%nbj(jb,kgrd) 
    865          zh  = SUM(pdta_read_dz(jb,1,:) ) 
    866          ! 
    867          ! Warnings to flag differences in the input and model topgraphy - is this useful/necessary? 
     858         ! 
     859         ! --- max jk where input data /= FillValue --- ! 
     860         ipkmax = 1 
     861         DO jkb = 2, ipkb 
     862            IF( pdta_read(jb,1,jkb) /= pfv )   ipkmax = MAX( ipkmax, jkb ) 
     863         END DO 
     864         ! 
     865         ! --- calculate depth at t,u,v points --- ! 
    868866         SELECT CASE( kgrd )                          
    869          CASE(1) 
    870             IF( ABS( (zh - ht_n(ji,jj)) / ht_n(ji,jj)) * tmask(ji,jj,1) > 0.01_wp ) THEN 
    871                WRITE(ctmp1,"(I10.10)") jb  
    872                CALL ctl_warn('fld_bdy_interp: T depths differ between grids at BDY point '//TRIM(ctmp1)//' by more than 1%') 
    873                !   IF(lwp) WRITE(numout,*) 'DEPTHT', zh, sum(e3t_n(ji,jj,:), mask=tmask(ji,jj,:)==1),  ht_n(ji,jj), jb, jb, ji, jj 
    874             ENDIF 
    875          CASE(2) 
    876             IF( ABS( (zh - hu_n(ji,jj)) * r1_hu_n(ji,jj)) * umask(ji,jj,1) > 0.01_wp ) THEN 
    877                WRITE(ctmp1,"(I10.10)") jb  
    878                CALL ctl_warn('fld_bdy_interp: U depths differ between grids at BDY point '//TRIM(ctmp1)//' by more than 1%') 
    879                !   IF(lwp) WRITE(numout,*) 'DEPTHU', zh, SUM(e3u_n(ji,jj,:), mask=umask(ji,jj,:)==1),  SUM(umask(ji,jj,:)), & 
    880                !      &                hu_n(ji,jj), jb, jb, ji, jj, narea-1, pdta_read(jb,1,:) 
    881             ENDIF 
    882          CASE(3) 
    883             IF( ABS( (zh - hv_n(ji,jj)) * r1_hv_n(ji,jj)) * vmask(ji,jj,1) > 0.01_wp ) THEN 
    884                WRITE(ctmp1,"(I10.10)") jb 
    885                CALL ctl_warn('fld_bdy_interp: V depths differ between grids at BDY point '//TRIM(ctmp1)//' by more than 1%') 
    886             ENDIF 
    887          END SELECT 
    888          ! 
    889          SELECT CASE( kgrd )                          
    890          CASE(1) 
    891             ! depth of T points: 
     867         CASE(1)            ! depth of T points: 
    892868            zdepth(:) = gdept_n(ji,jj,:) 
    893          CASE(2) 
    894             ! depth of U points: we must not use gdept_n as we don't want to do a communication 
    895             !   --> copy what is done for gdept_n in domvvl... 
     869         CASE(2)            ! depth of U points: we must not use gdept_n as we don't want to do a communication 
     870            !                 --> copy what is done for gdept_n in domvvl... 
    896871            zdhalf(1) = 0.0_wp 
    897872            zdepth(1) = 0.5_wp * e3uw_n(ji,jj,1) 
     
    903878               zcoef = ( umask(ji,jj,jk) - wumask(ji,jj,jk) ) 
    904879               zdhalf(jk) = zdhalf(jk-1) + e3u_n(ji,jj,jk-1) 
    905                zdepth(jk) =      zcoef  * ( zdhalf(jk  ) + 0.5 * e3uw_n(ji,jj,jk))  & 
    906                   &         + (1-zcoef) * ( zdepth(jk-1) +       e3uw_n(ji,jj,jk)) 
     880               zdepth(jk) =       zcoef  * ( zdhalf(jk  ) + 0.5 * e3uw_n(ji,jj,jk))  & 
     881                  &         + (1.-zcoef) * ( zdepth(jk-1) +       e3uw_n(ji,jj,jk)) 
    907882            END DO 
    908          CASE(3) 
    909             ! depth of V points: we must not use gdept_n as we don't want to do a communication 
    910             !   --> copy what is done for gdept_n in domvvl... 
     883         CASE(3)            ! depth of V points: we must not use gdept_n as we don't want to do a communication 
     884            !                 --> copy what is done for gdept_n in domvvl... 
    911885            zdhalf(1) = 0.0_wp 
    912886            zdepth(1) = 0.5_wp * e3vw_n(ji,jj,1) 
     
    918892               zcoef = ( vmask(ji,jj,jk) - wvmask(ji,jj,jk) ) 
    919893               zdhalf(jk) = zdhalf(jk-1) + e3v_n(ji,jj,jk-1) 
    920                zdepth(jk) =      zcoef  * ( zdhalf(jk  ) + 0.5 * e3vw_n(ji,jj,jk))  & 
    921                   &         + (1-zcoef) * ( zdepth(jk-1) +       e3vw_n(ji,jj,jk)) 
     894               zdepth(jk) =       zcoef  * ( zdhalf(jk  ) + 0.5 * e3vw_n(ji,jj,jk))  & 
     895                  &         + (1.-zcoef) * ( zdepth(jk-1) +       e3vw_n(ji,jj,jk)) 
    922896            END DO 
    923897         END SELECT 
    924          ! 
    925          DO jk = 1, jpk                       
    926             IF(     zdepth(jk) < pdta_read_z(jb,1,          1) ) THEN                ! above the first level of external data 
    927                pdta(jb,1,jk) =  pdta_read(jb,1,1) 
    928             ELSEIF( zdepth(jk) > pdta_read_z(jb,1,ipkb) ) THEN                       ! below the last level of external data  
    929                pdta(jb,1,jk) =  pdta_read(jb,1,MAXLOC(pdta_read_z(jb,1,:),1)) 
    930             ELSE                                                             ! inbetween: vertical interpolation between jkb & jkb+1 
    931                DO jkb = 1, ipkb-1                                            ! when  gdept_n(jkb) < zdepth(jk) < gdept_n(jkb+1) 
    932                   IF( ( ( zdepth(jk) - pdta_read_z(jb,1,jkb) ) * ( zdepth(jk) - pdta_read_z(jb,1,jkb+1) ) <= 0._wp ) & 
    933                      &    .AND. ( pdta_read_z(jb,1,jkb+1) /= zfv_alt) ) THEN  ! linear interpolation between 2 levels 
     898         !          
     899         ! --- interpolate bdy data on the model grid --- ! 
     900         DO jk = 1, jpk 
     901            IF(     zdepth(jk) <= pdta_read_z(jb,1,1)      ) THEN   ! above the first level of external data 
     902               pdta(jb,1,jk) = pdta_read(jb,1,1) 
     903            ELSEIF( zdepth(jk) >  pdta_read_z(jb,1,ipkmax) ) THEN   ! below the last level of external data /= FillValue 
     904               pdta(jb,1,jk) = pdta_read(jb,1,ipkmax) 
     905            ELSE                                                    ! inbetween: vertical interpolation between jkb & jkb+1 
     906               DO jkb = 1, ipkmax-1 
     907                  IF( ( ( zdepth(jk) - pdta_read_z(jb,1,jkb) ) * ( zdepth(jk) - pdta_read_z(jb,1,jkb+1) ) ) <= 0._wp ) THEN ! linear interpolation between 2 levels 
    934908                     zi = ( zdepth(jk) - pdta_read_z(jb,1,jkb) ) / ( pdta_read_z(jb,1,jkb+1) - pdta_read_z(jb,1,jkb) ) 
    935                      pdta(jb,1,jk) = pdta_read(jb,1,jkb) + ( pdta_read  (jb,1,jkb+1) - pdta_read  (jb,1,jkb) ) * zi 
     909                     pdta(jb,1,jk) = pdta_read(jb,1,jkb) + zi * ( pdta_read(jb,1,jkb+1) - pdta_read(jb,1,jkb) ) 
    936910                  ENDIF 
    937911               END DO 
    938912            ENDIF 
    939          END DO   ! jpk 
     913         END DO 
    940914         ! 
    941915      END DO   ! ipi 
    942        
    943       IF(kgrd == 2) THEN                                  ! do we need to adjust the transport term? 
     916 
     917      ! --- mask data and adjust transport --- ! 
     918      SELECT CASE( kgrd )                          
     919 
     920      CASE(1)                                 ! mask data (probably unecessary) 
    944921         DO jb = 1, ipi 
    945922            ji = idx_bdy(kbdy)%nbi(jb,kgrd) 
    946923            jj = idx_bdy(kbdy)%nbj(jb,kgrd) 
    947             zh  = SUM(pdta_read_dz(jb,1,:) ) 
     924            DO jk = 1, jpk                       
     925               pdta(jb,1,jk) = pdta(jb,1,jk) * tmask(ji,jj,jk) 
     926            END DO 
     927         END DO 
     928          
     929      CASE(2)                                 ! adjust the U-transport term 
     930         DO jb = 1, ipi 
     931            ji = idx_bdy(kbdy)%nbi(jb,kgrd) 
     932            jj = idx_bdy(kbdy)%nbj(jb,kgrd) 
    948933            ztrans = 0._wp 
     934            DO jkb = 1, ipkb                              ! calculate transport on input grid 
     935               IF( pdta_read(jb,1,jkb) /= pfv )   ztrans = ztrans + pdta_read(jb,1,jkb) * pdta_read_dz(jb,1,jkb) 
     936            ENDDO 
    949937            ztrans_new = 0._wp 
    950             DO jkb = 1, ipkb                              ! calculate transport on input grid 
    951                ztrans     = ztrans     + pdta_read(jb,1,jkb) * pdta_read_dz(jb, 1,jkb) 
    952             ENDDO 
    953938            DO jk = 1, jpk                                ! calculate transport on model grid 
    954                ztrans_new = ztrans_new +      pdta(jb,1,jk ) *        e3u_n(ji,jj,jk ) * umask(ji,jj,jk) 
     939               ztrans_new = ztrans_new + pdta(jb,1,jk ) * e3u_n(ji,jj,jk) * umask(ji,jj,jk) 
    955940            ENDDO 
    956941            DO jk = 1, jpk                                ! make transport correction 
    957942               IF(ldtotvel) THEN ! bdy data are total velocity so adjust bt transport term to match input data 
    958943                  pdta(jb,1,jk) = ( pdta(jb,1,jk) + ( ztrans - ztrans_new ) * r1_hu_n(ji,jj) ) * umask(ji,jj,jk) 
    959                ELSE ! we're just dealing with bc velocity so bt transport term should sum to zero 
    960                   pdta(jb,1,jk) =   pdta(jb,1,jk) + (  0._wp - ztrans_new ) * r1_hu_n(ji,jj)  * umask(ji,jj,jk) 
     944               ELSE              ! we're just dealing with bc velocity so bt transport term should sum to zero 
     945                  pdta(jb,1,jk) = ( pdta(jb,1,jk) + (  0._wp - ztrans_new ) * r1_hu_n(ji,jj) ) * umask(ji,jj,jk) 
    961946               ENDIF 
    962947            ENDDO 
    963948         ENDDO 
    964       ENDIF 
    965        
    966       IF(kgrd == 3) THEN                                  ! do we need to adjust the transport term? 
     949 
     950      CASE(3)                                 ! adjust the V-transport term 
    967951         DO jb = 1, ipi 
    968952            ji = idx_bdy(kbdy)%nbi(jb,kgrd) 
    969953            jj = idx_bdy(kbdy)%nbj(jb,kgrd) 
    970             zh  = SUM(pdta_read_dz(jb,1,:) ) 
    971954            ztrans = 0._wp 
     955            DO jkb = 1, ipkb                              ! calculate transport on input grid 
     956               IF( pdta_read(jb,1,jkb) /= pfv )   ztrans = ztrans + pdta_read(jb,1,jkb) * pdta_read_dz(jb,1,jkb) 
     957            ENDDO 
    972958            ztrans_new = 0._wp 
    973             DO jkb = 1, ipkb                              ! calculate transport on input grid 
    974                ztrans     = ztrans     + pdta_read(jb,1,jkb) * pdta_read_dz(jb, 1,jkb) 
    975             ENDDO 
    976959            DO jk = 1, jpk                                ! calculate transport on model grid 
    977                ztrans_new = ztrans_new +      pdta(jb,1,jk ) *        e3v_n(ji,jj,jk ) * vmask(ji,jj,jk) 
     960               ztrans_new = ztrans_new + pdta(jb,1,jk ) * e3v_n(ji,jj,jk) * vmask(ji,jj,jk) 
    978961            ENDDO 
    979962            DO jk = 1, jpk                                ! make transport correction 
    980963               IF(ldtotvel) THEN ! bdy data are total velocity so adjust bt transport term to match input data 
    981964                  pdta(jb,1,jk) = ( pdta(jb,1,jk) + ( ztrans - ztrans_new ) * r1_hv_n(ji,jj) ) * vmask(ji,jj,jk) 
    982                ELSE ! we're just dealing with bc velocity so bt transport term should sum to zero 
    983                   pdta(jb,1,jk) =   pdta(jb,1,jk) + (  0._wp - ztrans_new ) * r1_hv_n(ji,jj)  * vmask(ji,jj,jk) 
     965               ELSE              ! we're just dealing with bc velocity so bt transport term should sum to zero 
     966                  pdta(jb,1,jk) = ( pdta(jb,1,jk) + (  0._wp - ztrans_new ) * r1_hv_n(ji,jj) ) * vmask(ji,jj,jk) 
    984967               ENDIF 
    985968            ENDDO 
    986969         ENDDO 
    987       ENDIF 
    988        
     970      END SELECT 
     971 
    989972   END SUBROUTINE fld_bdy_interp 
    990973 
    991  
     974    
    992975   SUBROUTINE fld_rot( kt, sd ) 
    993976      !!--------------------------------------------------------------------- 
  • NEMO/branches/UKMO/NEMO_4.0.2_ENHANCE-02_ISF_nemo/src/OCE/SBC/sbc_oce.F90

    r11521 r12706  
    104104   !!              Ocean Surface Boundary Condition fields 
    105105   !!---------------------------------------------------------------------- 
    106    INTEGER , PUBLIC ::  ncpl_qsr_freq            !: qsr coupling frequency per days from atmosphere 
     106   INTEGER , PUBLIC ::  ncpl_qsr_freq = 0        !: qsr coupling frequency per days from atmosphere (used by top) 
    107107   ! 
    108108   LOGICAL , PUBLIC ::   lhftau = .FALSE.        !: HF tau used in TKE: mean(stress module) - module(mean stress) 
  • NEMO/branches/UKMO/NEMO_4.0.2_ENHANCE-02_ISF_nemo/src/OCE/SBC/sbcblk.F90

    r12143 r12706  
    801801      REAL(wp), DIMENSION(jpi,jpj)     ::   zevap, zsnw   ! evaporation and snw distribution after wind blowing (SI3) 
    802802      REAL(wp), DIMENSION(jpi,jpj)     ::   zrhoa 
     803      REAL(wp), DIMENSION(jpi,jpj)     ::   ztmp, ztmp2 
    803804      !!--------------------------------------------------------------------- 
    804805      ! 
     
    913914         qtr_ice_top(:,:,:) = 0._wp  
    914915      END WHERE 
     916      ! 
     917 
     918      IF( iom_use('evap_ao_cea') .OR. iom_use('hflx_evap_cea') ) THEN 
     919         ztmp(:,:) = zevap(:,:) * ( 1._wp - at_i_b(:,:) )  
     920         CALL iom_put( 'evap_ao_cea'  , ztmp(:,:) * tmask(:,:,1) )   ! ice-free oce evap (cell average) 
     921         CALL iom_put( 'hflx_evap_cea', ztmp(:,:) * sst_m(:,:) * rcp * tmask(:,:,1) )   ! heat flux from evap (cell average) 
     922      ENDIF 
     923      IF( iom_use('hflx_rain_cea') ) THEN 
     924         ztmp(:,:) = rcp * ( SUM( (ptsu-rt0) * a_i_b, dim=3 ) + sst_m(:,:) * ( 1._wp - at_i_b(:,:) ) ) 
     925         CALL iom_put( 'hflx_rain_cea', ( tprecip(:,:) - sprecip(:,:) ) * ztmp(:,:) )   ! heat flux from rain (cell average) 
     926      ENDIF 
     927      IF( iom_use('hflx_snow_cea') .OR. iom_use('hflx_snow_ao_cea') .OR. iom_use('hflx_snow_ai_cea')  )  THEN 
     928          WHERE( SUM( a_i_b, dim=3 ) > 1.e-10 ) ;   ztmp(:,:) = rcpi * SUM( (ptsu-rt0) * a_i_b, dim=3 ) / SUM( a_i_b, dim=3 ) 
     929          ELSEWHERE                             ;   ztmp(:,:) = rcp * sst_m(:,:)     
     930          ENDWHERE 
     931          ztmp2(:,:) = sprecip(:,:) * ( ztmp(:,:) - rLfus )  
     932          CALL iom_put('hflx_snow_cea'   , ztmp2(:,:) ) ! heat flux from snow (cell average) 
     933          CALL iom_put('hflx_snow_ao_cea', ztmp2(:,:) * ( 1._wp - zsnw(:,:) ) ) ! heat flux from snow (over ocean) 
     934          CALL iom_put('hflx_snow_ai_cea', ztmp2(:,:) *           zsnw(:,:)   ) ! heat flux from snow (over ice) 
     935      ENDIF 
    915936      ! 
    916937      IF(ln_ctl) THEN 
  • NEMO/branches/UKMO/NEMO_4.0.2_ENHANCE-02_ISF_nemo/src/OCE/SBC/sbccpl.F90

    r12143 r12706  
    574574      IF ( TRIM( sn_rcv_emp%clcat    ) == 'yes' )   srcv(jpr_ievp)%nct       = nn_cats_cpl 
    575575 
     576#if defined key_si3 
     577      IF( ln_cndflx .AND. .NOT.ln_cndemulate ) THEN  
     578         IF( .NOT.srcv(jpr_ts_ice)%laction )  & 
     579            &   CALL ctl_stop( 'sbc_cpl_init: srcv(jpr_ts_ice)%laction should be set to true when ln_cndflx=T' )      
     580      ENDIF 
     581#endif 
    576582      !                                                      ! ------------------------- ! 
    577583      !                                                      !      Wave breaking        !     
     
    863869      ELSE IF( sn_snd_crt%clvgrd /= 'T' ) THEN   
    864870         CALL ctl_stop( 'sn_snd_crt%clvgrd must be equal to T' ) 
    865          ssnd(jps_ocx1:jps_ivz1)%clgrid  = 'T'      ! all oce and ice components on the same unique grid 
    866871      ENDIF 
    867872      ssnd(jps_ocx1:jps_ivz1)%laction = .TRUE.   ! default: all are send 
     
    10411046      ENDIF 
    10421047      xcplmask(:,:,0) = 1. - SUM( xcplmask(:,:,1:nn_cplmodel), dim = 3 ) 
    1043       ! 
    1044       ncpl_qsr_freq = cpl_freq( 'O_QsrOce' ) + cpl_freq( 'O_QsrMix' ) + cpl_freq( 'I_QsrOce' ) + cpl_freq( 'I_QsrMix' ) 
    1045       IF( ln_dm2dc .AND. ln_cpl .AND. ncpl_qsr_freq /= 86400 )   & 
    1046          &   CALL ctl_stop( 'sbc_cpl_init: diurnal cycle reconstruction (ln_dm2dc) needs daily couping for solar radiation' ) 
    1047       IF( ln_dm2dc .AND. ln_cpl ) ncpl_qsr_freq = 86400 / ncpl_qsr_freq 
    10481048      ! 
    10491049   END SUBROUTINE sbc_cpl_init 
     
    11111111      REAL(wp), DIMENSION(jpi,jpj) ::   ztx, zty, zmsk, zemp, zqns, zqsr 
    11121112      !!---------------------------------------------------------------------- 
     1113      ! 
     1114      IF( kt == nit000 ) THEN 
     1115      !   cannot be done in the init phase when we use agrif as cpl_freq requires that oasis_enddef is done 
     1116         ncpl_qsr_freq = cpl_freq( 'O_QsrOce' ) + cpl_freq( 'O_QsrMix' ) + cpl_freq( 'I_QsrOce' ) + cpl_freq( 'I_QsrMix' ) 
     1117         IF( ln_dm2dc .AND. ncpl_qsr_freq /= 86400 )   & 
     1118            &   CALL ctl_stop( 'sbc_cpl_rcv: diurnal cycle reconstruction (ln_dm2dc) needs daily couping for solar radiation' ) 
     1119         ncpl_qsr_freq = 86400 / ncpl_qsr_freq   ! used by top 
     1120      ENDIF 
    11131121      ! 
    11141122      IF( ln_mixcpl )   zmsk(:,:) = 1. - xcplmask(:,:,0) 
     
    12441252      IF( srcv(jpr_co2)%laction )   atm_co2(:,:) = frcv(jpr_co2)%z3(:,:,1) 
    12451253      ! 
    1246       !                                                      ! ================== ! 
    1247       !                                                      !   ice skin temp.   ! 
    1248       !                                                      ! ================== ! 
    1249 #if defined key_si3 
    1250       ! needed by Met Office 
    1251       IF( srcv(jpr_ts_ice)%laction ) THEN  
    1252          WHERE    ( frcv(jpr_ts_ice)%z3(:,:,:) > 0.0  )   ;   tsfc_ice(:,:,:) = 0.0  
    1253          ELSEWHERE( frcv(jpr_ts_ice)%z3(:,:,:) < -60. )   ;   tsfc_ice(:,:,:) = -60. 
    1254          ELSEWHERE                                        ;   tsfc_ice(:,:,:) = frcv(jpr_ts_ice)%z3(:,:,:) 
    1255          END WHERE 
    1256       ENDIF  
    1257 #endif 
    12581254      !                                                      ! ========================= !  
    12591255      !                                                      ! Mean Sea Level Pressure   !   (taum)  
     
    16351631      !!                   sprecip           solid precipitation over the ocean   
    16361632      !!---------------------------------------------------------------------- 
    1637       REAL(wp), INTENT(in), DIMENSION(:,:)             ::   picefr     ! ice fraction                [0 to 1] 
    1638       !                                                !!           ! optional arguments, used only in 'mixed oce-ice' case 
    1639       REAL(wp), INTENT(in), DIMENSION(:,:,:), OPTIONAL ::   palbi      ! all skies ice albedo  
    1640       REAL(wp), INTENT(in), DIMENSION(:,:  ), OPTIONAL ::   psst       ! sea surface temperature     [Celsius] 
    1641       REAL(wp), INTENT(in), DIMENSION(:,:,:), OPTIONAL ::   pist       ! ice surface temperature     [Kelvin] 
    1642       REAL(wp), INTENT(in), DIMENSION(:,:,:), OPTIONAL ::   phs        ! snow depth                  [m] 
    1643       REAL(wp), INTENT(in), DIMENSION(:,:,:), OPTIONAL ::   phi        ! ice thickness               [m] 
     1633      REAL(wp), INTENT(in)   , DIMENSION(:,:)             ::   picefr     ! ice fraction                [0 to 1] 
     1634      !                                                   !!           ! optional arguments, used only in 'mixed oce-ice' case or for Met-Office coupling 
     1635      REAL(wp), INTENT(in)   , DIMENSION(:,:,:), OPTIONAL ::   palbi      ! all skies ice albedo  
     1636      REAL(wp), INTENT(in)   , DIMENSION(:,:  ), OPTIONAL ::   psst       ! sea surface temperature     [Celsius] 
     1637      REAL(wp), INTENT(inout), DIMENSION(:,:,:), OPTIONAL ::   pist       ! ice surface temperature     [Kelvin] => inout for Met-Office 
     1638      REAL(wp), INTENT(in)   , DIMENSION(:,:,:), OPTIONAL ::   phs        ! snow depth                  [m] 
     1639      REAL(wp), INTENT(in)   , DIMENSION(:,:,:), OPTIONAL ::   phi        ! ice thickness               [m] 
    16441640      ! 
    16451641      INTEGER  ::   ji, jj, jl   ! dummy loop index 
     
    16481644      REAL(wp), DIMENSION(jpi,jpj)     ::   zemp_tot, zemp_ice, zemp_oce, ztprecip, zsprecip  , zevap_oce, zdevap_ice 
    16491645      REAL(wp), DIMENSION(jpi,jpj)     ::   zqns_tot, zqns_oce, zqsr_tot, zqsr_oce, zqprec_ice, zqemp_oce, zqemp_ice 
    1650       REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zqns_ice, zqsr_ice, zdqns_ice, zqevap_ice, zevap_ice    !!gm , zfrqsr_tr_i 
     1646      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zqns_ice, zqsr_ice, zdqns_ice, zqevap_ice, zevap_ice, zqtr_ice_top, ztsu 
    16511647      !!---------------------------------------------------------------------- 
    16521648      ! 
     
    17741770      IF( srcv(jpr_cal)%laction )   CALL iom_put( 'calving_cea' , frcv(jpr_cal)%z3(:,:,1) * tmask(:,:,1)                )  ! calving 
    17751771      IF( srcv(jpr_icb)%laction )   CALL iom_put( 'iceberg_cea' , frcv(jpr_icb)%z3(:,:,1) * tmask(:,:,1)                )  ! icebergs 
    1776       IF( iom_use('snowpre') )      CALL iom_put( 'snowpre'     , sprecip(:,:)                                          )  ! Snow 
    1777       IF( iom_use('precip') )       CALL iom_put( 'precip'      , tprecip(:,:)                                          )  ! total  precipitation 
    1778       IF( iom_use('rain') )         CALL iom_put( 'rain'        , tprecip(:,:) - sprecip(:,:)                           )  ! liquid precipitation  
    1779       IF( iom_use('snow_ao_cea') )  CALL iom_put( 'snow_ao_cea' , sprecip(:,:) * ( 1._wp - zsnw(:,:) )                  )  ! Snow over ice-free ocean  (cell average) 
    1780       IF( iom_use('snow_ai_cea') )  CALL iom_put( 'snow_ai_cea' , sprecip(:,:) *           zsnw(:,:)                    )  ! Snow over sea-ice         (cell average) 
    1781       IF( iom_use('subl_ai_cea') )  CALL iom_put( 'subl_ai_cea' , frcv(jpr_ievp)%z3(:,:,1) * picefr(:,:) * tmask(:,:,1) )  ! Sublimation over sea-ice (cell average) 
    1782       IF( iom_use('evap_ao_cea') )  CALL iom_put( 'evap_ao_cea' , ( frcv(jpr_tevp)%z3(:,:,1)  & 
     1772      CALL iom_put( 'snowpre'     , sprecip(:,:)                                          )  ! Snow 
     1773      CALL iom_put( 'precip'      , tprecip(:,:)                                          )  ! total  precipitation 
     1774      IF ( iom_use('rain') ) CALL iom_put( 'rain'        , tprecip(:,:) - sprecip(:,:)                           )  ! liquid precipitation  
     1775      IF ( iom_use('snow_ao_cea') ) CALL iom_put( 'snow_ao_cea' , sprecip(:,:) * ( 1._wp - zsnw(:,:) )                  )  ! Snow over ice-free ocean  (cell average) 
     1776      IF ( iom_use('snow_ai_cea') ) CALL iom_put( 'snow_ai_cea' , sprecip(:,:) *           zsnw(:,:)                    )  ! Snow over sea-ice         (cell average) 
     1777      IF ( iom_use('rain_ao_cea') ) CALL iom_put( 'rain_ao_cea' , ( tprecip(:,:) - sprecip(:,:) ) * picefr(:,:)         )  ! liquid precipitation over ocean (cell average) 
     1778      IF ( iom_use('subl_ai_cea') ) CALL iom_put( 'subl_ai_cea' , frcv(jpr_ievp)%z3(:,:,1) * picefr(:,:) * tmask(:,:,1) )  ! Sublimation over sea-ice (cell average) 
     1779      IF ( iom_use('evap_ao_cea') ) CALL iom_put( 'evap_ao_cea' , ( frcv(jpr_tevp)%z3(:,:,1)  & 
    17831780         &                                                        - frcv(jpr_ievp)%z3(:,:,1) * picefr(:,:) ) * tmask(:,:,1) )  ! ice-free oce evap (cell average) 
    17841781      ! note: runoff output is done in sbcrnf (which includes icebergs too) and iceshelf output is done in sbcisf 
     
    18151812! ** NEED TO SORT OUT HOW THIS SHOULD WORK IN THE MULTI-CATEGORY CASE - CURRENTLY NOT ALLOWED WHEN INTERFACE INITIALISED ** 
    18161813         zqns_tot(:,:  ) = frcv(jpr_qnsmix)%z3(:,:,1) 
    1817          zqns_ice(:,:,1) = frcv(jpr_qnsmix)%z3(:,:,1)    & 
    1818             &            + frcv(jpr_dqnsdt)%z3(:,:,1) * ( pist(:,:,1) - ( (rt0 + psst(:,:  ) ) * ziceld(:,:)   & 
    1819             &                                           + pist(:,:,1) * picefr(:,:) ) ) 
     1814         IF ( TRIM(sn_rcv_qsr%clcat) == 'yes' ) THEN 
     1815            DO jl = 1, jpl 
     1816               zqns_ice(:,:,jl) = frcv(jpr_qnsmix)%z3(:,:,jl)    & 
     1817                  &             + frcv(jpr_dqnsdt)%z3(:,:,jl) * ( pist(:,:,jl) - ( ( rt0 + psst(:,:) ) * ziceld(:,:)   & 
     1818                  &                                                                + pist(:,:,jl) * picefr(:,:) ) ) 
     1819            END DO 
     1820         ELSE 
     1821            DO jl = 1, jpl 
     1822               zqns_ice(:,:,jl) = frcv(jpr_qnsmix)%z3(:,:, 1)    & 
     1823                  &             + frcv(jpr_dqnsdt)%z3(:,:, 1) * ( pist(:,:,jl) - ( ( rt0 + psst(:,:) ) * ziceld(:,:)   & 
     1824                  &                                                                + pist(:,:,jl) * picefr(:,:) ) ) 
     1825            END DO 
     1826         ENDIF 
    18201827      END SELECT 
    18211828      !                                      
     
    19021909#endif 
    19031910      ! outputs 
    1904       IF ( srcv(jpr_cal)%laction       ) CALL iom_put('hflx_cal_cea'    , - frcv(jpr_cal)%z3(:,:,1) * rLfus )                      ! latent heat from calving 
    1905       IF ( srcv(jpr_icb)%laction       ) CALL iom_put('hflx_icb_cea'    , - frcv(jpr_icb)%z3(:,:,1) * rLfus )                      ! latent heat from icebergs melting 
    1906       IF ( iom_use('hflx_rain_cea')    ) CALL iom_put('hflx_rain_cea'  , ( tprecip(:,:) - sprecip(:,:) ) * zcptrain(:,:) )        ! heat flux from rain (cell average) 
    1907       IF ( iom_use('hflx_evap_cea')    ) CALL iom_put('hflx_evap_cea'   , ( frcv(jpr_tevp)%z3(:,:,1) - frcv(jpr_ievp)%z3(:,:,1) & 
    1908            &                                                              * picefr(:,:) ) * zcptn(:,:) * tmask(:,:,1) )            ! heat flux from evap (cell average) 
    1909       IF ( iom_use('hflx_snow_cea')    ) CALL iom_put('hflx_snow_cea'   , sprecip(:,:) * ( zcptsnw(:,:) - rLfus )  )               ! heat flux from snow (cell average) 
    1910       IF ( iom_use('hflx_snow_ao_cea') ) CALL iom_put('hflx_snow_ao_cea', sprecip(:,:) * ( zcptsnw(:,:) - rLfus ) & 
    1911            &                                                              * ( 1._wp - zsnw(:,:) )                  )               ! heat flux from snow (over ocean) 
    1912       IF ( iom_use('hflx_snow_ai_cea') ) CALL iom_put('hflx_snow_ai_cea', sprecip(:,:) * ( zcptsnw(:,:) - rLfus ) &  
    1913            &                                                              *           zsnw(:,:)                    )               ! heat flux from snow (over ice) 
     1911      IF ( srcv(jpr_cal)%laction ) CALL iom_put('hflx_cal_cea'    , - frcv(jpr_cal)%z3(:,:,1) * rLfus )   ! latent heat from calving 
     1912      IF ( srcv(jpr_icb)%laction ) CALL iom_put('hflx_icb_cea'    , - frcv(jpr_icb)%z3(:,:,1) * rLfus )   ! latent heat from icebergs melting 
     1913      IF ( iom_use(   'hflx_rain_cea') ) CALL iom_put('hflx_rain_cea' , ( tprecip(:,:) - sprecip(:,:) ) * zcptrain(:,:) )        ! heat flux from rain (cell average) 
     1914      IF ( iom_use(   'hflx_evap_cea') ) CALL iom_put('hflx_evap_cea' , ( frcv(jpr_tevp)%z3(:,:,1) - frcv(jpr_ievp)%z3(:,:,1) * picefr(:,:) ) & 
     1915           &                         * zcptn(:,:) * tmask(:,:,1) )            ! heat flux from evap (cell average) 
     1916      IF ( iom_use(   'hflx_prec_cea') ) CALL iom_put('hflx_prec_cea' ,    sprecip(:,:) * ( zcptsnw(:,:) - rLfus )  &                    ! heat flux from all precip (cell avg) 
     1917         &                          + ( tprecip(:,:) - sprecip(:,:) ) * zcptrain(:,:) ) 
     1918      IF ( iom_use(   'hflx_snow_cea') ) CALL iom_put('hflx_snow_cea'   , sprecip(:,:) * ( zcptsnw(:,:) - rLfus )  )               ! heat flux from snow (cell average) 
     1919      IF ( iom_use('hflx_snow_ao_cea') ) CALL iom_put('hflx_snow_ao_cea', sprecip(:,:) * ( zcptsnw(:,:) - rLfus ) * ( 1._wp - zsnw(:,:) ) )   ! heat flux from snow (over ocean) 
     1920      IF ( iom_use('hflx_snow_ai_cea') ) CALL iom_put('hflx_snow_ai_cea', sprecip(:,:) * ( zcptsnw(:,:) - rLfus ) *  zsnw(:,:) )              ! heat flux from snow (over ice) 
    19141921      ! note: hflx for runoff and iceshelf are done in sbcrnf and sbcisf resp. 
    19151922      ! 
     
    19291936            END DO 
    19301937         ENDIF 
    1931          zqsr_tot(:,:  ) = frcv(jpr_qsrmix)%z3(:,:,1) 
    1932          zqsr_ice(:,:,1) = frcv(jpr_qsrice)%z3(:,:,1) 
    19331938      CASE( 'oce and ice' ) 
    19341939         zqsr_tot(:,:  ) =  ziceld(:,:) * frcv(jpr_qsroce)%z3(:,:,1) 
     
    19501955!       Create solar heat flux over ice using incoming solar heat flux and albedos 
    19511956!       ( see OASIS3 user guide, 5th edition, p39 ) 
    1952          zqsr_ice(:,:,1) = frcv(jpr_qsrmix)%z3(:,:,1) * ( 1.- palbi(:,:,1) )   & 
    1953             &            / (  1.- ( alb_oce_mix(:,:  ) * ziceld(:,:)       & 
    1954             &                     + palbi      (:,:,1) * picefr(:,:) ) ) 
     1957         IF ( TRIM(sn_rcv_qsr%clcat) == 'yes' ) THEN 
     1958            DO jl = 1, jpl 
     1959               zqsr_ice(:,:,jl) = frcv(jpr_qsrmix)%z3(:,:,jl) * ( 1.- palbi(:,:,jl) )   & 
     1960                  &            / (  1.- ( alb_oce_mix(:,:   ) * ziceld(:,:)       & 
     1961                  &                     + palbi      (:,:,jl) * picefr(:,:) ) ) 
     1962            END DO 
     1963         ELSE 
     1964            DO jl = 1, jpl 
     1965               zqsr_ice(:,:,jl) = frcv(jpr_qsrmix)%z3(:,:, 1) * ( 1.- palbi(:,:,jl) )   & 
     1966                  &            / (  1.- ( alb_oce_mix(:,:   ) * ziceld(:,:)       & 
     1967                  &                     + palbi      (:,:,jl) * picefr(:,:) ) ) 
     1968            END DO 
     1969         ENDIF 
    19551970      CASE( 'none'      )       ! Not available as for now: needs additional coding   
    19561971      !                         ! since fields received, here zqsr_tot,  are not defined with none option 
     
    20122027      !                                                      ! ========================= ! 
    20132028      CASE ('coupled') 
    2014          qml_ice(:,:,:) = frcv(jpr_topm)%z3(:,:,:) 
    2015          qcn_ice(:,:,:) = frcv(jpr_botm)%z3(:,:,:) 
     2029         IF( ln_mixcpl ) THEN 
     2030            DO jl=1,jpl 
     2031               qml_ice(:,:,jl) = qml_ice(:,:,jl) * xcplmask(:,:,0) + frcv(jpr_topm)%z3(:,:,jl) * zmsk(:,:) 
     2032               qcn_ice(:,:,jl) = qcn_ice(:,:,jl) * xcplmask(:,:,0) + frcv(jpr_botm)%z3(:,:,jl) * zmsk(:,:) 
     2033            ENDDO 
     2034         ELSE 
     2035            qml_ice(:,:,:) = frcv(jpr_topm)%z3(:,:,:) 
     2036            qcn_ice(:,:,:) = frcv(jpr_botm)%z3(:,:,:) 
     2037         ENDIF 
    20162038      END SELECT 
    2017       ! 
    20182039      !                                                      ! ========================= ! 
    20192040      !                                                      !      Transmitted Qsr      !   [W/m2] 
     
    20222043         ! 
    20232044         !                    ! ===> used prescribed cloud fraction representative for polar oceans in summer (0.81) 
    2024          ztri = 0.18 * ( 1.0 - cldf_ice ) + 0.35 * cldf_ice    ! surface transmission parameter (Grenfell Maykut 77) 
     2045         ztri = 0.18 * ( 1.0 - cldf_ice ) + 0.35 * cldf_ice    ! surface transmission when hi>10cm (Grenfell Maykut 77) 
    20252046         ! 
    2026          qtr_ice_top(:,:,:) = ztri * qsr_ice(:,:,:) 
    2027          WHERE( phs(:,:,:) >= 0.0_wp )   qtr_ice_top(:,:,:) = 0._wp            ! snow fully opaque 
    2028          WHERE( phi(:,:,:) <= 0.1_wp )   qtr_ice_top(:,:,:) = qsr_ice(:,:,:)   ! thin ice transmits all solar radiation 
     2047         WHERE    ( phs(:,:,:) <= 0._wp .AND. phi(:,:,:) <  0.1_wp )       ! linear decrease from hi=0 to 10cm   
     2048            zqtr_ice_top(:,:,:) = qsr_ice(:,:,:) * ( ztri + ( 1._wp - ztri ) * ( 1._wp - phi(:,:,:) * 10._wp ) ) 
     2049         ELSEWHERE( phs(:,:,:) <= 0._wp .AND. phi(:,:,:) >= 0.1_wp )       ! constant (ztri) when hi>10cm 
     2050            zqtr_ice_top(:,:,:) = qsr_ice(:,:,:) * ztri 
     2051         ELSEWHERE                                                         ! zero when hs>0 
     2052            zqtr_ice_top(:,:,:) = 0._wp 
     2053         END WHERE 
    20292054         !      
    20302055      ELSEIF( ln_cndflx .AND. .NOT.ln_cndemulate ) THEN      !==  conduction flux as surface forcing  ==! 
     
    20322057         !                    ! ===> here we must receive the qtr_ice_top array from the coupler 
    20332058         !                           for now just assume zero (fully opaque ice) 
    2034          qtr_ice_top(:,:,:) = 0._wp 
     2059         zqtr_ice_top(:,:,:) = 0._wp 
     2060         ! 
     2061      ENDIF 
     2062      ! 
     2063      IF( ln_mixcpl ) THEN 
     2064         DO jl=1,jpl 
     2065            qtr_ice_top(:,:,jl) = qtr_ice_top(:,:,jl) * xcplmask(:,:,0) + zqtr_ice_top(:,:,jl) * zmsk(:,:) 
     2066         ENDDO 
     2067      ELSE 
     2068         qtr_ice_top(:,:,:) = zqtr_ice_top(:,:,:) 
     2069      ENDIF 
     2070      !                                                      ! ================== ! 
     2071      !                                                      !   ice skin temp.   ! 
     2072      !                                                      ! ================== ! 
     2073      ! needed by Met Office 
     2074      IF( srcv(jpr_ts_ice)%laction ) THEN  
     2075         WHERE    ( frcv(jpr_ts_ice)%z3(:,:,:) > 0.0  )   ;   ztsu(:,:,:) =   0. + rt0  
     2076         ELSEWHERE( frcv(jpr_ts_ice)%z3(:,:,:) < -60. )   ;   ztsu(:,:,:) = -60. + rt0 
     2077         ELSEWHERE                                        ;   ztsu(:,:,:) = frcv(jpr_ts_ice)%z3(:,:,:) + rt0 
     2078         END WHERE 
     2079         ! 
     2080         IF( ln_mixcpl ) THEN 
     2081            DO jl=1,jpl 
     2082               pist(:,:,jl) = pist(:,:,jl) * xcplmask(:,:,0) + ztsu(:,:,jl) * zmsk(:,:) 
     2083            ENDDO 
     2084         ELSE 
     2085            pist(:,:,:) = ztsu(:,:,:) 
     2086         ENDIF 
    20352087         ! 
    20362088      ENDIF 
     
    21952247         CASE default    ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_thick%clcat' ) 
    21962248         END SELECT 
    2197          IF( ssnd(jps_fice)%laction )   CALL cpl_snd( jps_fice, isec, ztmp3, info ) 
     2249         CALL cpl_snd( jps_fice, isec, ztmp3, info ) 
    21982250      ENDIF 
    21992251 
     
    22552307      !                                                      !      Ice melt ponds       !  
    22562308      !                                                      ! ------------------------- ! 
    2257       ! needed by Met Office 
     2309      ! needed by Met Office: 1) fraction of ponded ice 2) local/actual pond depth  
    22582310      IF( ssnd(jps_a_p)%laction .OR. ssnd(jps_ht_p)%laction ) THEN  
    22592311         SELECT CASE( sn_snd_mpnd%cldes)   
     
    22612313            SELECT CASE( sn_snd_mpnd%clcat )   
    22622314            CASE( 'yes' )   
    2263                ztmp3(:,:,1:jpl) =  a_ip(:,:,1:jpl) 
    2264                ztmp4(:,:,1:jpl) =  v_ip(:,:,1:jpl)   
     2315               ztmp3(:,:,1:jpl) =  a_ip_frac(:,:,1:jpl) 
     2316               ztmp4(:,:,1:jpl) =  h_ip(:,:,1:jpl)   
    22652317            CASE( 'no' )   
    22662318               ztmp3(:,:,:) = 0.0   
    22672319               ztmp4(:,:,:) = 0.0   
    22682320               DO jl=1,jpl   
    2269                  ztmp3(:,:,1) = ztmp3(:,:,1) + a_ip(:,:,jpl)   
    2270                  ztmp4(:,:,1) = ztmp4(:,:,1) + v_ip(:,:,jpl)  
     2321                 ztmp3(:,:,1) = ztmp3(:,:,1) + a_ip_frac(:,:,jpl)   
     2322                 ztmp4(:,:,1) = ztmp4(:,:,1) + h_ip(:,:,jpl)  
    22712323               ENDDO   
    22722324            CASE default   ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_mpnd%clcat' )   
     
    23062358      !                                                      !  CO2 flux from PISCES     !  
    23072359      !                                                      ! ------------------------- ! 
    2308       IF( ssnd(jps_co2)%laction .AND. l_co2cpl )   CALL cpl_snd( jps_co2, isec, RESHAPE ( oce_co2, (/jpi,jpj,1/) ) , info ) 
     2360      IF( ssnd(jps_co2)%laction .AND. l_co2cpl )   THEN   
     2361         ztmp1(:,:) = oce_co2(:,:) * 1000.  ! conversion in molC/m2/s  
     2362         CALL cpl_snd( jps_co2, isec, RESHAPE ( ztmp1, (/jpi,jpj,1/) ) , info )  
     2363      ENDIF  
    23092364      ! 
    23102365      !                                                      ! ------------------------- ! 
  • NEMO/branches/UKMO/NEMO_4.0.2_ENHANCE-02_ISF_nemo/src/OCE/SBC/sbcmod.F90

    r12143 r12706  
    236236#endif 
    237237      ! 
     238      ! 
     239      IF( sbc_ssr_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'sbc_init : unable to allocate sbc_ssr arrays' ) 
     240      IF( .NOT.ln_ssr ) THEN               !* Initialize qrp and erp if no restoring  
     241         qrp(:,:) = 0._wp 
     242         erp(:,:) = 0._wp 
     243      ENDIF 
     244      ! 
     245 
    238246      IF( nn_ice == 0 ) THEN        !* No sea-ice in the domain : ice fraction is always zero 
    239247         IF( nn_components /= jp_iam_opa )   fr_i(:,:) = 0._wp    ! except for OPA in SAS-OPA coupled case 
     
    536544         CALL iom_put( "taum"  , taum       )                   ! wind stress module 
    537545         CALL iom_put( "wspd"  , wndm       )                   ! wind speed  module over free ocean or leads in presence of sea-ice 
     546         CALL iom_put( "qrp", qrp )                             ! heat flux damping 
     547         CALL iom_put( "erp", erp )                             ! freshwater flux damping 
    538548      ENDIF 
    539549      ! 
  • NEMO/branches/UKMO/NEMO_4.0.2_ENHANCE-02_ISF_nemo/src/OCE/SBC/sbcrnf.F90

    r12143 r12706  
    2020   USE sbc_oce        ! surface boundary condition variables 
    2121   USE eosbn2         ! Equation Of State 
    22    USE closea        ! closed seas 
     22   USE closea, ONLY: l_clo_rnf, clo_rnf ! closed seas 
    2323   ! 
    2424   USE in_out_manager ! I/O manager 
     
    4242   REAL(wp)                   ::      rn_dep_max        !: depth over which runoffs is spread       (ln_rnf_depth_ini =T) 
    4343   INTEGER                    ::      nn_rnf_depth_file !: create (=1) a runoff depth file or not (=0) 
     44   LOGICAL                    ::   ln_rnf_icb        !: iceberg flux is specified in a file 
    4445   LOGICAL                    ::   ln_rnf_tem        !: temperature river runoffs attribute specified in a file 
    4546   LOGICAL           , PUBLIC ::   ln_rnf_sal        !: salinity    river runoffs attribute specified in a file 
    4647   TYPE(FLD_N)       , PUBLIC ::   sn_rnf            !: information about the runoff file to be read 
    4748   TYPE(FLD_N)                ::   sn_cnf            !: information about the runoff mouth file to be read 
     49   TYPE(FLD_N)                ::   sn_i_rnf        !: information about the iceberg flux file to be read 
    4850   TYPE(FLD_N)                ::   sn_s_rnf          !: information about the salinities of runoff file to be read 
    4951   TYPE(FLD_N)                ::   sn_t_rnf          !: information about the temperatures of runoff file to be read 
     
    6466 
    6567   TYPE(FLD),        ALLOCATABLE, DIMENSION(:) ::   sf_rnf       ! structure: river runoff (file information, fields read) 
     68   TYPE(FLD),        ALLOCATABLE, DIMENSION(:) ::   sf_i_rnf     ! structure: iceberg flux (file information, fields read) 
    6669   TYPE(FLD),        ALLOCATABLE, DIMENSION(:) ::   sf_s_rnf     ! structure: river runoff salinity (file information, fields read)   
    6770   TYPE(FLD),        ALLOCATABLE, DIMENSION(:) ::   sf_t_rnf     ! structure: river runoff temperature (file information, fields read)   
     
    111114      !                                            !-------------------! 
    112115      ! 
    113       IF( .NOT. l_rnfcpl )   CALL fld_read ( kt, nn_fsbc, sf_rnf   )    ! Read Runoffs data and provide it at kt 
     116      ! 
     117      IF( .NOT. l_rnfcpl )  THEN 
     118                            CALL fld_read ( kt, nn_fsbc, sf_rnf   )    ! Read Runoffs data and provide it at kt ( runoffs + iceberg ) 
     119         IF( ln_rnf_icb )   CALL fld_read ( kt, nn_fsbc, sf_i_rnf )    ! idem for iceberg flux if required 
     120      ENDIF 
    114121      IF(   ln_rnf_tem   )   CALL fld_read ( kt, nn_fsbc, sf_t_rnf )    ! idem for runoffs temperature if required 
    115122      IF(   ln_rnf_sal   )   CALL fld_read ( kt, nn_fsbc, sf_s_rnf )    ! idem for runoffs salinity    if required 
     
    117124      IF( MOD( kt - 1, nn_fsbc ) == 0 ) THEN 
    118125         ! 
    119          IF( .NOT. l_rnfcpl )   rnf(:,:) = rn_rfact * ( sf_rnf(1)%fnow(:,:,1) ) * tmask(:,:,1)       ! updated runoff value at time step kt 
     126         IF( .NOT. l_rnfcpl ) THEN 
     127             rnf(:,:) = rn_rfact * ( sf_rnf(1)%fnow(:,:,1) ) * tmask(:,:,1)  ! updated runoff value at time step kt 
     128             IF( ln_rnf_icb ) THEN 
     129                fwficb(:,:) = rn_rfact * ( sf_i_rnf(1)%fnow(:,:,1) ) * tmask(:,:,1)  ! updated runoff value at time step kt 
     130                CALL iom_put( 'iceberg_cea'  , fwficb(:,:)  )         ! output iceberg flux 
     131                CALL iom_put( 'hflx_icb_cea' , fwficb(:,:) * rLfus )   ! output Heat Flux into Sea Water due to Iceberg Thermodynamics --> 
     132             ENDIF 
     133         ENDIF 
    120134         ! 
    121135         !                                                           ! set temperature & salinity content of runoffs 
     
    128142         ELSE                                                        ! use SST as runoffs temperature 
    129143            !CEOD River is fresh water so must at least be 0 unless we consider ice 
    130             rnf_tsc(:,:,jp_tem) = MAX(sst_m(:,:),0.0_wp) * rnf(:,:) * r1_rau0 
     144            rnf_tsc(:,:,jp_tem) = MAX( sst_m(:,:), 0.0_wp ) * rnf(:,:) * r1_rau0 
    131145         ENDIF 
    132146         !                                                           ! use runoffs salinity data 
    133147         IF( ln_rnf_sal )   rnf_tsc(:,:,jp_sal) = ( sf_s_rnf(1)%fnow(:,:,1) ) * rnf(:,:) * r1_rau0 
    134148         !                                                           ! else use S=0 for runoffs (done one for all in the init) 
    135          IF( iom_use('runoffs') )        CALL iom_put( 'runoffs'     , rnf(:,:)                         )   ! output runoff mass flux 
     149                                         CALL iom_put( 'runoffs'     , rnf(:,:)                         )   ! output runoff mass flux 
    136150         IF( iom_use('hflx_rnf_cea') )   CALL iom_put( 'hflx_rnf_cea', rnf_tsc(:,:,jp_tem) * rau0 * rcp )   ! output runoff sensible heat (W/m2) 
    137151      ENDIF 
     
    238252      REAL(wp), DIMENSION(jpi,jpj,2) :: zrnfcl     
    239253      !! 
    240       NAMELIST/namsbc_rnf/ cn_dir            , ln_rnf_depth, ln_rnf_tem, ln_rnf_sal,   & 
    241          &                 sn_rnf, sn_cnf    , sn_s_rnf    , sn_t_rnf  , sn_dep_rnf,   & 
     254      NAMELIST/namsbc_rnf/ cn_dir            , ln_rnf_depth, ln_rnf_tem, ln_rnf_sal, ln_rnf_icb,   & 
     255         &                 sn_rnf, sn_cnf    , sn_i_rnf, sn_s_rnf    , sn_t_rnf  , sn_dep_rnf,   & 
    242256         &                 ln_rnf_mouth      , rn_hrnf     , rn_avt_rnf, rn_rfact,     & 
    243257         &                 ln_rnf_depth_ini  , rn_dep_max  , rn_rnf_max, nn_rnf_depth_file 
     
    261275      !                                   ! ============ 
    262276      ! 
    263       REWIND( numnam_ref )              ! Namelist namsbc_rnf in reference namelist : Runoffs  
     277      REWIND( numnam_ref ) 
    264278      READ  ( numnam_ref, namsbc_rnf, IOSTAT = ios, ERR = 901) 
    265279901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namsbc_rnf in reference namelist' ) 
    266280 
    267       REWIND( numnam_cfg )              ! Namelist namsbc_rnf in configuration namelist : Runoffs 
     281      REWIND( numnam_cfg ) 
    268282      READ  ( numnam_cfg, namsbc_rnf, IOSTAT = ios, ERR = 902 ) 
    269283902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namsbc_rnf in configuration namelist' ) 
     
    295309         IF( sn_rnf%ln_tint ) ALLOCATE( sf_rnf(1)%fdta(jpi,jpj,1,2) ) 
    296310         CALL fld_fill( sf_rnf, (/ sn_rnf /), cn_dir, 'sbc_rnf_init', 'read runoffs data', 'namsbc_rnf', no_print ) 
     311         ! 
     312         IF( ln_rnf_icb ) THEN                      ! Create (if required) sf_i_rnf structure 
     313            IF(lwp) WRITE(numout,*) 
     314            IF(lwp) WRITE(numout,*) '          iceberg flux read in a file' 
     315            ALLOCATE( sf_i_rnf(1), STAT=ierror  ) 
     316            IF( ierror > 0 ) THEN 
     317               CALL ctl_stop( 'sbc_rnf_init: unable to allocate sf_i_rnf structure' )   ;   RETURN 
     318            ENDIF 
     319            ALLOCATE( sf_i_rnf(1)%fnow(jpi,jpj,1)   ) 
     320            IF( sn_i_rnf%ln_tint ) ALLOCATE( sf_i_rnf(1)%fdta(jpi,jpj,1,2) ) 
     321            CALL fld_fill (sf_i_rnf, (/ sn_i_rnf /), cn_dir, 'sbc_rnf_init', 'read iceberg flux data', 'namsbc_rnf' ) 
     322         ELSE 
     323            fwficb(:,:) = 0._wp 
     324         ENDIF 
     325 
    297326      ENDIF 
    298327      ! 
     
    337366               IF( h_rnf(ji,jj) > 0._wp ) THEN 
    338367                  jk = 2 
    339                   DO WHILE ( jk /= mbkt(ji,jj) .AND. gdept_0(ji,jj,jk) < h_rnf(ji,jj) ) ;  jk = jk + 1 
     368                  DO WHILE ( jk < mbkt(ji,jj) .AND. gdept_0(ji,jj,jk) < h_rnf(ji,jj) ) ;  jk = jk + 1 
    340369                  END DO 
    341370                  nk_rnf(ji,jj) = jk 
     
    394423               IF( zrnfcl(ji,jj,1) > 0._wp ) THEN 
    395424                  jk = 2 
    396                   DO WHILE ( jk /= mbkt(ji,jj) .AND. gdept_0(ji,jj,jk) < h_rnf(ji,jj) ) ;  jk = jk + 1 
     425                  DO WHILE ( jk < mbkt(ji,jj) .AND. gdept_0(ji,jj,jk) < h_rnf(ji,jj) ) ;  jk = jk + 1 
    397426                  END DO 
    398427                  nk_rnf(ji,jj) = jk 
     
    435464         !                                      !    - mixed upstream-centered (ln_traadv_cen2=T) 
    436465         ! 
    437          IF ( ln_rnf_depth )   CALL ctl_warn( 'sbc_rnf_init: increased mixing turned on but effects may already',   & 
     466         IF( ln_rnf_depth )   CALL ctl_warn( 'sbc_rnf_init: increased mixing turned on but effects may already',   & 
    438467            &                                              'be spread through depth by ln_rnf_depth'               ) 
    439468         ! 
  • NEMO/branches/UKMO/NEMO_4.0.2_ENHANCE-02_ISF_nemo/src/OCE/SBC/sbcssr.F90

    r12143 r12706  
    3030   PUBLIC   sbc_ssr        ! routine called in sbcmod 
    3131   PUBLIC   sbc_ssr_init   ! routine called in sbcmod 
     32   PUBLIC   sbc_ssr_alloc  ! routine called in sbcmod 
    3233 
    3334   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   erp   !: evaporation damping   [kg/m2/s] 
    3435   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   qrp   !: heat flux damping        [w/m2] 
     36   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   coefice   !: under ice relaxation coefficient 
    3537 
    3638   !                                   !!* Namelist namsbc_ssr * 
     
    4143   LOGICAL         ::   ln_sssr_bnd     ! flag to bound erp term  
    4244   REAL(wp)        ::   rn_sssr_bnd     ! ABS(Max./Min.) value of erp term [mm/day] 
     45   INTEGER         ::   nn_sssr_ice     ! Control of restoring under ice 
    4346 
    4447   REAL(wp) , ALLOCATABLE, DIMENSION(:) ::   buffer   ! Temporary buffer for exchange 
     
    97100                  END DO 
    98101               END DO 
    99                CALL iom_put( "qrp", qrp )                             ! heat flux damping 
     102            ENDIF 
     103            ! 
     104            IF( nn_sssr /= 0 .AND. nn_sssr_ice /= 1 ) THEN 
     105              ! use fraction of ice ( fr_i ) to adjust relaxation under ice if nn_sssr_ice .ne. 1 
     106              ! n.b. coefice is initialised and fixed to 1._wp if nn_sssr_ice = 1 
     107               DO jj = 1, jpj 
     108                  DO ji = 1, jpi 
     109                     SELECT CASE ( nn_sssr_ice ) 
     110                       CASE ( 0 )    ;  coefice(ji,jj) = 1._wp - fr_i(ji,jj)              ! no/reduced damping under ice 
     111                       CASE  DEFAULT ;  coefice(ji,jj) = 1._wp + ( nn_sssr_ice - 1 ) * fr_i(ji,jj) ! reinforced damping (x nn_sssr_ice) under ice ) 
     112                     END SELECT 
     113                  END DO 
     114               END DO 
    100115            ENDIF 
    101116            ! 
     
    105120                  DO ji = 1, jpi 
    106121                     zerp = zsrp * ( 1. - 2.*rnfmsk(ji,jj) )   &      ! No damping in vicinity of river mouths 
     122                        &        *   coefice(ji,jj)            &      ! Optional control of damping under sea-ice 
    107123                        &        * ( sss_m(ji,jj) - sf_sss(1)%fnow(ji,jj,1) ) * tmask(ji,jj,1) 
    108124                     sfx(ji,jj) = sfx(ji,jj) + zerp                 ! salt flux 
     
    110126                  END DO 
    111127               END DO 
    112                CALL iom_put( "erp", erp )                             ! freshwater flux damping 
    113128               ! 
    114129            ELSEIF( nn_sssr == 2 ) THEN                               !* Salinity damping term (volume flux (emp) and associated heat flux (qns) 
     
    118133                  DO ji = 1, jpi                             
    119134                     zerp = zsrp * ( 1. - 2.*rnfmsk(ji,jj) )   &      ! No damping in vicinity of river mouths 
     135                        &        *   coefice(ji,jj)            &      ! Optional control of damping under sea-ice 
    120136                        &        * ( sss_m(ji,jj) - sf_sss(1)%fnow(ji,jj,1) )   & 
    121137                        &        / MAX(  sss_m(ji,jj), 1.e-20   ) * tmask(ji,jj,1) 
     
    126142                  END DO 
    127143               END DO 
    128                CALL iom_put( "erp", erp )                             ! freshwater flux damping 
    129144            ENDIF 
    130145            ! 
     
    154169      CHARACTER(len=100) ::  cn_dir          ! Root directory for location of ssr files 
    155170      TYPE(FLD_N) ::   sn_sst, sn_sss        ! informations about the fields to be read 
    156       NAMELIST/namsbc_ssr/ cn_dir, nn_sstr, nn_sssr, rn_dqdt, rn_deds, sn_sst, sn_sss, ln_sssr_bnd, rn_sssr_bnd 
     171      NAMELIST/namsbc_ssr/ cn_dir, nn_sstr, nn_sssr, rn_dqdt, rn_deds, sn_sst, & 
     172              & sn_sss, ln_sssr_bnd, rn_sssr_bnd, nn_sssr_ice 
    157173      INTEGER     ::  ios 
    158174      !!---------------------------------------------------------------------- 
     
    182198         WRITE(numout,*) '         flag to bound erp term                 ln_sssr_bnd = ', ln_sssr_bnd 
    183199         WRITE(numout,*) '         ABS(Max./Min.) erp threshold           rn_sssr_bnd = ', rn_sssr_bnd, ' mm/day' 
    184       ENDIF 
    185       ! 
    186       !                            !* Allocate erp and qrp array 
    187       ALLOCATE( qrp(jpi,jpj), erp(jpi,jpj), STAT=ierror ) 
    188       IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_ssr: unable to allocate erp and qrp array' ) 
     200         WRITE(numout,*) '      Cntrl of surface restoration under ice nn_sssr_ice    = ', nn_sssr_ice 
     201         WRITE(numout,*) '          ( 0 = no restoration under ice)' 
     202         WRITE(numout,*) '          ( 1 = restoration everywhere  )' 
     203         WRITE(numout,*) '          (>1 = enhanced restoration under ice  )' 
     204      ENDIF 
    189205      ! 
    190206      IF( nn_sstr == 1 ) THEN      !* set sf_sst structure & allocate arrays 
     
    216232      ENDIF 
    217233      ! 
     234      coefice(:,:) = 1._wp         !  Initialise coefice to 1._wp ; will not need to be changed if nn_sssr_ice=1 
    218235      !                            !* Initialize qrp and erp if no restoring  
    219236      IF( nn_sstr /= 1                   )   qrp(:,:) = 0._wp 
     
    221238      ! 
    222239   END SUBROUTINE sbc_ssr_init 
     240          
     241   INTEGER FUNCTION sbc_ssr_alloc() 
     242      !!---------------------------------------------------------------------- 
     243      !!               ***  FUNCTION sbc_ssr_alloc  *** 
     244      !!---------------------------------------------------------------------- 
     245      sbc_ssr_alloc = 0       ! set to zero if no array to be allocated 
     246      IF( .NOT. ALLOCATED( erp ) ) THEN 
     247         ALLOCATE( qrp(jpi,jpj), erp(jpi,jpj), coefice(jpi,jpj), STAT= sbc_ssr_alloc ) 
     248         ! 
     249         IF( lk_mpp                  )   CALL mpp_sum ( 'sbcssr', sbc_ssr_alloc ) 
     250         IF( sbc_ssr_alloc /= 0 )   CALL ctl_warn('sbc_ssr_alloc: failed to allocate arrays.') 
     251         ! 
     252      ENDIF 
     253   END FUNCTION 
    223254       
    224255   !!====================================================================== 
Note: See TracChangeset for help on using the changeset viewer.