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 12470 for branches/UKMO – NEMO

Changeset 12470 for branches/UKMO


Ignore:
Timestamp:
2020-02-26T13:23:10+01:00 (4 years ago)
Author:
jcastill
Message:

Correct the calculation of wind velocity for shelf flux direct forcing - clean up the shelf flux direct forcing code

File:
1 edited

Legend:

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

    r11277 r12470  
    135135         slf_i(jp_emp ) = sn_emp 
    136136         ! 
    137             ALLOCATE( sf(jpfld), STAT=ierror )        ! set sf structure 
    138             IF( ln_shelf_flx ) slf_i(jp_press) = sn_press 
    139     
    140             ! define local jpfld depending on shelf_flx logical 
    141             IF( ln_shelf_flx ) THEN 
    142                jpfld_local = jpfld 
    143             ELSE 
    144                jpfld_local = jpfld-1 
    145             ENDIF 
    146             ! 
     137         ALLOCATE( sf(jpfld), STAT=ierror )        ! set sf structure 
     138 
     139         ! define local jpfld depending on shelf_flx logical 
     140         IF( ln_shelf_flx ) THEN 
     141            slf_i(jp_press) = sn_press 
     142            jpfld_local = jpfld 
     143         ELSE 
     144            jpfld_local = jpfld-1 
     145         ENDIF 
     146         ! 
    147147         IF( ierror > 0 ) THEN    
    148148            CALL ctl_stop( 'sbc_flx: unable to allocate sf structure' )   ;   RETURN   
     
    163163      IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN                        ! update ocean fluxes at each SBC frequency 
    164164 
    165          !!UKMO SHELF wind speed relative to surface currents - put here to allow merging with coupling branch 
     165         IF( ln_dm2dc ) THEN   ;   qsr(:,:) = sbc_dcy( sf(jp_qsr)%fnow(:,:,1) )   ! modify now Qsr to include the diurnal cycle 
     166         ELSE                  ;   qsr(:,:) =          sf(jp_qsr)%fnow(:,:,1) 
     167         ENDIF 
     168!CDIR COLLAPSE 
     169         !!UKMO SHELF read wind instead of momentum components (in the U,V grid)  
    166170         IF( ln_shelf_flx ) THEN 
    167171            CALL wrk_alloc( jpi,jpj, zwnd_i, zwnd_j ) 
    168172 
     173            ! wind speed relative to surface currents 
    169174            IF( ln_rel_wind ) THEN 
    170175               DO jj = 1, jpj 
     
    178183               zwnd_j(:,:) = sf(jp_vtau)%fnow(:,:,1) 
    179184            ENDIF 
    180          ENDIF 
    181  
    182          IF( ln_dm2dc ) THEN   ;   qsr(:,:) = sbc_dcy( sf(jp_qsr)%fnow(:,:,1) )   ! modify now Qsr to include the diurnal cycle 
    183          ELSE                  ;   qsr(:,:) =          sf(jp_qsr)%fnow(:,:,1) 
    184          ENDIF 
    185 !CDIR COLLAPSE 
    186             !!UKMO SHELF effect of atmospheric pressure on SSH 
    187             ! If using ln_apr_dyn, this is done there so don't repeat here. 
    188             IF( ln_shelf_flx .AND. .NOT. ln_apr_dyn) THEN 
    189                DO jj = 1, jpjm1 
    190                   DO ji = 1, jpim1 
    191                      apgu(ji,jj) = (-1.0/rau0)*(sf(jp_press)%fnow(ji+1,jj,1)-sf(jp_press)%fnow(ji,jj,1))/e1u(ji,jj) 
    192                      apgv(ji,jj) = (-1.0/rau0)*(sf(jp_press)%fnow(ji,jj+1,1)-sf(jp_press)%fnow(ji,jj,1))/e2v(ji,jj) 
    193                   END DO 
    194                END DO 
    195             ENDIF ! ln_shelf_flx 
    196  
    197          DO jj = 1, jpj                                           ! set the ocean fluxes from read fields 
    198             DO ji = 1, jpi 
    199                    IF( ln_shelf_flx ) THEN 
    200                       !! UKMO SHELF - need atmospheric pressure to calculate Haney forcing 
    201                       pressnow(ji,jj) = sf(jp_press)%fnow(ji,jj,1) 
    202                       !! UKMO SHELF flux files contain wind speed not wind stress 
    203                       totwindspd = sqrt(zwnd_i(ji,jj)*zwnd_i(ji,jj) + zwnd_j(ji,jj)*zwnd_j(ji,jj)) 
    204                       cs = 0.63 + (0.066 * totwindspd) 
    205                       utau(ji,jj) = cs * (rhoa/rau0) * zwnd_i(ji,jj) * totwindspd 
    206                       vtau(ji,jj) = cs * (rhoa/rau0) * zwnd_j(ji,jj) * totwindspd 
    207                    ELSE 
    208                       utau(ji,jj) = sf(jp_utau)%fnow(ji,jj,1) 
    209                       vtau(ji,jj) = sf(jp_vtau)%fnow(ji,jj,1) 
    210                    ENDIF 
    211                    qsr (ji,jj) = sf(jp_qsr )%fnow(ji,jj,1) 
    212                    IF( ln_foam_flx .OR. ln_shelf_flx ) THEN 
    213                       !! UKMO FOAM flux files contain non-solar heat flux (qns) rather than total heat flux (qtot) 
    214                       qns (ji,jj) = sf(jp_qtot)%fnow(ji,jj,1) 
    215                       !! UKMO FOAM flux files contain the net DOWNWARD freshwater flux P-E rather then E-P 
    216                       emp (ji,jj) = -1. * sf(jp_emp )%fnow(ji,jj,1) 
    217                    ELSE 
    218                       qns (ji,jj) = sf(jp_qtot)%fnow(ji,jj,1) - sf(jp_qsr)%fnow(ji,jj,1) 
    219                       emp (ji,jj) = sf(jp_emp )%fnow(ji,jj,1) 
    220                    ENDIF 
    221             END DO 
    222          END DO 
    223          !                                                        ! add modification due to drag coefficient read from wave forcing  
    224          !                                                        ! this code is inefficient but put here to allow merging with another UKMO branch  
    225          IF( ln_shelf_flx ) THEN  
    226             ! calculate first the wind module, as it will be used later  
     185 
     186            ! calculate the wind speed module 
    227187            DO jj = 2, jpjm1  
    228188               DO ji = fs_2, fs_jpim1    ! vect. opt.  
     
    233193            END DO  
    234194            CALL lbc_lnk( wndm(:,:), 'T', 1. )  
    235          
     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)  
    236216            IF( ln_cdgw .AND. nn_drag == jp_std ) THEN  
    237217               IF( cpl_wdrag ) THEN   
     
    263243               END DO  
    264244               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 
    265253            ENDIF  
    266          ENDIF 
     254         ELSE 
     255            DO jj = 1, jpj                                           ! set the ocean fluxes from read fields 
     256               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 
    267273         !                                                        ! add to qns the heat due to e-p 
    268274         qns(:,:) = qns(:,:) - emp(:,:) * sst_m(:,:) * rcp        ! mass flux is at SST 
    269275         ! 
    270     
    271             !! UKMO FOAM wind fluxes need lbc_lnk calls owing to a bug in interp.exe 
    272             IF( ln_foam_flx ) THEN 
    273                CALL lbc_lnk( utau(:,:), 'U', -1. ) 
    274                CALL lbc_lnk( vtau(:,:), 'V', -1. ) 
    275             ENDIF 
    276      
     276         !! 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 
    277281         !                                                        ! module of wind stress and wind speed at T-point 
    278282         zcoef = 1. / ( zrhoa * zcdrag ) 
     
    285289               zmod = 0.5 * SQRT( ztx * ztx + zty * zty ) 
    286290               taum(ji,jj) = zmod 
    287                IF ( .NOT. (ln_shelf_flx .AND. ln_cpl)) THEN 
    288                wndm(ji,jj) = SQRT( zmod * zcoef ) 
     291               IF ( .NOT. ln_shelf_flx ) THEN 
     292                  wndm(ji,jj) = SQRT( zmod * zcoef ) 
    289293               ENDIF 
    290294            END DO 
    291295         END DO 
    292296         taum(:,:) = taum(:,:) * tmask(:,:,1) ; wndm(:,:) = wndm(:,:) * tmask(:,:,1) 
    293          CALL lbc_lnk( taum(:,:), 'T', 1. )   ;   CALL lbc_lnk( wndm(:,:), 'T', 1. ) 
     297         CALL lbc_lnk( taum(:,:), 'T', 1. )  ;   CALL lbc_lnk( wndm(:,:), 'T', 1. ) 
    294298 
    295299         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.