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 5357 for branches/2015 – NEMO

Changeset 5357 for branches/2015


Ignore:
Timestamp:
2015-06-04T20:39:20+02:00 (9 years ago)
Author:
clem
Message:

LIM3: change the interface between the ice and atm for both coupled and forced modes. Some work still needs to be done to deal with sublimation in coupled mode.

Location:
branches/2015/dev_r5218_CNRS17_coupling/NEMOGCM/NEMO
Files:
15 edited

Legend:

Unmodified
Added
Removed

  • branches/2015/dev_r5218_CNRS17_coupling/NEMOGCM/NEMO/LIM_SRC_3/limitd_th.F90

    r5202 r5357  
    9191      !!------------------------------------------------------------------ 
    9292 
    93       CALL wrk_alloc( jpi,jpj, zremap_flag )    ! integer 
    94       CALL wrk_alloc( jpi,jpj,jpl-1, zdonor )   ! integer 
     93      CALL wrk_alloc( jpi,jpj, zremap_flag ) 
     94      CALL wrk_alloc( jpi,jpj,jpl-1, zdonor ) 
    9595      CALL wrk_alloc( jpi,jpj,jpl, zdhice, g0, g1, hL, hR, zht_i_b, dummy_es ) 
    9696      CALL wrk_alloc( jpi,jpj,jpl-1, zdaice, zdvice )    
    9797      CALL wrk_alloc( jpi,jpj,jpl+1, zhbnew, kkstart = 0 )    
    9898      CALL wrk_alloc( (jpi+1)*(jpj+1), zvetamin, zvetamax )    
    99       CALL wrk_alloc( (jpi+1)*(jpj+1), nind_i, nind_j )   ! integer  
     99      CALL wrk_alloc( (jpi+1)*(jpj+1), nind_i, nind_j )  
    100100      CALL wrk_alloc( jpi,jpj, zhb0,zhb1,vt_i_init,vt_i_final,vt_s_init,vt_s_final,et_i_init,et_i_final,et_s_init,et_s_final ) 
    101101 
    … …  
    128128               rswitch           = MAX( 0.0, SIGN( 1.0, a_i(ji,jj,jl) - epsi10 ) )     !0 if no ice and 1 if yes 
    129129               ht_i(ji,jj,jl)    = v_i(ji,jj,jl) / MAX( a_i(ji,jj,jl), epsi10 ) * rswitch 
    130                rswitch           = MAX( 0.0, SIGN( 1.0, a_i_b(ji,jj,jl) - epsi10) )    !0 if no ice and 1 if yes 
     130               rswitch           = MAX( 0.0, SIGN( 1.0, a_i_b(ji,jj,jl) - epsi10) ) 
    131131               zht_i_b(ji,jj,jl) = v_i_b(ji,jj,jl) / MAX( a_i_b(ji,jj,jl), epsi10 ) * rswitch 
    132                zdhice(ji,jj,jl) = ht_i(ji,jj,jl) - zht_i_b(ji,jj,jl)  
     132               IF( a_i(ji,jj,jl) > epsi10 )   zdhice(ji,jj,jl) = ht_i(ji,jj,jl) - zht_i_b(ji,jj,jl) ! clem: useless IF statement?  
    133133            END DO 
    134134         END DO 
    … …  
    172172            ! 
    173173            zhbnew(ii,ij,jl) = hi_max(jl) 
    174             IF ( a_i_b(ii,ij,jl) > epsi10 .AND. a_i_b(ii,ij,jl+1) > epsi10 ) THEN 
     174            IF    ( a_i_b(ii,ij,jl) > epsi10 .AND. a_i_b(ii,ij,jl+1) > epsi10 ) THEN 
    175175               !interpolate between adjacent category growth rates 
    176176               zslope           = ( zdhice(ii,ij,jl+1) - zdhice(ii,ij,jl) ) / ( zht_i_b(ii,ij,jl+1) - zht_i_b(ii,ij,jl) ) 
    177177               zhbnew(ii,ij,jl) = hi_max(jl) + zdhice(ii,ij,jl) + zslope * ( hi_max(jl) - zht_i_b(ii,ij,jl) ) 
    178             ELSEIF ( a_i_b(ii,ij,jl) > epsi10) THEN 
     178            ELSEIF( a_i_b(ii,ij,jl) > epsi10) THEN 
    179179               zhbnew(ii,ij,jl) = hi_max(jl) + zdhice(ii,ij,jl) 
    180             ELSEIF ( a_i_b(ii,ij,jl+1) > epsi10) THEN 
     180            ELSEIF( a_i_b(ii,ij,jl+1) > epsi10) THEN 
    181181               zhbnew(ii,ij,jl) = hi_max(jl) + zdhice(ii,ij,jl+1) 
    182182            ENDIF 
    … …  
    187187            ii = nind_i(ji) 
    188188            ij = nind_j(ji) 
    189             IF( a_i(ii,ij,jl) > epsi10 .AND. ht_i(ii,ij,jl) >= zhbnew(ii,ij,jl) ) THEN 
     189 
     190            ! clem: we do not want ht_i to be too close to either HR or HL otherwise a division by nearly 0 is possible  
     191            ! in lim_itd_fitline in the case (HR-HL) = 3(Hice - HL) or = 3(HR - Hice) 
     192            IF    ( a_i(ii,ij,jl  ) > epsi10 .AND. ht_i(ii,ij,jl  ) > ( zhbnew(ii,ij,jl) - epsi10 ) ) THEN 
    190193               zremap_flag(ii,ij) = 0 
    191             ELSEIF( a_i(ii,ij,jl+1) > epsi10 .AND. ht_i(ii,ij,jl+1) <= zhbnew(ii,ij,jl) ) THEN 
     194            ELSEIF( a_i(ii,ij,jl+1) > epsi10 .AND. ht_i(ii,ij,jl+1) < ( zhbnew(ii,ij,jl) + epsi10 ) ) THEN 
    192195               zremap_flag(ii,ij) = 0 
    193196            ENDIF 
    194197 
    195198            !- 4.3 Check that each zhbnew does not exceed maximal values hi_max   
     199            IF( zhbnew(ii,ij,jl) < hi_max(jl-1) ) zremap_flag(ii,ij) = 0 
    196200            IF( zhbnew(ii,ij,jl) > hi_max(jl+1) ) zremap_flag(ii,ij) = 0 
    197             IF( zhbnew(ii,ij,jl) < hi_max(jl-1) ) zremap_flag(ii,ij) = 0 
     201            ! clem bug: why is not the following instead? 
     202            !!IF( zhbnew(ii,ij,jl) < hi_max(jl-1) ) zremap_flag(ii,ij) = 0 
     203            !!IF( zhbnew(ii,ij,jl) > hi_max(jl  ) ) zremap_flag(ii,ij) = 0 
     204  
    198205         END DO 
    199206 
    … …  
    219226      DO jj = 1, jpj 
    220227         DO ji = 1, jpi 
    221             zhb0(ji,jj) = hi_max(0) ! 0eme 
    222             zhb1(ji,jj) = hi_max(1) ! 1er 
    223  
    224             zhbnew(ji,jj,klbnd-1) = 0._wp 
     228            zhb0(ji,jj) = hi_max(0) 
     229            zhb1(ji,jj) = hi_max(1) 
    225230 
    226231            IF( a_i(ji,jj,kubnd) > epsi10 ) THEN 
    227232               zhbnew(ji,jj,kubnd) = MAX( hi_max(kubnd-1), 3._wp * ht_i(ji,jj,kubnd) - 2._wp * zhbnew(ji,jj,kubnd-1) ) 
    228233            ELSE 
    229                zhbnew(ji,jj,kubnd) = hi_max(kubnd)   
    230                !!? clem bug: since hi_max(jpl)=99, this limit is very high  
    231                !!? but I think it is erased in fitline subroutine  
     234!clem bug               zhbnew(ji,jj,kubnd) = hi_max(kubnd)   
     235               zhbnew(ji,jj,kubnd) = hi_max(kubnd-1) ! not used anyway 
     236            ENDIF 
     237 
     238            ! clem: we do not want ht_i_b to be too close to either HR or HL otherwise a division by nearly 0 is possible  
     239            ! in lim_itd_fitline in the case (HR-HL) = 3(Hice - HL) or = 3(HR - Hice) 
     240            IF    ( zht_i_b(ji,jj,klbnd) < ( zhb0(ji,jj) + epsi10 ) )  THEN 
     241               zremap_flag(ji,jj) = 0 
     242            ELSEIF( zht_i_b(ji,jj,klbnd) > ( zhb1(ji,jj) - epsi10 ) )  THEN 
     243               zremap_flag(ji,jj) = 0 
    232244            ENDIF 
    233245 
    … …  
    248260 
    249261         IF( a_i(ii,ij,klbnd) > epsi10 ) THEN 
     262 
    250263            zdh0 = zdhice(ii,ij,klbnd) !decrease of ice thickness in the lower category 
    251             IF( zdh0 < 0.0 ) THEN !remove area from category 1 
     264 
     265            IF( zdh0 < 0.0 ) THEN      !remove area from category 1 
    252266               zdh0 = MIN( -zdh0, hi_max(klbnd) ) 
    253  
    254267               !Integrate g(1) from 0 to dh0 to estimate area melted 
    255268               zetamax = MIN( zdh0, hR(ii,ij,klbnd) ) - hL(ii,ij,klbnd) 
     269 
    256270               IF( zetamax > 0.0 ) THEN 
    257                   zx1  = zetamax 
    258                   zx2  = 0.5 * zetamax * zetamax  
    259                   zda0 = g1(ii,ij,klbnd) * zx2 + g0(ii,ij,klbnd) * zx1 !ice area removed 
    260                   ! Constrain new thickness <= ht_i 
    261                   zdamax = a_i(ii,ij,klbnd) * (1.0 - ht_i(ii,ij,klbnd) / zht_i_b(ii,ij,klbnd) ) ! zdamax > 0 
    262                   !ice area lost due to melting of thin ice 
    263                   zda0   = MIN( zda0, zdamax ) 
    264  
     271                  zx1    = zetamax 
     272                  zx2    = 0.5 * zetamax * zetamax  
     273                  zda0   = g1(ii,ij,klbnd) * zx2 + g0(ii,ij,klbnd) * zx1                        ! ice area removed 
     274                  zdamax = a_i(ii,ij,klbnd) * (1.0 - ht_i(ii,ij,klbnd) / zht_i_b(ii,ij,klbnd) ) ! Constrain new thickness <= ht_i                 
     275                  zda0   = MIN( zda0, zdamax )                                                  ! ice area lost due to melting  
     276                                                                                                !     of thin ice (zdamax > 0) 
    265277                  ! Remove area, conserving volume 
    266278                  ht_i(ii,ij,klbnd) = ht_i(ii,ij,klbnd) * a_i(ii,ij,klbnd) / ( a_i(ii,ij,klbnd) - zda0 ) 
    … …  
    269281               ENDIF 
    270282 
    271             ELSE ! if ice accretion ! a_i > epsi10; zdh0 > 0 
     283            ELSE ! if ice accretion zdh0 > 0 
     284               ! zhbnew was 0, and is shifted to the right to account for thin ice growth in openwater (F0 = f1) 
    272285               zhbnew(ii,ij,klbnd-1) = MIN( zdh0, hi_max(klbnd) )  
    273                ! zhbnew was 0, and is shifted to the right to account for thin ice 
    274                ! growth in openwater (F0 = f1) 
    275             ENDIF ! zdh0  
    276  
    277          ENDIF ! a_i > epsi10 
     286            ENDIF 
     287 
     288         ENDIF 
    278289 
    279290      END DO 
    … …  
    303314 
    304315            IF (zhbnew(ii,ij,jl) > hi_max(jl)) THEN ! transfer from jl to jl+1 
    305  
    306316               ! left and right integration limits in eta space 
    307317               zvetamin(ji) = MAX( hi_max(jl), hL(ii,ij,jl) ) - hL(ii,ij,jl) 
    308                zvetamax(ji) = MIN (zhbnew(ii,ij,jl), hR(ii,ij,jl) ) - hL(ii,ij,jl) 
     318               zvetamax(ji) = MIN( zhbnew(ii,ij,jl), hR(ii,ij,jl) ) - hL(ii,ij,jl) 
    309319               zdonor(ii,ij,jl) = jl 
    310320 
    311             ELSE  ! zhbnew(jl) <= hi_max(jl) ; transfer from jl+1 to jl 
    312  
     321            ELSE                                    ! zhbnew(jl) <= hi_max(jl) ; transfer from jl+1 to jl 
    313322               ! left and right integration limits in eta space 
    314323               zvetamin(ji) = 0.0 
    … …  
    316325               zdonor(ii,ij,jl) = jl + 1 
    317326 
    318             ENDIF  ! zhbnew(jl) > hi_max(jl) 
     327            ENDIF 
    319328 
    320329            zetamax = MAX( zvetamax(ji), zvetamin(ji) ) ! no transfer if etamax < etamin 
    … …  
    333342 
    334343         END DO 
    335       END DO ! jl klbnd -> kubnd - 1 
     344      END DO 
    336345 
    337346      !!---------------------------------------------------------------------------------------------- 
    … …  
    375384      ENDIF 
    376385 
    377       CALL wrk_dealloc( jpi,jpj, zremap_flag )    ! integer 
    378       CALL wrk_dealloc( jpi,jpj,jpl-1, zdonor )   ! integer 
     386      CALL wrk_dealloc( jpi,jpj, zremap_flag ) 
     387      CALL wrk_dealloc( jpi,jpj,jpl-1, zdonor ) 
    379388      CALL wrk_dealloc( jpi,jpj,jpl, zdhice, g0, g1, hL, hR, zht_i_b, dummy_es ) 
    380389      CALL wrk_dealloc( jpi,jpj,jpl-1, zdaice, zdvice )    
    381390      CALL wrk_dealloc( jpi,jpj,jpl+1, zhbnew, kkstart = 0 )    
    382391      CALL wrk_dealloc( (jpi+1)*(jpj+1), zvetamin, zvetamax )    
    383       CALL wrk_dealloc( (jpi+1)*(jpj+1), nind_i, nind_j )   ! integer  
     392      CALL wrk_dealloc( (jpi+1)*(jpj+1), nind_i, nind_j )  
    384393      CALL wrk_dealloc( jpi,jpj, zhb0,zhb1,vt_i_init,vt_i_final,vt_s_init,vt_s_final,et_i_init,et_i_final,et_s_init,et_s_final ) 
    385394 
    … …  
    406415      INTEGER , DIMENSION(jpi,jpj), INTENT(in   ) ::   zremap_flag  ! 
    407416      ! 
    408       INTEGER ::   ji,jj           ! horizontal indices 
     417      INTEGER  ::   ji,jj        ! horizontal indices 
    409418      REAL(wp) ::   zh13         ! HbL + 1/3 * (HbR - HbL) 
    410419      REAL(wp) ::   zh23         ! HbL + 2/3 * (HbR - HbL) 
    … …  
    413422      !!------------------------------------------------------------------ 
    414423      ! 
    415       ! 
    416424      DO jj = 1, jpj 
    417425         DO ji = 1, jpi 
    418426            ! 
    419427            IF( zremap_flag(ji,jj) == 1 .AND. a_i(ji,jj,num_cat) > epsi10   & 
    420                &                        .AND. hice(ji,jj)        > 0._wp     ) THEN 
     428               &                        .AND. hice(ji,jj)        > 0._wp ) THEN 
    421429 
    422430               ! Initialize hL and hR 
    423  
    424431               hL(ji,jj) = HbL(ji,jj) 
    425432               hR(ji,jj) = HbR(ji,jj) 
    426433 
    427434               ! Change hL or hR if hice falls outside central third of range 
    428  
    429435               zh13 = 1.0 / 3.0 * ( 2.0 * hL(ji,jj) + hR(ji,jj) ) 
    430436               zh23 = 1.0 / 3.0 * ( hL(ji,jj) + 2.0 * hR(ji,jj) ) 
    … …  
    435441 
    436442               ! Compute coefficients of g(eta) = g0 + g1*eta 
    437  
    438443               zdhr = 1._wp / (hR(ji,jj) - hL(ji,jj)) 
    439444               zwk1 = 6._wp * a_i(ji,jj,num_cat) * zdhr 
    … …  
    442447               g1(ji,jj) = 2._wp * zdhr * zwk1 * ( zwk2 - 0.5 ) 
    443448               ! 
    444             ELSE                   ! remap_flag = .false. or a_i < epsi10  
     449            ELSE  ! remap_flag = .false. or a_i < epsi10  
    445450               hL(ji,jj) = 0._wp 
    446451               hR(ji,jj) = 0._wp 
    447452               g0(ji,jj) = 0._wp 
    448453               g1(ji,jj) = 0._wp 
    449             ENDIF                  ! a_i > epsi10 
     454            ENDIF 
    450455            ! 
    451456         END DO 
    … …  
    471476 
    472477      INTEGER ::   ji, jj, jl, jl2, jl1, jk   ! dummy loop indices 
    473       INTEGER ::   ii, ij          ! indices when changing from 2D-1D is done 
     478      INTEGER ::   ii, ij                     ! indices when changing from 2D-1D is done 
    474479 
    475480      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zaTsfn 
    … …  
    484489      INTEGER, POINTER, DIMENSION(:) ::   nind_i, nind_j   ! compressed indices for i/j directions 
    485490 
    486       INTEGER ::   nbrem             ! number of cells with ice to transfer 
    487  
    488       LOGICAL ::   zdaice_negative         ! true if daice < -puny 
    489       LOGICAL ::   zdvice_negative         ! true if dvice < -puny 
    490       LOGICAL ::   zdaice_greater_aicen    ! true if daice > aicen 
    491       LOGICAL ::   zdvice_greater_vicen    ! true if dvice > vicen 
     491      INTEGER  ::   nbrem             ! number of cells with ice to transfer 
    492492      !!------------------------------------------------------------------ 
    493493 
    494494      CALL wrk_alloc( jpi,jpj,jpl, zaTsfn ) 
    495495      CALL wrk_alloc( jpi,jpj, zworka ) 
    496       CALL wrk_alloc( (jpi+1)*(jpj+1), nind_i, nind_j )   ! integer 
     496      CALL wrk_alloc( (jpi+1)*(jpj+1), nind_i, nind_j ) 
    497497 
    498498      !---------------------------------------------------------------------------------------------- 
    … …  
    504504      END DO 
    505505 
    506 !clem: I think the following is wrong (if enabled, it creates negative concentration/volume around -epsi10) 
    507 !      !---------------------------------------------------------------------------------------------- 
    508 !      ! 2) Check for daice or dvice out of range, allowing for roundoff error 
    509 !      !---------------------------------------------------------------------------------------------- 
    510 !      ! Note: zdaice < 0 or zdvice < 0 usually happens when category jl 
    511 !      ! has a small area, with h(n) very close to a boundary.  Then 
    512 !      ! the coefficients of g(h) are large, and the computed daice and 
    513 !      ! dvice can be in error. If this happens, it is best to transfer 
    514 !      ! either the entire category or nothing at all, depending on which 
    515 !      ! side of the boundary hice(n) lies. 
    516 !      !----------------------------------------------------------------- 
    517 !      DO jl = klbnd, kubnd-1 
    518 ! 
    519 !         zdaice_negative = .false. 
    520 !         zdvice_negative = .false. 
    521 !         zdaice_greater_aicen = .false. 
    522 !         zdvice_greater_vicen = .false. 
    523 ! 
    524 !         DO jj = 1, jpj 
    525 !            DO ji = 1, jpi 
    526 ! 
    527 !               IF (zdonor(ji,jj,jl) > 0) THEN 
    528 !                  jl1 = zdonor(ji,jj,jl) 
    529 ! 
    530 !                  IF (zdaice(ji,jj,jl) < 0.0) THEN 
    531 !                     IF (zdaice(ji,jj,jl) > -epsi10) THEN 
    532 !                        IF ( ( jl1 == jl   .AND. ht_i(ji,jj,jl1) >  hi_max(jl) ) .OR.  & 
    533 !                             ( jl1 == jl+1 .AND. ht_i(ji,jj,jl1) <= hi_max(jl) ) ) THEN 
    534 !                           zdaice(ji,jj,jl) = a_i(ji,jj,jl1)  ! shift entire category 
    535 !                           zdvice(ji,jj,jl) = v_i(ji,jj,jl1) 
    536 !                        ELSE 
    537 !                           zdaice(ji,jj,jl) = 0.0 ! shift no ice 
    538 !                           zdvice(ji,jj,jl) = 0.0 
    539 !                        ENDIF 
    540 !                     ELSE 
    541 !                        zdaice_negative = .true. 
    542 !                     ENDIF 
    543 !                  ENDIF 
    544 ! 
    545 !                  IF (zdvice(ji,jj,jl) < 0.0) THEN 
    546 !                     IF (zdvice(ji,jj,jl) > -epsi10 ) THEN 
    547 !                        IF ( ( jl1 == jl   .AND. ht_i(ji,jj,jl1) >  hi_max(jl) ) .OR.  & 
    548 !                             ( jl1 == jl+1 .AND. ht_i(ji,jj,jl1) <= hi_max(jl) ) ) THEN 
    549 !                           zdaice(ji,jj,jl) = a_i(ji,jj,jl1) ! shift entire category 
    550 !                           zdvice(ji,jj,jl) = v_i(ji,jj,jl1)  
    551 !                        ELSE 
    552 !                           zdaice(ji,jj,jl) = 0.0    ! shift no ice 
    553 !                           zdvice(ji,jj,jl) = 0.0 
    554 !                        ENDIF 
    555 !                    ELSE 
    556 !                       zdvice_negative = .true. 
    557 !                    ENDIF 
    558 !                 ENDIF 
    559 ! 
    560 !                  ! If daice is close to aicen, set daice = aicen. 
    561 !                  IF (zdaice(ji,jj,jl) > a_i(ji,jj,jl1) - epsi10 ) THEN 
    562 !                     IF (zdaice(ji,jj,jl) < a_i(ji,jj,jl1)+epsi10) THEN 
    563 !                        zdaice(ji,jj,jl) = a_i(ji,jj,jl1) 
    564 !                        zdvice(ji,jj,jl) = v_i(ji,jj,jl1)  
    565 !                    ELSE 
    566 !                       zdaice_greater_aicen = .true. 
    567 !                    ENDIF 
    568 !                  ENDIF 
    569 ! 
    570 !                  IF (zdvice(ji,jj,jl) > v_i(ji,jj,jl1)-epsi10) THEN 
    571 !                     IF (zdvice(ji,jj,jl) < v_i(ji,jj,jl1)+epsi10) THEN 
    572 !                        zdaice(ji,jj,jl) = a_i(ji,jj,jl1) 
    573 !                        zdvice(ji,jj,jl) = v_i(ji,jj,jl1)  
    574 !                     ELSE 
    575 !                        zdvice_greater_vicen = .true. 
    576 !                     ENDIF 
    577 !                  ENDIF 
    578 ! 
    579 !               ENDIF               ! donor > 0 
    580 !            END DO 
    581 !         END DO 
    582 ! 
    583 !      END DO 
    584 !clem 
    585506      !------------------------------------------------------------------------------- 
    586       ! 3) Transfer volume and energy between categories 
     507      ! 2) Transfer volume and energy between categories 
    587508      !------------------------------------------------------------------------------- 
    588509 
    … …  
    604525 
    605526            jl1 = zdonor(ii,ij,jl) 
    606             rswitch       = MAX( 0._wp , SIGN( 1._wp , v_i(ii,ij,jl1) - epsi20 ) ) 
    607             zworka(ii,ij) = zdvice(ii,ij,jl) / MAX( v_i(ii,ij,jl1), epsi20 ) * rswitch 
     527            rswitch       = MAX( 0._wp , SIGN( 1._wp , v_i(ii,ij,jl1) - epsi10 ) ) 
     528            zworka(ii,ij) = zdvice(ii,ij,jl) / MAX( v_i(ii,ij,jl1), epsi10 ) * rswitch 
    608529            IF( jl1 == jl) THEN   ;   jl2 = jl1+1 
    609530            ELSE                  ;   jl2 = jl  
    … …  
    613534            ! Ice areas 
    614535            !-------------- 
    615  
    616536            a_i(ii,ij,jl1) = a_i(ii,ij,jl1) - zdaice(ii,ij,jl) 
    617537            a_i(ii,ij,jl2) = a_i(ii,ij,jl2) + zdaice(ii,ij,jl) 
    … …  
    620540            ! Ice volumes 
    621541            !-------------- 
    622  
    623542            v_i(ii,ij,jl1) = v_i(ii,ij,jl1) - zdvice(ii,ij,jl)  
    624543            v_i(ii,ij,jl2) = v_i(ii,ij,jl2) + zdvice(ii,ij,jl) 
    … …  
    627546            ! Snow volumes 
    628547            !-------------- 
    629  
    630548            zdvsnow        = v_s(ii,ij,jl1) * zworka(ii,ij) 
    631549            v_s(ii,ij,jl1) = v_s(ii,ij,jl1) - zdvsnow 
    … …  
    635553            ! Snow heat content   
    636554            !-------------------- 
    637  
    638555            zdesnow            = e_s(ii,ij,1,jl1) * zworka(ii,ij) 
    639556            e_s(ii,ij,1,jl1)   = e_s(ii,ij,1,jl1) - zdesnow 
    … …  
    643560            ! Ice age  
    644561            !-------------- 
    645  
    646562            zdo_aice           = oa_i(ii,ij,jl1) * zdaice(ii,ij,jl) 
    647563            oa_i(ii,ij,jl1)    = oa_i(ii,ij,jl1) - zdo_aice 
    … …  
    651567            ! Ice salinity 
    652568            !-------------- 
    653  
    654569            zdsm_vice          = smv_i(ii,ij,jl1) * zworka(ii,ij) 
    655570            smv_i(ii,ij,jl1)   = smv_i(ii,ij,jl1) - zdsm_vice 
    … …  
    659574            ! Surface temperature 
    660575            !--------------------- 
    661  
    662576            zdaTsf             = t_su(ii,ij,jl1) * zdaice(ii,ij,jl) 
    663577            zaTsfn(ii,ij,jl1)  = zaTsfn(ii,ij,jl1) - zdaTsf 
    … …  
    710624      CALL wrk_dealloc( jpi,jpj,jpl, zaTsfn ) 
    711625      CALL wrk_dealloc( jpi,jpj, zworka ) 
    712       CALL wrk_dealloc( (jpi+1)*(jpj+1), nind_i, nind_j )   ! integer 
     626      CALL wrk_dealloc( (jpi+1)*(jpj+1), nind_i, nind_j ) 
    713627      ! 
    714628   END SUBROUTINE lim_itd_shiftice 
    … …  
    859773      ENDIF 
    860774      ! 
    861       CALL wrk_dealloc( jpi,jpj,jpl, zdonor )   ! interger 
     775      CALL wrk_dealloc( jpi,jpj,jpl, zdonor ) 
    862776      CALL wrk_dealloc( jpi,jpj,jpl, zdaice, zdvice ) 
    863777      CALL wrk_dealloc( jpi,jpj, vt_i_init, vt_i_final, vt_s_init, vt_s_final ) 

  • branches/2015/dev_r5218_CNRS17_coupling/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90

    r5352 r5357  
    9494      !!              - fr_i    : ice fraction 
    9595      !!              - tn_ice  : sea-ice surface temperature 
    96       !!              - alb_ice : sea-ice albedo (ln_cpl=T) 
     96      !!              - alb_ice : sea-ice albedo (only useful in coupled mode) 
    9797      !! 
    9898      !! References : Goosse, H. et al. 1996, Bul. Soc. Roy. Sc. Liege, 65, 87-90. 
    … …  
    101101      !!              The ref should be Rousset et al., 2015 
    102102      !!--------------------------------------------------------------------- 
    103       INTEGER, INTENT(in) ::   kt                                   ! number of iteration 
    104       INTEGER  ::   ji, jj, jl, jk                                  ! dummy loop indices 
    105       REAL(wp) ::   zemp                                            ! local scalars 
    106       REAL(wp) ::   zf_mass                                         ! Heat flux associated with mass exchange ice->ocean (W.m-2) 
    107       REAL(wp) ::   zfcm1                                           ! New solar flux received by the ocean 
     103      INTEGER, INTENT(in) ::   kt                                  ! number of iteration 
     104      INTEGER  ::   ji, jj, jl, jk                                 ! dummy loop indices 
     105      REAL(wp) ::   zqmass                                         ! Heat flux associated with mass exchange ice->ocean (W.m-2) 
     106      REAL(wp) ::   zqsr                                           ! New solar flux received by the ocean 
    108107      ! 
    109108      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zalb_cs, zalb_os     ! 2D/3D workspace 
    … …  
    126125            !      heat flux at the ocean surface      ! 
    127126            !------------------------------------------! 
    128             ! Solar heat flux reaching the ocean = zfcm1 (W.m-2)  
     127            ! Solar heat flux reaching the ocean = zqsr (W.m-2)  
    129128            !--------------------------------------------------- 
    130             IF( ln_cpl ) THEN  
    131                !!! LIM2 version zqsr = qsr_tot(ji,jj) + ( fstric(ji,jj) - qsr_ice(ji,jj,1) ) * ( 1.0 - pfrld(ji,jj) ) 
    132                zfcm1 = qsr_tot(ji,jj) 
    133                DO jl = 1, jpl 
    134                   zfcm1 = zfcm1 + ( ftr_ice(ji,jj,jl) - qsr_ice(ji,jj,jl) ) * a_i_b(ji,jj,jl) 
    135                END DO 
    136             ELSE 
    137                !!! LIM2 version zqsr = pfrld(ji,jj) * qsr(ji,jj)  + ( 1.  - pfrld(ji,jj) ) * fstric(ji,jj) 
    138                zfcm1   = pfrld(ji,jj) * qsr(ji,jj) 
    139                DO jl = 1, jpl 
    140                   zfcm1   = zfcm1 + a_i_b(ji,jj,jl) * ftr_ice(ji,jj,jl) 
    141                END DO 
    142             ENDIF 
     129            zqsr = qsr_tot(ji,jj) 
     130            DO jl = 1, jpl 
     131               zqsr = zqsr - a_i_b(ji,jj,jl) * (  qsr_ice(ji,jj,jl) - ftr_ice(ji,jj,jl) )  
     132            END DO 
    143133 
    144134            ! Total heat flux reaching the ocean = hfx_out (W.m-2)  
    145135            !--------------------------------------------------- 
    146             zf_mass        = hfx_thd(ji,jj) + hfx_dyn(ji,jj) + hfx_res(ji,jj) ! heat flux from snow is 0 (T=0 degC) 
    147             hfx_out(ji,jj) = hfx_out(ji,jj) + zf_mass + zfcm1 
     136            zqmass         = hfx_thd(ji,jj) + hfx_dyn(ji,jj) + hfx_res(ji,jj) ! heat flux from snow is 0 (T=0 degC) 
     137            hfx_out(ji,jj) = hfx_out(ji,jj) + zqmass + zqsr 
    148138 
    149139            ! Add the residual from heat diffusion equation (W.m-2) 
    … …  
    153143            ! New qsr and qns used to compute the oceanic heat flux at the next time step 
    154144            !--------------------------------------------------- 
    155             qsr(ji,jj) = zfcm1                                       
    156             qns(ji,jj) = hfx_out(ji,jj) - zfcm1               
     145            qsr(ji,jj) = zqsr                                       
     146            qns(ji,jj) = hfx_out(ji,jj) - zqsr               
    157147 
    158148            !------------------------------------------! 
    … …  
    167157            !                     Even if i see Ice melting as a FW and SALT flux 
    168158            !         
    169             !  computing freshwater exchanges at the ice/ocean interface 
    170             IF( ln_cpl ) THEN  
    171                 zemp =   emp_tot(ji,jj)                                    &   ! net mass flux over grid cell 
    172                    &   - emp_ice(ji,jj) * ( 1._wp - pfrld(ji,jj) )         &   ! minus the mass flux intercepted by sea ice 
    173                    &   + sprecip(ji,jj) * ( pfrld(ji,jj) - pfrld(ji,jj)**rn_betas )   ! 
    174             ELSE 
    175                zemp =   emp(ji,jj)     *           pfrld(ji,jj)            &   ! evaporation over oceanic fraction 
    176                   &   - tprecip(ji,jj) * ( 1._wp - pfrld(ji,jj) )          &   ! all precipitation reach the ocean 
    177                   &   + sprecip(ji,jj) * ( 1._wp - pfrld(ji,jj)**rn_betas )       ! except solid precip intercepted by sea-ice 
    178             ENDIF 
    179  
    180159            ! mass flux from ice/ocean 
    181160            wfx_ice(ji,jj) = wfx_bog(ji,jj) + wfx_bom(ji,jj) + wfx_sum(ji,jj) + wfx_sni(ji,jj)   & 
    … …  
    184163            ! mass flux at the ocean/ice interface 
    185164            fmmflx(ji,jj) = - ( wfx_ice(ji,jj) + wfx_snw(ji,jj) ) * r1_rdtice  ! F/M mass flux save at least for biogeochemical model 
    186             emp(ji,jj)    = zemp - wfx_ice(ji,jj) - wfx_snw(ji,jj)             ! mass flux + F/M mass flux (always ice/ocean mass exchange) 
     165            emp(ji,jj)    = emp_oce(ji,jj) - wfx_ice(ji,jj) - wfx_snw(ji,jj)   ! mass flux + F/M mass flux (always ice/ocean mass exchange) 
    187166             
    188167         END DO 
    … …  
    213192      tn_ice(:,:,:) = t_su(:,:,:)           ! Ice surface temperature                       
    214193 
    215       !------------------------------------------------! 
    216       !    Snow/ice albedo (only if sent to coupler)   ! 
    217       !------------------------------------------------! 
    218       IF( ln_cpl ) THEN          ! coupled case 
    219  
    220             CALL wrk_alloc( jpi, jpj, jpl, zalb_cs, zalb_os ) 
    221  
    222             CALL albedo_ice( t_su, ht_i, ht_s, zalb_cs, zalb_os )  ! cloud-sky and overcast-sky ice albedos 
    223  
    224             alb_ice(:,:,:) = ( 1. - cldf_ice ) * zalb_cs(:,:,:) + cldf_ice * zalb_os(:,:,:) 
    225  
    226             CALL wrk_dealloc( jpi, jpj, jpl, zalb_cs, zalb_os ) 
    227  
    228       ENDIF 
     194      !------------------------------------------------------------------------! 
     195      !    Snow/ice albedo (only if sent to coupler, useless in forced mode)   ! 
     196      !------------------------------------------------------------------------! 
     197      CALL wrk_alloc( jpi, jpj, jpl, zalb_cs, zalb_os )     
     198      CALL albedo_ice( t_su, ht_i, ht_s, zalb_cs, zalb_os )  ! cloud-sky and overcast-sky ice albedos 
     199      alb_ice(:,:,:) = ( 1. - cldf_ice ) * zalb_cs(:,:,:) + cldf_ice * zalb_os(:,:,:) 
     200      CALL wrk_dealloc( jpi, jpj, jpl, zalb_cs, zalb_os ) 
    229201 
    230202      ! conservation test 

  • branches/2015/dev_r5218_CNRS17_coupling/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90

    r5352 r5357  
    4848   PRIVATE 
    4949 
    50    PUBLIC   lim_thd        ! called by limstp module 
    51    PUBLIC   lim_thd_init   ! called by sbc_lim_init 
     50   PUBLIC   lim_thd         ! called by limstp module 
     51   PUBLIC   lim_thd_init    ! called by sbc_lim_init 
    5252 
    5353   !! * Substitutions 
    … …  
    9090      REAL(wp), PARAMETER :: zch        = 0.0057_wp       ! heat transfer coefficient 
    9191      ! 
    92       REAL(wp), POINTER, DIMENSION(:,:) ::  zqsr, zqns 
    9392      !!------------------------------------------------------------------- 
    94       CALL wrk_alloc( jpi,jpj, zqsr, zqns ) 
    9593 
    9694      IF( nn_timing == 1 )  CALL timing_start('limthd') 
    … …  
    134132      ! 2) Partial computation of forcing for the thermodynamic sea ice model.      ! 
    135133      !-----------------------------------------------------------------------------! 
    136  
    137       !--- Ocean solar and non solar fluxes to be used in zqld 
    138       IF ( .NOT. ln_cpl ) THEN   ! --- forced case, fluxes to the lead are the same as over the ocean 
    139          ! 
    140          zqsr(:,:) = qsr(:,:)      ; zqns(:,:) = qns(:,:) 
    141          ! 
    142       ELSE                       ! --- coupled case, fluxes to the lead are total - intercepted 
    143          ! 
    144          zqsr(:,:) = qsr_tot(:,:)  ; zqns(:,:) = qns_tot(:,:) 
    145          ! 
    146          DO jl = 1, jpl 
    147             DO jj = 1, jpj 
    148                DO ji = 1, jpi 
    149                   zqsr(ji,jj) = zqsr(ji,jj) - qsr_ice(ji,jj,jl) * a_i_b(ji,jj,jl) 
    150                   zqns(ji,jj) = zqns(ji,jj) - qns_ice(ji,jj,jl) * a_i_b(ji,jj,jl) 
    151                END DO 
    152             END DO 
    153          END DO 
    154          ! 
    155       ENDIF 
    156  
    157134      DO jj = 1, jpj 
    158135         DO ji = 1, jpi 
    … …  
    165142            !           !  temperature and turbulent mixing (McPhee, 1992) 
    166143            ! 
    167  
    168144            ! --- Energy received in the lead, zqld is defined everywhere (J.m-2) --- ! 
    169             ! REMARK valid at least in forced mode from clem 
    170             ! precip is included in qns but not in qns_ice 
    171             IF ( ln_cpl ) THEN 
    172                zqld =  tmask(ji,jj,1) * rdt_ice *  & 
    173                   &    (   zqsr(ji,jj) * frq_m(ji,jj) + zqns(ji,jj)                     &   ! pfrld already included in coupled mode 
    174                   &    + ( pfrld(ji,jj)**rn_betas - pfrld(ji,jj) ) * sprecip(ji,jj) *   &   ! heat content of precip 
    175                   &      ( cpic * ( MIN( tatm_ice(ji,jj), rt0_snow ) - rt0 ) - lfus )   & 
    176                   &    + ( 1._wp - pfrld(ji,jj) ) * ( tprecip(ji,jj) - sprecip(ji,jj) ) * rcp * ( tatm_ice(ji,jj) - rt0 ) ) 
    177             ELSE 
    178                zqld =  tmask(ji,jj,1) * rdt_ice *  & 
    179                   &      ( pfrld(ji,jj) * ( zqsr(ji,jj) * frq_m(ji,jj) + zqns(ji,jj) )   & 
    180                   &    + ( pfrld(ji,jj)**rn_betas - pfrld(ji,jj) ) * sprecip(ji,jj) *    &  ! heat content of precip 
    181                   &      ( cpic * ( MIN( tatm_ice(ji,jj), rt0_snow ) - rt0 ) - lfus )    & 
    182                   &    + ( 1._wp - pfrld(ji,jj) ) * ( tprecip(ji,jj) - sprecip(ji,jj) ) * rcp * ( tatm_ice(ji,jj) - rt0 ) ) 
    183             ENDIF 
     145            zqld =  tmask(ji,jj,1) * rdt_ice *  & 
     146               &    ( pfrld(ji,jj) * qsr_oce(ji,jj) * frq_m(ji,jj) + pfrld(ji,jj) * qns_oce(ji,jj) + qemp_oce(ji,jj) ) 
    184147 
    185148            ! --- Energy needed to bring ocean surface layer until its freezing (<0, J.m-2) --- ! 
    … …  
    208171            ! Net heat flux on top of ice-ocean [W.m-2] 
    209172            ! ----------------------------------------- 
    210             !     heat flux at the ocean surface + precip 
    211             !   + heat flux at the ice   surface  
    212             hfx_in(ji,jj) = hfx_in(ji,jj)                                                                                         &  
    213                ! heat flux above the ocean 
    214                &    +             pfrld(ji,jj)   * ( zqns(ji,jj) + zqsr(ji,jj) )                                                  & 
    215                ! latent heat of precip (note that precip is included in qns but not in qns_ice) 
    216                &    +   ( 1._wp - pfrld(ji,jj) ) * sprecip(ji,jj) * ( cpic * ( MIN( tatm_ice(ji,jj), rt0_snow ) - rt0 ) - lfus )  & 
    217                &    +   ( 1._wp - pfrld(ji,jj) ) * ( tprecip(ji,jj) - sprecip(ji,jj) ) * rcp * ( tatm_ice(ji,jj) - rt0 )          & 
    218                ! heat flux above the ice 
    219                &    +   SUM(    a_i_b(ji,jj,:)   * ( qns_ice(ji,jj,:) + qsr_ice(ji,jj,:) ) ) 
     173            hfx_in(ji,jj) = qns_tot(ji,jj) + qsr_tot(ji,jj)  
    220174 
    221175            ! ----------------------------------------------------------------------------- 
    222             ! Net heat flux that is retroceded to the ocean or taken from the ocean [W.m-2] 
     176            ! Net heat flux on top of the ocean after ice thermo (1st step) [W.m-2] 
    223177            ! ----------------------------------------------------------------------------- 
    224178            !     First  step here              :  non solar + precip - qlead - qturb 
    225179            !     Second step in limthd_dh      :  heat remaining if total melt (zq_rema)  
    226180            !     Third  step in limsbc         :  heat from ice-ocean mass exchange (zf_mass) + solar 
    227             hfx_out(ji,jj) = hfx_out(ji,jj)                                                                                       &  
    228                ! Non solar heat flux received by the ocean 
    229                &    +        pfrld(ji,jj) * zqns(ji,jj)                                                                            & 
    230                ! latent heat of precip (note that precip is included in qns but not in qns_ice) 
    231                &    +      ( pfrld(ji,jj)**rn_betas - pfrld(ji,jj) ) * sprecip(ji,jj)       & 
    232                &         * ( cpic * ( MIN( tatm_ice(ji,jj), rt0_snow ) - rt0 ) - lfus )  & 
    233                &    +      ( 1._wp - pfrld(ji,jj) ) * ( tprecip(ji,jj) - sprecip(ji,jj) ) * rcp * ( tatm_ice(ji,jj) - rt0 )       & 
    234                ! heat flux taken from the ocean where there is open water ice formation 
    235                &    -      qlead(ji,jj) * r1_rdtice                                                                               & 
    236                ! heat flux taken from the ocean during bottom growth/melt (fhld should be 0 while bott growth) 
    237                &    -      at_i(ji,jj) * fhtur(ji,jj)                                                                             & 
    238                &    -      at_i(ji,jj) *  fhld(ji,jj) 
    239  
     181            hfx_out(ji,jj) =   pfrld(ji,jj) * qns_oce(ji,jj) + qemp_oce(ji,jj)  &  ! Non solar heat flux received by the ocean                
     182               &             - qlead(ji,jj) * r1_rdtice                         &  ! heat flux taken from the ocean where there is open water ice formation 
     183               &             - at_i(ji,jj) * fhtur(ji,jj)                       &  ! heat flux taken by turbulence 
     184               &             - at_i(ji,jj) *  fhld(ji,jj)                          ! heat flux taken during bottom growth/melt  
     185                                                                                   !    (fhld should be 0 while bott growth) 
    240186         END DO 
    241187      END DO 
    … …  
    410356      IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limthd', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 
    411357 
    412       CALL wrk_dealloc( jpi,jpj, zqsr, zqns ) 
    413  
    414358      !------------------------------------------------------------------------------| 
    415359      !  6) Transport of ice between thickness categories.                           | 
    … …  
    470414   END SUBROUTINE lim_thd  
    471415 
     416  
    472417   SUBROUTINE lim_thd_temp( kideb, kiut ) 
    473418      !!----------------------------------------------------------------------- 
    … …  
    568513         END DO 
    569514          
    570          CALL tab_2d_1d( nbpb, tatm_ice_1d(1:nbpb), tatm_ice(:,:)  , jpi, jpj, npb(1:nbpb) ) 
     515         CALL tab_2d_1d( nbpb, qprec_ice_1d(1:nbpb), qprec_ice(:,:) , jpi, jpj, npb(1:nbpb) ) 
    571516         CALL tab_2d_1d( nbpb, qsr_ice_1d (1:nbpb), qsr_ice(:,:,jl) , jpi, jpj, npb(1:nbpb) ) 
    572517         CALL tab_2d_1d( nbpb, fr1_i0_1d  (1:nbpb), fr1_i0          , jpi, jpj, npb(1:nbpb) ) 
    … …  
    574519         CALL tab_2d_1d( nbpb, qns_ice_1d (1:nbpb), qns_ice(:,:,jl) , jpi, jpj, npb(1:nbpb) ) 
    575520         CALL tab_2d_1d( nbpb, ftr_ice_1d (1:nbpb), ftr_ice(:,:,jl) , jpi, jpj, npb(1:nbpb) ) 
    576          IF( .NOT. ln_cpl ) THEN 
    577             CALL tab_2d_1d( nbpb, qla_ice_1d (1:nbpb), qla_ice(:,:,jl) , jpi, jpj, npb(1:nbpb) ) 
    578             CALL tab_2d_1d( nbpb, dqla_ice_1d(1:nbpb), dqla_ice(:,:,jl), jpi, jpj, npb(1:nbpb) ) 
    579          ENDIF 
     521         CALL tab_2d_1d( nbpb, evap_ice_1d (1:nbpb), evap_ice(:,:,jl), jpi, jpj, npb(1:nbpb) ) 
    580522         CALL tab_2d_1d( nbpb, dqns_ice_1d(1:nbpb), dqns_ice(:,:,jl), jpi, jpj, npb(1:nbpb) ) 
    581523         CALL tab_2d_1d( nbpb, t_bo_1d     (1:nbpb), t_bo            , jpi, jpj, npb(1:nbpb) ) 

  • branches/2015/dev_r5218_CNRS17_coupling/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90

    r5220 r5357  
    2929   PRIVATE 
    3030 
    31    PUBLIC   lim_thd_dh   ! called by lim_thd 
     31   PUBLIC   lim_thd_dh      ! called by lim_thd 
     32   PUBLIC   lim_thd_snwblow ! called in sbcblk/sbccpl and here 
     33 
     34   INTERFACE lim_thd_snwblow 
     35      MODULE PROCEDURE lim_thd_snwblow_1d, lim_thd_snwblow_2d 
     36   END INTERFACE 
    3237 
    3338   !!---------------------------------------------------------------------- 
    … …  
    7176      REAL(wp) ::   zfdum        
    7277      REAL(wp) ::   zfracs       ! fractionation coefficient for bottom salt entrapment 
    73       REAL(wp) ::   zcoeff       ! dummy argument for snowfall partitioning over ice and leads 
    74       REAL(wp) ::   zs_snic  ! snow-ice salinity 
     78      REAL(wp) ::   zs_snic      ! snow-ice salinity 
    7579      REAL(wp) ::   zswi1        ! switch for computation of bottom salinity 
    7680      REAL(wp) ::   zswi12       ! switch for computation of bottom salinity 
    … …  
    103107      REAL(wp), POINTER, DIMENSION(:) ::   zqh_s       ! total snow heat content (J.m-2) 
    104108      REAL(wp), POINTER, DIMENSION(:) ::   zq_s        ! total snow enthalpy     (J.m-3) 
     109      REAL(wp), POINTER, DIMENSION(:) ::   zsnw        ! distribution of snow after wind blowing 
    105110 
    106111      REAL(wp) :: zswitch_sal 
    … …  
    118123 
    119124      CALL wrk_alloc( jpij, zqprec, zq_su, zq_bo, zf_tt, zq_rema ) 
    120       CALL wrk_alloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, zqh_i, zqh_s, zq_s ) 
     125      CALL wrk_alloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, zqh_i, zqh_s, zq_s, zsnw ) 
    121126      CALL wrk_alloc( jpij, nlay_i, zdeltah, zh_i ) 
    122127      CALL wrk_alloc( jpij, nlay_i, icount ) 
    … …  
    219224 
    220225      zdeltah(:,:) = 0._wp 
     226      CALL lim_thd_snwblow( 1. - at_i_1d, zsnw ) ! snow distribution over ice after wind blowing 
    221227      DO ji = kideb, kiut 
    222228         !----------- 
    … …  
    224230         !----------- 
    225231         ! thickness change 
    226          zcoeff = ( 1._wp - ( 1._wp - at_i_1d(ji) )**rn_betas ) / at_i_1d(ji)  
    227          zdh_s_pre(ji) = zcoeff * sprecip_1d(ji) * rdt_ice * r1_rhosn 
    228          ! enthalpy of the precip (>0, J.m-3) (tatm_ice is now in K) 
    229          zqprec   (ji) = rhosn * ( cpic * ( rt0 - MIN( tatm_ice_1d(ji), rt0_snow) ) + lfus )    
     232         zdh_s_pre(ji) = zsnw(ji) * sprecip_1d(ji) * rdt_ice * r1_rhosn / at_i_1d(ji) 
     233         ! enthalpy of the precip (>0, J.m-3) 
     234         zqprec   (ji) = - qprec_ice_1d(ji)    
    230235         IF( sprecip_1d(ji) == 0._wp ) zqprec(ji) = 0._wp 
    231236         ! heat flux from snow precip (>0, W.m-2) 
    … …  
    280285      ! clem comment: ice should also sublimate 
    281286      zdeltah(:,:) = 0._wp 
    282       IF( ln_cpl ) THEN 
    283          ! coupled mode: sublimation already included in emp_ice (to do in limsbc_ice) 
    284          zdh_s_sub(:)      =  0._wp  
    285       ELSE 
    286          ! forced  mode: snow thickness change due to sublimation 
    287          DO ji = kideb, kiut 
    288             zdh_s_sub(ji)  =  MAX( - ht_s_1d(ji) , - qla_ice_1d(ji) / ( rhosn * lsub ) * rdt_ice ) 
    289             ! Heat flux by sublimation [W.m-2], < 0 
    290             !      sublimate first snow that had fallen, then pre-existing snow 
    291             zdeltah(ji,1)  = MAX( zdh_s_sub(ji), - zdh_s_pre(ji) ) 
    292             hfx_sub_1d(ji) = hfx_sub_1d(ji) + ( zdeltah(ji,1) * zqprec(ji) + ( zdh_s_sub(ji) - zdeltah(ji,1) ) * q_s_1d(ji,1)  & 
    293                &                              ) * a_i_1d(ji) * r1_rdtice 
    294             ! Mass flux by sublimation 
    295             wfx_sub_1d(ji) =  wfx_sub_1d(ji) - rhosn * a_i_1d(ji) * zdh_s_sub(ji) * r1_rdtice 
    296             ! new snow thickness 
    297             ht_s_1d(ji)    =  MAX( 0._wp , ht_s_1d(ji) + zdh_s_sub(ji) ) 
    298             ! update precipitations after sublimation and correct sublimation 
    299             zdh_s_pre(ji) = zdh_s_pre(ji) + zdeltah(ji,1) 
    300             zdh_s_sub(ji) = zdh_s_sub(ji) - zdeltah(ji,1) 
    301          END DO 
    302       ENDIF 
    303  
     287      ! coupled mode: sublimation is set to 0 (evap_ice = 0) until further notice 
     288      ! forced  mode: snow thickness change due to sublimation 
     289      DO ji = kideb, kiut 
     290         zdh_s_sub(ji)  =  MAX( - ht_s_1d(ji) , - evap_ice_1d(ji) * r1_rhosn * rdt_ice ) 
     291         ! Heat flux by sublimation [W.m-2], < 0 
     292         !      sublimate first snow that had fallen, then pre-existing snow 
     293         zdeltah(ji,1)  = MAX( zdh_s_sub(ji), - zdh_s_pre(ji) ) 
     294         hfx_sub_1d(ji) = hfx_sub_1d(ji) + ( zdeltah(ji,1) * zqprec(ji) + ( zdh_s_sub(ji) - zdeltah(ji,1) ) * q_s_1d(ji,1)  & 
     295            &                              ) * a_i_1d(ji) * r1_rdtice 
     296         ! Mass flux by sublimation 
     297         wfx_sub_1d(ji) =  wfx_sub_1d(ji) - rhosn * a_i_1d(ji) * zdh_s_sub(ji) * r1_rdtice 
     298         ! new snow thickness 
     299         ht_s_1d(ji)    =  MAX( 0._wp , ht_s_1d(ji) + zdh_s_sub(ji) ) 
     300         ! update precipitations after sublimation and correct sublimation 
     301         zdh_s_pre(ji) = zdh_s_pre(ji) + zdeltah(ji,1) 
     302         zdh_s_sub(ji) = zdh_s_sub(ji) - zdeltah(ji,1) 
     303      END DO 
     304       
    304305      ! --- Update snow diags --- ! 
    305306      DO ji = kideb, kiut 
    … …  
    689690       
    690691      CALL wrk_dealloc( jpij, zqprec, zq_su, zq_bo, zf_tt, zq_rema ) 
    691       CALL wrk_dealloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, zqh_i, zqh_s, zq_s ) 
     692      CALL wrk_dealloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, zqh_i, zqh_s, zq_s, zsnw ) 
    692693      CALL wrk_dealloc( jpij, nlay_i, zdeltah, zh_i ) 
    693694      CALL wrk_dealloc( jpij, nlay_i, icount ) 
    … …  
    695696      ! 
    696697   END SUBROUTINE lim_thd_dh 
     698 
     699 
     700   !!-------------------------------------------------------------------------- 
     701   !! INTERFACE lim_thd_snwblow 
     702   !! ** Purpose :   Compute distribution of precip over the ice 
     703   !!-------------------------------------------------------------------------- 
     704   SUBROUTINE lim_thd_snwblow_2d( pin, pout ) 
     705      REAL(wp), DIMENSION(:,:), INTENT(in)  :: pin   ! previous fraction lead ( pfrld or (1. - a_i_b) ) 
     706      REAL(wp), DIMENSION(:,:), INTENT(out) :: pout 
     707      pout = ( 1._wp - ( pin )**rn_betas ) 
     708   END SUBROUTINE lim_thd_snwblow_2d 
     709 
     710   SUBROUTINE lim_thd_snwblow_1d( pin, pout ) 
     711      REAL(wp), DIMENSION(:), INTENT(in)  :: pin 
     712      REAL(wp), DIMENSION(:), INTENT(out) :: pout 
     713      pout = ( 1._wp - ( pin )**rn_betas ) 
     714   END SUBROUTINE lim_thd_snwblow_1d 
     715 
    697716    
    698717#else 

  • branches/2015/dev_r5218_CNRS17_coupling/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90

    r5220 r5357  
    2424   USE wrk_nemo       ! work arrays 
    2525   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)   
    26    USE sbc_oce, ONLY : ln_cpl 
    2726 
    2827   IMPLICIT NONE 
    … …  
    746745      !-------------------------------------------------------------------------! 
    747746      DO ji = kideb, kiut 
    748          ! forced mode only : update of latent heat fluxes (sublimation) (always >=0, upward flux)  
    749          IF( .NOT. ln_cpl) qla_ice_1d (ji) = MAX( 0._wp, qla_ice_1d (ji) + dqla_ice_1d(ji) * ( t_su_1d(ji) - ztsub(ji) ) ) 
    750747         !                                ! surface ice conduction flux 
    751748         isnow(ji)       = 1._wp - MAX( 0._wp, SIGN( 1._wp, -ht_s_1d(ji) ) ) 

  • branches/2015/dev_r5218_CNRS17_coupling/NEMOGCM/NEMO/LIM_SRC_3/limwri.F90

    r5202 r5357  
    176176      CALL iom_put( "utau_ice"    , utau_ice            )        ! wind stress over ice along i-axis at I-point 
    177177      CALL iom_put( "vtau_ice"    , vtau_ice            )        ! wind stress over ice along j-axis at I-point 
    178       CALL iom_put( "snowpre"     , sprecip             )        ! snow precipitation  
     178      CALL iom_put( "snowpre"     , sprecip * 86400.    )        ! snow precipitation  
    179179      CALL iom_put( "micesalt"    , smt_i               )        ! mean ice salinity 
    180180 
    … …  
    232232      CALL iom_put ('hfxdif'     , hfx_dif(:,:)         )   !   
    233233      CALL iom_put ('hfxopw'     , hfx_opw(:,:)         )   !   
    234       CALL iom_put ('hfxtur'     , fhtur(:,:) * at_i(:,:) ) ! turbulent heat flux at ice base  
     234      CALL iom_put ('hfxtur'     , fhtur(:,:) * SUM(a_i_b(:,:,:), dim=3) ) ! turbulent heat flux at ice base  
    235235      CALL iom_put ('hfxdhc'     , diag_heat(:,:)       )   ! Heat content variation in snow and ice  
    236236      CALL iom_put ('hfxspr'     , hfx_spr(:,:)         )   ! Heat content of snow precip  

  • branches/2015/dev_r5218_CNRS17_coupling/NEMOGCM/NEMO/LIM_SRC_3/thd_ice.F90

    r5167 r5357  
    8989   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   fhld_1d       !: <==> the 2D  fhld 
    9090   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   dqns_ice_1d   !: <==> the 2D  dqns_ice 
    91    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   qla_ice_1d    !: <==> the 2D  qla_ice 
    92    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   dqla_ice_1d   !: <==> the 2D  dqla_ice 
    93    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   tatm_ice_1d   !: <==> the 2D  tatm_ice 
     91   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   evap_ice_1d   !: <==> the 2D  evap_ice 
     92   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   qprec_ice_1d  !: <==> the 2D  qprec_ice 
    9493   !                                                     ! to reintegrate longwave flux inside the ice thermodynamics 
    9594   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   i0            !: fraction of radiation transmitted to the ice 
    … …  
    153152         &      fhld_1d    (jpij) , wfx_sub_1d (jpij) , wfx_bog_1d (jpij) , wfx_bom_1d(jpij) ,  & 
    154153         &      wfx_sum_1d(jpij)  , wfx_sni_1d (jpij) , wfx_opw_1d (jpij) , wfx_res_1d(jpij) ,  & 
    155          &      dqns_ice_1d(jpij) , qla_ice_1d (jpij) , dqla_ice_1d(jpij) ,                     & 
    156          &      tatm_ice_1d(jpij) , i0         (jpij) ,                                         &   
     154         &      dqns_ice_1d(jpij) , evap_ice_1d (jpij),                                         & 
     155         &      qprec_ice_1d(jpij), i0         (jpij) ,                                         &   
    157156         &      sfx_bri_1d (jpij) , sfx_bog_1d (jpij) , sfx_bom_1d (jpij) , sfx_sum_1d (jpij),  & 
    158157         &      sfx_sni_1d (jpij) , sfx_opw_1d (jpij) , sfx_res_1d (jpij) ,                     & 

  • branches/2015/dev_r5218_CNRS17_coupling/NEMOGCM/NEMO/OPA_SRC/IOM/restart.F90

    r5352 r5357  
    213213         CALL iom_get( numror, jpdom_autoglo, 'sshb'   , sshb    ) 
    214214! EM Attention Ceci doit etre reimplemente correctement 
    215 !EM         IF( lk_lim3 .OR. ( nn_components == jp_iam_opa ) )   CALL iom_get( numror, jpdom_autoglo, 'fse3t_b', fse3t_b(:,:,:) ) 
    216          CALL iom_get( numror, jpdom_autoglo, 'fse3t_b', fse3t_b(:,:,:) ) 
     215         IF( lk_lim3 .OR. ( nn_components == jp_iam_opa ) )   CALL iom_get( numror, jpdom_autoglo, 'fse3t_b', fse3t_b(:,:,:) ) 
     216!clem         CALL iom_get( numror, jpdom_autoglo, 'fse3t_b', fse3t_b(:,:,:) ) 
    217217      ELSE 
    218218         neuler = 0 

  • branches/2015/dev_r5218_CNRS17_coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_ice.F90

    r5220 r5357  
    6969   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   fr1_i0         !: Solar surface transmission parameter, thick ice  [-] 
    7070   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   fr2_i0         !: Solar surface transmission parameter, thin ice   [-] 
    71    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   emp_ice        !: sublimation - precip over sea ice            [kg/m2] 
     71   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   emp_ice        !: sublimation - precip over sea ice          [kg/m2/s] 
    7272 
    7373   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   topmelt            !: category topmelt 
    7474   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   botmelt            !: category botmelt 
     75 
     76#if defined  key_lim3 
     77   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   evap_ice       !: sublimation                              [kg/m2/s] 
     78   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   devap_ice      !: sublimation sensitivity                [kg/m2/s/K] 
     79   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   qns_oce        !: non solar heat flux over ocean              [W/m2] 
     80   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   qsr_oce        !: non solar heat flux over ocean              [W/m2] 
     81   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   qemp_oce       !: heat flux of precip and evap over ocean     [W/m2] 
     82   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   qemp_ice       !: heat flux of precip and evap over ice       [W/m2] 
     83   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   qprec_ice      !: heat flux of precip over ice                [J/m3] 
     84   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   emp_oce        !: evap - precip over ocean                 [kg/m2/s] 
     85#endif 
     86#if defined key_lim3 || defined key_lim2 
     87   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   wndm_ice       !: wind speed module at T-point                 [m/s] 
     88#endif 
    7589 
    7690#if defined key_cice 
    … …  
    100114#endif 
    101115 
    102 #if defined key_lim3 || defined key_cice 
    103    ! not used with LIM2 
     116#if defined key_cice 
    104117   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   tatm_ice       !: air temperature [K] 
    105118#endif 
    … …  
    125138      ALLOCATE( qns_ice (jpi,jpj,jpl) , qsr_ice (jpi,jpj,jpl) ,     & 
    126139         &      qla_ice (jpi,jpj,jpl) , dqla_ice(jpi,jpj,jpl) ,     & 
    127          &      dqns_ice(jpi,jpj,jpl) , tn_ice  (jpi,jpj,jpl) ,     & 
    128          &      alb_ice (jpi,jpj,jpl) ,                             & 
    129          &      utau_ice(jpi,jpj)     , vtau_ice(jpi,jpj)     ,     & 
     140         &      dqns_ice(jpi,jpj,jpl) , tn_ice  (jpi,jpj,jpl) , alb_ice (jpi,jpj,jpl) ,   & 
     141         &      utau_ice(jpi,jpj)     , vtau_ice(jpi,jpj)     , wndm_ice(jpi,jpj)     ,   & 
    130142         &      fr1_i0  (jpi,jpj)     , fr2_i0  (jpi,jpj)     ,     & 
    131 #if defined key_lim3 
    132          &      tatm_ice(jpi,jpj)     ,                             & 
    133 #endif 
    134143#if defined key_lim2 
    135144         &      a_i(jpi,jpj,jpl)      ,                             & 
     145#endif 
     146#if defined key_lim3 
     147         &      evap_ice(jpi,jpj,jpl) , devap_ice(jpi,jpj,jpl) , qprec_ice(jpi,jpj) ,  & 
     148         &      qemp_ice(jpi,jpj)     , qemp_oce(jpi,jpj)      ,                       & 
     149         &      qns_oce (jpi,jpj)     , qsr_oce (jpi,jpj)      , emp_oce (jpi,jpj)  ,  & 
    136150#endif 
    137151         &      emp_ice(jpi,jpj)      ,  STAT= ierr(1) ) 

  • branches/2015/dev_r5218_CNRS17_coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_clio.F90

    r5126 r5357  
    4646   PUBLIC sbc_blk_clio        ! routine called by sbcmod.F90  
    4747   PUBLIC blk_ice_clio        ! routine called by sbcice_lim.F90  
     48   PUBLIC blk_ice_clio_tau    ! routine called by sbcice_lim.F90  
     49   PUBLIC blk_ice_clio_flx    ! routine called by sbcice_lim.F90  
    4850 
    4951   INTEGER , PARAMETER ::   jpfld   = 7           ! maximum number of files to read  
    … …  
    399401   END SUBROUTINE blk_oce_clio 
    400402 
    401  
    402403   SUBROUTINE blk_ice_clio(  pst   , palb_cs, palb_os, palb,  & 
    403404      &                      p_taui, p_tauj, p_qns , p_qsr,   & 
    … …  
    405406      &                      p_tpr , p_spr ,                  & 
    406407      &                      p_fr1 , p_fr2 , cd_grid, pdim  ) 
     408 
    407409      !!--------------------------------------------------------------------------- 
    408410      !!                     ***  ROUTINE blk_ice_clio  *** 
    409       !!                  
     411      !! 
    410412      !!  ** Purpose :   Computation of the heat fluxes at ocean and snow/ice 
    411413      !!       surface the solar heat at ocean and snow/ice surfaces and the  
    412414      !!       sensitivity of total heat fluxes to the SST variations 
    413415      !!          
    414       !!  ** Method  :   The flux of heat at the ice and ocean surfaces are derived 
    415       !!       from semi-empirical ( or bulk ) formulae which relate the flux to  
    416       !!       the properties of the surface and of the lower atmosphere. Here, we 
    417       !!       follow the work of Oberhuber, 1988    
    418       !! 
    419       !!  ** Action  :   call albedo_oce/albedo_ice to compute ocean/ice albedo  
    420       !!               - snow precipitation 
    421       !!               - solar flux at the ocean and ice surfaces 
    422       !!               - the long-wave radiation for the ocean and sea/ice 
    423       !!               - turbulent heat fluxes over water and ice 
    424       !!               - evaporation over water 
    425       !!               - total heat fluxes sensitivity over ice (dQ/dT) 
    426       !!               - latent heat flux sensitivity over ice (dQla/dT) 
    427       !!               - qns  :  modified the non solar heat flux over the ocean 
    428       !!                         to take into account solid precip latent heat flux 
     416      !!  ** Action  :   Call of blk_ice_clio_tau and blk_ice_clio_flx 
     417      !! 
    429418      !!---------------------------------------------------------------------- 
     419 
    430420      REAL(wp), INTENT(in   ), DIMENSION(:,:,:)   ::   pst      ! ice surface temperature                   [Kelvin] 
    431421      REAL(wp), INTENT(in   ), DIMENSION(:,:,:)   ::   palb_cs  ! ice albedo (clear    sky) (alb_ice_cs)         [-] 
    … …  
    445435      CHARACTER(len=1), INTENT(in   )             ::   cd_grid  ! type of sea-ice grid ("C" or "B" grid) 
    446436      INTEGER, INTENT(in   )                      ::   pdim     ! number of ice categories 
     437 
     438      CALL blk_ice_clio_tau( p_taui, p_tauj, cd_grid ) 
     439      CALL blk_ice_clio_flx(  pst   , palb_cs, palb_os, palb,        & 
     440         &                    p_qns , p_qsr, p_qla , p_dqns, p_dqla, & 
     441         &                    p_tpr , p_spr ,p_fr1 , p_fr2 , pdim  ) 
     442 
     443   END SUBROUTINE blk_ice_clio 
     444 
     445   SUBROUTINE blk_ice_clio_tau( p_taui, p_tauj, cd_grid ) 
     446      !!--------------------------------------------------------------------------- 
     447      !!                     ***  ROUTINE blk_ice_clio_tau  *** 
     448      !!                  
     449      !!  ** Purpose :   Computation momentum flux at the ice-atm interface   
     450      !!          
     451      !!  ** Method  :   Read utau from a forcing file. Rearrange if C-grid 
     452      !! 
     453      !!---------------------------------------------------------------------- 
     454      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   p_taui   ! surface ice stress at I-point (i-component) [N/m2] 
     455      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   p_tauj   ! surface ice stress at I-point (j-component) [N/m2] 
     456      CHARACTER(len=1), INTENT(in   )             ::   cd_grid  ! type of sea-ice grid ("C" or "B" grid) 
     457      !! 
     458      INTEGER  ::   ji, jj    ! dummy loop indices 
     459      REAL(wp) ::   zcoef 
     460      !! 
     461      !!--------------------------------------------------------------------- 
     462      ! 
     463 
     464      IF( nn_timing == 1 )  CALL timing_start('blk_ice_clio_tau') 
     465 
     466      SELECT CASE( cd_grid ) 
     467 
     468      CASE( 'C' )                          ! C-grid ice dynamics 
     469 
     470         zcoef  = cai / cao                         ! Change from air-sea stress to air-ice stress 
     471         p_taui(:,:) = zcoef * utau(:,:) 
     472         p_tauj(:,:) = zcoef * vtau(:,:) 
     473 
     474      CASE( 'I' )                          ! I-grid ice dynamics:  I-point (i.e. F-point lower-left corner) 
     475 
     476         zcoef  = 0.5_wp * cai / cao                ! Change from air-sea stress to air-ice stress 
     477         DO jj = 2, jpj         ! stress from ocean U- and V-points to ice U,V point 
     478            DO ji = 2, jpi   ! I-grid : no vector opt. 
     479               p_taui(ji,jj) = zcoef * ( utau(ji-1,jj  ) + utau(ji-1,jj-1) ) 
     480               p_tauj(ji,jj) = zcoef * ( vtau(ji  ,jj-1) + vtau(ji-1,jj-1) ) 
     481            END DO 
     482         END DO 
     483 
     484         CALL lbc_lnk( p_taui(:,:), 'I', -1. )   ;   CALL lbc_lnk( p_tauj(:,:), 'I', -1. )   ! I-point 
     485 
     486      END SELECT 
     487 
     488      IF(ln_ctl) THEN 
     489         CALL prt_ctl(tab2d_1=p_taui , clinfo1=' blk_ice_clio: p_taui : ', tab2d_2=p_tauj , clinfo2=' p_tauj : ') 
     490      ENDIF 
     491 
     492      IF( nn_timing == 1 )  CALL timing_stop('blk_ice_clio_tau') 
     493 
     494   END SUBROUTINE blk_ice_clio_tau 
     495 
     496   SUBROUTINE blk_ice_clio_flx(  pst   , palb_cs, palb_os, palb,  & 
     497      &                          p_qns , p_qsr, p_qla , p_dqns, p_dqla, & 
     498      &                          p_tpr , p_spr ,p_fr1 , p_fr2 , pdim  ) 
     499      !!--------------------------------------------------------------------------- 
     500      !!                     ***  ROUTINE blk_ice_clio_flx *** 
     501      !!                  
     502      !!  ** Purpose :   Computation of the heat fluxes at ocean and snow/ice 
     503      !!       surface the solar heat at ocean and snow/ice surfaces and the  
     504      !!       sensitivity of total heat fluxes to the SST variations 
     505      !!          
     506      !!  ** Method  :   The flux of heat at the ice and ocean surfaces are derived 
     507      !!       from semi-empirical ( or bulk ) formulae which relate the flux to  
     508      !!       the properties of the surface and of the lower atmosphere. Here, we 
     509      !!       follow the work of Oberhuber, 1988    
     510      !! 
     511      !!  ** Action  :   call albedo_oce/albedo_ice to compute ocean/ice albedo  
     512      !!               - snow precipitation 
     513      !!               - solar flux at the ocean and ice surfaces 
     514      !!               - the long-wave radiation for the ocean and sea/ice 
     515      !!               - turbulent heat fluxes over water and ice 
     516      !!               - evaporation over water 
     517      !!               - total heat fluxes sensitivity over ice (dQ/dT) 
     518      !!               - latent heat flux sensitivity over ice (dQla/dT) 
     519      !!               - qns  :  modified the non solar heat flux over the ocean 
     520      !!                         to take into account solid precip latent heat flux 
     521      !!---------------------------------------------------------------------- 
     522      REAL(wp), INTENT(in   ), DIMENSION(:,:,:)   ::   pst      ! ice surface temperature                   [Kelvin] 
     523      REAL(wp), INTENT(in   ), DIMENSION(:,:,:)   ::   palb_cs  ! ice albedo (clear    sky) (alb_ice_cs)         [-] 
     524      REAL(wp), INTENT(in   ), DIMENSION(:,:,:)   ::   palb_os  ! ice albedo (overcast sky) (alb_ice_os)         [-] 
     525      REAL(wp), INTENT(  out), DIMENSION(:,:,:)   ::   palb     ! ice albedo (actual value)                      [-] 
     526      REAL(wp), INTENT(  out), DIMENSION(:,:,:)   ::   p_qns    ! non solar heat flux over ice (T-point)      [W/m2] 
     527      REAL(wp), INTENT(  out), DIMENSION(:,:,:)   ::   p_qsr    !     solar heat flux over ice (T-point)      [W/m2] 
     528      REAL(wp), INTENT(  out), DIMENSION(:,:,:)   ::   p_qla    ! latent    heat flux over ice (T-point)      [W/m2] 
     529      REAL(wp), INTENT(  out), DIMENSION(:,:,:)   ::   p_dqns   ! non solar heat sensistivity  (T-point)      [W/m2] 
     530      REAL(wp), INTENT(  out), DIMENSION(:,:,:)   ::   p_dqla   ! latent    heat sensistivity  (T-point)      [W/m2] 
     531      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   p_tpr    ! total precipitation          (T-point)   [Kg/m2/s] 
     532      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   p_spr    ! solid precipitation          (T-point)   [Kg/m2/s] 
     533      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   p_fr1    ! 1sr fraction of qsr penetration in ice         [-] 
     534      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   p_fr2    ! 2nd fraction of qsr penetration in ice         [-] 
     535      INTEGER, INTENT(in   )                      ::   pdim     ! number of ice categories 
    447536      !! 
    448537      INTEGER  ::   ji, jj, jl    ! dummy loop indices 
    449538      INTEGER  ::   ijpl          ! number of ice categories (size of 3rd dim of input arrays) 
    450539      !! 
    451       REAL(wp) ::   zcoef, zmt1, zmt2, zmt3, ztatm3     ! temporary scalars 
     540      REAL(wp) ::   zmt1, zmt2, zmt3, ztatm3                    ! temporary scalars 
    452541      REAL(wp) ::   ztaevbk, zind1, zind2, zind3, ztamr         !    -         - 
    453542      REAL(wp) ::   zesi, zqsati, zdesidt                       !    -         - 
    … …  
    463552      !!--------------------------------------------------------------------- 
    464553      ! 
    465       IF( nn_timing == 1 )  CALL timing_start('blk_ice_clio') 
     554      IF( nn_timing == 1 )  CALL timing_start('blk_ice_clio_flx') 
    466555      ! 
    467556      CALL wrk_alloc( jpi,jpj, ztatm, zqatm, zevsqr, zrhoa ) 
    … …  
    471560      zpatm = 101000.                        ! atmospheric pressure  (assumed constant  here) 
    472561 
    473 #if defined key_lim3       
    474       tatm_ice(:,:) = sf(jp_tair)%fnow(:,:,1)   ! LIM3: make Tair available in sea-ice. WARNING allocated after call to ice_init 
    475 #endif 
    476       !                                                        ! surface ocean fluxes computed with CLIO bulk formulea 
    477       !------------------------------------! 
    478       !   momentum fluxes  (utau, vtau )   ! 
    479       !------------------------------------! 
    480  
    481       SELECT CASE( cd_grid ) 
    482       CASE( 'C' )                          ! C-grid ice dynamics 
    483          zcoef  = cai / cao                         ! Change from air-sea stress to air-ice stress 
    484          p_taui(:,:) = zcoef * utau(:,:) 
    485          p_tauj(:,:) = zcoef * vtau(:,:) 
    486       CASE( 'I' )                          ! I-grid ice dynamics:  I-point (i.e. F-point lower-left corner) 
    487          zcoef  = 0.5_wp * cai / cao                ! Change from air-sea stress to air-ice stress 
    488          DO jj = 2, jpj         ! stress from ocean U- and V-points to ice U,V point 
    489             DO ji = 2, jpi   ! I-grid : no vector opt. 
    490                p_taui(ji,jj) = zcoef * ( utau(ji-1,jj  ) + utau(ji-1,jj-1) ) 
    491                p_tauj(ji,jj) = zcoef * ( vtau(ji  ,jj-1) + vtau(ji-1,jj-1) ) 
    492             END DO 
    493          END DO 
    494          CALL lbc_lnk( p_taui(:,:), 'I', -1. )   ;   CALL lbc_lnk( p_tauj(:,:), 'I', -1. )   ! I-point 
    495       END SELECT 
    496  
    497  
     562      !-------------------------------------------------------------------------------- 
    498563      !  Determine cloud optical depths as a function of latitude (Chou et al., 1981). 
    499564      !  and the correction factor for taking into account  the effect of clouds  
    500       !------------------------------------------------------ 
     565      !-------------------------------------------------------------------------------- 
     566 
    501567!CDIR NOVERRCHK 
    502568!CDIR COLLAPSE 
    … …  
    653719         CALL prt_ctl(tab3d_1=p_dqla , clinfo1=' blk_ice_clio: p_dqla : ', tab3d_2=pst    , clinfo2=' pst    : ', kdim=ijpl) 
    654720         CALL prt_ctl(tab2d_1=p_tpr  , clinfo1=' blk_ice_clio: p_tpr  : ', tab2d_2=p_spr  , clinfo2=' p_spr  : ') 
    655          CALL prt_ctl(tab2d_1=p_taui , clinfo1=' blk_ice_clio: p_taui : ', tab2d_2=p_tauj , clinfo2=' p_tauj : ') 
    656721      ENDIF 
    657722 
    … …  
    659724      CALL wrk_dealloc( jpi,jpj,pdim, z_qlw, z_qsb ) 
    660725      ! 
    661       IF( nn_timing == 1 )  CALL timing_stop('blk_ice_clio') 
    662       ! 
    663    END SUBROUTINE blk_ice_clio 
    664  
     726      IF( nn_timing == 1 )  CALL timing_stop('blk_ice_clio_flx') 
     727      ! 
     728   END SUBROUTINE blk_ice_clio_flx 
    665729 
    666730   SUBROUTINE blk_clio_qsr_oce( pqsr_oce ) 

  • branches/2015/dev_r5218_CNRS17_coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90

    r5065 r5357  
    4646   USE sbc_ice         ! Surface boundary condition: ice fields 
    4747   USE lib_fortran     ! to use key_nosignedzero 
     48#if defined key_lim3 
     49   USE ice, ONLY       : u_ice, v_ice, jpl, pfrld, a_i_b 
     50   USE limthd_dh       ! for CALL lim_thd_snwblow 
     51#elif defined key_lim2 
     52   USE ice_2, ONLY     : u_ice, v_ice 
     53   USE par_ice_2 
     54#endif 
    4855 
    4956   IMPLICIT NONE 
    … …  
    5158 
    5259   PUBLIC   sbc_blk_core         ! routine called in sbcmod module 
    53    PUBLIC   blk_ice_core         ! routine called in sbc_ice_lim module 
     60#if defined key_lim2 || defined key_lim3 
     61   PUBLIC   blk_ice_core_tau     ! routine called in sbc_ice_lim module 
     62   PUBLIC   blk_ice_core_flx     ! routine called in sbc_ice_lim module 
     63#endif 
    5464   PUBLIC   blk_ice_meanqsr      ! routine called in sbc_ice_lim module 
    5565   PUBLIC   turb_core_2z         ! routine calles in sbcblk_mfs module 
    … …  
    376386      emp (:,:) = (  zevap(:,:)                                          &   ! mass flux (evap. - precip.) 
    377387         &         - sf(jp_prec)%fnow(:,:,1) * rn_pfac  ) * tmask(:,:,1) 
    378       qns(:,:) = zqlw(:,:) - zqsb(:,:) - zqla(:,:)                                &   ! Downward Non Solar flux 
     388      ! 
     389      qns(:,:) = zqlw(:,:) - zqsb(:,:) - zqla(:,:)                                &   ! Downward Non Solar  
    379390         &     - sf(jp_snow)%fnow(:,:,1) * rn_pfac * lfus                         &   ! remove latent melting heat for solid precip 
    380391         &     - zevap(:,:) * pst(:,:) * rcp                                      &   ! remove evap heat content at SST 
    … …  
    383394         &     + sf(jp_snow)%fnow(:,:,1) * rn_pfac                                &   ! add solid  precip heat content at min(Tair,Tsnow) 
    384395         &     * ( MIN( sf(jp_tair)%fnow(:,:,1), rt0_snow ) - rt0 ) * cpic * tmask(:,:,1) 
     396      ! 
     397#if defined key_lim3 
     398      qns_oce(:,:) = zqlw(:,:) - zqsb(:,:) - zqla(:,:)                                ! non solar without emp (only needed by LIM3) 
     399      qsr_oce(:,:) = qsr(:,:) 
     400#endif 
    385401      ! 
    386402      CALL iom_put( "qlw_oce",   zqlw )                 ! output downward longwave heat over the ocean 
    … …  
    406422  
    407423    
    408    SUBROUTINE blk_ice_core(  pst   , pui   , pvi   , palb ,   & 
    409       &                      p_taui, p_tauj, p_qns , p_qsr,   & 
    410       &                      p_qla , p_dqns, p_dqla,          & 
    411       &                      p_tpr , p_spr ,                  & 
    412       &                      p_fr1 , p_fr2 , cd_grid, pdim  )  
    413       !!--------------------------------------------------------------------- 
    414       !!                     ***  ROUTINE blk_ice_core  *** 
     424#if defined key_lim2 || defined key_lim3 
     425   SUBROUTINE blk_ice_core_tau 
     426      !!--------------------------------------------------------------------- 
     427      !!                     ***  ROUTINE blk_ice_core_tau  *** 
    415428      !! 
    416429      !! ** Purpose :   provide the surface boundary condition over sea-ice 
    417430      !! 
    418       !! ** Method  :   compute momentum, heat and freshwater exchanged 
    419       !!                between atmosphere and sea-ice using CORE bulk 
    420       !!                formulea, ice variables and read atmmospheric fields. 
     431      !! ** Method  :   compute momentum using CORE bulk 
     432      !!                formulea, ice variables and read atmospheric fields. 
    421433      !!                NB: ice drag coefficient is assumed to be a constant 
    422       !!  
    423       !! caution : the net upward water flux has with mm/day unit 
    424       !!--------------------------------------------------------------------- 
    425       REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   pst      ! ice surface temperature (>0, =rt0 over land) [Kelvin] 
    426       REAL(wp), DIMENSION(:,:)  , INTENT(in   ) ::   pui      ! ice surface velocity (i- and i- components      [m/s] 
    427       REAL(wp), DIMENSION(:,:)  , INTENT(in   ) ::   pvi      !    at I-point (B-grid) or U & V-point (C-grid) 
    428       REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   palb     ! ice albedo (all skies)                            [%] 
    429       REAL(wp), DIMENSION(:,:)  , INTENT(  out) ::   p_taui   ! i- & j-components of surface ice stress        [N/m2] 
    430       REAL(wp), DIMENSION(:,:)  , INTENT(  out) ::   p_tauj   !   at I-point (B-grid) or U & V-point (C-grid) 
    431       REAL(wp), DIMENSION(:,:,:), INTENT(  out) ::   p_qns    ! non solar heat flux over ice (T-point)         [W/m2] 
    432       REAL(wp), DIMENSION(:,:,:), INTENT(  out) ::   p_qsr    !     solar heat flux over ice (T-point)         [W/m2] 
    433       REAL(wp), DIMENSION(:,:,:), INTENT(  out) ::   p_qla    ! latent    heat flux over ice (T-point)         [W/m2] 
    434       REAL(wp), DIMENSION(:,:,:), INTENT(  out) ::   p_dqns   ! non solar heat sensistivity  (T-point)         [W/m2] 
    435       REAL(wp), DIMENSION(:,:,:), INTENT(  out) ::   p_dqla   ! latent    heat sensistivity  (T-point)         [W/m2] 
    436       REAL(wp), DIMENSION(:,:)  , INTENT(  out) ::   p_tpr    ! total precipitation          (T-point)      [Kg/m2/s] 
    437       REAL(wp), DIMENSION(:,:)  , INTENT(  out) ::   p_spr    ! solid precipitation          (T-point)      [Kg/m2/s] 
    438       REAL(wp), DIMENSION(:,:)  , INTENT(  out) ::   p_fr1    ! 1sr fraction of qsr penetration in ice (T-point)  [%] 
    439       REAL(wp), DIMENSION(:,:)  , INTENT(  out) ::   p_fr2    ! 2nd fraction of qsr penetration in ice (T-point)  [%] 
    440       CHARACTER(len=1)          , INTENT(in   ) ::   cd_grid  ! ice grid ( C or B-grid) 
    441       INTEGER                   , INTENT(in   ) ::   pdim     ! number of ice categories 
    442       !! 
    443       INTEGER  ::   ji, jj, jl    ! dummy loop indices 
    444       INTEGER  ::   ijpl          ! number of ice categories (size of 3rd dim of input arrays) 
    445       REAL(wp) ::   zst2, zst3 
    446       REAL(wp) ::   zcoef_wnorm, zcoef_wnorm2, zcoef_dqlw, zcoef_dqla, zcoef_dqsb 
    447       REAL(wp) ::   zztmp                                        ! temporary variable 
    448       REAL(wp) ::   zwnorm_f, zwndi_f , zwndj_f                  ! relative wind module and components at F-point 
    449       REAL(wp) ::             zwndi_t , zwndj_t                  ! relative wind components at T-point 
    450       !! 
    451       REAL(wp), DIMENSION(:,:)  , POINTER ::   z_wnds_t          ! wind speed ( = | U10m - U_ice | ) at T-point 
    452       REAL(wp), DIMENSION(:,:,:), POINTER ::   z_qlw             ! long wave heat flux over ice 
    453       REAL(wp), DIMENSION(:,:,:), POINTER ::   z_qsb             ! sensible  heat flux over ice 
    454       REAL(wp), DIMENSION(:,:,:), POINTER ::   z_dqlw            ! long wave heat sensitivity over ice 
    455       REAL(wp), DIMENSION(:,:,:), POINTER ::   z_dqsb            ! sensible  heat sensitivity over ice 
    456       !!--------------------------------------------------------------------- 
    457       ! 
    458       IF( nn_timing == 1 )  CALL timing_start('blk_ice_core') 
    459       ! 
    460       CALL wrk_alloc( jpi,jpj, z_wnds_t ) 
    461       CALL wrk_alloc( jpi,jpj,pdim, z_qlw, z_qsb, z_dqlw, z_dqsb )  
    462  
    463       ijpl  = pdim                            ! number of ice categories 
    464  
     434      !!--------------------------------------------------------------------- 
     435      INTEGER  ::   ji, jj    ! dummy loop indices 
     436      REAL(wp) ::   zcoef_wnorm, zcoef_wnorm2 
     437      REAL(wp) ::   zwnorm_f, zwndi_f , zwndj_f               ! relative wind module and components at F-point 
     438      REAL(wp) ::             zwndi_t , zwndj_t               ! relative wind components at T-point 
     439      !!--------------------------------------------------------------------- 
     440      ! 
     441      IF( nn_timing == 1 )  CALL timing_start('blk_ice_core_tau') 
     442      ! 
    465443      ! local scalars ( place there for vector optimisation purposes) 
    466444      zcoef_wnorm  = rhoa * Cice 
    467445      zcoef_wnorm2 = rhoa * Cice * 0.5 
    468       zcoef_dqlw   = 4.0 * 0.95 * Stef 
    469       zcoef_dqla   = -Ls * Cice * 11637800. * (-5897.8) 
    470       zcoef_dqsb   = rhoa * cpa * Cice 
    471446 
    472447!!gm brutal.... 
    473       z_wnds_t(:,:) = 0.e0 
    474       p_taui  (:,:) = 0.e0 
    475       p_tauj  (:,:) = 0.e0 
     448      utau_ice  (:,:) = 0._wp 
     449      vtau_ice  (:,:) = 0._wp 
     450      wndm_ice  (:,:) = 0._wp 
    476451!!gm end 
    477452 
    478 #if defined key_lim3 
    479       tatm_ice(:,:) = sf(jp_tair)%fnow(:,:,1)   ! LIM3: make Tair available in sea-ice. WARNING allocated after call to ice_init 
    480 #endif 
    481453      ! ----------------------------------------------------------------------------- ! 
    482454      !    Wind components and module relative to the moving ocean ( U10m - U_ice )   ! 
    483455      ! ----------------------------------------------------------------------------- ! 
    484       SELECT CASE( cd_grid ) 
     456      SELECT CASE( cp_ice_msh ) 
    485457      CASE( 'I' )                  ! B-grid ice dynamics :   I-point (i.e. F-point with sea-ice indexation) 
    486458         !                           and scalar wind at T-point ( = | U10m - U_ice | ) (masked) 
    … …  
    489461               ! ... scalar wind at I-point (fld being at T-point) 
    490462               zwndi_f = 0.25 * (  sf(jp_wndi)%fnow(ji-1,jj  ,1) + sf(jp_wndi)%fnow(ji  ,jj  ,1)   & 
    491                   &              + sf(jp_wndi)%fnow(ji-1,jj-1,1) + sf(jp_wndi)%fnow(ji  ,jj-1,1)  ) - rn_vfac * pui(ji,jj) 
     463                  &              + sf(jp_wndi)%fnow(ji-1,jj-1,1) + sf(jp_wndi)%fnow(ji  ,jj-1,1)  ) - rn_vfac * u_ice(ji,jj) 
    492464               zwndj_f = 0.25 * (  sf(jp_wndj)%fnow(ji-1,jj  ,1) + sf(jp_wndj)%fnow(ji  ,jj  ,1)   & 
    493                   &              + sf(jp_wndj)%fnow(ji-1,jj-1,1) + sf(jp_wndj)%fnow(ji  ,jj-1,1)  ) - rn_vfac * pvi(ji,jj) 
     465                  &              + sf(jp_wndj)%fnow(ji-1,jj-1,1) + sf(jp_wndj)%fnow(ji  ,jj-1,1)  ) - rn_vfac * v_ice(ji,jj) 
    494466               zwnorm_f = zcoef_wnorm * SQRT( zwndi_f * zwndi_f + zwndj_f * zwndj_f ) 
    495467               ! ... ice stress at I-point 
    496                p_taui(ji,jj) = zwnorm_f * zwndi_f 
    497                p_tauj(ji,jj) = zwnorm_f * zwndj_f 
     468               utau_ice(ji,jj) = zwnorm_f * zwndi_f 
     469               vtau_ice(ji,jj) = zwnorm_f * zwndj_f 
    498470               ! ... scalar wind at T-point (fld being at T-point) 
    499                zwndi_t = sf(jp_wndi)%fnow(ji,jj,1) - rn_vfac * 0.25 * (  pui(ji,jj+1) + pui(ji+1,jj+1)   & 
    500                   &                                                    + pui(ji,jj  ) + pui(ji+1,jj  )  ) 
    501                zwndj_t = sf(jp_wndj)%fnow(ji,jj,1) - rn_vfac * 0.25 * (  pvi(ji,jj+1) + pvi(ji+1,jj+1)   & 
    502                   &                                                    + pvi(ji,jj  ) + pvi(ji+1,jj  )  ) 
    503                z_wnds_t(ji,jj)  = SQRT( zwndi_t * zwndi_t + zwndj_t * zwndj_t ) * tmask(ji,jj,1) 
     471               zwndi_t = sf(jp_wndi)%fnow(ji,jj,1) - rn_vfac * 0.25 * (  u_ice(ji,jj+1) + u_ice(ji+1,jj+1)   & 
     472                  &                                                    + u_ice(ji,jj  ) + u_ice(ji+1,jj  )  ) 
     473               zwndj_t = sf(jp_wndj)%fnow(ji,jj,1) - rn_vfac * 0.25 * (  v_ice(ji,jj+1) + v_ice(ji+1,jj+1)   & 
     474                  &                                                    + v_ice(ji,jj  ) + v_ice(ji+1,jj  )  ) 
     475               wndm_ice(ji,jj)  = SQRT( zwndi_t * zwndi_t + zwndj_t * zwndj_t ) * tmask(ji,jj,1) 
    504476            END DO 
    505477         END DO 
    506          CALL lbc_lnk( p_taui  , 'I', -1. ) 
    507          CALL lbc_lnk( p_tauj  , 'I', -1. ) 
    508          CALL lbc_lnk( z_wnds_t, 'T',  1. ) 
     478         CALL lbc_lnk( utau_ice, 'I', -1. ) 
     479         CALL lbc_lnk( vtau_ice, 'I', -1. ) 
     480         CALL lbc_lnk( wndm_ice, 'T',  1. ) 
    509481         ! 
    510482      CASE( 'C' )                  ! C-grid ice dynamics :   U & V-points (same as ocean) 
    511483         DO jj = 2, jpj 
    512484            DO ji = fs_2, jpi   ! vect. opt. 
    513                zwndi_t = (  sf(jp_wndi)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( pui(ji-1,jj  ) + pui(ji,jj) )  ) 
    514                zwndj_t = (  sf(jp_wndj)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( pvi(ji  ,jj-1) + pvi(ji,jj) )  ) 
    515                z_wnds_t(ji,jj) = SQRT( zwndi_t * zwndi_t + zwndj_t * zwndj_t ) * tmask(ji,jj,1) 
     485               zwndi_t = (  sf(jp_wndi)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( u_ice(ji-1,jj  ) + u_ice(ji,jj) )  ) 
     486               zwndj_t = (  sf(jp_wndj)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( v_ice(ji  ,jj-1) + v_ice(ji,jj) )  ) 
     487               wndm_ice(ji,jj) = SQRT( zwndi_t * zwndi_t + zwndj_t * zwndj_t ) * tmask(ji,jj,1) 
    516488            END DO 
    517489         END DO 
    518490         DO jj = 2, jpjm1 
    519491            DO ji = fs_2, fs_jpim1   ! vect. opt. 
    520                p_taui(ji,jj) = zcoef_wnorm2 * ( z_wnds_t(ji+1,jj  ) + z_wnds_t(ji,jj) )                          & 
    521                   &          * ( 0.5 * (sf(jp_wndi)%fnow(ji+1,jj,1) + sf(jp_wndi)%fnow(ji,jj,1) ) - rn_vfac * pui(ji,jj) ) 
    522                p_tauj(ji,jj) = zcoef_wnorm2 * ( z_wnds_t(ji,jj+1  ) + z_wnds_t(ji,jj) )                          & 
    523                   &          * ( 0.5 * (sf(jp_wndj)%fnow(ji,jj+1,1) + sf(jp_wndj)%fnow(ji,jj,1) ) - rn_vfac * pvi(ji,jj) ) 
     492               utau_ice(ji,jj) = zcoef_wnorm2 * ( wndm_ice(ji+1,jj  ) + wndm_ice(ji,jj) )                          & 
     493                  &          * ( 0.5 * (sf(jp_wndi)%fnow(ji+1,jj,1) + sf(jp_wndi)%fnow(ji,jj,1) ) - rn_vfac * u_ice(ji,jj) ) 
     494               vtau_ice(ji,jj) = zcoef_wnorm2 * ( wndm_ice(ji,jj+1  ) + wndm_ice(ji,jj) )                          & 
     495                  &          * ( 0.5 * (sf(jp_wndj)%fnow(ji,jj+1,1) + sf(jp_wndj)%fnow(ji,jj,1) ) - rn_vfac * v_ice(ji,jj) ) 
    524496            END DO 
    525497         END DO 
    526          CALL lbc_lnk( p_taui  , 'U', -1. ) 
    527          CALL lbc_lnk( p_tauj  , 'V', -1. ) 
    528          CALL lbc_lnk( z_wnds_t, 'T',  1. ) 
     498         CALL lbc_lnk( utau_ice, 'U', -1. ) 
     499         CALL lbc_lnk( vtau_ice, 'V', -1. ) 
     500         CALL lbc_lnk( wndm_ice, 'T',  1. ) 
    529501         ! 
    530502      END SELECT 
     503 
     504      IF(ln_ctl) THEN 
     505         CALL prt_ctl(tab2d_1=utau_ice  , clinfo1=' blk_ice_core: utau_ice : ', tab2d_2=vtau_ice  , clinfo2=' vtau_ice : ') 
     506         CALL prt_ctl(tab2d_1=wndm_ice  , clinfo1=' blk_ice_core: wndm_ice : ') 
     507      ENDIF 
     508 
     509      IF( nn_timing == 1 )  CALL timing_stop('blk_ice_core_tau') 
     510       
     511   END SUBROUTINE blk_ice_core_tau 
     512 
     513 
     514   SUBROUTINE blk_ice_core_flx( ptsu, palb ) 
     515      !!--------------------------------------------------------------------- 
     516      !!                     ***  ROUTINE blk_ice_core_flx  *** 
     517      !! 
     518      !! ** Purpose :   provide the surface boundary condition over sea-ice 
     519      !! 
     520      !! ** Method  :   compute heat and freshwater exchanged 
     521      !!                between atmosphere and sea-ice using CORE bulk 
     522      !!                formulea, ice variables and read atmmospheric fields. 
     523      !!  
     524      !! caution : the net upward water flux has with mm/day unit 
     525      !!--------------------------------------------------------------------- 
     526      REAL(wp), DIMENSION(:,:,:), INTENT(in)  ::   ptsu          ! sea ice surface temperature 
     527      REAL(wp), DIMENSION(:,:,:), INTENT(in)  ::   palb          ! ice albedo (all skies) 
     528      !! 
     529      INTEGER  ::   ji, jj, jl    ! dummy loop indices 
     530      REAL(wp) ::   zst2, zst3 
     531      REAL(wp) ::   zcoef_dqlw, zcoef_dqla, zcoef_dqsb 
     532      REAL(wp) ::   zztmp, z1_lsub                               ! temporary variable 
     533      !! 
     534      REAL(wp), DIMENSION(:,:,:), POINTER ::   z_qlw             ! long wave heat flux over ice 
     535      REAL(wp), DIMENSION(:,:,:), POINTER ::   z_qsb             ! sensible  heat flux over ice 
     536      REAL(wp), DIMENSION(:,:,:), POINTER ::   z_dqlw            ! long wave heat sensitivity over ice 
     537      REAL(wp), DIMENSION(:,:,:), POINTER ::   z_dqsb            ! sensible  heat sensitivity over ice 
     538      REAL(wp), DIMENSION(:,:)  , POINTER ::   zevap, zsnw       ! evaporation and snw distribution after wind blowing (LIM3) 
     539      !!--------------------------------------------------------------------- 
     540      ! 
     541      IF( nn_timing == 1 )  CALL timing_start('blk_ice_core_flx') 
     542      ! 
     543      CALL wrk_alloc( jpi,jpj,jpl, z_qlw, z_qsb, z_dqlw, z_dqsb )  
     544 
     545      ! local scalars ( place there for vector optimisation purposes) 
     546      zcoef_dqlw   = 4.0 * 0.95 * Stef 
     547      zcoef_dqla   = -Ls * Cice * 11637800. * (-5897.8) 
     548      zcoef_dqsb   = rhoa * cpa * Cice 
    531549 
    532550      zztmp = 1. / ( 1. - albo ) 
    533551      !                                     ! ========================== ! 
    534       DO jl = 1, ijpl                       !  Loop over ice categories  ! 
     552      DO jl = 1, jpl                        !  Loop over ice categories  ! 
    535553         !                                  ! ========================== ! 
    536554         DO jj = 1 , jpj 
    … …  
    539557               !      I   Radiative FLUXES   ! 
    540558               ! ----------------------------! 
    541                zst2 = pst(ji,jj,jl) * pst(ji,jj,jl) 
    542                zst3 = pst(ji,jj,jl) * zst2 
     559               zst2 = ptsu(ji,jj,jl) * ptsu(ji,jj,jl) 
     560               zst3 = ptsu(ji,jj,jl) * zst2 
    543561               ! Short Wave (sw) 
    544                p_qsr(ji,jj,jl) = zztmp * ( 1. - palb(ji,jj,jl) ) * qsr(ji,jj) 
     562               qsr_ice(ji,jj,jl) = zztmp * ( 1. - palb(ji,jj,jl) ) * qsr(ji,jj) 
    545563               ! Long  Wave (lw) 
    546                z_qlw(ji,jj,jl) = 0.95 * ( sf(jp_qlw)%fnow(ji,jj,1) - Stef * pst(ji,jj,jl) * zst3 ) * tmask(ji,jj,1) 
     564               z_qlw(ji,jj,jl) = 0.95 * ( sf(jp_qlw)%fnow(ji,jj,1) - Stef * ptsu(ji,jj,jl) * zst3 ) * tmask(ji,jj,1) 
    547565               ! lw sensitivity 
    548566               z_dqlw(ji,jj,jl) = zcoef_dqlw * zst3                                                
    … …  
    554572               ! ... turbulent heat fluxes 
    555573               ! Sensible Heat 
    556                z_qsb(ji,jj,jl) = rhoa * cpa * Cice * z_wnds_t(ji,jj) * ( pst(ji,jj,jl) - sf(jp_tair)%fnow(ji,jj,1) ) 
     574               z_qsb(ji,jj,jl) = rhoa * cpa * Cice * wndm_ice(ji,jj) * ( ptsu(ji,jj,jl) - sf(jp_tair)%fnow(ji,jj,1) ) 
    557575               ! Latent Heat 
    558                p_qla(ji,jj,jl) = rn_efac * MAX( 0.e0, rhoa * Ls  * Cice * z_wnds_t(ji,jj)   &                            
    559                   &                         * (  11637800. * EXP( -5897.8 / pst(ji,jj,jl) ) / rhoa - sf(jp_humi)%fnow(ji,jj,1)  ) ) 
    560                ! Latent heat sensitivity for ice (Dqla/Dt) 
    561                IF( p_qla(ji,jj,jl) > 0._wp ) THEN 
    562                   p_dqla(ji,jj,jl) = rn_efac * zcoef_dqla * z_wnds_t(ji,jj) / ( zst2 ) * EXP( -5897.8 / pst(ji,jj,jl) ) 
     576               qla_ice(ji,jj,jl) = rn_efac * MAX( 0.e0, rhoa * Ls  * Cice * wndm_ice(ji,jj)   &                            
     577                  &                         * (  11637800. * EXP( -5897.8 / ptsu(ji,jj,jl) ) / rhoa - sf(jp_humi)%fnow(ji,jj,1)  ) ) 
     578              ! Latent heat sensitivity for ice (Dqla/Dt) 
     579               IF( qla_ice(ji,jj,jl) > 0._wp ) THEN 
     580                  dqla_ice(ji,jj,jl) = rn_efac * zcoef_dqla * wndm_ice(ji,jj) / ( zst2 ) * EXP( -5897.8 / ptsu(ji,jj,jl) ) 
    563581               ELSE 
    564                   p_dqla(ji,jj,jl) = 0._wp 
     582                  dqla_ice(ji,jj,jl) = 0._wp 
    565583               ENDIF 
    566584 
    567585               ! Sensible heat sensitivity (Dqsb_ice/Dtn_ice) 
    568                z_dqsb(ji,jj,jl) = zcoef_dqsb * z_wnds_t(ji,jj) 
     586               z_dqsb(ji,jj,jl) = zcoef_dqsb * wndm_ice(ji,jj) 
    569587 
    570588               ! ----------------------------! 
    … …  
    572590               ! ----------------------------! 
    573591               ! Downward Non Solar flux 
    574                p_qns (ji,jj,jl) =     z_qlw (ji,jj,jl) - z_qsb (ji,jj,jl) - p_qla (ji,jj,jl) 
     592               qns_ice (ji,jj,jl) =     z_qlw (ji,jj,jl) - z_qsb (ji,jj,jl) - qla_ice (ji,jj,jl) 
    575593               ! Total non solar heat flux sensitivity for ice 
    576                p_dqns(ji,jj,jl) = - ( z_dqlw(ji,jj,jl) + z_dqsb(ji,jj,jl) + p_dqla(ji,jj,jl) ) 
     594               dqns_ice(ji,jj,jl) = - ( z_dqlw(ji,jj,jl) + z_dqsb(ji,jj,jl) + dqla_ice(ji,jj,jl) ) 
    577595            END DO 
    578596            ! 
    … …  
    581599      END DO 
    582600      ! 
     601      tprecip(:,:) = sf(jp_prec)%fnow(:,:,1) * rn_pfac      ! total precipitation [kg/m2/s] 
     602      sprecip(:,:) = sf(jp_snow)%fnow(:,:,1) * rn_pfac      ! solid precipitation [kg/m2/s] 
     603      CALL iom_put( 'snowpre', sprecip * 86400. )                  ! Snow precipitation 
     604      CALL iom_put( 'precip' , tprecip * 86400. )                  ! Total precipitation 
     605 
     606#if defined  key_lim3 
     607      CALL wrk_alloc( jpi,jpj, zevap, zsnw )  
     608 
     609      ! --- evaporation --- ! 
     610      z1_lsub = 1._wp / Lsub 
     611      evap_ice (:,:,:) = qla_ice (:,:,:) * z1_lsub ! sublimation 
     612      devap_ice(:,:,:) = dqla_ice(:,:,:) * z1_lsub 
     613      zevap    (:,:)   = emp(:,:) + tprecip(:,:)   ! evaporation over ocean 
     614 
     615      ! --- evaporation minus precipitation --- ! 
     616      CALL lim_thd_snwblow( pfrld, zsnw )  ! snow distribution over ice after wind blowing  
     617      emp_oce(:,:) = pfrld(:,:) * zevap(:,:) - ( tprecip(:,:) - sprecip(:,:) ) - sprecip(:,:) * (1._wp - zsnw ) 
     618      emp_ice(:,:) = SUM( a_i_b(:,:,:) * evap_ice(:,:,:), dim=3 ) - sprecip(:,:) * zsnw 
     619      emp_tot(:,:) = emp_oce(:,:) + emp_ice(:,:) 
     620 
     621      ! --- heat flux associated with emp --- ! 
     622      qemp_oce(:,:) = - pfrld(:,:) * zevap(:,:) * sst_m(:,:) * rcp                               & ! evap 
     623         &          + ( tprecip(:,:) - sprecip(:,:) ) * ( sf(jp_tair)%fnow(:,:,1) - rt0 ) * rcp  & ! liquid precip 
     624         &          +   sprecip(:,:) * ( 1._wp - zsnw ) *                                        & ! solid precip 
     625         &              ( ( MIN( sf(jp_tair)%fnow(:,:,1), rt0_snow ) - rt0 ) * cpic * tmask(:,:,1) - lfus ) 
     626      qemp_ice(:,:) =   sprecip(:,:) * zsnw *                                                    & ! solid precip (only) 
     627         &              ( ( MIN( sf(jp_tair)%fnow(:,:,1), rt0_snow ) - rt0 ) * cpic * tmask(:,:,1) - lfus ) 
     628 
     629      ! --- total solar and non solar fluxes --- ! 
     630      qns_tot(:,:) = pfrld(:,:) * qns_oce(:,:) + SUM( a_i_b(:,:,:) * qns_ice(:,:,:), dim=3 ) + qemp_ice(:,:) + qemp_oce(:,:) 
     631      qsr_tot(:,:) = pfrld(:,:) * qsr_oce(:,:) + SUM( a_i_b(:,:,:) * qsr_ice(:,:,:), dim=3 ) 
     632 
     633      ! --- heat content of precip over ice in J/m3 (to be used in 1D-thermo) --- ! 
     634      qprec_ice(:,:) = rhosn * ( ( MIN( sf(jp_tair)%fnow(:,:,1), rt0_snow ) - rt0 ) * cpic * tmask(:,:,1) - lfus ) 
     635 
     636      CALL wrk_dealloc( jpi,jpj, zevap, zsnw )  
     637#endif 
     638 
    583639      !-------------------------------------------------------------------- 
    584640      ! FRACTIONs of net shortwave radiation which is not absorbed in the 
    … …  
    586642      ! ( Maykut and Untersteiner, 1971 ; Ebert and Curry, 1993 ) 
    587643      ! 
    588       p_fr1(:,:) = ( 0.18 * ( 1.0 - cldf_ice ) + 0.35 * cldf_ice ) 
    589       p_fr2(:,:) = ( 0.82 * ( 1.0 - cldf_ice ) + 0.65 * cldf_ice ) 
    590       ! 
    591       p_tpr(:,:) = sf(jp_prec)%fnow(:,:,1) * rn_pfac      ! total precipitation [kg/m2/s] 
    592       p_spr(:,:) = sf(jp_snow)%fnow(:,:,1) * rn_pfac      ! solid precipitation [kg/m2/s] 
    593       CALL iom_put( 'snowpre', p_spr * 86400. )                  ! Snow precipitation 
    594       CALL iom_put( 'precip' , p_tpr * 86400. )                  ! Total precipitation 
     644      fr1_i0(:,:) = ( 0.18 * ( 1.0 - cldf_ice ) + 0.35 * cldf_ice ) 
     645      fr2_i0(:,:) = ( 0.82 * ( 1.0 - cldf_ice ) + 0.65 * cldf_ice ) 
     646      ! 
    595647      ! 
    596648      IF(ln_ctl) THEN 
    597          CALL prt_ctl(tab3d_1=p_qla   , clinfo1=' blk_ice_core: p_qla  : ', tab3d_2=z_qsb   , clinfo2=' z_qsb  : ', kdim=ijpl) 
    598          CALL prt_ctl(tab3d_1=z_qlw   , clinfo1=' blk_ice_core: z_qlw  : ', tab3d_2=p_dqla  , clinfo2=' p_dqla : ', kdim=ijpl) 
    599          CALL prt_ctl(tab3d_1=z_dqsb  , clinfo1=' blk_ice_core: z_dqsb : ', tab3d_2=z_dqlw  , clinfo2=' z_dqlw : ', kdim=ijpl) 
    600          CALL prt_ctl(tab3d_1=p_dqns  , clinfo1=' blk_ice_core: p_dqns : ', tab3d_2=p_qsr   , clinfo2=' p_qsr  : ', kdim=ijpl) 
    601          CALL prt_ctl(tab3d_1=pst     , clinfo1=' blk_ice_core: pst    : ', tab3d_2=p_qns   , clinfo2=' p_qns  : ', kdim=ijpl) 
    602          CALL prt_ctl(tab2d_1=p_tpr   , clinfo1=' blk_ice_core: p_tpr  : ', tab2d_2=p_spr   , clinfo2=' p_spr  : ') 
    603          CALL prt_ctl(tab2d_1=p_taui  , clinfo1=' blk_ice_core: p_taui : ', tab2d_2=p_tauj  , clinfo2=' p_tauj : ') 
    604          CALL prt_ctl(tab2d_1=z_wnds_t, clinfo1=' blk_ice_core: z_wnds_t : ') 
     649         CALL prt_ctl(tab3d_1=qla_ice , clinfo1=' blk_ice_core: qla_ice  : ', tab3d_2=z_qsb   , clinfo2=' z_qsb    : ', kdim=jpl) 
     650         CALL prt_ctl(tab3d_1=z_qlw   , clinfo1=' blk_ice_core: z_qlw    : ', tab3d_2=dqla_ice, clinfo2=' dqla_ice : ', kdim=jpl) 
     651         CALL prt_ctl(tab3d_1=z_dqsb  , clinfo1=' blk_ice_core: z_dqsb   : ', tab3d_2=z_dqlw  , clinfo2=' z_dqlw   : ', kdim=jpl) 
     652         CALL prt_ctl(tab3d_1=dqns_ice, clinfo1=' blk_ice_core: dqns_ice : ', tab3d_2=qsr_ice , clinfo2=' qsr_ice  : ', kdim=jpl) 
     653         CALL prt_ctl(tab3d_1=ptsu    , clinfo1=' blk_ice_core: ptsu     : ', tab3d_2=qns_ice , clinfo2=' qns_ice  : ', kdim=jpl) 
     654         CALL prt_ctl(tab2d_1=tprecip , clinfo1=' blk_ice_core: tprecip  : ', tab2d_2=sprecip , clinfo2=' sprecip  : ') 
    605655      ENDIF 
    606656 
    607       CALL wrk_dealloc( jpi,jpj,   z_wnds_t ) 
    608       CALL wrk_dealloc( jpi,jpj,   pdim, z_qlw, z_qsb, z_dqlw, z_dqsb ) 
    609       ! 
    610       IF( nn_timing == 1 )  CALL timing_stop('blk_ice_core') 
    611       ! 
    612    END SUBROUTINE blk_ice_core 
     657      CALL wrk_dealloc( jpi,jpj,jpl, z_qlw, z_qsb, z_dqlw, z_dqsb ) 
     658      ! 
     659      IF( nn_timing == 1 )  CALL timing_stop('blk_ice_core_flx') 
     660       
     661   END SUBROUTINE blk_ice_core_flx 
     662#endif 
    613663 
    614664 

  • branches/2015/dev_r5218_CNRS17_coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90

    r5352 r5357  
    4848   USE ice_domain_size, only: ncat 
    4949#endif 
     50#if defined key_lim3 
     51   USE limthd_dh       ! for CALL lim_thd_snwblow 
     52#endif 
     53 
    5054   IMPLICIT NONE 
    5155   PRIVATE 
    … …  
    13591363      REAL(wp), POINTER, DIMENSION(:,:  ) ::   zemp_tot, zemp_ice, zsprecip, ztprecip, zqns_tot, zqsr_tot 
    13601364      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zqns_ice, zqsr_ice, zdqns_ice 
     1365      REAL(wp), POINTER, DIMENSION(:,:  ) ::   zevap, zsnw, zqns_oce, zqprec_ice, zqemp_oce ! for LIM3 
    13611366      !!---------------------------------------------------------------------- 
    13621367      ! 
    … …  
    14691474            &                                                   +          pist(:,:,1)   * zicefr(:,:) ) ) 
    14701475      END SELECT 
    1471       ztmp(:,:) = p_frld(:,:) * zsprecip(:,:) * lfus 
    1472       zqns_tot(:,:) = zqns_tot(:,:)                       &            ! zqns_tot update over free ocean with: 
    1473          &          - ztmp(:,:)                           &            ! remove the latent heat flux of solid precip. melting 
    1474          &          - (  zemp_tot(:,:)                    &            ! remove the heat content of mass flux (assumed to be at SST) 
    1475          &             - zemp_ice(:,:) * zicefr(:,:)  ) * zcptn(:,:)  
    1476       IF( iom_use('hflx_snow_cea') )   & 
    1477          CALL iom_put( 'hflx_snow_cea', ztmp + sprecip(:,:) * zcptn(:,:) )   ! heat flux from snow (cell average) 
    14781476!!gm 
    14791477!!    currently it is taken into account in leads budget but not in the zqns_tot, and thus not in  
    … …  
    14911489      ENDIF 
    14921490 
     1491      ztmp(:,:) = p_frld(:,:) * zsprecip(:,:) * lfus 
     1492      IF( iom_use('hflx_snow_cea') )    CALL iom_put( 'hflx_snow_cea', ztmp + sprecip(:,:) * zcptn(:,:) )   ! heat flux from snow (cell average) 
     1493 
     1494#if defined key_lim3 
     1495      CALL wrk_alloc( jpi,jpj, zevap, zsnw, zqns_oce, zqprec_ice, zqemp_oce )  
     1496 
     1497      ! --- evaporation --- ! 
     1498      ! clem: evap_ice is set to 0 for LIM3 since we still do not know what to do with sublimation 
     1499      ! the problem is: the atm. imposes both mass evaporation and heat removed from the snow/ice 
     1500      !                 but it is incoherent WITH the ice model   
     1501      DO jl=1,jpl 
     1502         evap_ice(:,:,jl) = 0._wp  ! should be: frcv(jpr_ievp)%z3(:,:,1) 
     1503      ENDDO 
     1504      zevap(:,:) = zemp_tot(:,:) + ztprecip(:,:) ! evaporation over ocean 
     1505 
     1506      ! --- evaporation minus precipitation --- ! 
     1507      emp_oce(:,:) = emp_tot(:,:) - emp_ice(:,:) 
     1508 
     1509      ! --- non solar flux over ocean --- ! 
     1510      !         note: p_frld cannot be = 0 since we limit the ice concentration to amax 
     1511      zqns_oce = 0._wp 
     1512      WHERE( p_frld /= 0._wp )  zqns_oce(:,:) = ( zqns_tot(:,:) - SUM( a_i * zqns_ice, dim=3 ) ) / p_frld(:,:) 
     1513 
     1514      ! --- heat flux associated with emp --- ! 
     1515      CALL lim_thd_snwblow( p_frld, zsnw )  ! snow distribution over ice after wind blowing 
     1516      zqemp_oce(:,:) = -      zevap(:,:)                   * p_frld(:,:)      *   zcptn(:,:)   &      ! evap 
     1517         &             + ( ztprecip(:,:) - zsprecip(:,:) )                    *   zcptn(:,:)   &      ! liquid precip 
     1518         &             +   zsprecip(:,:)                   * ( 1._wp - zsnw ) * ( zcptn(:,:) - lfus ) ! solid precip over ocean 
     1519      qemp_ice(:,:)  = -   frcv(jpr_ievp)%z3(:,:,1)        * zicefr(:,:)      *   zcptn(:,:)   &      ! ice evap 
     1520         &             +   zsprecip(:,:)                   * zsnw             * ( zcptn(:,:) - lfus ) ! solid precip over ice 
     1521 
     1522      ! --- heat content of precip over ice in J/m3 (to be used in 1D-thermo) --- ! 
     1523      zqprec_ice(:,:) = rhosn * ( zcptn(:,:) - lfus ) 
     1524 
     1525      ! --- total non solar flux --- ! 
     1526      zqns_tot(:,:) = zqns_tot(:,:) + qemp_ice(:,:) + zqemp_oce(:,:) 
     1527 
     1528      ! --- in case both coupled/forced are active, we must mix values --- !  
    14931529      IF( ln_mixcpl ) THEN 
     1530         qns_tot(:,:) = qns_tot(:,:) * xcplmask(:,:,0) + zqns_tot(:,:)* zmsk(:,:) 
     1531         qns_oce(:,:) = qns_oce(:,:) * xcplmask(:,:,0) + zqns_oce(:,:)* zmsk(:,:) 
     1532         DO jl=1,jpl 
     1533            qns_ice(:,:,jl) = qns_ice(:,:,jl) * xcplmask(:,:,0) +  zqns_ice(:,:,jl)* zmsk(:,:) 
     1534         ENDDO 
     1535         qprec_ice(:,:) = qprec_ice(:,:) * xcplmask(:,:,0) + zqprec_ice(:,:)* zmsk(:,:) 
     1536         qemp_oce(:,:) = qemp_oce(:,:) * xcplmask(:,:,0) + zqemp_oce(:,:)* zmsk(:,:) 
     1537!!clem         evap_ice(:,:) = evap_ice(:,:) * xcplmask(:,:,0) 
     1538      ELSE 
     1539         qns_tot  (:,:  ) = zqns_tot  (:,:  ) 
     1540         qns_oce  (:,:  ) = zqns_oce  (:,:  ) 
     1541         qns_ice  (:,:,:) = zqns_ice  (:,:,:) 
     1542         qprec_ice(:,:)   = zqprec_ice(:,:) 
     1543         qemp_oce (:,:)   = zqemp_oce (:,:) 
     1544      ENDIF 
     1545 
     1546      CALL wrk_dealloc( jpi,jpj, zevap, zsnw, zqns_oce, zqprec_ice, zqemp_oce )  
     1547 
     1548#else 
     1549 
     1550      ! clem: this formulation is certainly wrong. 
     1551      zqns_tot(:,:) = zqns_tot(:,:)                       &            ! zqns_tot update over free ocean with: 
     1552         &          - ztmp(:,:)                           &            ! remove the latent heat flux of solid precip. melting 
     1553         &          - (  zemp_tot(:,:)                    &            ! remove the heat content of mass flux (assumed to be at SST) 
     1554         &             - zemp_ice(:,:) * zicefr(:,:)  ) * zcptn(:,:)  
     1555 
     1556     IF( ln_mixcpl ) THEN 
    14941557         qns_tot(:,:) = qns(:,:) * p_frld(:,:) + SUM( qns_ice(:,:,:) * a_i(:,:,:), dim=3 )   ! total flux from blk 
    14951558         qns_tot(:,:) = qns_tot(:,:) * xcplmask(:,:,0) +  zqns_tot(:,:)* zmsk(:,:) 
    … …  
    15011564         qns_ice(:,:,:) = zqns_ice(:,:,:) 
    15021565      ENDIF 
     1566 
     1567#endif 
    15031568 
    15041569      !                                                      ! ========================= ! 

  • branches/2015/dev_r5218_CNRS17_coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90

    r5220 r5357  
    116116 
    117117      IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN     !  Ice time-step only 
     118 
    118119         !-----------------------!                                            
    119120         ! --- Bulk Formulae --- !                                            
    … …  
    125126         t_bo(:,:) = ( eos_fzp( sss_m ) + rt0 ) * tmask(:,:,1) + rt0 * ( 1._wp - tmask(:,:,1) )   
    126127         !                                                                                       
    127          ! Ice albedo 
    128          CALL wrk_alloc( jpi,jpj    , zutau_ice, zvtau_ice) 
    129          CALL wrk_alloc( jpi,jpj,jpl, zalb_os, zalb_cs, zalb_ice ) 
    130          CALL albedo_ice( t_su, ht_i, ht_s, zalb_cs, zalb_os )  ! cloud-sky and overcast-sky ice albedos 
    131  
    132          ! CORE and COUPLED bulk formulations 
    133          SELECT CASE( kblk ) 
    134          CASE( jp_core , jp_cpl ) 
    135  
    136             ! albedo depends on cloud fraction because of non-linear spectral effects 
    137             zalb_ice(:,:,:) = ( 1. - cldf_ice ) * zalb_cs(:,:,:) + cldf_ice * zalb_os(:,:,:) 
    138             ! In CLIO the cloud fraction is read in the climatology and the all-sky albedo  
    139             ! (zalb_ice) is computed within the bulk routine 
    140              
    141          END SELECT 
     128!!clem         ! Ice albedo 
     129!!clem         CALL wrk_alloc( jpi,jpj,jpl, zalb_os, zalb_cs, zalb_ice ) 
     130!!clem         CALL albedo_ice( t_su, ht_i, ht_s, zalb_cs, zalb_os )  ! cloud-sky and overcast-sky ice albedos 
     131!! 
     132!!         ! CORE and COUPLED bulk formulations 
     133!!         SELECT CASE( kblk ) 
     134!!         CASE( jp_core , jp_cpl ) 
     135!!            ! albedo depends on cloud fraction because of non-linear spectral effects 
     136!!            zalb_ice(:,:,:) = ( 1. - cldf_ice ) * zalb_cs(:,:,:) + cldf_ice * zalb_os(:,:,:) 
     137!!            ! In CLIO the cloud fraction is read in the climatology and the all-sky albedo  
     138!!            ! (zalb_ice) is computed within the bulk routine 
     139!!clem         END SELECT 
    142140          
    143141         ! Mask sea ice surface temperature (set to rt0 over land) 
    … …  
    156154         SELECT CASE( kblk ) 
    157155         CASE( jp_clio )                                       ! CLIO bulk formulation 
    158             CALL blk_ice_clio( t_su , zalb_cs    , zalb_os    , zalb_ice  ,               & 
    159                &                      utau_ice   , vtau_ice   , qns_ice   , qsr_ice   ,   & 
    160                &                      qla_ice    , dqns_ice   , dqla_ice  ,               & 
    161                &                      tprecip    , sprecip    ,                           & 
    162                &                      fr1_i0     , fr2_i0     , cp_ice_msh, jpl  ) 
    163             !          
    164             IF( nn_limflx /= 2 )   CALL ice_lim_flx( t_su, zalb_ice, qns_ice, qsr_ice ,   & 
    165                &                                           dqns_ice, qla_ice, dqla_ice, nn_limflx ) 
     156!!clem            CALL blk_ice_clio( t_su , zalb_cs    , zalb_os    , zalb_ice  ,               & 
     157!!               &                      utau_ice   , vtau_ice   , qns_ice   , qsr_ice   ,   & 
     158!!               &                      qla_ice    , dqns_ice   , dqla_ice  ,               & 
     159!!               &                      tprecip    , sprecip    ,                           & 
     160!!               &                      fr1_i0     , fr2_i0     , cp_ice_msh, jpl  ) 
     161!!            !          
     162!!            IF( nn_limflx /= 2 )   CALL ice_lim_flx( t_su, zalb_ice, qns_ice, qsr_ice ,   & 
     163!!               &                                           dqns_ice, qla_ice, dqla_ice, nn_limflx ) 
     164            CALL blk_ice_clio_tau( utau_ice, vtau_ice, cp_ice_msh ) 
    166165 
    167166         CASE( jp_core )                                       ! CORE bulk formulation 
    168             CALL blk_ice_core( t_su , u_ice     , v_ice     , zalb_ice   ,               & 
    169                &                      utau_ice  , vtau_ice  , qns_ice    , qsr_ice   ,   & 
    170                &                      qla_ice   , dqns_ice  , dqla_ice   ,               & 
    171                &                      tprecip   , sprecip   ,                            & 
    172                &                      fr1_i0    , fr2_i0    , cp_ice_msh, jpl  ) 
    173                ! 
    174             IF( nn_limflx /= 2 )   CALL ice_lim_flx( t_su, zalb_ice, qns_ice, qsr_ice ,   & 
    175                &                                           dqns_ice, qla_ice, dqla_ice, nn_limflx ) 
     167!!clem            CALL blk_ice_core( t_su , u_ice     , v_ice     , zalb_ice   ,               & 
     168!!clem               &                      utau_ice  , vtau_ice  , qns_ice    , qsr_ice   ,   & 
     169!!clem               &                      qla_ice   , dqns_ice  , dqla_ice   ,               & 
     170!!clem               &                      tprecip   , sprecip   ,                            & 
     171!!clem               &                      fr1_i0    , fr2_i0    , cp_ice_msh, jpl  ) 
     172!!clem            IF( nn_limflx /= 2 )   CALL ice_lim_flx( t_su, zalb_ice, qns_ice, qsr_ice ,   & 
     173!!clem               &                                           dqns_ice, qla_ice, dqla_ice, nn_limflx ) 
     174            CALL blk_ice_core_tau 
    176175            ! 
    177176         CASE ( jp_cpl ) 
    … …  
    182181          
    183182         IF( ln_mixcpl) THEN 
     183            CALL wrk_alloc( jpi,jpj    , zutau_ice, zvtau_ice) 
    184184            CALL sbc_cpl_ice_tau( zutau_ice , zvtau_ice ) 
    185185            utau_ice(:,:) = utau_ice(:,:) * xcplmask(:,:,0) + zutau_ice(:,:) * ( 1. - xcplmask(:,:,0) ) 
    186186            vtau_ice(:,:) = vtau_ice(:,:) * xcplmask(:,:,0) + zvtau_ice(:,:) * ( 1. - xcplmask(:,:,0) ) 
     187            CALL wrk_dealloc( jpi,jpj  , zutau_ice, zvtau_ice) 
    187188         ENDIF 
    188189 
    … …  
    229230         phicif(:,:)  = vt_i(:,:) 
    230231          
     232         ! Ice albedo 
     233         CALL wrk_alloc( jpi,jpj,jpl, zalb_os, zalb_cs, zalb_ice ) 
     234         CALL albedo_ice( t_su, ht_i, ht_s, zalb_cs, zalb_os )  ! cloud-sky and overcast-sky ice albedos 
     235  
    231236         SELECT CASE( kblk ) 
     237         CASE( jp_clio )                                       ! CLIO bulk formulation 
     238            ! In CLIO the cloud fraction is read in the climatology and the all-sky albedo  
     239            ! (zalb_ice) is computed within the bulk routine 
     240            CALL blk_ice_clio_flx( t_su , zalb_cs, zalb_os  , zalb_ice, qns_ice   , qsr_ice   ,    & 
     241               &                      qla_ice, dqns_ice   , dqla_ice  , tprecip, sprecip    ,  & 
     242               &                      fr1_i0     , fr2_i0     , jpl  ) 
     243            !          
     244            IF( nn_limflx /= 2 )   CALL ice_lim_flx( t_su, zalb_ice, qns_ice, qsr_ice ,   & 
     245               &                                           dqns_ice, evap_ice, devap_ice, nn_limflx ) 
     246 
     247         CASE( jp_core )                                       ! CORE bulk formulation 
     248            ! albedo depends on cloud fraction because of non-linear spectral effects 
     249            zalb_ice(:,:,:) = ( 1. - cldf_ice ) * zalb_cs(:,:,:) + cldf_ice * zalb_os(:,:,:) 
     250            CALL blk_ice_core_flx( t_su, zalb_ice ) 
     251            IF( nn_limflx /= 2 )   CALL ice_lim_flx( t_su, zalb_ice, qns_ice, qsr_ice ,   & 
     252               &                                           dqns_ice, evap_ice, devap_ice, nn_limflx ) 
     253 
    232254         CASE ( jp_cpl ) 
     255            ! albedo depends on cloud fraction because of non-linear spectral effects 
     256            zalb_ice(:,:,:) = ( 1. - cldf_ice ) * zalb_cs(:,:,:) + cldf_ice * zalb_os(:,:,:) 
    233257            CALL sbc_cpl_ice_flx( p_frld=pfrld, palbi=zalb_ice, psst=sst_m, pist=t_su    ) 
    234258            IF( nn_limflx == 2 )   CALL ice_lim_flx( t_su, zalb_ice, qns_ice, qsr_ice ,   & 
    235                &                                           dqns_ice, qla_ice, dqla_ice, nn_limflx ) 
     259               &                                           dqns_ice, evap_ice, devap_ice, nn_limflx ) 
    236260            ! Latent heat flux is forced to 0 in coupled: it is included in qns (non-solar heat flux) 
    237             qla_ice  (:,:,:) = 0._wp 
    238             dqla_ice (:,:,:) = 0._wp 
     261            evap_ice  (:,:,:) = 0._wp 
     262            devap_ice (:,:,:) = 0._wp 
     263 
    239264         END SELECT 
     265         CALL wrk_dealloc( jpi,jpj,jpl, zalb_os, zalb_cs, zalb_ice ) 
     266 
    240267         ! 
    241268         CALL lim_thd( kt )                         ! Ice thermodynamics  
    … …  
    256283         IF( ln_icectl )  CALL lim_ctl( kt )        ! alerts in case of model crash 
    257284         ! 
    258          CALL wrk_dealloc( jpi,jpj  , zutau_ice, zvtau_ice) 
    259          CALL wrk_dealloc( jpi,jpj,jpl, zalb_os, zalb_cs, zalb_ice ) 
    260285         ! 
    261286      ENDIF   ! End sea-ice time step only 
    … …  
    486511    
    487512   SUBROUTINE ice_lim_flx( ptn_ice, palb_ice, pqns_ice, pqsr_ice,   & 
    488          &                          pdqn_ice, pqla_ice, pdql_ice, k_limflx ) 
     513         &                          pdqn_ice, pevap_ice, pdevap_ice, k_limflx ) 
    489514      !!--------------------------------------------------------------------- 
    490515      !!                  ***  ROUTINE ice_lim_flx  *** 
    … …  
    504529      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pqsr_ice   ! net solar flux 
    505530      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pdqn_ice   ! non solar flux sensitivity 
    506       REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pqla_ice   ! latent heat flux 
    507       REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pdql_ice   ! latent heat flux sensitivity 
     531      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pevap_ice  ! sublimation 
     532      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pdevap_ice ! sublimation sensitivity 
    508533      ! 
    509534      INTEGER  ::   jl      ! dummy loop index 
    … …  
    514539      REAL(wp), POINTER, DIMENSION(:,:) :: z_qsr_m   ! Mean solar heat flux over all categories 
    515540      REAL(wp), POINTER, DIMENSION(:,:) :: z_qns_m   ! Mean non solar heat flux over all categories 
    516       REAL(wp), POINTER, DIMENSION(:,:) :: z_qla_m   ! Mean latent heat flux over all categories 
     541      REAL(wp), POINTER, DIMENSION(:,:) :: z_evap_m  ! Mean sublimation over all categories 
    517542      REAL(wp), POINTER, DIMENSION(:,:) :: z_dqn_m   ! Mean d(qns)/dT over all categories 
    518       REAL(wp), POINTER, DIMENSION(:,:) :: z_dql_m   ! Mean d(qla)/dT over all categories 
     543      REAL(wp), POINTER, DIMENSION(:,:) :: z_devap_m ! Mean d(evap)/dT over all categories 
    519544      !!---------------------------------------------------------------------- 
    520545 
    … …  
    524549      SELECT CASE( k_limflx )                              !==  averaged on all ice categories  ==! 
    525550      CASE( 0 , 1 ) 
    526          CALL wrk_alloc( jpi,jpj, z_qsr_m, z_qns_m, z_qla_m, z_dqn_m, z_dql_m) 
     551         CALL wrk_alloc( jpi,jpj, z_qsr_m, z_qns_m, z_evap_m, z_dqn_m, z_devap_m) 
    527552         ! 
    528553         z_qns_m(:,:) = fice_ice_ave ( pqns_ice (:,:,:) ) 
    529554         z_qsr_m(:,:) = fice_ice_ave ( pqsr_ice (:,:,:) ) 
    530555         z_dqn_m(:,:) = fice_ice_ave ( pdqn_ice (:,:,:) ) 
    531          z_qla_m(:,:) = fice_ice_ave ( pqla_ice (:,:,:) ) 
    532          z_dql_m(:,:) = fice_ice_ave ( pdql_ice (:,:,:) ) 
     556         z_evap_m(:,:) = fice_ice_ave ( pevap_ice (:,:,:) ) 
     557         z_devap_m(:,:) = fice_ice_ave ( pdevap_ice (:,:,:) ) 
    533558         DO jl = 1, jpl 
    534559            pdqn_ice(:,:,jl) = z_dqn_m(:,:) 
    535             pdql_ice(:,:,jl) = z_dql_m(:,:) 
     560            pdevap_ice(:,:,jl) = z_devap_m(:,:) 
    536561         END DO 
    537562         ! 
    … …  
    539564            pqns_ice(:,:,jl) = z_qns_m(:,:) 
    540565            pqsr_ice(:,:,jl) = z_qsr_m(:,:) 
    541             pqla_ice(:,:,jl) = z_qla_m(:,:) 
     566            pevap_ice(:,:,jl) = z_evap_m(:,:) 
    542567         END DO 
    543568         ! 
    544          CALL wrk_dealloc( jpi,jpj, z_qsr_m, z_qns_m, z_qla_m, z_dqn_m, z_dql_m) 
     569         CALL wrk_dealloc( jpi,jpj, z_qsr_m, z_qns_m, z_evap_m, z_dqn_m, z_devap_m) 
    545570      END SELECT 
    546571 
    … …  
    553578         DO jl = 1, jpl 
    554579            pqns_ice(:,:,jl) = pqns_ice(:,:,jl) + pdqn_ice(:,:,jl) * (ptn_ice(:,:,jl) - ztem_m(:,:)) 
    555             pqla_ice(:,:,jl) = pqla_ice(:,:,jl) + pdql_ice(:,:,jl) * (ptn_ice(:,:,jl) - ztem_m(:,:)) 
     580            pevap_ice(:,:,jl) = pevap_ice(:,:,jl) + pdevap_ice(:,:,jl) * (ptn_ice(:,:,jl) - ztem_m(:,:)) 
    556581            pqsr_ice(:,:,jl) = pqsr_ice(:,:,jl) * ( 1._wp - palb_ice(:,:,jl) ) / ( 1._wp - zalb_m(:,:) )  
    557582         END DO 
    … …  
    603628      wfx_spr(:,:) = 0._wp   ;    
    604629       
    605       hfx_in (:,:) = 0._wp   ;   hfx_out(:,:) = 0._wp 
    606630      hfx_thd(:,:) = 0._wp   ;    
    607631      hfx_snw(:,:) = 0._wp   ;   hfx_opw(:,:) = 0._wp 
    … …  
    620644       
    621645   END SUBROUTINE sbc_lim_diag0 
    622        
     646 
     647      
    623648   FUNCTION fice_cell_ave ( ptab ) 
    624649      !!-------------------------------------------------------------------------- 

  • branches/2015/dev_r5218_CNRS17_coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim_2.F90

    r5220 r5357  
    192192 
    193193         CASE( jp_core )           ! CORE bulk formulation 
    194             CALL blk_ice_core( zsist, u_ice      , v_ice      , zalb_ice   ,            & 
    195                &                      utau_ice   , vtau_ice   , qns_ice    , qsr_ice,   & 
    196                &                      qla_ice    , dqns_ice   , dqla_ice   ,            & 
    197                &                      tprecip    , sprecip    ,                         & 
    198                &                      fr1_i0     , fr2_i0     , cp_ice_msh , jpl  ) 
     194            CALL blk_ice_core_tau 
     195            CALL blk_ice_core_flx( zsist, zalb_ice ) 
    199196            IF( ltrcdm2dc_ice )   CALL blk_ice_meanqsr( zalb_ice, qsr_ice_mean, jpl ) 
    200197 

  • branches/2015/dev_r5218_CNRS17_coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90

    r5299 r5357  
    417417      !                                                           ! (update freshwater fluxes) 
    418418!RBbug do not understand why see ticket 667 
    419       !clem-bugsal CALL lbc_lnk( emp, 'T', 1. ) 
     419!clem: it looks like it is necessary for the north fold (in certain circumstances). Don't know why. 
     420      CALL lbc_lnk( emp, 'T', 1. ) 
    420421      ! 
    421422      IF( kt == nit000 ) THEN                          !   set the forcing field at nit000 - 1    ! 
Note: See TracChangeset for help on using the changeset viewer.