Changeset 12567


Ignore:
Timestamp:
2020-03-18T11:37:00+01:00 (7 months ago)
Author:
jcastill
Message:

Review the previous changes and make sure that wind and momentum components and module are in the right grid; reduce the number of loops to the minimum

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/UKMO/AMM15_v3_6_STABLE_package_collate_coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbcflx.F90

    r12470 r12567  
    135135         slf_i(jp_emp ) = sn_emp 
    136136         ! 
    137          ALLOCATE( sf(jpfld), STAT=ierror )        ! set sf structure 
    138  
    139137         ! define local jpfld depending on shelf_flx logical 
    140138         IF( ln_shelf_flx ) THEN 
     
    145143         ENDIF 
    146144         ! 
     145         ALLOCATE( sf(jpfld), STAT=ierror )        ! set sf structure 
     146         ! 
    147147         IF( ierror > 0 ) THEN    
    148148            CALL ctl_stop( 'sbc_flx: unable to allocate sf structure' )   ;   RETURN   
     
    166166         ELSE                  ;   qsr(:,:) =          sf(jp_qsr)%fnow(:,:,1) 
    167167         ENDIF 
    168 !CDIR COLLAPSE 
    169          !!UKMO SHELF read wind instead of momentum components (in the U,V grid)  
     168 
     169         !!UKMO SHELF wind speed relative to surface currents - put here to allow merging with coupling branch 
    170170         IF( ln_shelf_flx ) THEN 
    171171            CALL wrk_alloc( jpi,jpj, zwnd_i, zwnd_j ) 
    172172 
    173             ! wind speed relative to surface currents 
    174             IF( ln_rel_wind ) THEN 
    175                DO jj = 1, jpj 
    176                   DO ji = 1, jpi 
     173            ! set the ocean fluxes from read fields 
     174            DO jj = 1, jpj 
     175               DO ji = 1, jpi 
     176                  !! UKMO SHELF - need atmospheric pressure to calculate Haney forcing 
     177                  pressnow(ji,jj) = sf(jp_press)%fnow(ji,jj,1) 
     178                  !! UKMO SHELF flux files contain wind speed not wind stress 
     179                  IF( ln_rel_wind ) THEN 
    177180                     zwnd_i(ji,jj) = sf(jp_utau)%fnow(ji,jj,1) - rn_wfac * ssu_m(ji,jj) 
    178181                     zwnd_j(ji,jj) = sf(jp_vtau)%fnow(ji,jj,1) - rn_wfac * ssv_m(ji,jj) 
    179                   END DO 
    180                END DO 
    181             ELSE 
    182                zwnd_i(:,:) = sf(jp_utau)%fnow(:,:,1) 
    183                zwnd_j(:,:) = sf(jp_vtau)%fnow(:,:,1) 
    184             ENDIF 
    185  
    186             ! calculate the wind speed module 
    187             DO jj = 2, jpjm1  
    188                DO ji = fs_2, fs_jpim1    ! vect. opt.  
    189                   ztx = zwnd_i(ji-1,jj  ) + zwnd_i(ji,jj)  
    190                   zty = zwnd_j(ji  ,jj-1) + zwnd_j(ji,jj)  
    191                   wndm(ji,jj) = 0.5 * SQRT( ztx * ztx + zty * zty )  
    192                END DO  
    193             END DO  
    194             CALL lbc_lnk( wndm(:,:), 'T', 1. )  
    195  
    196             qsr(:,:) = sf(jp_qsr )%fnow(:,:,1) 
    197             !! UKMO FOAM flux files contain non-solar heat flux (qns) rather than total heat flux (qtot) 
    198             qns(:,:) = sf(jp_qtot)%fnow(:,:,1) 
    199             !! UKMO FOAM flux files contain the net DOWNWARD freshwater flux P-E rather then E-P 
    200             emp(:,:) = -1. * sf(jp_emp )%fnow(:,:,1) 
    201             ! need atmospheric pressure to calculate Haney forcing 
    202             pressnow(:,:) = sf(jp_press)%fnow(:,:,1) 
    203  
    204             !!Effect of atmospheric pressure on SSH 
    205             ! If using ln_apr_dyn, this is done there so don't repeat here. 
    206             IF( .NOT. ln_apr_dyn ) THEN 
    207                DO jj = 1, jpjm1 
    208                   DO ji = 1, jpim1 
    209                      apgu(ji,jj) = (-1.0/rau0)*(pressnow(ji+1,jj)-pressnow(ji,jj))/e1u(ji,jj) 
    210                      apgv(ji,jj) = (-1.0/rau0)*(pressnow(ji,jj+1)-pressnow(ji,jj))/e2v(ji,jj) 
    211                   END DO 
    212                END DO 
    213             ENDIF 
    214  
    215             ! Calculate momentum from wind components (several options)  
    216             IF( ln_cdgw .AND. nn_drag == jp_std ) THEN  
    217                IF( cpl_wdrag ) THEN   
    218                   ! reset utau and vtau to the wind components: the momentum will  
    219                   ! be calculated from the coupled value of the drag coefficient  
    220                   DO jj = 1, jpj  
    221                      DO ji = 1, jpi  
     182                  ELSE 
     183                     zwnd_i(ji,jj) = sf(jp_utau)%fnow(ji,jj,1) 
     184                     zwnd_j(ji,jj) = sf(jp_vtau)%fnow(ji,jj,1) 
     185                  ENDIF 
     186                  wndm(ji,jj) = sqrt(zwnd_i(ji,jj)*zwnd_i(ji,jj) + zwnd_j(ji,jj)*zwnd_j(ji,jj)) 
     187 
     188                  ! add modification due to drag coefficient read from wave forcing  
     189                  IF( ln_cdgw .AND. nn_drag == jp_std ) THEN  
     190                     IF( cpl_wdrag ) THEN   
     191                        ! reset utau and vtau to the wind components: the momentum will  
     192                        ! be calculated from the coupled value of the drag coefficient  
    222193                        utau(ji,jj) = zwnd_i(ji,jj)  
    223194                        vtau(ji,jj) = zwnd_j(ji,jj)  
    224                      END DO  
    225                   END DO  
    226                ELSE  
    227                   DO jj = 1, jpjm1  
    228                      DO ji = 1, jpim1  
    229                         utau(ji,jj) = zrhoa * 0.5 * ( cdn_wave(ji,jj) + cdn_wave(ji+1,jj) ) * zwnd_i(ji,jj) * &  
    230                                                                         0.5 * ( wndm(ji,jj) + wndm(ji+1,jj) )  
    231                         vtau(ji,jj) = zrhoa * 0.5 * ( cdn_wave(ji,jj) + cdn_wave(ji,jj+1) ) * zwnd_j(ji,jj) * &  
    232                                                                         0.5 * ( wndm(ji,jj) + wndm(ji,jj+1) )  
    233                      END DO  
    234                   END DO  
    235                   CALL lbc_lnk_multi( utau(:,:), 'U', -1., vtau(:,:), 'V', -1. )  
    236                ENDIF  
    237             ELSE IF( nn_drag == jp_const ) THEN  
    238                DO jj = 1, jpjm1  
    239                   DO ji = 1, jpim1  
    240                      utau(ji,jj) = zrhoa * zcdrag * zwnd_i(ji,jj) * 0.5 * ( wndm(ji,jj) + wndm(ji+1,jj) )  
    241                      vtau(ji,jj) = zrhoa * zcdrag * zwnd_j(ji,jj) * 0.5 * ( wndm(ji,jj) + wndm(ji,jj+1) )  
    242                   END DO  
    243                END DO  
    244                CALL lbc_lnk_multi( utau(:,:), 'U', -1., vtau(:,:), 'V', -1. )  
    245             ELSE IF( nn_drag == jp_ukmo ) THEN  
    246                DO jj = 1, jpj 
    247                   DO ji = 1, jpi 
    248                       cs = 0.63 + (0.066 * wndm(ji,jj)) 
    249                       utau(ji,jj) = cs * (rhoa/rau0) * zwnd_i(ji,jj) * 0.5 * ( wndm(ji,jj) + wndm(ji+1,jj) ) 
    250                       vtau(ji,jj) = cs * (rhoa/rau0) * zwnd_j(ji,jj) * 0.5 * ( wndm(ji,jj) + wndm(ji,jj+1) ) 
    251                   END DO 
    252                END DO 
    253             ENDIF  
     195                     ELSE  
     196                        utau(ji,jj) = zrhoa * cdn_wave(ji,jj) * zwnd_i(ji,jj) * wndm(ji,jj) 
     197                        vtau(ji,jj) = zrhoa * cdn_wave(ji,jj) * zwnd_j(ji,jj) * wndm(ji,jj) 
     198                     ENDIF  
     199                  ELSE IF( nn_drag == jp_const ) THEN  
     200                     utau(ji,jj) = zrhoa * zcdrag * zwnd_i(ji,jj) * wndm(ji,jj) 
     201                     vtau(ji,jj) = zrhoa * zcdrag * zwnd_j(ji,jj) * wndm(ji,jj) 
     202                  ELSE IF( nn_drag == jp_ukmo ) THEN  
     203                     cs = 0.63 + (0.066 * wndm(ji,jj)) 
     204                     utau(ji,jj) = cs * (rhoa/rau0) * zwnd_i(ji,jj) * wndm(ji,jj) 
     205                     vtau(ji,jj) = cs * (rhoa/rau0) * zwnd_j(ji,jj) * wndm(ji,jj) 
     206                  ENDIF  
     207                  taum(ji,jj) = sqrt(utau(ji,jj)*utau(ji,jj) + vtau(ji,jj)*vtau(ji,jj)) 
     208 
     209                  qsr (ji,jj) = sf(jp_qsr )%fnow(ji,jj,1) 
     210                  !! UKMO FOAM flux files contain non-solar heat flux (qns) rather than total heat flux (qtot) 
     211                  qns (ji,jj) = sf(jp_qtot)%fnow(ji,jj,1) 
     212                  !! UKMO FOAM flux files contain the net DOWNWARD freshwater flux P-E rather then E-P 
     213                  emp (ji,jj) = -1. * sf(jp_emp )%fnow(ji,jj,1) 
     214               END DO 
     215            END DO 
     216            ! Move tau components to the right grid 
     217            DO jj = 1, jpjm1 
     218               DO ji = 1, jpim1 
     219                  utau(ji,jj) = 0.5 * ( utau(ji,jj) + utau(ji+1,jj) ) 
     220                  vtau(ji,jj) = 0.5 * ( vtau(ji,jj) + vtau(ji,jj+1) ) 
     221 
     222                  !!UKMO SHELF effect of atmospheric pressure on SSH 
     223                  ! If using ln_apr_dyn, this is done there so don't repeat here. 
     224                  IF( .NOT. ln_apr_dyn) THEN 
     225                     apgu(ji,jj) = (-1.0/rau0)*(sf(jp_press)%fnow(ji+1,jj,1)-sf(jp_press)%fnow(ji,jj,1))/e1u(ji,jj) 
     226                     apgv(ji,jj) = (-1.0/rau0)*(sf(jp_press)%fnow(ji,jj+1,1)-sf(jp_press)%fnow(ji,jj,1))/e2v(ji,jj) 
     227                  ENDIF 
     228               END DO 
     229            END DO 
    254230         ELSE 
    255231            DO jj = 1, jpj                                           ! set the ocean fluxes from read fields 
    256232               DO ji = 1, jpi 
    257                       utau(ji,jj) = sf(jp_utau)%fnow(ji,jj,1) 
    258                       vtau(ji,jj) = sf(jp_vtau)%fnow(ji,jj,1) 
    259  
    260                       qsr (ji,jj) = sf(jp_qsr )%fnow(ji,jj,1) 
    261                       IF( ln_foam_flx ) THEN 
    262                          !! UKMO FOAM flux files contain non-solar heat flux (qns) rather than total heat flux (qtot) 
    263                          qns (ji,jj) = sf(jp_qtot)%fnow(ji,jj,1) 
    264                          !! UKMO FOAM flux files contain the net DOWNWARD freshwater flux P-E rather then E-P 
    265                          emp (ji,jj) = -1. * sf(jp_emp )%fnow(ji,jj,1) 
    266                       ELSE 
    267                          qns (ji,jj) = sf(jp_qtot)%fnow(ji,jj,1) - sf(jp_qsr)%fnow(ji,jj,1) 
    268                          emp (ji,jj) = sf(jp_emp )%fnow(ji,jj,1) 
    269                       ENDIF 
    270                END DO 
    271             END DO 
    272          ENDIF ! ln_shelf_flx 
     233                  utau(ji,jj) = sf(jp_utau)%fnow(ji,jj,1) 
     234                  vtau(ji,jj) = sf(jp_vtau)%fnow(ji,jj,1) 
     235                  qsr (ji,jj) = sf(jp_qsr )%fnow(ji,jj,1) 
     236                  IF( ln_foam_flx ) THEN 
     237                     !! UKMO FOAM flux files contain non-solar heat flux (qns) rather than total heat flux (qtot) 
     238                     qns (ji,jj) = sf(jp_qtot)%fnow(ji,jj,1) 
     239                     !! UKMO FOAM flux files contain the net DOWNWARD freshwater flux P-E rather then E-P 
     240                     emp (ji,jj) = -1. * sf(jp_emp )%fnow(ji,jj,1) 
     241                  ELSE 
     242                     qns (ji,jj) = sf(jp_qtot)%fnow(ji,jj,1) - sf(jp_qsr)%fnow(ji,jj,1) 
     243                     emp (ji,jj) = sf(jp_emp )%fnow(ji,jj,1) 
     244                  ENDIF 
     245               END DO 
     246            END DO 
     247            !                                                        ! module of wind stress and wind speed at T-point 
     248            zcoef = 1. / ( zrhoa * zcdrag ) 
     249!CDIR NOVERRCHK 
     250            DO jj = 2, jpjm1 
     251!CDIR NOVERRCHK 
     252               DO ji = fs_2, fs_jpim1   ! vect. opt. 
     253                  ztx = utau(ji-1,jj  ) + utau(ji,jj)  
     254                  zty = vtau(ji  ,jj-1) + vtau(ji,jj)  
     255                  zmod = 0.5 * SQRT( ztx * ztx + zty * zty ) 
     256                  taum(ji,jj) = zmod 
     257                  wndm(ji,jj) = SQRT( zmod * zcoef ) 
     258               END DO 
     259            END DO 
     260         ENDIF 
    273261         !                                                        ! add to qns the heat due to e-p 
    274262         qns(:,:) = qns(:,:) - emp(:,:) * sst_m(:,:) * rcp        ! mass flux is at SST 
    275263         ! 
    276264         !! UKMO FOAM wind fluxes need lbc_lnk calls owing to a bug in interp.exe 
    277          IF( ln_foam_flx ) THEN 
    278             CALL lbc_lnk( utau(:,:), 'U', -1. ) 
    279             CALL lbc_lnk( vtau(:,:), 'V', -1. ) 
    280          ENDIF 
    281          !                                                        ! module of wind stress and wind speed at T-point 
    282          zcoef = 1. / ( zrhoa * zcdrag ) 
    283 !CDIR NOVERRCHK 
    284          DO jj = 2, jpjm1 
    285 !CDIR NOVERRCHK 
    286             DO ji = fs_2, fs_jpim1   ! vect. opt. 
    287                ztx = utau(ji-1,jj  ) + utau(ji,jj)  
    288                zty = vtau(ji  ,jj-1) + vtau(ji,jj)  
    289                zmod = 0.5 * SQRT( ztx * ztx + zty * zty ) 
    290                taum(ji,jj) = zmod 
    291                IF ( .NOT. ln_shelf_flx ) THEN 
    292                   wndm(ji,jj) = SQRT( zmod * zcoef ) 
    293                ENDIF 
    294             END DO 
    295          END DO 
     265         IF( ln_foam_flx .OR. ln_shelf_flx ) THEN 
     266            CALL lbc_lnk_multi( utau(:,:), 'U', -1., vtau(:,:), 'V', -1. )  
     267         ENDIF 
     268     
    296269         taum(:,:) = taum(:,:) * tmask(:,:,1) ; wndm(:,:) = wndm(:,:) * tmask(:,:,1) 
    297          CALL lbc_lnk( taum(:,:), 'T', 1. )  ;   CALL lbc_lnk( wndm(:,:), 'T', 1. ) 
     270         CALL lbc_lnk_multi( taum(:,:), 'T', 1., wndm(:,:), 'T', 1. ) 
    298271 
    299272         IF( nitend-nit000 <= 100 .AND. lwp ) THEN                ! control print (if less than 100 time-step asked) 
Note: See TracChangeset for help on using the changeset viewer.