ROMS-TOPAZ column_physics 접합

HCJung-jbnu commented 1 year ago

BIO source/sink 이슈

FENNEL :: OCEAN%t(i, j, k, nnew, itrc) 에서 nnew에 increment*Hz 로 전달 -> TOPAZ도 똑같이...

    DO itrc=1,NBT
      ibio=idbio(itrc)
      DO k=1,N(ng)
        DO i=Istr,Iend
          Bio_old(i,k,ibio)=MAX(0.0_r8,t(i,j,k,nstp,ibio))
          Bio(i,k,ibio)=Bio_old(i,k,ibio)
        END DO
      END DO
    END DO
.....
!-----------------------------------------------------------------------
!  Update global tracer variables: Add increment due to BGC processes
!  to tracer array in time index "nnew". Index "nnew" is solution after
!  advection and mixing and has transport units (m Tunits) hence the
!  increment is multiplied by Hz.  Notice that we need to subtract
!  original values "Bio_old" at the top of the routine to just account
!  for the concentractions affected by BGC processes. This also takes
!  into account any constraints (non-negative concentrations, carbon
!  concentration range) specified before entering BGC kernel. If "Bio"
!  were unchanged by BGC processes, the increment would be exactly
!  zero. Notice that final tracer values, t(:,:,:,nnew,:) are not
!  bounded >=0 so that we can preserve total inventory of N and
!  C even when advection causes tracer concentration to go negative.
!  (J. Wilkin and H. Arango, Apr 27, 2012)  :: nnew 는 tendency*dt
!-----------------------------------------------------------------------
    DO itrc=1,NBT
      ibio=idbio(itrc)
      DO k=1,N(ng)
        DO i=Istr,Iend
          cff=Bio(i,k,ibio)-Bio_old(i,k,ibio)
#ifdef MASKING
          cff=cff*rmask(i,j)
# ifdef WET_DRY
          cff=cff*rmask_wet(i,j)
# endif
#endif
          t(i,j,k,nnew,ibio)=t(i,j,k,nnew,ibio)+cff*Hz(i,j,k) <-----TOPAZ에 적용
        END DO
      END DO
    END DO
  END DO J_LOOP

cobalt source/sink 이후 "Set value for boundaries" 는 무엇?

! Set value for boundaries <-- what's this?
DO it=1,NOBGC
 CALL bc_r3d_tile (ng, tile,                                     &
&                      LBi, UBi, LBj, UBj, 1, UBk,                 &
&                      obgc(:,:,:,nnew,it))
END DO

HCJung-jbnu commented 1 year ago

topaz.h : SUB. biology


       DO k = 0, N(ng) !!In TOPAZ, 0=TOP, N(ng)=BOTTOM, as opposed to ROMS
       DO j = Jstr,Jend
       DO i = Istr,Iend
         TP(ng)%depth_w(i,j,k)=GRID(ng)%z_w(i,j,N(ng)-k)
       END DO
       END DO
       END DO

       DO k = 1, N(ng)
       DO j = Jstr,Jend
       DO i = Istr,Iend
         TP(ng)%temp(i,j,k)=OCEAN(ng)%t(i,j,N(ng)-k+1,nstp(ng),itemp)
         TP(ng)%salt(i,j,k)=OCEAN(ng)%t(i,j,N(ng)-k+1,nstp(ng),isalt)
         TP(ng)%rho(i,j,k)=OCEAN(ng)%rho(i,j,N(ng)-k+1)+1000.0 !! rho anomaly + 1000.0
         TP(ng)%dzt(i,j,k)=GRID(ng)%Hz(i,j,N(ng)-k+1) !!rho point vertical thickness(m)
!         TP(ng)%dzt(i,j,k)=GRID(ng)%z_w(i,j,N(ng)-k+1)-GRID(ng)%z_w(i,j,N(ng)-k)
       END DO
       END DO
       END DO

       DO j = Jstr,Jend
       DO i = Istr,Iend
         TP(ng)%dzwt(i,j,0)=GRID(ng)%Hz(i,j,N(ng))/2.0
         TP(ng)%dzwt(i,j,N(ng))=GRID(ng)%Hz(i,j,1)/2.0
!         TP(ng)%dzwt(i,j,0)=GRID(ng)%z_w(i,j,N(ng))-GRID(ng)%z_r(i,j,N(ng))
!         TP(ng)%dzwt(i,j,N(ng))=(GRID(ng)%z_w(i,j,N(ng))-GRID(ng)%z_w(i,j,0))-(GRID(ng)%z_w(i,j,N(ng))-GRID(ng)%z_r(i,j,1))
       DO k = 1, N(ng)-1
          TP(ng)%dzwt(i,j,k)=(GRID(ng)%Hz(i,j,N(ng)-k+1)/2.0)+(GRID(ng)%Hz(i,j,N(ng)-k)/2.0)
!         TP(ng)%dzwt(i,j,k)=(GRID(ng)%z_w(i,j,N(ng))-GRID(ng)%z_r(i,j,N(ng)-k)) - (GRID(ng)%z_w(i,j,N(ng))-GRID(ng)%z_r(i,j,N(ng)-k+1))
       END DO
       END DO
       END DO

       cff1 = rho0 * 550.0_r8
       DO j=Jstr,Jend
       DO i=Istr,Iend
         TP(ng)%qsr(i,j)=FORCES(ng)%srflx(i,j)*rho0*Cp !!degC m/s -> Watts/m2
         TP(ng)%visfrac(i,j)=1.0-lmd_r1(INT(MIXING(ng)%Jwtype(i,j))) ! define visible and IR fraction according to variable water type
         TP(ng)%k_ir(i,j)=1.0/lmd_mu1(INT(MIXING(ng)%Jwtype(i,j))) ! define InfraRed decay scale according to variable water type
         TP(ng)%hmlpt(i,j)=mxl_depth(i,j)
         TP(ng)%apr(i,j)=FORCES(ng)%Pair(i,j)*100.0 !!mb*100.0 = Pa
#ifdef BULK_FLUXES
         TP(ng)%wndm(i,j)=SQRT(FORCES(ng)%Uwind(i,j)**2+FORCES(ng)%Vwind(i,j)**2) !!m/s
#else
         TP(ng)%wndm(i,j)=SQRT(cff1*SQRT((0.5_r8*(FORCES(ng)%sustr(i,j)+FORCES(ng)%sustr(i+1,j)))**2+  &
                               (0.5_r8*(FORCES(ng)%svstr(i,j)+FORCES(ng)%svstr(i,j+1)))**2))
#endif
         TP(ng)%flux_fed(i,j)=FORCES(ng)%soluble_fe(i,j)
         TP(ng)%flux_lith(i,j)=FORCES(ng)%mineral_fe(i,j)
         TP(ng)%flux_no3_wet(i,j)=FORCES(ng)%wetdep_no3(i,j)
         TP(ng)%flux_no3_dry(i,j)=FORCES(ng)%drydep_no3(i,j)
         TP(ng)%flux_nh4_wet(i,j)=FORCES(ng)%wetdep_nh4(i,j)
         TP(ng)%flux_nh4_dry(i,j)=FORCES(ng)%drydep_nh4(i,j)
       END DO
       END DO

       DO ntrc = 3, NT(ng) !!ntrc=1=temp, =2=salt
       DO k = 1, N(ng)
       DO j = Jstr,Jend
       DO i = Istr,Iend
         TP(ng)%prog_trc(i,j,k,ntrc-2)=OCEAN(ng)%t(i,j,N(ng)-k+1,nstp(ng),ntrc)
         TP(ng)%prog_trc_old(i,j,k,ntrc-2)=OCEAN(ng)%t(i,j,N(ng)-k+1,nstp(ng),ntrc)
       END DO
       END DO
       END DO
       END DO
       DEALLOCATE(mxl_depth)

       CALL topaz_optic(ng, rho0)

       CALL generic_TOPAZ_column_physics(ng, tile, Istr, Iend, Jstr, Jend, nstep_bio, dt(ng))

       DO ntrc = 3, NT(ng) !!ntrc=1=temp, =2=salt
       DO k = 1, N(ng)
       DO j = Jstr,Jend
       DO i = Istr,Iend
         cff = TP(ng)%prog_trc(i,j,k,ntrc-2) - TP(ng)%prog_trc_old(i,j,k,ntrc-2)
#ifdef MASKING
         cff=cff*GRID(ng)%rmask(i,j)
# ifdef WET_DRY
         cff=cff*GRID(ng)%rmask_wet(i,j)
# endif
#endif
         OCEAN(ng)%t(i,j,N(ng)-k+1,nnew(ng),ntrc)=OCEAN(ng)%t(i,j,N(ng)-k+1,nnew(ng),ntrc) + cff * GRID(ng)%Hz(i,j,N(ng)-k+1)
       END DO
       END DO
       END DO
       END DO

       DO ntrc = 1, NDbio3d
       DO k = 1, N(ng)
       DO j = Jstr,Jend
       DO i = Istr,Iend
#ifdef MASKING
         DIAGS(ng)%DiaBio3d(i,j,N(ng)-k+1,ntrc)=TP(ng)%diag_trc(i,j,k,ntrc) * GRID(ng)%rmask(i,j)
# ifdef WET_DRY
         DIAGS(ng)%DiaBio3d(i,j,N(ng)-k+1,ntrc)=TP(ng)%diag_trc(i,j,k,ntrc) * GRID(ng)%rmask_full(i,j)
# endif
#else
         DIAGS(ng)%DiaBio3d(i,j,N(ng)-k+1,ntrc)=TP(ng)%diag_trc(i,j,k,ntrc)
#endif
       END DO
       END DO
       END DO
       END DO
       DO j = Jstr,Jend
       DO i = Istr,Iend
         DIAGS(ng)%DiaBio2d(i,j,1)=TP(ng)%co2_flux_alpha(i,j)
         DIAGS(ng)%DiaBio2d(i,j,2)=TP(ng)%co2_flux_sc_no(i,j)
         DIAGS(ng)%DiaBio2d(i,j,3)=TP(ng)%co2_flux_csurf(i,j)
         DIAGS(ng)%DiaBio2d(i,j,4)=TP(ng)%o2_flux_alpha(i,j)
         DIAGS(ng)%DiaBio2d(i,j,5)=TP(ng)%o2_flux_sc_no(i,j)
         DIAGS(ng)%DiaBio2d(i,j,6)=TP(ng)%o2_flux_csurf(i,j)
       END DO
       END DO

HCJung-jbnu commented 1 year ago

generic_TOPAZ_mode.F90 : SUB. topaz_optic


      do i = isc, iec; do j = jsc, jec
            f_vis(i,j)=TP(ng)%visfrac(i,j)
            f_ir(i,j)=1.0-TP(ng)%visfrac(i,j)
      end do; end do

      do n= 1, num_diag_tracers
        if ( T_diag(n)%name .eq. 'chl' ) then
          do i = isc, iec; do j = jsc, jec; do k = 1, Grids%nk
            chl(i,j,k)=TP(ng)%diag_trc(i,j,k,n-1)*rho0/1000.0 !!ug/kg -> mg/m3
          end do; end do; end do
        endif
      end do

      do i = isc, iec
      do j = jsc, jec

        TP(ng)%decayW(i,j,0)=1.0 !! k=1=TOP=surface W, k=Grids%nk=BOTTOM W

        k=1
        vis_irr(i,j)=TP(ng)%qsr(i,j)*f_vis(i,j)
        red(i,j,k)=0.5*exp(-TP(ng)%dzt(i,j,k)*(0.225+(0.037*chl(i,j,k)**0.629)))
        blue(i,j,k)=0.5*exp(-TP(ng)%dzt(i,j,k)*(0.0232+(0.074*chl(i,j,k)**0.674)))

        sw_frac_zw(i,j,k)=f_vis(i,j)*(red(i,j,k)+blue(i,j,k))
        sw_frac_zt(i,j,k)=0.5*(f_vis(i,j)+sw_frac_zw(i,j,k))
        INTP(ng)%opacity(i,j,k)=-log(sw_frac_zw(i,j,k)/f_vis(i,j))/TP(ng)%dzt(i,j,k)

        INTP(ng)%par(i,j,k)=vis_irr(i,j)*exp(-INTP(ng)%opacity(i,j,k)*TP(ng)%dzt(i,j,k)*0.5)
        vis_irr(i,j)=vis_irr(i,j)*exp(-INTP(ng)%opacity(i,j,k)*TP(ng)%dzt(i,j,k))

        TP(ng)%decayW(i,j,k)=(f_ir(i,j)*exp(-TP(ng)%k_ir(i,j)*ABS(TP(ng)%depth_w(i,j,k))))+red(i,j,k)+blue(i,j,k)

        do k = 2, Grids%nk
          red(i,j,k)=red(i,j,k-1)*exp(-TP(ng)%dzt(i,j,k)*(0.2250+(0.037*chl(i,j,k)**0.629)))
          blue(i,j,k)=blue(i,j,k-1)*exp(-TP(ng)%dzt(i,j,k)*(0.0232+(0.074*chl(i,j,k)**0.674)))

          sw_frac_zw(i,j,k)=f_vis(i,j)*(red(i,j,k)+blue(i,j,k))  !!<---- ROMS-optic TOTAL decayW(i,j,k,1)
          sw_frac_zt(i,j,k)=0.5*(sw_frac_zw(i,j,k-1)+sw_frac_zw(i,j,k))
          INTP(ng)%opacity(i,j,k)=-log(sw_frac_zw(i,j,k)/sw_frac_zw(i,j,k-1))/TP(ng)%dzt(i,j,k)

          if ( ABS(TP(ng)%depth_w(i,j,k-1)) <= 200.0 ) then !!mom:200m, nemo:400m, aphotic zone depth
            INTP(ng)%par(i,j,k)=vis_irr(i,j)*exp(-INTP(ng)%opacity(i,j,k)*TP(ng)%dzt(i,j,k)*0.5)
            vis_irr(i,j)=vis_irr(i,j)*exp(-INTP(ng)%opacity(i,j,k)*TP(ng)%dzt(i,j,k))
            TP(ng)%decayW(i,j,k)=(f_ir(i,j)*exp(-TP(ng)%k_ir(i,j)*ABS(TP(ng)%depth_w(i,j,k))))+red(i,j,k)+blue(i,j,k)
          else
            INTP(ng)%par(i,j,k)=0.0
            vis_irr(i,j)=0.0
            TP(ng)%decayW(i,j,k)=0.0
          end if

      end do
      end do
      end do

      deallocate(chl, sw_frac_zt, sw_frac_zw, red, blue, f_vis, f_ir, vis_irr)

HCJung-jbnu commented 1 year ago

pre_step3d.F : SUB. pre_step3d_tile

.......
#   ifdef TOPAZ
!
! Use the fraction of solar shortwave radiation, as computed by optic_manizza from TOPAZ,
! (at vertical W-points) penetrating water column.
!
      DO k=0,N(ng)
        DO j=Jstr,Jend
          DO i=Istr,Iend
             swdk(i,j,k) = TP(ng)%decayW(i,j,N(ng)-k)
          END DO
        END DO
      ENDDO
#   else
.......

HCJung-jbnu / ROMS-TOPAZ

ROMS-TOPAZ column_physics 접합 #15

BIO source/sink 이슈

topaz.h : SUB. biology

generic_TOPAZ_mode.F90 : SUB. topaz_optic

pre_step3d.F : SUB. pre_step3d_tile