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 14701 – NEMO

Changeset 14701


Ignore:
Timestamp:
2021-04-13T15:59:31+02:00 (3 years ago)
Author:
davestorkey
Message:

2021/ENHANCE-01_davestorkey_fix_3D_momentum_trends : science changes

Location:
NEMO/branches/2021/ENHANCE-01_davestorkey_fix_3D_momentum_trends
Files:
7 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2021/ENHANCE-01_davestorkey_fix_3D_momentum_trends/cfgs/SHARED/field_def_nemo-oce.xml

    r14224 r14701  
    10521052  </field_group> 
    10531053 
    1054   <field_group id="trendU" grid_ref="grid_U_3D"> 
    1055     <!-- variables available with ln_dyn_trd --> 
    1056     <field id="utrd_hpg"       long_name="i-trend: hydrostatic pressure gradient"          unit="m/s^2"                        /> 
    1057     <field id="utrd_spg"       long_name="i-trend: surface     pressure gradient"          unit="m/s^2"                        /> 
    1058     <field id="utrd_spgexp"    long_name="i-trend: surface pressure gradient (explicit)"   unit="m/s^2"                        /> 
    1059     <field id="utrd_spgflt"    long_name="i-trend: surface pressure gradient (filtered)"   unit="m/s^2"                        /> 
    1060     <field id="utrd_keg"       long_name="i-trend: KE gradient         or hor. adv."       unit="m/s^2"                        /> 
    1061     <field id="utrd_rvo"       long_name="i-trend: relative  vorticity or metric term"     unit="m/s^2"                        /> 
    1062     <field id="utrd_pvo"       long_name="i-trend: planetary vorticity"                    unit="m/s^2"                        /> 
    1063     <field id="utrd_zad"       long_name="i-trend: vertical  advection"                    unit="m/s^2"                        /> 
    1064     <field id="utrd_udx"       long_name="i-trend: U.dx[U]"                                unit="m/s^2"                        /> 
    1065     <field id="utrd_ldf"       long_name="i-trend: lateral   diffusion"                    unit="m/s^2"                        /> 
    1066     <field id="utrd_zdf"       long_name="i-trend: vertical  diffusion"                    unit="m/s^2"                        /> 
    1067     <field id="utrd_tau"       long_name="i-trend: wind stress "                           unit="m/s^2"   grid_ref="grid_U_2D" /> 
    1068     <field id="utrd_bfr"       long_name="i-trend: bottom friction (explicit)"             unit="m/s^2"                        /> 
    1069     <field id="utrd_bfri"      long_name="i-trend: bottom friction (implicit)"             unit="m/s^2"                        /> 
    1070     <field id="utrd_tot"       long_name="i-trend: total momentum trend before atf"        unit="m/s^2"                        /> 
    1071     <field id="utrd_atf"       long_name="i-trend: asselin time filter trend"              unit="m/s^2"                        /> 
    1072   </field_group> 
    1073  
    1074   <field_group id="trendV" grid_ref="grid_V_3D"> 
    1075     <!-- variables available with ln_dyn_trd --> 
    1076     <field id="vtrd_hpg"       long_name="j-trend: hydrostatic pressure gradient"          unit="m/s^2"                        /> 
    1077     <field id="vtrd_spg"       long_name="j-trend: surface     pressure gradient"          unit="m/s^2"                        /> 
    1078     <field id="vtrd_spgexp"    long_name="j-trend: surface pressure gradient (explicit)"   unit="m/s^2"                        /> 
    1079     <field id="vtrd_spgflt"    long_name="j-trend: surface pressure gradient (filtered)"   unit="m/s^2"                        /> 
    1080     <field id="vtrd_keg"       long_name="j-trend: KE gradient         or hor. adv."       unit="m/s^2"                        /> 
    1081     <field id="vtrd_rvo"       long_name="j-trend: relative  vorticity or metric term"     unit="m/s^2"                        /> 
    1082     <field id="vtrd_pvo"       long_name="j-trend: planetary vorticity"                    unit="m/s^2"                        /> 
    1083     <field id="vtrd_zad"       long_name="j-trend: vertical  advection"                    unit="m/s^2"                        /> 
    1084     <field id="vtrd_vdy"       long_name="i-trend: V.dx[V]"                                unit="m/s^2"                        /> 
    1085     <field id="vtrd_ldf"       long_name="j-trend: lateral   diffusion"                    unit="m/s^2"                        /> 
    1086     <field id="vtrd_zdf"       long_name="j-trend: vertical  diffusion"                    unit="m/s^2"                        /> 
    1087     <field id="vtrd_tau"       long_name="j-trend: wind stress "                           unit="m/s^2"   grid_ref="grid_V_2D" /> 
    1088     <field id="vtrd_bfr"       long_name="j-trend: bottom friction (explicit)"             unit="m/s^2"                        /> 
    1089     <field id="vtrd_bfri"      long_name="j-trend: bottom friction (implicit)"             unit="m/s^2"                        /> 
    1090     <field id="vtrd_tot"       long_name="j-trend: total momentum trend before atf"        unit="m/s^2"                        /> 
    1091     <field id="vtrd_atf"       long_name="j-trend: asselin time filter trend"              unit="m/s^2"                        /> 
    1092   </field_group> 
     1054   <field_group id="trendU" grid_ref="grid_U_3D"> 
     1055     <!-- variables available with ln_dyn_trd --> 
     1056     <field id="utrd_hpg"       long_name="i-trend: hydrostatic pressure gradient"          unit="m/s^2"                        /> 
     1057     <field id="utrd_spg2d"     long_name="i-trend: surface pressure gradient: true trend"       unit="m/s^2"    grid_ref="grid_U_2D" /> 
     1058     <field id="utrd_spgexp"    long_name="i-trend: surface pressure gradient (explicit)"   unit="m/s^2"                        /> 
     1059     <field id="utrd_keg"       long_name="i-trend: KE gradient         or hor. adv."       unit="m/s^2"                        /> 
     1060     <field id="utrd_rvo"       long_name="i-trend: relative  vorticity or metric term"     unit="m/s^2"                        /> 
     1061     <field id="utrd_pvo"       long_name="i-trend: planetary vorticity: 3D component"      unit="m/s^2"                        /> 
     1062     <field id="utrd_pvo2d"     long_name="i-trend: planetary vorticity: 2D component"      unit="m/s^2"    grid_ref="grid_U_2D" /> 
     1063     <field id="utrd_pvo_corr"     long_name="i-trend: planetary vorticity: correction"      unit="m/s^2"    grid_ref="grid_U_2D" /> 
     1064     <field id="utrd_zad"       long_name="i-trend: vertical  advection"                    unit="m/s^2"                        /> 
     1065     <field id="utrd_udx"       long_name="i-trend: U.dx[U]"                                unit="m/s^2"                        /> 
     1066     <field id="utrd_ldf"       long_name="i-trend: lateral   diffusion"                    unit="m/s^2"                        /> 
     1067     <field id="utrd_zdf"       long_name="i-trend: vertical  diffusion"                    unit="m/s^2"                        /> 
     1068     <field id="utrd_tau"       long_name="i-trend: wind stress "                           unit="m/s^2"   grid_ref="grid_U_2D" /> 
     1069     <field id="utrd_bfr"       long_name="i-trend: bottom friction (explicit)"             unit="m/s^2"                        />    
     1070     <field id="utrd_bfri"      long_name="i-trend: bottom friction (implicit)"             unit="m/s^2"                        />    
     1071     <field id="utrd_bfr2d"     long_name="i-trend: bottom friction: 2D component"      unit="m/s^2"    grid_ref="grid_U_2D" /> 
     1072     <field id="utrd_tot"       long_name="i-trend: total momentum trend before atf"        unit="m/s^2"                        />    
     1073     <field id="utrd_atf"       long_name="i-trend: asselin time filter trend"              unit="m/s^2"                        />    
     1074     <!-- thickness weighted versions --> 
     1075     <field id="utrd_hpg_e3u"       unit="m2/s^2" > utrd_hpg * e3u </field> 
     1076     <field id="utrd_spg_e3u"       unit="m2/s^2" > utrd_spg * e3u </field> 
     1077     <field id="utrd_spgexp_e3u"    unit="m2/s^2" > utrd_spgexp * e3u </field> 
     1078     <field id="utrd_keg_e3u"       unit="m2/s^2" > utrd_keg * e3u </field> 
     1079     <field id="utrd_rvo_e3u"       unit="m2/s^2" > utrd_rvo * e3u </field> 
     1080     <field id="utrd_pvo_e3u"       unit="m2/s^2" > utrd_pvo * e3u </field> 
     1081     <field id="utrd_zad_e3u"       unit="m2/s^2" > utrd_zad * e3u </field> 
     1082     <field id="utrd_ldf_e3u"       unit="m2/s^2" > utrd_ldf * e3u </field> 
     1083     <field id="utrd_zdf_e3u"       unit="m2/s^2" > utrd_zdf * e3u </field> 
     1084     <field id="utrd_tot_e3u"       unit="m2/s^2" > utrd_tot * e3u </field> 
     1085     <field id="utrd_atf_e3u"       unit="m2/s^2" > utrd_atf * e3u </field> 
     1086     <field id="utrd_bta_e3u"       unit="m2/s^2" > utrd_bta * e3u </field>   
     1087     <field id="utrd_bfr_e3u"       unit="m2/s^2" > utrd_bfr * e3u </field>   
     1088     <field id="utrd_bfri_e3u"      unit="m2/s^2" > utrd_bfri * e3u </field>   
     1089   </field_group> 
     1090 
     1091   <field_group id="trendV" grid_ref="grid_V_3D"> 
     1092     <!-- variables available with ln_dyn_trd --> 
     1093     <field id="vtrd_hpg"       long_name="j-trend: hydrostatic pressure gradient"          unit="m/s^2"                        /> 
     1094     <field id="vtrd_spg2d"     long_name="j-trend: surface pressure gradient: true trend"       unit="m/s^2"    grid_ref="grid_V_2D" /> 
     1095     <field id="vtrd_spgexp"    long_name="j-trend: surface pressure gradient (explicit)"   unit="m/s^2"                        /> 
     1096     <field id="vtrd_keg"       long_name="j-trend: KE gradient         or hor. adv."       unit="m/s^2"                        /> 
     1097     <field id="vtrd_rvo"       long_name="j-trend: relative  vorticity or metric term"     unit="m/s^2"                        /> 
     1098     <field id="vtrd_pvo"       long_name="j-trend: planetary vorticity: 3D component"      unit="m/s^2"                        /> 
     1099     <field id="vtrd_pvo2d"     long_name="j-trend: planetary vorticity: 2D component"      unit="m/s^2"    grid_ref="grid_V_2D" /> 
     1100     <field id="vtrd_pvo_corr"     long_name="j-trend: planetary vorticity: correction"      unit="m/s^2"    grid_ref="grid_V_2D" /> 
     1101     <field id="vtrd_zad"       long_name="j-trend: vertical  advection"                    unit="m/s^2"                        /> 
     1102     <field id="vtrd_vdy"       long_name="j-trend: V.dx[V]"                                unit="m/s^2"                        /> 
     1103     <field id="vtrd_ldf"       long_name="j-trend: lateral   diffusion"                    unit="m/s^2"                        /> 
     1104     <field id="vtrd_zdf"       long_name="j-trend: vertical  diffusion"                    unit="m/s^2"                        /> 
     1105     <field id="vtrd_tau"       long_name="j-trend: wind stress "                           unit="m/s^2"   grid_ref="grid_V_2D" /> 
     1106     <field id="vtrd_bfr"       long_name="j-trend: bottom friction (explicit)"             unit="m/s^2"                        />    
     1107     <field id="vtrd_bfri"      long_name="j-trend: bottom friction (implicit)"             unit="m/s^2"                        />    
     1108     <field id="vtrd_bfr2d"     long_name="j-trend: bottom friction: 2D component"      unit="m/s^2"    grid_ref="grid_V_2D" /> 
     1109     <field id="vtrd_tot"       long_name="j-trend: total momentum trend before atf"        unit="m/s^2"                        />    
     1110     <field id="vtrd_atf"       long_name="j-trend: asselin time filter trend"              unit="m/s^2"                        />    
     1111     <!-- thickness weighted versions --> 
     1112     <field id="vtrd_hpg_e3v"       unit="m2/s^2" > vtrd_hpg * e3v </field> 
     1113     <field id="vtrd_spg_e3v"       unit="m2/s^2" > vtrd_spg * e3v </field> 
     1114     <field id="vtrd_spgexp_e3v"    unit="m2/s^2" > vtrd_spgexp * e3v </field> 
     1115     <field id="vtrd_keg_e3v"       unit="m2/s^2" > vtrd_keg * e3v </field> 
     1116     <field id="vtrd_rvo_e3v"       unit="m2/s^2" > vtrd_rvo * e3v </field> 
     1117     <field id="vtrd_pvo_e3v"       unit="m2/s^2" > vtrd_pvo * e3v </field> 
     1118     <field id="vtrd_zad_e3v"       unit="m2/s^2" > vtrd_zad * e3v </field> 
     1119     <field id="vtrd_ldf_e3v"       unit="m2/s^2" > vtrd_ldf * e3v </field> 
     1120     <field id="vtrd_zdf_e3v"       unit="m2/s^2" > vtrd_zdf * e3v </field> 
     1121     <field id="vtrd_tot_e3v"       unit="m2/s^2" > vtrd_tot * e3v </field> 
     1122     <field id="vtrd_atf_e3v"       unit="m2/s^2" > vtrd_atf * e3v </field> 
     1123     <field id="vtrd_bta_e3v"       unit="m2/s^2" > vtrd_bta * e3v </field>   
     1124     <field id="vtrd_bfr_e3v"       unit="m2/s^2" > vtrd_bfr * e3v </field>   
     1125     <field id="vtrd_bfri_e3v"      unit="m2/s^2" > vtrd_bfri * e3v </field>   
     1126   </field_group> 
    10931127 
    10941128 
  • NEMO/branches/2021/ENHANCE-01_davestorkey_fix_3D_momentum_trends/src/OCE/DYN/dynspg.F90

    r14225 r14701  
    9090      IF( ln_timing )   CALL timing_start('dyn_spg') 
    9191      ! 
    92       IF( l_trddyn )   THEN                      ! temporary save of ta and sa trends 
     92      IF( l_trddyn .AND. nspg == np_EXP )   THEN   ! temporary save of ta and sa trends 
    9393         ALLOCATE( ztrdu(jpi,jpj,jpk) , ztrdv(jpi,jpj,jpk) ) 
    9494         ztrdu(:,:,:) = puu(:,:,:,Krhs) 
     
    170170      END SELECT 
    171171      ! 
    172       IF( l_trddyn )   THEN                  ! save the surface pressure gradient trends for further diagnostics 
     172      IF( l_trddyn .AND. nspg == np_EXP )   THEN   ! save the surface pressure gradient trends for further diagnostics 
    173173         ztrdu(:,:,:) = puu(:,:,:,Krhs) - ztrdu(:,:,:) 
    174174         ztrdv(:,:,:) = pvv(:,:,:,Krhs) - ztrdv(:,:,:) 
    175175         CALL trd_dyn( ztrdu, ztrdv, jpdyn_spg, kt, Kmm ) 
     176         ! In this case we also need to finalise the write-out of the PVO 3D diagnostic with zero correction 
     177         ztrdu(:,:,1) = 0._wp 
     178         ztrdv(:,:,1) = 0._wp 
     179         CALL trd_dyn( ztrdu(:,:,1), ztrdv(:,:,1), jpdyn_pvo, kt, Kmm ) 
    176180         DEALLOCATE( ztrdu , ztrdv ) 
    177181      ENDIF 
  • NEMO/branches/2021/ENHANCE-01_davestorkey_fix_3D_momentum_trends/src/OCE/DYN/dynspg_ts.F90

    r14433 r14701  
    6060   USE iom             ! IOM library 
    6161   USE restart         ! only for lrst_oce 
     62   USE trd_oce         ! trends: ocean variables 
     63   USE trddyn          ! trend manager: dynamics 
    6264 
    6365   USE iom   ! to remove 
     
    154156      REAL(wp) ::   r1_Dt_b, z1_hu, z1_hv          ! local scalars 
    155157      REAL(wp) ::   za0, za1, za2, za3              !   -      - 
    156       REAL(wp) ::   zztmp, zldg               !   -      - 
     158      REAL(wp) ::   zztmp, zldg                     !   -      - 
    157159      REAL(wp) ::   zhu_bck, zhv_bck, zhdiv         !   -      - 
    158160      REAL(wp) ::   zun_save, zvn_save              !   -      - 
     
    175177      REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: ztwdmask, zuwdmask, zvwdmask ! ROMS wetting and drying masks at t,u,v points 
    176178      REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zuwdav2, zvwdav2    ! averages over the sub-steps of zuwdmask and zvwdmask 
     179      REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zspgtrdu, zspgtrdv, zpvotrdu, zpvotrdv  ! SPG and PVO trends (if l_trddyn) 
     180      REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: ztautrdu, ztautrdv, zbfrtrdu, zbfrtrdv  ! TAU and BFR trends (if l_trddyn) 
    177181      REAL(wp) ::   zt0substep !   Time of day at the beginning of the time substep 
    178182      !!---------------------------------------------------------------------- 
     
    181185      !                                         !* Allocate temporary arrays 
    182186      IF( ln_wd_dl ) ALLOCATE( ztwdmask(jpi,jpj), zuwdmask(jpi,jpj), zvwdmask(jpi,jpj), zuwdav2(jpi,jpj), zvwdav2(jpi,jpj)) 
     187      ! 
     188      IF( l_trddyn ) THEN 
     189          ALLOCATE( zspgtrdu(jpi,jpj), zspgtrdv(jpi,jpj), zpvotrdu(jpi,jpj), zpvotrdv(jpi,jpj), & 
     190         &          ztautrdu(jpi,jpj), ztautrdv(jpi,jpj), zbfrtrdu(jpi,jpj), zbfrtrdv(jpi,jpj) ) 
     191          zspgtrdu(:,:) = 0._wp 
     192          zspgtrdv(:,:) = 0._wp 
     193          zpvotrdu(:,:) = 0._wp 
     194          zpvotrdv(:,:) = 0._wp 
     195          ztautrdu(:,:) = 0._wp 
     196          ztautrdv(:,:) = 0._wp 
     197          zbfrtrdu(:,:) = 0._wp 
     198          zbfrtrdv(:,:) = 0._wp 
     199      ENDIF 
     200      ! 
     201      zu_trd(:,:) = 0._wp 
     202      zv_trd(:,:) = 0._wp 
     203      zu_spg(:,:) = 0._wp 
     204      zv_spg(:,:) = 0._wp 
     205      ! 
     206      zmdi=1.e+20                               !  missing data indicator for masking 
    183207      ! 
    184208      zwdramp = r_rn_wdmin1               ! simplest ramp  
     
    263287            &                                                                          zu_trd, zv_trd   )   ! ==>> out 
    264288         ! 
    265          DO_2D( 0, 0, 0, 0 )                          ! Remove coriolis term (and possibly spg) from barotropic trend 
     289         IF( l_trddyn ) THEN 
     290            ! send correction to baroclinic planetary vorticity trend to trd_dyn 
     291            CALL trd_dyn( zu_trd, zv_trd, jpdyn_pvo_corr, kt, Kmm ) 
     292         ENDIF 
     293         ! 
     294         DO_2D( 0, 0, 0, 0 )                          ! Remove coriolis term from barotropic trend 
    266295            zu_frc(ji,jj) = zu_frc(ji,jj) - zu_trd(ji,jj) * ssumask(ji,jj) 
    267296            zv_frc(ji,jj) = zv_frc(ji,jj) - zv_trd(ji,jj) * ssvmask(ji,jj) 
     
    277306         END_2D 
    278307      ELSE                !* remove baroclinic drag AND provide the barotropic drag coefficients 
     308         IF( l_trddyn ) THEN 
     309            zbfrtrdu(:,:) = zu_frc(:,:) 
     310            zbfrtrdv(:,:) = zv_frc(:,:) 
     311         ENDIF 
     312         ! 
    279313         CALL dyn_drg_init( Kbb, Kmm, puu, pvv, puu_b, pvv_b, zu_frc, zv_frc, zCdU_u, zCdU_v ) 
     314         ! 
     315         IF( l_trddyn ) THEN 
     316            ! bottom friction trend diagnostic: bottom friction due to baroclinic currents 
     317            zbfrtrdu(:,:) = zu_frc(:,:) - zbfrtrdu(:,:) 
     318            zbfrtrdv(:,:) = zv_frc(:,:) - zbfrtrdv(:,:)  
     319         ENDIF 
    280320      ENDIF 
    281321      ! 
     
    301341      !                                   !=  Add wind forcing  =! 
    302342      !                                   !  ------------------  ! 
     343      IF( l_trddyn ) THEN 
     344         ztautrdu(:,:) = zu_frc(:,:) 
     345         ztautrdv(:,:) = zv_frc(:,:) 
     346      ENDIF 
     347      ! 
    303348      IF( ln_bt_fw ) THEN 
    304349         DO_2D( 0, 0, 0, 0 ) 
     
    313358         END_2D 
    314359      ENDIF   
     360      ! 
     361      IF( l_trddyn ) THEN 
     362         ! wind stress trend diagnostic 
     363         ztautrdu(:,:) = zu_frc(:,:) - ztautrdu(:,:) 
     364         ztautrdv(:,:) = zv_frc(:,:) - ztautrdv(:,:)  
     365      ENDIF 
    315366      ! 
    316367      !              !----------------! 
     
    583634         ENDIF 
    584635         ! 
     636         IF( l_trddyn ) THEN 
     637            za2 = wgtbtp2(jn) 
     638            zspgtrdu(:,:) = zspgtrdu(:,:) + za2 * zu_spg(:,:) * ssumask(:,:) 
     639            zspgtrdv(:,:) = zspgtrdv(:,:) + za2 * zv_spg(:,:) * ssvmask(:,:) 
     640         ENDIF 
     641         ! 
    585642         ! Add Coriolis trend: 
    586643         ! zwz array below or triads normally depend on sea level with ln_linssh=F and should be updated 
     
    588645         ! Recall that zhU and zhV hold fluxes at jn+0.5 (extrapolated not backward interpolated) 
    589646         CALL dyn_cor_2d( zhtp2_e, zhup2_e, zhvp2_e, ua_e, va_e, zhU, zhV,    zu_trd, zv_trd   ) 
     647         ! 
     648         IF( l_trddyn ) THEN 
     649            za2 = wgtbtp2(jn) 
     650            zpvotrdu(:,:) = zpvotrdu(:,:) + za2 * zu_trd(:,:) * ssumask(:,:) 
     651            zpvotrdv(:,:) = zpvotrdv(:,:) + za2 * zv_trd(:,:) * ssvmask(:,:) 
     652         ENDIF 
    590653         ! 
    591654         ! Add tidal astronomical forcing if defined 
     
    604667               zv_trd(ji,jj) = zv_trd(ji,jj) + zCdU_v(ji,jj) * vn_e(ji,jj) * hvr_e(ji,jj) 
    605668            END_2D 
     669            IF( l_trddyn ) THEN 
     670               za2 = wgtbtp2(jn) 
     671               zbfrtrdu(:,:) = zbfrtrdu(:,:) + za2 * zCdU_u(:,:) * un_e(:,:) * hur_e(:,:) 
     672               zbfrtrdv(:,:) = zbfrtrdv(:,:) + za2 * zCdU_v(:,:) * vn_e(:,:) * hvr_e(:,:) 
     673            ENDIF 
    606674         ENDIF 
    607675         ! 
     
    828896      IF( ln_wd_il )   DEALLOCATE( zcpx, zcpy ) 
    829897      IF( ln_wd_dl )   DEALLOCATE( ztwdmask, zuwdmask, zvwdmask, zuwdav2, zvwdav2 ) 
     898      ! 
     899      IF( l_trddyn ) THEN 
     900         CALL trd_dyn( zspgtrdu, zspgtrdv, jpdyn_spg, kt ) 
     901         CALL trd_dyn( zpvotrdu, zpvotrdv, jpdyn_pvo, kt ) 
     902         CALL trd_dyn( ztautrdu, ztautrdv, jpdyn_tau, kt ) 
     903         CALL trd_dyn( zbfrtrdu, zbfrtrdv, jpdyn_bfr, kt ) 
     904         DEALLOCATE( zspgtrdu, zspgtrdv, zpvotrdu, zpvotrdv, ztautrdu, ztautrdv, zbfrtrdu, zbfrtrdv ) 
     905      ENDIF 
    830906      ! 
    831907      CALL iom_put( "baro_u" , puu_b(:,:,Kmm) )  ! Barotropic  U Velocity 
  • NEMO/branches/2021/ENHANCE-01_davestorkey_fix_3D_momentum_trends/src/OCE/DYN/dynvor.F90

    r14433 r14701  
    138138         ztrdu(:,:,:) = puu(:,:,:,Krhs) - ztrdu(:,:,:) 
    139139         ztrdv(:,:,:) = pvv(:,:,:,Krhs) - ztrdv(:,:,:) 
    140          CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt, Kmm ) 
     140         CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo_save, kt, Kmm ) 
    141141         ! 
    142142         IF( n_dynadv /= np_LIN_dyn ) THEN   !* relative vorticity or metric trend (only in non-linear case) 
  • NEMO/branches/2021/ENHANCE-01_davestorkey_fix_3D_momentum_trends/src/OCE/DYN/dynzdf.F90

    r13497 r14701  
    7272      INTEGER  ::   ji, jj, jk         ! dummy loop indices 
    7373      INTEGER  ::   iku, ikv           ! local integers 
     74      INTEGER  ::   ikbu, ikbv         ! local integers 
    7475      REAL(wp) ::   zzwi, ze3ua, zdt   ! local scalars 
    7576      REAL(wp) ::   zzws, ze3va        !   -      - 
     
    9596      !                             !* explicit top/bottom drag case 
    9697      IF( .NOT.ln_drgimp )   CALL zdf_drg_exp( kt, Kmm, puu(:,:,:,Kbb), pvv(:,:,:,Kbb), puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )  ! add top/bottom friction trend to (puu(Kaa),pvv(Kaa)) 
    97       ! 
    98       ! 
    99       IF( l_trddyn )   THEN         !* temporary save of ta and sa trends 
    100          ALLOCATE( ztrdu(jpi,jpj,jpk), ztrdv(jpi,jpj,jpk) )  
    101          ztrdu(:,:,:) = puu(:,:,:,Krhs) 
    102          ztrdv(:,:,:) = pvv(:,:,:,Krhs) 
    103       ENDIF 
    10498      ! 
    10599      !              !==  RHS: Leap-Frog time stepping on all trends but the vertical mixing  ==!   (put in puu(:,:,:,Kaa),pvv(:,:,:,Kaa)) 
     
    131125            pvv(ji,jj,jk,Kaa) = ( pvv(ji,jj,jk,Kaa) - vv_b(ji,jj,Kaa) ) * vmask(ji,jj,jk) 
    132126         END_3D 
     127         IF( l_trddyn )   THEN         !* temporary save of ta and sa trends 
     128            ALLOCATE( ztrdu(jpi,jpj,jpk), ztrdv(jpi,jpj,jpk) )  
     129            ztrdu(:,:,:) = puu(:,:,:,Kaa) 
     130            ztrdv(:,:,:) = pvv(:,:,:,Kaa) 
     131         ENDIF 
    133132         DO_2D( 0, 0, 0, 0 )      ! Add bottom/top stress due to barotropic component only 
    134133            iku = mbku(ji,jj)         ! ocean bottom level at u- and v-points  
     
    141140            pvv(ji,jj,ikv,Kaa) = pvv(ji,jj,ikv,Kaa) + rDt * 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) * vv_b(ji,jj,Kaa) / ze3va 
    142141         END_2D 
     142         IF( l_trddyn )   THEN                      ! save explicit part of bottom friction trends 
     143            ztrdu(:,:,:) = ( puu(:,:,:,Kaa) - ztrdu(:,:,:) ) / rDt  
     144            ztrdv(:,:,:) = ( pvv(:,:,:,Kaa) - ztrdv(:,:,:) ) / rDt  
     145            CALL trd_dyn( ztrdu, ztrdv, jpdyn_bfr, kt, Kmm ) 
     146         ENDIF 
    143147         IF( ln_isfcav.OR.ln_drgice_imp ) THEN    ! Ocean cavities (ISF) 
    144148            DO_2D( 0, 0, 0, 0 ) 
     
    155159      ENDIF 
    156160      ! 
     161      IF( l_trddyn )   THEN         !* temporary save of ta and sa trends 
     162         IF( .NOT. ALLOCATED(ztrdu) ) ALLOCATE( ztrdu(jpi,jpj,jpk), ztrdv(jpi,jpj,jpk) )  
     163         ztrdu(:,:,:) = puu(:,:,:,Kaa) 
     164         ztrdv(:,:,:) = pvv(:,:,:,Kaa) 
     165      ENDIF 
     166      ! 
    157167      !              !==  Vertical diffusion on u  ==! 
    158168      ! 
     
    432442      ! 
    433443      IF( l_trddyn )   THEN                      ! save the vertical diffusive trends for further diagnostics 
    434          ztrdu(:,:,:) = ( puu(:,:,:,Kaa) - puu(:,:,:,Kbb) ) / rDt - ztrdu(:,:,:) 
    435          ztrdv(:,:,:) = ( pvv(:,:,:,Kaa) - pvv(:,:,:,Kbb) ) / rDt - ztrdv(:,:,:) 
     444         ztrdu(:,:,:) = ( puu(:,:,:,Kaa) - ztrdu(:,:,:) ) / rDt  
     445         ztrdv(:,:,:) = ( pvv(:,:,:,Kaa) - ztrdv(:,:,:) ) / rDt 
    436446         CALL trd_dyn( ztrdu, ztrdv, jpdyn_zdf, kt, Kmm ) 
     447         ! 
     448         IF( ln_drgimp ) THEN 
     449            ztrdu(:,:,:) = 0._wp    ;   ztrdv(:,:,:) = 0._wp 
     450            DO_2D( 2, jpim1, 2, jpjm1 ) 
     451               ikbu = mbku(ji,jj)          ! deepest ocean u- & v-levels 
     452               ikbv = mbkv(ji,jj) 
     453               ztrdu(ji,jj,ikbu) = 0.5 * ( rCdU_bot(ji+1,jj) + rCdU_bot(ji,jj) )  &  
     454     &                         * puu(ji,jj,ikbu,Kmm) / e3u(ji,jj,ikbu,Kmm) 
     455               ztrdv(ji,jj,ikbv) = 0.5 * ( rCdU_bot(ji,jj+1) + rCdU_bot(ji,jj) )  & 
     456     &                         * pvv(ji,jj,ikbv,Kmm) / e3v(ji,jj,ikbv,Kmm) 
     457            END_2D 
     458            CALL trd_dyn( ztrdu, ztrdv, jpdyn_bfri, kt, Kmm ) 
     459         ENDIF 
     460         ! 
    437461         DEALLOCATE( ztrdu, ztrdv )  
    438462      ENDIF 
  • NEMO/branches/2021/ENHANCE-01_davestorkey_fix_3D_momentum_trends/src/OCE/TRD/trd_oce.F90

    r14239 r14701  
    6060   ! 
    6161   !                                                  !!!* Momentum trends indices 
    62    INTEGER, PUBLIC, PARAMETER ::   jptot_dyn  = 13     !: Total trend nb: change it when adding/removing one indice below 
     62   INTEGER, PUBLIC, PARAMETER ::   jptot_dyn  = 17     !: Total trend nb: change it when adding/removing one indice below 
    6363   !                               ===============     !   
    64    INTEGER, PUBLIC, PARAMETER ::   jpdyn_hpg  =  1     !: hydrostatic pressure gradient  
    65    INTEGER, PUBLIC, PARAMETER ::   jpdyn_spg  =  2     !: surface     pressure gradient 
    66    INTEGER, PUBLIC, PARAMETER ::   jpdyn_keg  =  3     !: kinetic energy gradient  or horizontal advection 
    67    INTEGER, PUBLIC, PARAMETER ::   jpdyn_rvo  =  4     !: relative  vorticity      or metric term 
    68    INTEGER, PUBLIC, PARAMETER ::   jpdyn_pvo  =  5     !: planetary vorticity 
    69    INTEGER, PUBLIC, PARAMETER ::   jpdyn_zad  =  6     !: vertical advection 
    70    INTEGER, PUBLIC, PARAMETER ::   jpdyn_ldf  =  7     !: horizontal diffusion    
    71    INTEGER, PUBLIC, PARAMETER ::   jpdyn_zdf  =  8     !: vertical   diffusion 
    72    INTEGER, PUBLIC, PARAMETER ::   jpdyn_bfr  =  9     !: bottom  stress  
    73    INTEGER, PUBLIC, PARAMETER ::   jpdyn_atf  = 10     !: Asselin time filter 
    74    INTEGER, PUBLIC, PARAMETER ::   jpdyn_tau  = 11     !: surface stress 
    75    INTEGER, PUBLIC, PARAMETER ::   jpdyn_bfri = 12     !: implicit bottom friction (ln_drgimp=.TRUE.) 
    76    INTEGER, PUBLIC, PARAMETER ::   jpdyn_ken  = 13     !: use for calculation of KE 
     64   INTEGER, PUBLIC, PARAMETER ::   jpdyn_hpg   =  1     !: hydrostatic pressure gradient  
     65   INTEGER, PUBLIC, PARAMETER ::   jpdyn_hpg_save =  2  !: hydrostatic pressure gradient (saved value) 
     66   INTEGER, PUBLIC, PARAMETER ::   jpdyn_hpg_corr =  3  !: hydrostatic pressure gradient (initial correction) 
     67   INTEGER, PUBLIC, PARAMETER ::   jpdyn_spg   =  4     !: surface     pressure gradient 
     68   INTEGER, PUBLIC, PARAMETER ::   jpdyn_keg   =  5     !: kinetic energy gradient  or horizontal advection 
     69   INTEGER, PUBLIC, PARAMETER ::   jpdyn_rvo   =  6     !: relative  vorticity      or metric term 
     70   INTEGER, PUBLIC, PARAMETER ::   jpdyn_pvo   =  7     !: planetary vorticity 
     71   INTEGER, PUBLIC, PARAMETER ::   jpdyn_pvo_save =  8  !: planetary vorticity (saved value) 
     72   INTEGER, PUBLIC, PARAMETER ::   jpdyn_pvo_corr =  9  !: planetary vorticity (initial correction) 
     73   INTEGER, PUBLIC, PARAMETER ::   jpdyn_zad   = 10     !: vertical advection 
     74   INTEGER, PUBLIC, PARAMETER ::   jpdyn_ldf   = 11     !: horizontal diffusion    
     75   INTEGER, PUBLIC, PARAMETER ::   jpdyn_zdf   = 12     !: vertical   diffusion 
     76   INTEGER, PUBLIC, PARAMETER ::   jpdyn_bfr   = 13     !: bottom  stress  
     77   INTEGER, PUBLIC, PARAMETER ::   jpdyn_atf   = 14     !: Asselin time filter 
     78   INTEGER, PUBLIC, PARAMETER ::   jpdyn_tau   = 15     !: surface stress 
     79   INTEGER, PUBLIC, PARAMETER ::   jpdyn_bfri  = 16     !: implicit bottom friction (ln_drgimp=.TRUE.) 
     80   INTEGER, PUBLIC, PARAMETER ::   jpdyn_ken   = 17     !: use for calculation of KE 
    7781   ! 
    7882   !!---------------------------------------------------------------------- 
  • NEMO/branches/2021/ENHANCE-01_davestorkey_fix_3D_momentum_trends/src/OCE/TRD/trddyn.F90

    r14433 r14701  
    1818   USE sbc_oce        ! surface boundary condition: ocean 
    1919   USE zdf_oce        ! ocean vertical physics: variables 
    20 !!gm   USE zdfdrg         ! ocean vertical physics: bottom friction 
    2120   USE trd_oce        ! trends: ocean variables 
    2221   USE trdken         ! trends: Kinetic ENergy  
     
    3534   PUBLIC trd_dyn        ! called by all dynXXX modules 
    3635 
     36   INTERFACE trd_dyn 
     37      module procedure trd_dyn_3d, trd_dyn_2d 
     38   END INTERFACE 
     39 
     40   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), SAVE :: zutrd_hpg, zvtrd_hpg 
     41   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), SAVE :: zutrd_pvo, zvtrd_pvo 
     42   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), SAVE :: zutrd_bfr, zvtrd_bfr 
     43   REAL(wp), ALLOCATABLE, DIMENSION(:,:)  , SAVE :: zutrd_tau, zvtrd_tau 
     44 
    3745   !! * Substitutions 
    3846#  include "do_loop_substitute.h90" 
     
    4553CONTAINS 
    4654 
    47    SUBROUTINE trd_dyn( putrd, pvtrd, ktrd, kt, Kmm ) 
     55   SUBROUTINE trd_dyn_3d( putrd, pvtrd, ktrd, kt, Kmm ) 
    4856      !!--------------------------------------------------------------------- 
    49       !!                  ***  ROUTINE trd_mod  *** 
     57      !!                  ***  ROUTINE trd_dyn_3d  *** 
    5058      !!  
    5159      !! ** Purpose :   Dispatch momentum trend computation, e.g. 3D output,  
     
    5765      INTEGER                   , INTENT(in   ) ::   kt             ! time step 
    5866      INTEGER                   , INTENT(in   ) ::   Kmm            ! time level index 
     67      REAL(wp), ALLOCATABLE, DIMENSION(:,:)     ::   zue, zve       ! temporary 2D arrays 
     68      INTEGER                                   ::   jk 
    5969      !!---------------------------------------------------------------------- 
    6070      ! 
     
    6575!!gm NB : here a lbc_lnk should probably be added 
    6676 
     77      SELECT CASE( ktrd ) 
     78      CASE( jpdyn_pvo_save )  
     79         ! 
     80         ! save 3D coriolis trends to possibly have barotropic part corrected later before writing out 
     81         ALLOCATE( zutrd_pvo(jpi,jpj,jpk), zvtrd_pvo(jpi,jpj,jpk) ) 
     82         zutrd_pvo(:,:,:) = putrd(:,:,:) 
     83         zvtrd_pvo(:,:,:) = pvtrd(:,:,:) 
     84 
     85      CASE( jpdyn_spg )  
     86         ! For explicit scheme SPG trends come here as 3D fields 
     87         ! Add SPG trend to 3D HPG trend and also output as 2D diagnostic in own right. 
     88         CALL trd_dyn_iom_2d( putrd(:,:,1), pvtrd(:,:,1), jpdyn_spg, kt, Kmm )  
     89         putrd(:,:,:) = putrd(:,:,:) + zutrd_hpg(:,:,:)  
     90         pvtrd(:,:,:) = pvtrd(:,:,:) + zvtrd_hpg(:,:,:)  
     91         DEALLOCATE( zutrd_hpg, zvtrd_hpg ) 
     92 
     93      CASE( jpdyn_bfr ) 
     94         ! 
     95         ! Add 3D part of BFR trend minus its depth-mean part to depth-mean trend already saved. 
     96         ALLOCATE( zue(jpi,jpj), zve(jpi,jpj) ) 
     97         zue(:,:) = e3u_a(:,:,1) * putrd(:,:,1) * umask(:,:,1) 
     98         zve(:,:) = e3v_a(:,:,1) * pvtrd(:,:,1) * vmask(:,:,1) 
     99         DO jk = 2, jpkm1 
     100            zue(:,:) = zue(:,:) + e3u_a(:,:,jk) * putrd(:,:,jk) * umask(:,:,jk) 
     101            zve(:,:) = zve(:,:) + e3v_a(:,:,jk) * pvtrd(:,:,jk) * vmask(:,:,jk) 
     102         END DO 
     103         DO jk = 1, jpkm1 
     104            putrd(:,:,jk) = zutrd_bfr(:,:,jk) + putrd(:,:,jk) - zue(:,:) * r1_hu_a(:,:) 
     105            pvtrd(:,:,jk) = zvtrd_bfr(:,:,jk) + pvtrd(:,:,jk) - zve(:,:) * r1_hv_a(:,:) 
     106         END DO 
     107         ! Update locally saved BFR trends to add to ZDF trend. 
     108         zutrd_bfr(:,:,:) = putrd(:,:,:)  
     109         zvtrd_bfr(:,:,:) = pvtrd(:,:,:) 
     110 
     111      CASE( jpdyn_zdf )  
     112         ! ZDF trend: Remove barotropic component and add wind stress and bottom friction 
     113         !            trends from dynspg_ts. Also adding on the bottom stress for the  
     114         !            baroclinic solution in the case of explicit bottom friction.  
     115         ALLOCATE( zue(jpi,jpj), zve(jpi,jpj) ) 
     116         zue(:,:) = e3u_a(:,:,1) * putrd(:,:,1) * umask(:,:,1) 
     117         zve(:,:) = e3v_a(:,:,1) * pvtrd(:,:,1) * vmask(:,:,1) 
     118         DO jk = 2, jpkm1 
     119            zue(:,:) = zue(:,:) + e3u_a(:,:,jk) * putrd(:,:,jk) * umask(:,:,jk) 
     120            zve(:,:) = zve(:,:) + e3v_a(:,:,jk) * pvtrd(:,:,jk) * vmask(:,:,jk) 
     121         END DO 
     122         DO jk = 1, jpkm1 
     123            putrd(:,:,jk) = zutrd_tau(:,:) + zutrd_bfr(:,:,jk) + putrd(:,:,jk) - zue(:,:) * r1_hu_a(:,:) 
     124            pvtrd(:,:,jk) = zvtrd_tau(:,:) + zvtrd_bfr(:,:,jk) + pvtrd(:,:,jk) - zve(:,:) * r1_hv_a(:,:) 
     125         END DO 
     126         DEALLOCATE( zue, zve, zutrd_tau, zvtrd_tau, zutrd_bfr, zvtrd_bfr ) 
     127 
     128      END SELECT 
     129 
     130      IF ( ktrd /= jpdyn_pvo_save ) THEN 
     131         ! 
     132         !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
     133         !   3D output of momentum and/or tracers trends using IOM interface 
     134         !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
     135         IF( ln_dyn_trd )   CALL trd_dyn_iom_3d( putrd, pvtrd, ktrd, kt, Kmm ) 
     136 
     137         !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
     138         !  Integral Constraints Properties for momentum and/or tracers trends 
     139         !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
     140         IF( ln_glo_trd )   CALL trd_glo( putrd, pvtrd, ktrd, 'DYN', kt, Kmm ) 
     141 
     142         !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
     143         !  Kinetic Energy trends 
     144         !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
     145         IF( ln_KE_trd  )   CALL trd_ken( putrd, pvtrd, ktrd, kt, Kmm ) 
     146 
     147         !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
     148         !  Vorticity trends 
     149         !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
     150         IF( ln_vor_trd )   CALL trd_vor( putrd, pvtrd, ktrd, kt, Kmm ) 
     151 
     152         !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
     153         !  Mixed layer trends for active tracers 
     154         !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
     155         !!gm      IF( ln_dyn_mxl )   CALL trd_mxl_dyn    
     156         ! 
     157      ENDIF 
     158      ! 
     159   END SUBROUTINE trd_dyn_3d 
     160 
     161 
     162   SUBROUTINE trd_dyn_2d( putrd, pvtrd, ktrd, kt, Kmm ) 
     163      !!--------------------------------------------------------------------- 
     164      !!                  ***  ROUTINE trd_mod  *** 
     165      !!  
     166      !! ** Purpose :   Dispatch momentum trend computation, e.g. 2D output,  
     167      !!              integral constraints, barotropic vorticity, kinetic enrgy,  
     168      !!              and/or mixed layer budget. 
     169      !!---------------------------------------------------------------------- 
     170      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   putrd, pvtrd   ! U and V trends  
     171      INTEGER                 , INTENT(in   ) ::   ktrd           ! trend index 
     172      INTEGER                 , INTENT(in   ) ::   kt             ! time step 
     173      INTEGER                 , INTENT(in   ) ::   Kmm            ! time level index 
     174      INTEGER                                 ::   jk 
     175      !!---------------------------------------------------------------------- 
     176      ! 
     177      putrd(:,:) = putrd(:,:) * umask(:,:,1)                       ! mask the trends 
     178      pvtrd(:,:) = pvtrd(:,:) * vmask(:,:,1) 
     179      ! 
     180 
     181!!gm NB : here a lbc_lnk should probably be added 
     182 
     183      SELECT CASE(ktrd) 
     184 
     185      CASE( jpdyn_pvo_corr ) 
     186         ! 
     187         ! Remove "first-guess" barotropic coriolis trend from 3D PVO trend.  
     188         DO jk = 1, jpkm1 
     189            zutrd_pvo(:,:,jk) = zutrd_pvo(:,:,jk) - putrd(:,:) 
     190            zvtrd_pvo(:,:,jk) = zvtrd_pvo(:,:,jk) - pvtrd(:,:) 
     191         ENDDO 
     192 
     193      CASE( jpdyn_spg ) 
     194          ! 
     195          ! For split-explicit scheme SPG trends come here as 2D fields 
     196          ! Add SPG trend to 3D HPG trend and also output as 2D diagnostic in own right. 
     197          DO jk = 1, jpkm1 
     198             zutrd_hpg(:,:,jk) = zutrd_hpg(:,:,jk) + putrd(:,:) 
     199             zvtrd_hpg(:,:,jk) = zvtrd_hpg(:,:,jk) + pvtrd(:,:) 
     200          ENDDO 
     201          CALL trd_dyn_3d( zutrd_hpg, zvtrd_hpg, jpdyn_hpg, kt ) 
     202          DEALLOCATE( zutrd_hpg, zvtrd_hpg ) 
     203 
     204      CASE( jpdyn_pvo ) 
     205          ! 
     206          ! Add 2D PVO trend to 3D PVO trend and also output as diagnostic in own right. 
     207          DO jk = 1, jpkm1 
     208             zutrd_pvo(:,:,jk) = zutrd_pvo(:,:,jk) + putrd(:,:) 
     209             zvtrd_pvo(:,:,jk) = zvtrd_pvo(:,:,jk) + pvtrd(:,:) 
     210          ENDDO 
     211          CALL trd_dyn_3d( zutrd_pvo, zvtrd_pvo, jpdyn_pvo, kt ) 
     212          DEALLOCATE( zutrd_pvo, zvtrd_pvo ) 
     213 
     214      CASE( jpdyn_tau ) 
     215          ! 
     216          ! Save 2D wind forcing trend locally (to be added to ZDF trend) 
     217          ! and output as a trend in its own right. 
     218          ALLOCATE( zutrd_tau(jpi,jpj), zvtrd_tau(jpi,jpj) ) 
     219          zutrd_tau(:,:) = putrd(:,:) 
     220          zvtrd_tau(:,:) = pvtrd(:,:) 
     221 
     222      CASE( jpdyn_bfr ) 
     223          ! 
     224          ! Create 3D BFR trend from 2D field and also output 2D field as diagnostic in own right. 
     225          ALLOCATE( zutrd_bfr(jpi,jpj,jpk), zvtrd_bfr(jpi,jpj,jpk) ) 
     226          zutrd_bfr(:,:,:) = 0.0 
     227          zvtrd_bfr(:,:,:) = 0.0 
     228          DO jk = 1, jpkm1 
     229             zutrd_bfr(:,:,jk) = putrd(:,:) * umask(:,:,jk) 
     230             zvtrd_bfr(:,:,jk) = pvtrd(:,:) * vmask(:,:,jk) 
     231          ENDDO 
     232 
     233      END SELECT 
     234 
    67235      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    68       !   3D output of momentum and/or tracers trends using IOM interface 
     236      !   2D output of momentum and/or tracers trends using IOM interface 
    69237      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    70       IF( ln_dyn_trd )   CALL trd_dyn_iom( putrd, pvtrd, ktrd, kt, Kmm ) 
     238      IF( ln_dyn_trd )   CALL trd_dyn_iom_2d( putrd, pvtrd, ktrd, kt, Kmm ) 
    71239          
    72       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    73       !  Integral Constraints Properties for momentum and/or tracers trends 
    74       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    75       IF( ln_glo_trd )   CALL trd_glo( putrd, pvtrd, ktrd, 'DYN', kt, Kmm ) 
    76  
    77       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    78       !  Kinetic Energy trends 
    79       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    80       IF( ln_KE_trd  )   CALL trd_ken( putrd, pvtrd, ktrd, kt, Kmm ) 
    81  
    82       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    83       !  Vorticity trends 
    84       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    85       IF( ln_vor_trd )   CALL trd_vor( putrd, pvtrd, ktrd, kt, Kmm ) 
     240 
     241!!$   CALLS TO THESE ROUTINES FOR 2D DIAGOSTICS NOT CODED YET 
     242!!$      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
     243!!$      !  Integral Constraints Properties for momentum and/or tracers trends 
     244!!$      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
     245!!$      IF( ln_glo_trd )   CALL trd_glo( putrd, pvtrd, ktrd, 'DYN', kt ) 
     246!!$ 
     247!!$      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
     248!!$      !  Kinetic Energy trends 
     249!!$      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
     250!!$      IF( ln_KE_trd  )   CALL trd_ken( putrd, pvtrd, ktrd, kt ) 
     251!!$ 
     252!!$      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
     253!!$      !  Vorticity trends 
     254!!$      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
     255!!$      IF( ln_vor_trd )   CALL trd_vor( putrd, pvtrd, ktrd, kt ) 
    86256 
    87257      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
     
    90260!!gm      IF( ln_dyn_mxl )   CALL trd_mxl_dyn    
    91261      ! 
    92    END SUBROUTINE trd_dyn 
    93  
    94  
    95    SUBROUTINE trd_dyn_iom( putrd, pvtrd, ktrd, kt, Kmm ) 
     262   END SUBROUTINE trd_dyn_2d 
     263 
     264 
     265   SUBROUTINE trd_dyn_iom_3d( putrd, pvtrd, ktrd, kt, Kmm ) 
    96266      !!--------------------------------------------------------------------- 
    97267      !!                  ***  ROUTINE trd_dyn_iom  *** 
     
    113283      CASE( jpdyn_hpg )   ;   CALL iom_put( "utrd_hpg", putrd )    ! hydrostatic pressure gradient 
    114284                              CALL iom_put( "vtrd_hpg", pvtrd ) 
    115       CASE( jpdyn_spg )   ;   CALL iom_put( "utrd_spg", putrd )    ! surface pressure gradient 
    116                               CALL iom_put( "vtrd_spg", pvtrd ) 
    117285      CASE( jpdyn_pvo )   ;   CALL iom_put( "utrd_pvo", putrd )    ! planetary vorticity 
    118286                              CALL iom_put( "vtrd_pvo", pvtrd ) 
     
    128296                                 z3dy(ji,jj,jk) = vv(ji,jj,jk,Kmm) * ( vv(ji,jj+1,jk,Kmm) - vv(ji,jj-1,jk,Kmm) ) / ( 2._wp * e2v(ji,jj) ) 
    129297                              END_3D 
    130                               CALL lbc_lnk( 'trddyn', z3dx, 'U', -1.0_wp, z3dy, 'V', -1.0_wp ) 
     298                              CALL lbc_lnk_multi( 'trddyn', z3dx, 'U', -1., z3dy, 'V', -1. ) 
    131299                              CALL iom_put( "utrd_udx", z3dx  ) 
    132300                              CALL iom_put( "vtrd_vdy", z3dy  ) 
     
    141309                              !                                    ! wind stress trends 
    142310                              ALLOCATE( z2dx(jpi,jpj) , z2dy(jpi,jpj) ) 
    143                               z2dx(:,:) = ( utau_b(:,:) + utau(:,:) ) / ( e3u(:,:,1,Kmm) * rho0 ) 
    144                               z2dy(:,:) = ( vtau_b(:,:) + vtau(:,:) ) / ( e3v(:,:,1,Kmm) * rho0 ) 
     311                              z2dx(:,:) = ( utau_b(:,:) + utau(:,:) ) / ( e3u_n(:,:,1) * rau0 ) 
     312                              z2dy(:,:) = ( vtau_b(:,:) + vtau(:,:) ) / ( e3v_n(:,:,1) * rau0 ) 
    145313                              CALL iom_put( "utrd_tau", z2dx ) 
    146314                              CALL iom_put( "vtrd_tau", z2dy ) 
    147315                              DEALLOCATE( z2dx , z2dy ) 
    148 !!gm  to be changed : computation should be done in dynzdf.F90 
    149 !!gm                + missing the top friction  
    150 !                              !                                    ! bottom stress tends (implicit case) 
    151 !                              IF( ln_drgimp ) THEN 
    152 !                                 ALLOCATE( z3dx(jpi,jpj,jpk) , z3dy(jpi,jpj,jpk) ) 
    153 !                             z3dx(:,:,:) = 0._wp   ;   z3dy(:,:,:) = 0._wp  ! after velocity known (now filed at this stage) 
    154 !                            DO jk = 1, jpkm1 
    155 !                                    DO jj = 2, jpjm1 
    156 !                                       DO ji = 2, jpim1 
    157 !                                      ikbu = mbku(ji,jj)          ! deepest ocean u- & v-levels 
    158 !                                          ikbv = mbkv(ji,jj) 
    159 !                                          z3dx(ji,jj,jk) = 0.5 * ( rCdU_bot(ji+1,jj) + rCdU_bot(ji,jj) ) &  
    160 !                                               &         * uu(ji,jj,ikbu,Kmm) / e3u(ji,jj,ikbu,Kmm) 
    161 !                                          z3dy(ji,jj,jk) = 0.5 * ( rCdU_bot(ji,jj+1) + rCdU_bot(ji,jj) ) & 
    162 !                                               &         * vv(ji,jj,ikbv,Kmm) / e3v(ji,jj,ikbv,Kmm) 
    163 !                                    END DO 
    164 !                                 END DO 
    165 !                              END DO 
    166 !                              CALL lbc_lnk( 'trddyn', z3dx, 'U', -1.0_wp, z3dy, 'V', -1.0_wp ) 
    167 !                              CALL iom_put( "utrd_bfr", z3dx ) 
    168 !                              CALL iom_put( "vtrd_bfr", z3dy ) 
    169 !                                 DEALLOCATE( z3dx , z3dy ) 
    170 !                              ENDIF 
    171 !!gm end 
    172       CASE( jpdyn_bfr )       ! called if ln_drgimp=F 
    173                               CALL iom_put( "utrd_bfr", putrd )    ! bottom friction (explicit case) 
     316      CASE( jpdyn_bfr )   ;   CALL iom_put( "utrd_bfr", putrd )    ! bottom friction (explicit case) 
    174317                              CALL iom_put( "vtrd_bfr", pvtrd ) 
     318      CASE( jpdyn_bfri)   ;   CALL iom_put( "utrd_bfri", putrd )    ! bottom friction (implicit case) 
     319                              CALL iom_put( "vtrd_bfri", pvtrd ) 
    175320      CASE( jpdyn_atf )   ;   CALL iom_put( "utrd_atf", putrd )        ! asselin filter trends  
    176321                              CALL iom_put( "vtrd_atf", pvtrd ) 
    177322      END SELECT 
    178323      ! 
    179    END SUBROUTINE trd_dyn_iom 
     324   END SUBROUTINE trd_dyn_iom_3d 
     325 
     326 
     327   SUBROUTINE trd_dyn_iom_2d( putrd, pvtrd, ktrd, kt, Kmm ) 
     328      !!--------------------------------------------------------------------- 
     329      !!                  ***  ROUTINE trd_dyn_iom  *** 
     330      !!  
     331      !! ** Purpose :   output 2D trends using IOM 
     332      !!---------------------------------------------------------------------- 
     333      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   putrd, pvtrd   ! U and V trends 
     334      INTEGER                 , INTENT(in   ) ::   ktrd           ! trend index 
     335      INTEGER                 , INTENT(in   ) ::   kt             ! time step 
     336      INTEGER                 , INTENT(in   ) ::   Kmm            ! time level index 
     337      ! 
     338      INTEGER ::   ji, jj, jk   ! dummy loop indices 
     339      INTEGER ::   ikbu, ikbv   ! local integers 
     340      REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   z2dx, z2dy   ! 2D workspace  
     341      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   z3dx, z3dy   ! 3D workspace  
     342      !!---------------------------------------------------------------------- 
     343      ! 
     344      SELECT CASE( ktrd ) 
     345      CASE( jpdyn_spg )      ;   CALL iom_put( "utrd_spg2d", putrd )      ! surface pressure gradient 
     346                                 CALL iom_put( "vtrd_spg2d", pvtrd ) 
     347      CASE( jpdyn_pvo )      ;   CALL iom_put( "utrd_pvo2d", putrd )      ! planetary vorticity (barotropic part) 
     348                                 CALL iom_put( "vtrd_pvo2d", pvtrd ) 
     349      CASE( jpdyn_hpg_corr ) ;   CALL iom_put( "utrd_hpg_corr", putrd )   ! horizontal pressure gradient correction 
     350                                 CALL iom_put( "vtrd_hpg_corr", pvtrd ) 
     351      CASE( jpdyn_pvo_corr ) ;   CALL iom_put( "utrd_pvo_corr", putrd )   ! planetary vorticity correction 
     352                                 CALL iom_put( "vtrd_pvo_corr", pvtrd ) 
     353      CASE( jpdyn_bfr )      ;   CALL iom_put( "utrd_bfr2d", putrd )      ! bottom friction due to barotropic currents 
     354                                 CALL iom_put( "vtrd_bfr2d", pvtrd ) 
     355      END SELECT 
     356      ! 
     357   END SUBROUTINE trd_dyn_iom_2d 
    180358 
    181359   !!====================================================================== 
Note: See TracChangeset for help on using the changeset viewer.