Changeset 8112


Ignore:
Timestamp:
2017-06-01T13:24:36+02:00 (3 years ago)
Author:
davestorkey
Message:

UKMO/branch/dev_r5518_GO6_package_FVPS: FVPS implementation

Location:
branches/UKMO/dev_r5518_GO6_package_FVPS/NEMOGCM
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • branches/UKMO/dev_r5518_GO6_package_FVPS/NEMOGCM/CONFIG/SHARED/field_def.xml

    r7393 r8112  
    8484         <field id="beta"         long_name="haline contraction"                                                        unit="1e3"    grid_ref="grid_T_3D" /> 
    8585         <field id="rhop"         long_name="potential density (sigma0)"        standard_name="sea_water_sigma_theta"   unit="kg/m3"  grid_ref="grid_T_3D" /> 
     86         <field id="rho"          long_name="in sity density"                                                           unit="kg/m3"  grid_ref="grid_T_3D" /> 
    8687 
    8788         <!-- Energy - horizontal divergence --> 
     
    415416         <field id="uocet"        long_name="ocean transport along i-axis times temperature (CRS)"                                               unit="degree_C*m/s"   grid_ref="grid_U_3D" /> 
    416417         <field id="uoces"        long_name="ocean transport along i-axis times salinity (CRS)"                                                  unit="0.001*m/s"   grid_ref="grid_U_3D" /> 
     418         <field id="hpgi"         long_name="horizontal pressure gradient along i-axis"                                                          unit="N/kg"        grid_ref="grid_U_3D" /> 
    417419 
    418420         <!-- variables available with MLE --> 
     
    459461         <field id="vocet"        long_name="ocean transport along j-axis times temperature (CRS)"                                               unit="degree_C*m/s"   grid_ref="grid_V_3D" /> 
    460462         <field id="voces"        long_name="ocean transport along j-axis times salinity (CRS)"                                                  unit="0.001*m/s"   grid_ref="grid_V_3D" /> 
     463         <field id="hpgj"         long_name="horizontal pressure gradient along j-axis"                                                          unit="N/kg"        grid_ref="grid_V_3D" /> 
    461464 
    462465         <!-- variables available with MLE --> 
     
    492495        <field id="woce"         long_name="ocean vertical velocity"              standard_name="upward_sea_water_velocity"   unit="m/s"  /> 
    493496        <field id="wocetr_eff"   long_name="effective ocean vertical transport"                                               unit="m3/s" /> 
     497        <field id="hpgdepth"     long_name="depth from HPG module" unit="m"  /> 
     498        <field id="pressure"     long_name="in situ pressure"      unit="N/m2"  /> 
    494499 
    495500        <!-- woce_eiv: available with key_traldf_eiv and key_diaeiv --> 
  • branches/UKMO/dev_r5518_GO6_package_FVPS/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90

    r6487 r8112  
    11331133         e3vw_0(:,:,jk) = e3w_1d(jk) 
    11341134      END DO 
    1135       DO jk = 1,jpk                         ! Computed as the minimum of neighbooring scale factors 
     1135      DO jk = 1,jpk                         ! Computed as the average of neighbouring scale factors 
    11361136         DO jj = 1, jpjm1 
    11371137            DO ji = 1, fs_jpim1   ! vector opt. 
    1138                e3u_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji+1,jj,jk) ) 
    1139                e3v_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji,jj+1,jk) ) 
    1140                e3uw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji+1,jj,jk) ) 
    1141                e3vw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji,jj+1,jk) ) 
     1138               e3u_0 (ji,jj,jk) = 0.5 * ( e3t_0(ji,jj,jk) + e3t_0(ji+1,jj,jk) ) 
     1139               e3v_0 (ji,jj,jk) = 0.5 * ( e3t_0(ji,jj,jk) + e3t_0(ji,jj+1,jk) ) 
     1140               e3uw_0(ji,jj,jk) = 0.5 * ( e3w_0(ji,jj,jk) + e3w_0(ji+1,jj,jk) ) 
     1141               e3vw_0(ji,jj,jk) = 0.5 * ( e3w_0(ji,jj,jk) + e3w_0(ji,jj+1,jk) ) 
    11421142            END DO 
    11431143         END DO 
     
    11731173         e3f_0(:,:,jk) = e3t_1d(jk) 
    11741174      END DO 
    1175       DO jk = 1, jpk                        ! Computed as the minimum of neighbooring V-scale factors 
     1175      DO jk = 1, jpk                        ! Computed as the average of neighbouring V-scale factors 
    11761176         DO jj = 1, jpjm1 
    11771177            DO ji = 1, fs_jpim1   ! vector opt. 
    1178                e3f_0(ji,jj,jk) = MIN( e3v_0(ji,jj,jk), e3v_0(ji+1,jj,jk) ) 
     1178               e3f_0(ji,jj,jk) = 0.5 * ( e3v_0(ji,jj,jk) + e3v_0(ji+1,jj,jk) ) 
    11791179            END DO 
    11801180         END DO 
  • branches/UKMO/dev_r5518_GO6_package_FVPS/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90

    r6486 r8112  
    4444   USE wrk_nemo        ! Memory Allocation 
    4545   USE timing          ! Timing 
     46   USE iom 
    4647 
    4748   IMPLICIT NONE 
     
    379380      !! 
    380381      INTEGER  ::   ji, jj, jk                 ! dummy loop indices 
    381       REAL(wp) ::   zcoef0, zuap, zvap, znad   ! temporary scalars 
    382       REAL(wp), POINTER, DIMENSION(:,:,:) ::  zhpi, zhpj 
    383       !!---------------------------------------------------------------------- 
    384       ! 
    385       CALL wrk_alloc( jpi,jpj,jpk, zhpi, zhpj ) 
     382      REAL(wp) ::   zcoef0, zuap, zvap, znad, zvnum   ! temporary scalars 
     383      REAL(wp), POINTER, DIMENSION(:,:)   ::  zunum, zudenom    ! temporary numerators for u and v velocities (for debugging)  
     384      REAL(wp), POINTER, DIMENSION(:,:)   ::  zpb, zpt   ! pressure at top and bottom of cells 
     385      REAL(wp), POINTER, DIMENSION(:,:)   ::  zzb, zzt   ! geopotential height (gz) at top and bottom of cells 
     386      REAL(wp), POINTER, DIMENSION(:,:,:) ::  zhpi, zhpj, zp3d, zz3d 
     387      !!---------------------------------------------------------------------- 
     388      ! 
     389      CALL wrk_alloc( jpi,jpj,jpk, zhpi, zhpj, zp3d, zz3d ) 
     390      CALL wrk_alloc( jpi,jpj, zunum, zudenom, zpb, zpt, zzb, zzt ) 
    386391      ! 
    387392      IF( kt == nit000 ) THEN 
     
    398403      ENDIF 
    399404 
    400       ! Surface value 
    401       DO jj = 2, jpjm1 
    402          DO ji = fs_2, fs_jpim1   ! vector opt. 
    403             ! hydrostatic pressure gradient along s-surfaces 
    404             zhpi(ji,jj,1) = zcoef0 / e1u(ji,jj) * ( fse3w(ji+1,jj  ,1) * ( znad + rhd(ji+1,jj  ,1) )   & 
    405                &                                  - fse3w(ji  ,jj  ,1) * ( znad + rhd(ji  ,jj  ,1) ) ) 
    406             zhpj(ji,jj,1) = zcoef0 / e2v(ji,jj) * ( fse3w(ji  ,jj+1,1) * ( znad + rhd(ji  ,jj+1,1) )   & 
    407                &                                  - fse3w(ji  ,jj  ,1) * ( znad + rhd(ji  ,jj  ,1) ) ) 
    408             ! s-coordinate pressure gradient correction 
    409             zuap = -zcoef0 * ( rhd   (ji+1,jj,1) + rhd   (ji,jj,1) + 2._wp * znad )   & 
    410                &           * ( fsde3w(ji+1,jj,1) - fsde3w(ji,jj,1) ) / e1u(ji,jj) 
    411             zvap = -zcoef0 * ( rhd   (ji,jj+1,1) + rhd   (ji,jj,1) + 2._wp * znad )   & 
    412                &           * ( fsde3w(ji,jj+1,1) - fsde3w(ji,jj,1) ) / e2v(ji,jj) 
    413             ! add to the general momentum trend 
    414             ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) + zuap 
    415             va(ji,jj,1) = va(ji,jj,1) + zhpj(ji,jj,1) + zvap 
    416          END DO 
    417       END DO 
    418  
    419       ! interior value (2=<jk=<jpkm1) 
    420       DO jk = 2, jpkm1 
     405! set the pressure values to zero at the top of the model (the atmospheric surface pressure is included in the barotropic step ?) 
     406      zpb(:,:) = 0._wp 
     407      zzb(:,:) = grav * sshn(:,:)     ! surface height had been incorrectly set to zero; corrected on 24 Apr 2017 
     408 
     409      if(lwp) write(numout,*) 'hpg_sco: grav is ',grav 
     410 
     411      DO jk = 1, jpkm1 
     412 
     413         zpt(:,:) = zpb(:,:) 
     414         zzt(:,:) = zzb(:,:)  
     415 
     416         zp3d(:,:,jk) = zpt(:,:) 
     417         zz3d(:,:,jk) = zzt(:,:) 
     418 
     419         ! integrate the pressure from the top of the cell to the bottom of the cell; a factor rho_0 has been removed from the pressure   
     420         ! corrected to calculate for loops from 1 to jpj and 1 to jpi 25 Apr 2017  
     421         DO jj = 1, jpj 
     422            DO ji = 1, jpi   ! vector opt. 
     423!!$               zpb(ji,jj) = zpt(ji,jj) + grav * fse3w_n(ji,jj,jk) * ( znad + rhd(ji,jj,jk) )     
     424               zpb(ji,jj) = zpt(ji,jj) + grav * fse3t_n(ji,jj,jk) * ( znad + rhd(ji,jj,jk) ) 
     425               zzb(ji,jj) = grav * ( sshn(ji,jj) - fsdepw_n(ji,jj,jk+1) )   ! surface height had been set to zero; corrected on 24 Apr 2017 
     426            END DO 
     427         END DO 
     428 
     429         ! calculate the Jacobian of the pressure forces 
    421430         DO jj = 2, jpjm1 
    422431            DO ji = fs_2, fs_jpim1   ! vector opt. 
    423                ! hydrostatic pressure gradient along s-surfaces 
    424                zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + zcoef0 / e1u(ji,jj)   & 
    425                   &           * (  fse3w(ji+1,jj,jk) * ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) + 2*znad )   & 
    426                   &              - fse3w(ji  ,jj,jk) * ( rhd(ji  ,jj,jk) + rhd(ji  ,jj,jk-1) + 2*znad )  ) 
    427                zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) + zcoef0 / e2v(ji,jj)   & 
    428                   &           * (  fse3w(ji,jj+1,jk) * ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) + 2*znad )   & 
    429                   &              - fse3w(ji,jj  ,jk) * ( rhd(ji,jj,  jk) + rhd(ji,jj  ,jk-1) + 2*znad )  ) 
    430                ! s-coordinate pressure gradient correction 
    431                zuap = -zcoef0 * ( rhd   (ji+1,jj  ,jk) + rhd   (ji,jj,jk) + 2._wp * znad )   & 
    432                   &           * ( fsde3w(ji+1,jj  ,jk) - fsde3w(ji,jj,jk) ) / e1u(ji,jj) 
    433                zvap = -zcoef0 * ( rhd   (ji  ,jj+1,jk) + rhd   (ji,jj,jk) + 2._wp * znad )   & 
    434                   &           * ( fsde3w(ji  ,jj+1,jk) - fsde3w(ji,jj,jk) ) / e2v(ji,jj) 
    435                ! add to the general momentum trend 
    436                ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk) + zuap 
    437                va(ji,jj,jk) = va(ji,jj,jk) + zhpj(ji,jj,jk) + zvap 
     432               zunum(ji, jj) = ( zpb(ji+1,jj) - zpt(ji  , jj) )*( zzb(ji  , jj) - zzt(ji+1, jj) ) - ( zpb(ji  , jj) - zpt(ji+1, jj) )*( zzb(ji+1, jj) - zzt(ji , jj) )    
     433               zudenom(ji, jj) =  zpb(ji+1, jj) - zpt(ji  , jj) + zpb(ji  , jj) - zpt(ji+1, jj)  
     434               zhpi(ji, jj, jk) = zunum(ji, jj) * r1_e1u(ji, jj) / zudenom(ji, jj)    
     435               ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk) 
     436 
     437               zvnum = ( zpb(ji,jj+1) - zpt(ji, jj  ) )*( zzb(ji, jj  ) - zzt(ji, jj+1) ) - ( zpb(ji, jj  ) - zpt(ji, jj+1) )*( zzb(ji, jj+1) - zzt(ji, jj  ) )    
     438               zhpj(ji, jj, jk) = zvnum * r1_e2v(ji, jj) / ( zpb(ji, jj+1) - zpt(ji, jj  ) + zpb(ji, jj  ) - zpt(ji, jj+1) )    
     439               va(ji,jj,jk) = va(ji,jj,jk) + zhpj(ji,jj,jk) 
    438440            END DO 
    439441         END DO 
    440       END DO 
    441       ! 
    442       CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zhpj ) 
     442 
     443! temporary diagnostic outputs (to be used on just a few timesteps on a small configuration)  
     444!!$         DO jj = 1, jpj 
     445!!$            write(numout, *) 'fse3w_n,   jk, jj = ', jk, jj, (fse3w_n(ji,jj,jk), ji = 1, jpi) 
     446!!$            write(numout, *) 'rhd,     jk, jj = ', jk, jj, (rhd(ji,jj,jk), ji = 1, jpi) 
     447!!$            write(numout, *) 'gdepw_n, jk+1, jj = ', jk+1, jj, (gdepw_n(ji,jj,jk+1), ji = 1, jpim1) 
     448!!$            write(numout, *) 'r1_e1u,  jk, jj = ', jk, jj, (r1_e1u(ji,jj), ji = 1, jpim1) 
     449!!$            write(numout, *) 'zpt,     jk, jj = ', jk, jj, (zpt(ji,jj), ji = 1, jpi) 
     450!!$            write(numout, *) 'zpb,     jk, jj = ', jk, jj, (zpb(ji,jj), ji = 1, jpi) 
     451!!$            write(numout, *) 'zzt,     jk, jj = ', jk, jj, (zzt(ji,jj), ji = 1, jpi) 
     452!!$            write(numout, *) 'zzb,     jk, jj = ', jk, jj, (zzb(ji,jj), ji = 1, jpi) 
     453!!$            write(numout, *) 'zunum,   jk, jj = ', jk, jj, (zunum(ji,jj), ji = 1, jpi) 
     454!!$            write(numout, *) 'zudenom, jk, jj = ', jk, jj, (zudenom(ji,jj), ji = 1, jpi) 
     455!!$            write(numout, *) 'tmask,   jk, jj = ', jk, jj, (tmask(ji,jj,jk), ji = 1, jpi) 
     456!!$            write(numout, *) 'umask,   jk, jj = ', jk, jj, (umask(ji,jj,jk), ji = 1, jpi) 
     457!!$            write(numout, *) 'vmask,   jk, jj = ', jk, jj, (vmask(ji,jj,jk), ji = 1, jpi) 
     458!!$ 
     459!!$            write(numout, *) 'zhpi,    jk, jj = ', jk, jj, (zhpi(ji,jj,jk), ji = 1, jpi) 
     460!!$            write(numout, *) 'zhpj,    jk, jj = ', jk, jj, (zhpj(ji,jj,jk), ji = 1, jpi) 
     461!!$         END DO  
     462 
     463      END DO 
     464      ! 
     465      CALL iom_put("pressure" , zp3d(:,:,:) ) 
     466      CALL iom_put("hpgdepth" , zz3d(:,:,:) ) 
     467      CALL iom_put("rho" , rhd(:,:,:) ) 
     468      CALL iom_put("hpgi", zhpi(:,:,:) ) 
     469      CALL iom_put("hpgj", zhpj(:,:,:) ) 
     470      ! 
     471      CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zhpj, zp3d, zz3d ) 
     472      CALL wrk_dealloc( jpi,jpj, zpb, zpt, zzb, zzt ) 
    443473      ! 
    444474   END SUBROUTINE hpg_sco 
Note: See TracChangeset for help on using the changeset viewer.