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 12565 for NEMO/branches/2020/dev_r12472_ASINTER-05_Masson_CurrentFeedback/src/OCE/SBC/sbcblk.F90 – NEMO

Ignore:
Timestamp:
2020-03-17T15:51:51+01:00 (4 years ago)
Author:
smasson
Message:

dev_r12472_ASINTER-05_Masson_CurrentFeedback: code inmplemenation, see #2156

File:
1 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2020/dev_r12472_ASINTER-05_Masson_CurrentFeedback/src/OCE/SBC/sbcblk.F90

    r12495 r12565  
    7474#endif 
    7575 
    76    INTEGER , PUBLIC            ::   jpfld         ! maximum number of files to read 
    77    INTEGER , PUBLIC, PARAMETER ::   jp_wndi = 1   ! index of 10m wind velocity (i-component) (m/s)    at T-point 
    78    INTEGER , PUBLIC, PARAMETER ::   jp_wndj = 2   ! index of 10m wind velocity (j-component) (m/s)    at T-point 
    79    INTEGER , PUBLIC, PARAMETER ::   jp_tair = 3   ! index of 10m air temperature             (Kelvin) 
    80    INTEGER , PUBLIC, PARAMETER ::   jp_humi = 4   ! index of specific humidity               ( % ) 
    81    INTEGER , PUBLIC, PARAMETER ::   jp_qsr  = 5   ! index of solar heat                      (W/m2) 
    82    INTEGER , PUBLIC, PARAMETER ::   jp_qlw  = 6   ! index of Long wave                       (W/m2) 
    83    INTEGER , PUBLIC, PARAMETER ::   jp_prec = 7   ! index of total precipitation (rain+snow) (Kg/m2/s) 
    84    INTEGER , PUBLIC, PARAMETER ::   jp_snow = 8   ! index of snow (solid prcipitation)       (kg/m2/s) 
    85    INTEGER , PUBLIC, PARAMETER ::   jp_slp  = 9   ! index of sea level pressure              (Pa) 
    86    INTEGER , PUBLIC, PARAMETER ::   jp_hpgi =10   ! index of ABL geostrophic wind or hpg (i-component) (m/s) at T-point 
    87    INTEGER , PUBLIC, PARAMETER ::   jp_hpgj =11   ! index of ABL geostrophic wind or hpg (j-component) (m/s) at T-point 
    88  
     76   INTEGER , PUBLIC, PARAMETER ::   jp_wndi  =  1   ! index of 10m wind velocity (i-component) (m/s)    at T-point 
     77   INTEGER , PUBLIC, PARAMETER ::   jp_wndj  =  2   ! index of 10m wind velocity (j-component) (m/s)    at T-point 
     78   INTEGER , PUBLIC, PARAMETER ::   jp_tair  =  3   ! index of 10m air temperature             (Kelvin) 
     79   INTEGER , PUBLIC, PARAMETER ::   jp_humi  =  4   ! index of specific humidity               ( % ) 
     80   INTEGER , PUBLIC, PARAMETER ::   jp_qsr   =  5   ! index of solar heat                      (W/m2) 
     81   INTEGER , PUBLIC, PARAMETER ::   jp_qlw   =  6   ! index of Long wave                       (W/m2) 
     82   INTEGER , PUBLIC, PARAMETER ::   jp_prec  =  7   ! index of total precipitation (rain+snow) (Kg/m2/s) 
     83   INTEGER , PUBLIC, PARAMETER ::   jp_snow  =  8   ! index of snow (solid prcipitation)       (kg/m2/s) 
     84   INTEGER , PUBLIC, PARAMETER ::   jp_slp   =  9   ! index of sea level pressure              (Pa) 
     85   INTEGER , PUBLIC, PARAMETER ::   jp_uoatm = 10   ! index of surface current (i-component) 
     86   !                                                !          seen by the atmospheric forcing (m/s) at T-point 
     87   INTEGER , PUBLIC, PARAMETER ::   jp_voatm = 11   ! index of surface current (j-component) 
     88   !                                                !          seen by the atmospheric forcing (m/s) at T-point 
     89   INTEGER , PUBLIC, PARAMETER ::   jp_hpgi  = 12   ! index of ABL geostrophic wind or hpg (i-component) (m/s) at T-point 
     90   INTEGER , PUBLIC, PARAMETER ::   jp_hpgj  = 13   ! index of ABL geostrophic wind or hpg (j-component) (m/s) at T-point 
     91   INTEGER , PUBLIC, PARAMETER ::   jpfld    = 13   ! maximum number of files to read 
     92 
     93   ! Warning: keep this structure allocatable for Agrif... 
    8994   TYPE(FLD), PUBLIC, ALLOCATABLE, DIMENSION(:) ::   sf   ! structure of input atmospheric fields (file informations, fields read) 
    9095 
     
    98103   LOGICAL  ::   ln_Cd_L15      ! ice-atm drag = F( ice concentration, atmospheric stability ) (Lupkes et al. JGR2015) 
    99104   ! 
     105   LOGICAL  ::   ln_crt_fbk     ! Add surface current feedback to the wind stress computation  (Renault et al. 2020) 
     106   REAL(wp) ::   rn_stau_a      ! Alpha and Beta coefficients of Renault et al. 2020, eq. 10: Stau = Alpha * Wnd + Beta 
     107   REAL(wp) ::   rn_stau_b      !  
     108   ! 
    100109   REAL(wp)         ::   rn_pfac   ! multiplication factor for precipitation 
    101110   REAL(wp), PUBLIC ::   rn_efac   ! multiplication factor for evaporation 
    102    REAL(wp), PUBLIC ::   rn_vfac   ! multiplication factor for ice/ocean velocity in the calculation of wind stress 
    103111   REAL(wp)         ::   rn_zqt    ! z(q,t) : height of humidity and temperature measurements 
    104112   REAL(wp)         ::   rn_zu     ! z(u)   : height of wind measurements 
     
    162170      !! 
    163171      CHARACTER(len=100)            ::   cn_dir                ! Root directory for location of atmospheric forcing files 
    164       TYPE(FLD_N), ALLOCATABLE, DIMENSION(:) ::   slf_i        ! array of namelist informations on the fields to read 
    165       TYPE(FLD_N) ::   sn_wndi, sn_wndj, sn_humi, sn_qsr       ! informations about the fields to be read 
    166       TYPE(FLD_N) ::   sn_qlw , sn_tair, sn_prec, sn_snow      !       "                        " 
    167       TYPE(FLD_N) ::   sn_slp , sn_hpgi, sn_hpgj               !       "                        " 
     172      TYPE(FLD_N), DIMENSION(jpfld) ::   slf_i                 ! array of namelist informations on the fields to read 
     173      TYPE(FLD_N) ::   sn_wndi, sn_wndj , sn_humi, sn_qsr      ! informations about the fields to be read 
     174      TYPE(FLD_N) ::   sn_qlw , sn_tair , sn_prec, sn_snow     !       "                        " 
     175      TYPE(FLD_N) ::   sn_slp , sn_uoatm, sn_voatm             !       "                        " 
     176      TYPE(FLD_N) ::   sn_hpgi, sn_hpgj                        !       "                        " 
     177      INTEGER     ::   ipka                                    ! number of levels in the atmospheric variable 
    168178      NAMELIST/namsbc_blk/ sn_wndi, sn_wndj, sn_humi, sn_qsr, sn_qlw ,                &   ! input fields 
    169          &                 sn_tair, sn_prec, sn_snow, sn_slp, sn_hpgi, sn_hpgj,       & 
     179         &                 sn_tair, sn_prec, sn_snow, sn_slp, sn_uoatm, sn_voatm,     & 
     180         &                 sn_hpgi, sn_hpgj,                                          & 
    170181         &                 ln_NCAR, ln_COARE_3p0, ln_COARE_3p6, ln_ECMWF,             &   ! bulk algorithm 
    171182         &                 cn_dir , rn_zqt, rn_zu,                                    & 
    172          &                 rn_pfac, rn_efac, rn_vfac, ln_Cd_L12, ln_Cd_L15,           & 
     183         &                 rn_pfac, rn_efac, ln_Cd_L12, ln_Cd_L15,                    & 
     184         &                 ln_crt_fbk, rn_stau_a, rn_stau_b,                          &   ! current feedback 
    173185         &                 ln_skin_cs, ln_skin_wl, ln_humi_sph, ln_humi_dpt, ln_humi_rlh  ! cool-skin / warm-layer !LB 
    174186      !!--------------------------------------------------------------------- 
     
    242254      !                                   !* set the bulk structure 
    243255      !                                      !- store namelist information in an array 
    244       IF( ln_blk ) jpfld = 9 
    245       IF( ln_abl ) jpfld = 11 
    246       ALLOCATE( slf_i(jpfld) ) 
    247       ! 
    248       slf_i(jp_wndi) = sn_wndi   ;   slf_i(jp_wndj) = sn_wndj 
    249       slf_i(jp_qsr ) = sn_qsr    ;   slf_i(jp_qlw ) = sn_qlw 
    250       slf_i(jp_tair) = sn_tair   ;   slf_i(jp_humi) = sn_humi 
    251       slf_i(jp_prec) = sn_prec   ;   slf_i(jp_snow) = sn_snow 
    252       slf_i(jp_slp ) = sn_slp 
    253       IF( ln_abl ) THEN 
    254          slf_i(jp_hpgi) = sn_hpgi   ;   slf_i(jp_hpgj) = sn_hpgj 
    255       END IF 
     256      ! 
     257      slf_i(jp_wndi ) = sn_wndi    ;   slf_i(jp_wndj ) = sn_wndj 
     258      slf_i(jp_qsr  ) = sn_qsr     ;   slf_i(jp_qlw  ) = sn_qlw 
     259      slf_i(jp_tair ) = sn_tair    ;   slf_i(jp_humi ) = sn_humi 
     260      slf_i(jp_prec ) = sn_prec    ;   slf_i(jp_snow ) = sn_snow 
     261      slf_i(jp_slp  ) = sn_slp 
     262      slf_i(jp_uoatm) = sn_uoatm   ;   slf_i(jp_voatm) = sn_voatm 
     263      slf_i(jp_hpgi ) = sn_hpgi    ;   slf_i(jp_hpgj ) = sn_hpgj 
     264      ! 
     265      IF( .NOT. ln_abl ) THEN   ! force to not use jp_hpgi and jp_hpgj, should already be done in namelist_* but we never know... 
     266         slf_i(jp_hpgi)%clname = 'NOT USED' 
     267         slf_i(jp_hpgj)%clname = 'NOT USED' 
     268      ENDIF 
    256269      ! 
    257270      !                                      !- allocate the bulk structure 
     
    264277      DO jfpr= 1, jpfld 
    265278         ! 
    266          IF( TRIM(sf(jfpr)%clrootname) == 'NOT USED' ) THEN    !--  not used field  --!   (only now allocated and set to zero) 
    267             ALLOCATE( sf(jfpr)%fnow(jpi,jpj,1) ) 
    268             sf(jfpr)%fnow(:,:,1) = 0._wp 
     279         IF(   ln_abl    .AND.                                                      & 
     280            &    ( jfpr == jp_wndi .OR. jfpr == jp_wndj .OR. jfpr == jp_humi .OR.   & 
     281            &      jfpr == jp_hpgi .OR. jfpr == jp_hpgj .OR. jfpr == jp_tair     )  ) THEN 
     282            ipka = jpka   ! ABL: some fields are 3D input 
     283         ELSE 
     284            ipka = 1 
     285         ENDIF 
     286         ! 
     287         ALLOCATE( sf(jfpr)%fnow(jpi,jpj,ipka) ) 
     288         ! 
     289         IF( TRIM(sf(jfpr)%clrootname) == 'NOT USED' ) THEN    !--  not used field  --!   (only now allocated and set to default) 
     290            IF(     jfpr == jp_slp  ) THEN 
     291               sf(jfpr)%fnow(:,:,1:ipka) = 101325._wp   ! use standard pressure in Pa 
     292            ELSEIF( jfpr == jp_prec .OR. jfpr == jp_snow .OR. jfpr == jp_uoatm .OR. jfpr == jp_voatm ) THEN 
     293               sf(jfpr)%fnow(:,:,1:ipka) = 0._wp        ! no precip or no snow or no surface currents 
     294            ELSEIF( ( jfpr == jp_hpgi .OR. jfpr == jp_hpgj ) .AND. .NOT. ln_abl ) THEN 
     295               DEALLOCATE( sf(jfpr)%fnow )              ! deallocate as not used in this case 
     296            ELSE 
     297               WRITE(ctmp1,*) 'sbc_blk_init: no default value defined for field number', jfpr 
     298               CALL ctl_stop( ctmp1 ) 
     299            ENDIF 
    269300         ELSE                                                  !-- used field  --! 
    270             IF(   ln_abl    .AND.                                                      & 
    271                &    ( jfpr == jp_wndi .OR. jfpr == jp_wndj .OR. jfpr == jp_humi .OR.   & 
    272                &      jfpr == jp_hpgi .OR. jfpr == jp_hpgj .OR. jfpr == jp_tair     )  ) THEN   ! ABL: some fields are 3D input 
    273                ALLOCATE( sf(jfpr)%fnow(jpi,jpj,jpka) ) 
    274                IF( sf(jfpr)%ln_tint )   ALLOCATE( sf(jfpr)%fdta(jpi,jpj,jpka,2) ) 
    275             ELSE                                                                                ! others or Bulk fields are 2D fiels 
    276                ALLOCATE( sf(jfpr)%fnow(jpi,jpj,1) ) 
    277                IF( sf(jfpr)%ln_tint )   ALLOCATE( sf(jfpr)%fdta(jpi,jpj,1,2) ) 
    278             ENDIF 
     301            IF( sf(jfpr)%ln_tint )   ALLOCATE( sf(jfpr)%fdta(jpi,jpj,ipka,2) )   ! allocate array for temporal interpolation 
    279302            ! 
    280303            IF( sf(jfpr)%freqh > 0. .AND. MOD( NINT(3600. * sf(jfpr)%freqh), nn_fsbc * NINT(rn_Dt) ) /= 0 )   & 
     
    327350         WRITE(numout,*) '      factor applied on precipitation (total & snow)      rn_pfac      = ', rn_pfac 
    328351         WRITE(numout,*) '      factor applied on evaporation                       rn_efac      = ', rn_efac 
    329          WRITE(numout,*) '      factor applied on ocean/ice velocity                rn_vfac      = ', rn_vfac 
    330352         WRITE(numout,*) '         (form absolute (=0) to relative winds(=1))' 
    331353         WRITE(numout,*) '      use ice-atm drag from Lupkes2012                    ln_Cd_L12    = ', ln_Cd_L12 
    332354         WRITE(numout,*) '      use ice-atm drag from Lupkes2015                    ln_Cd_L15    = ', ln_Cd_L15 
     355         WRITE(numout,*) '      use surface current feedback on wind stress         ln_crt_fbk   = ', ln_crt_fbk 
     356         IF(ln_crt_fbk) THEN 
     357         WRITE(numout,*) '         Renault et al. 2020, eq. 10: Stau = Alpha * Wnd + Beta' 
     358         WRITE(numout,*) '            Alpha                                         rn_stau_a    = ', rn_stau_a 
     359         WRITE(numout,*) '            Beta                                          rn_stau_b    = ', rn_stau_b 
     360         ENDIF 
    333361         ! 
    334362         WRITE(numout,*) 
     
    429457      !                                            ! compute the surface ocean fluxes using bulk formulea 
    430458      IF( MOD( kt - 1, nn_fsbc ) == 0 ) THEN 
    431          CALL blk_oce_1( kt, sf(jp_wndi)%fnow(:,:,1), sf(jp_wndj)%fnow(:,:,1),   &   !   <<= in 
    432             &                sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1),   &   !   <<= in 
    433             &                sf(jp_slp )%fnow(:,:,1), sst_m, ssu_m, ssv_m,       &   !   <<= in 
    434             &                sf(jp_qsr )%fnow(:,:,1), sf(jp_qlw )%fnow(:,:,1),   &   !   <<= in (wl/cs) 
    435             &                tsk_m, zssq, zcd_du, zsen, zevp )                       !   =>> out 
    436  
    437          CALL blk_oce_2(     sf(jp_tair)%fnow(:,:,1), sf(jp_qsr )%fnow(:,:,1),   &   !   <<= in 
    438             &                sf(jp_qlw )%fnow(:,:,1), sf(jp_prec)%fnow(:,:,1),   &   !   <<= in 
    439             &                sf(jp_snow)%fnow(:,:,1), tsk_m,                     &   !   <<= in 
    440             &                zsen, zevp )                                            !   <=> in out 
     459         CALL blk_oce_1( kt, sf(jp_wndi )%fnow(:,:,1), sf(jp_wndj )%fnow(:,:,1),   &   !   <<= in 
     460            &                sf(jp_tair )%fnow(:,:,1), sf(jp_humi )%fnow(:,:,1),   &   !   <<= in 
     461            &                sf(jp_slp  )%fnow(:,:,1), sst_m, ssu_m, ssv_m,        &   !   <<= in 
     462            &                sf(jp_uoatm)%fnow(:,:,1), sf(jp_voatm)%fnow(:,:,1),   &   !   <<= in 
     463            &                sf(jp_qsr  )%fnow(:,:,1), sf(jp_qlw  )%fnow(:,:,1),   &   !   <<= in (wl/cs) 
     464            &                tsk_m, zssq, zcd_du, zsen, zevp )                         !   =>> out 
     465 
     466         CALL blk_oce_2(     sf(jp_tair )%fnow(:,:,1), sf(jp_qsr  )%fnow(:,:,1),   &   !   <<= in 
     467            &                sf(jp_qlw  )%fnow(:,:,1), sf(jp_prec )%fnow(:,:,1),   &   !   <<= in 
     468            &                sf(jp_snow )%fnow(:,:,1), tsk_m,                      &   !   <<= in 
     469            &                zsen, zevp )                                              !   <=> in out 
    441470      ENDIF 
    442471      ! 
     
    470499 
    471500 
    472    SUBROUTINE blk_oce_1( kt, pwndi, pwndj , ptair, phumi, &  ! inp 
    473       &                  pslp , pst   , pu   , pv,        &  ! inp 
    474       &                  pqsr , pqlw  ,                   &  ! inp 
    475       &                  ptsk, pssq , pcd_du, psen , pevp   )  ! out 
     501   SUBROUTINE blk_oce_1( kt, pwndi, pwndj, ptair, phumi,        &  ! inp 
     502      &                      pslp , pst  , pu   , pv,            &  ! inp 
     503      &                      puatm, pvatm, pqsr , pqlw ,         &  ! inp 
     504      &                      ptsk , pssq , pcd_du, psen, pevp   )   ! out 
    476505      !!--------------------------------------------------------------------- 
    477506      !!                     ***  ROUTINE blk_oce_1  *** 
     
    498527      REAL(wp), INTENT(in   ), DIMENSION(:,:) ::   pu     ! surface current at U-point (i-component) [m/s] 
    499528      REAL(wp), INTENT(in   ), DIMENSION(:,:) ::   pv     ! surface current at V-point (j-component) [m/s] 
     529      REAL(wp), INTENT(in   ), DIMENSION(:,:) ::   puatm  ! surface current seen by the atm at T-point (i-component) [m/s] 
     530      REAL(wp), INTENT(in   ), DIMENSION(:,:) ::   pvatm  ! surface current seen by the atm at T-point (j-component) [m/s] 
    500531      REAL(wp), INTENT(in   ), DIMENSION(:,:) ::   pqsr   ! 
    501532      REAL(wp), INTENT(in   ), DIMENSION(:,:) ::   pqlw   ! 
     
    508539      INTEGER  ::   ji, jj               ! dummy loop indices 
    509540      REAL(wp) ::   zztmp                ! local variable 
     541      REAL(wp) ::   zstmax, zstau 
     542#if defined key_cyclone 
    510543      REAL(wp), DIMENSION(jpi,jpj) ::   zwnd_i, zwnd_j    ! wind speed components at T-point 
     544#endif 
     545      REAL(wp), DIMENSION(jpi,jpj) ::   ztau_i, ztau_j    ! wind stress components at T-point 
    511546      REAL(wp), DIMENSION(jpi,jpj) ::   zU_zu             ! bulk wind speed at height zu  [m/s] 
    512547      REAL(wp), DIMENSION(jpi,jpj) ::   ztpot             ! potential temperature of air at z=rn_zqt [K] 
     
    532567      zwnd_j(:,:) = 0._wp 
    533568      CALL wnd_cyc( kt, zwnd_i, zwnd_j )    ! add analytical tropical cyclone (Vincent et al. JGR 2012) 
    534       DO_2D_00_00 
    535          pwndi(ji,jj) = pwndi(ji,jj) + zwnd_i(ji,jj) 
    536          pwndj(ji,jj) = pwndj(ji,jj) + zwnd_j(ji,jj) 
     569      DO_2D_11_11 
     570         zwnd_i(ji,jj) = pwndi(ji,jj) + zwnd_i(ji,jj) 
     571         zwnd_j(ji,jj) = pwndj(ji,jj) + zwnd_j(ji,jj) 
     572         ! ... scalar wind at T-point (not masked) 
     573         wndm(ji,jj) = SQRT( zwnd_i(ji,jj) * zwnd_i(ji,jj) + zwnd_j(ji,jj) * zwnd_j(ji,jj) ) 
     574      END_2D 
     575#else 
     576      ! ... scalar wind module at T-point (not masked) 
     577      DO_2D_11_11 
     578         wndm(ji,jj) = SQRT(  pwndi(ji,jj) * pwndi(ji,jj) + pwndj(ji,jj) * pwndj(ji,jj)  ) 
    537579      END_2D 
    538580#endif 
    539       DO_2D_00_00 
    540          zwnd_i(ji,jj) = (  pwndi(ji,jj) - rn_vfac * 0.5 * ( pu(ji-1,jj  ) + pu(ji,jj) )  ) 
    541          zwnd_j(ji,jj) = (  pwndj(ji,jj) - rn_vfac * 0.5 * ( pv(ji  ,jj-1) + pv(ji,jj) )  ) 
    542       END_2D 
    543       CALL lbc_lnk_multi( 'sbcblk', zwnd_i, 'T', -1., zwnd_j, 'T', -1. ) 
    544       ! ... scalar wind ( = | U10m - U_oce | ) at T-point (masked) 
    545       wndm(:,:) = SQRT(  zwnd_i(:,:) * zwnd_i(:,:)   & 
    546          &             + zwnd_j(:,:) * zwnd_j(:,:)  ) * tmask(:,:,1) 
    547  
    548581      ! ----------------------------------------------------------------------------- ! 
    549582      !      I   Solar FLUX                                                           ! 
     
    597630         END_2D 
    598631      ENDIF 
    599  
    600  
    601632 
    602633      !! Time to call the user-selected bulk parameterization for 
     
    673704         pevp(:,:) = pevp(:,:) * tmask(:,:,1) 
    674705 
    675          ! Tau i and j component on T-grid points, using array "zcd_oce" as a temporary array... 
    676          zcd_oce = 0._wp 
    677          WHERE ( wndm > 0._wp ) zcd_oce = taum / wndm 
    678          zwnd_i = zcd_oce * zwnd_i 
    679          zwnd_j = zcd_oce * zwnd_j 
    680  
    681          CALL iom_put( "taum_oce", taum )   ! output wind stress module 
     706         DO_2D_11_11 
     707            IF( wndm(ji,jj) > 0._wp ) THEN 
     708               zztmp = taum(ji,jj) / wndm(ji,jj) 
     709#if defined key_cyclone 
     710               ztau_i(ji,jj) = zztmp * zwnd_i(ji,jj) 
     711               ztau_j(ji,jj) = zztmp * zwnd_j(ji,jj) 
     712#else 
     713               ztau_i(ji,jj) = zztmp * pwndi(ji,jj) 
     714               ztau_j(ji,jj) = zztmp * pwndj(ji,jj) 
     715#endif 
     716            ELSE 
     717               ztau_i(ji,jj) = 0._wp 
     718               ztau_j(ji,jj) = 0._wp                  
     719            ENDIF 
     720         END_2D 
     721 
     722         IF( ln_crt_fbk ) THEN   ! aply eq. 10 and 11 of Renault et al. 2020 (doi: 10.1029/2019MS001715) 
     723            zstmax = MIN( rn_stau_a * 3._wp + rn_stau_b, 0._wp )   ! set the max value of Stau corresponding to a wind of 3 m/s (<0) 
     724            DO_2D_01_01   ! end at jpj and jpi, as ztau_j(ji,jj+1) ztau_i(ji+1,jj) used in the next loop 
     725               zstau = MIN( rn_stau_a * wndm(ji,jj) + rn_stau_b, zstmax )   ! stau (<0) must be smaller than zstmax 
     726               ztau_i(ji,jj) = ztau_i(ji,jj) + zstau * ( 0.5_wp * ( pu(ji-1,jj  ) + pu(ji,jj) ) - puatm(ji,jj) ) 
     727               ztau_j(ji,jj) = ztau_j(ji,jj) + zstau * ( 0.5_wp * ( pv(ji  ,jj-1) + pv(ji,jj) ) - pvatm(ji,jj) ) 
     728               taum(ji,jj) = SQRT( ztau_i(ji,jj) * ztau_i(ji,jj) + ztau_j(ji,jj) * ztau_j(ji,jj) ) 
     729            END_2D 
     730         ENDIF 
    682731 
    683732         ! ... utau, vtau at U- and V_points, resp. 
    684733         !     Note the use of 0.5*(2-umask) in order to unmask the stress along coastlines 
    685734         !     Note the use of MAX(tmask(i,j),tmask(i+1,j) is to mask tau over ice shelves 
    686          DO_2D_10_10 
    687             utau(ji,jj) = 0.5 * ( 2. - umask(ji,jj,1) ) * ( zwnd_i(ji,jj) + zwnd_i(ji+1,jj  ) ) & 
    688                &          * MAX(tmask(ji,jj,1),tmask(ji+1,jj,1)) 
    689             vtau(ji,jj) = 0.5 * ( 2. - vmask(ji,jj,1) ) * ( zwnd_j(ji,jj) + zwnd_j(ji  ,jj+1) ) & 
    690                &          * MAX(tmask(ji,jj,1),tmask(ji,jj+1,1)) 
     735         DO_2D_00_00              ! start loop at 2, in case ln_crt_fbk = T 
     736            utau(ji,jj) = 0.5 * ( 2. - umask(ji,jj,1) ) * ( ztau_i(ji,jj) + ztau_i(ji+1,jj  ) ) & 
     737               &              * MAX(tmask(ji,jj,1),tmask(ji+1,jj,1)) 
     738            vtau(ji,jj) = 0.5 * ( 2. - vmask(ji,jj,1) ) * ( ztau_j(ji,jj) + ztau_j(ji  ,jj+1) ) & 
     739               &              * MAX(tmask(ji,jj,1),tmask(ji,jj+1,1)) 
    691740         END_2D 
    692          CALL lbc_lnk_multi( 'sbcblk', utau, 'U', -1., vtau, 'V', -1. ) 
     741 
     742         IF( ln_crt_fbk ) THEN 
     743            CALL lbc_lnk_multi( 'sbcblk', utau, 'U', -1., vtau, 'V', -1., taum, 'T', -1. ) 
     744         ELSE 
     745            CALL lbc_lnk_multi( 'sbcblk', utau, 'U', -1., vtau, 'V', -1. ) 
     746         ENDIF 
     747 
     748         CALL iom_put( "taum_oce", taum )   ! output wind stress module 
    693749 
    694750         IF(sn_cfctl%l_prtctl) THEN 
     
    861917      ! 
    862918      INTEGER  ::   ji, jj    ! dummy loop indices 
    863       REAL(wp) ::   zwndi_t , zwndj_t             ! relative wind components at T-point 
    864919      REAL(wp) ::   zootm_su                      ! sea-ice surface mean temperature 
    865920      REAL(wp) ::   zztmp1, zztmp2                ! temporary arrays 
     
    872927      ! ------------------------------------------------------------ ! 
    873928      ! C-grid ice dynamics :   U & V-points (same as ocean) 
    874       DO_2D_00_00 
    875          zwndi_t = (  pwndi(ji,jj) - rn_vfac * 0.5_wp * ( puice(ji-1,jj  ) + puice(ji,jj) )  ) 
    876          zwndj_t = (  pwndj(ji,jj) - rn_vfac * 0.5_wp * ( pvice(ji  ,jj-1) + pvice(ji,jj) )  ) 
    877          wndm_ice(ji,jj) = SQRT( zwndi_t * zwndi_t + zwndj_t * zwndj_t ) * tmask(ji,jj,1) 
     929      DO_2D_11_11 
     930         wndm_ice(ji,jj) = SQRT( pwndi(ji,jj) * pwndi(ji,jj) + pwndj(ji,jj) * pwndj(ji,jj) ) 
    878931      END_2D 
    879       CALL lbc_lnk( 'sbcblk', wndm_ice, 'T',  1. ) 
    880932      ! 
    881933      ! Make ice-atm. drag dependent on ice concentration 
     
    902954         ! C-grid ice dynamics :   U & V-points (same as ocean) 
    903955         DO_2D_00_00 
    904             putaui(ji,jj) = 0.5_wp * (  rhoa(ji+1,jj) * zcd_dui(ji+1,jj)             & 
    905                &                      + rhoa(ji  ,jj) * zcd_dui(ji  ,jj)  )          & 
    906                &         * ( 0.5_wp * ( pwndi(ji+1,jj) + pwndi(ji,jj) ) - rn_vfac * puice(ji,jj) ) 
    907             pvtaui(ji,jj) = 0.5_wp * (  rhoa(ji,jj+1) * zcd_dui(ji,jj+1)             & 
    908                &                      + rhoa(ji,jj  ) * zcd_dui(ji,jj  )  )          & 
    909                &         * ( 0.5_wp * ( pwndj(ji,jj+1) + pwndj(ji,jj) ) - rn_vfac * pvice(ji,jj) ) 
     956            putaui(ji,jj) = 0.25_wp * ( rhoa(ji+1,jj) * zcd_dui(ji+1,jj) + rhoa(ji,jj) * zcd_dui(ji,jj) )   & 
     957               &                    * (                   pwndi(ji+1,jj) +                 pwndi(ji,jj) ) 
     958            pvtaui(ji,jj) = 0.25_wp * ( rhoa(ji,jj+1) * zcd_dui(ji,jj+1) + rhoa(ji,jj) * zcd_dui(ji,jj) )   & 
     959               &                    * (                   pwndj(ji,jj+1) +                 pwndj(ji,jj) ) 
    910960         END_2D 
    911961         CALL lbc_lnk_multi( 'sbcblk', putaui, 'U', -1., pvtaui, 'V', -1. ) 
Note: See TracChangeset for help on using the changeset viewer.