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 10812 for branches/UKMO/r8395_India_uncoupled/NEMOGCM/NEMO/OPA_SRC/DOM/dtatsd.F90 – NEMO

Ignore:
Timestamp:
2019-03-29T15:46:28+01:00 (5 years ago)
Author:
jcastill
Message:

Add again the changes that were reverted before: they are needed after all

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/UKMO/r8395_India_uncoupled/NEMOGCM/NEMO/OPA_SRC/DOM/dtatsd.F90

    r10724 r10812  
    2323   USE wrk_nemo        ! Memory allocation 
    2424   USE timing          ! Timing 
     25   USE iom 
    2526 
    2627   IMPLICIT NONE 
     
    3132 
    3233   LOGICAL , PUBLIC ::   ln_tsd_init      !: T & S data flag 
     34   LOGICAL , PUBLIC ::   ln_tsd_interp    !: vertical interpolation flag 
    3335   LOGICAL , PUBLIC ::   ln_tsd_tradmp    !: internal damping toward input data flag 
    3436 
    3537   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_tsd   ! structure of input SST (file informations, fields read) 
     38   INTEGER                                 ::   jpk_init , inum_dta 
     39   INTEGER                                 ::   id ,linum   ! local integers 
     40   INTEGER                                 ::   zdim(4) 
    3641 
    3742   !!---------------------------------------------------------------------- 
     
    5358      LOGICAL, INTENT(in), OPTIONAL ::   ld_tradmp   ! force the initialization when tradp is used 
    5459      ! 
    55       INTEGER ::   ios, ierr0, ierr1, ierr2, ierr3   ! local integers 
     60      INTEGER ::   ios, ierr0, ierr1, ierr2, ierr3, ierr4, ierr5   ! local integers 
    5661      !! 
    5762      CHARACTER(len=100)            ::   cn_dir          ! Root directory for location of ssr files 
    58       TYPE(FLD_N), DIMENSION( jpts) ::   slf_i           ! array of namelist informations on the fields to read 
    59       TYPE(FLD_N)                   ::   sn_tem, sn_sal 
     63      TYPE(FLD_N), DIMENSION(jpts+2)::   slf_i           ! array of namelist informations on the fields to read 
     64      TYPE(FLD_N)                   ::   sn_tem, sn_sal, sn_dep, sn_msk 
     65       
    6066      !! 
    61       NAMELIST/namtsd/   ln_tsd_init, ln_tsd_tradmp, cn_dir, sn_tem, sn_sal 
     67      NAMELIST/namtsd/   ln_tsd_init, ln_tsd_interp, ln_tsd_tradmp, cn_dir, sn_tem, sn_sal, sn_dep, sn_msk 
    6268      !!---------------------------------------------------------------------- 
    6369      ! 
     
    6571      ! 
    6672      !  Initialisation 
    67       ierr0 = 0  ;  ierr1 = 0  ;  ierr2 = 0  ;  ierr3 = 0 
     73      ierr0 = 0  ;  ierr1 = 0  ;  ierr2 = 0  ;  ierr3 = 0  ; ierr4 = 0  ;  ierr5 = 0  
    6874      ! 
    6975      REWIND( numnam_ref )              ! Namelist namtsd in reference namelist :  
     
    8490         WRITE(numout,*) '   Namelist namtsd' 
    8591         WRITE(numout,*) '      Initialisation of ocean T & S with T &S input data   ln_tsd_init   = ', ln_tsd_init 
     92         WRITE(numout,*) '      iInterpolation of initial conditions in the vertical ln_tsd_interp = ', ln_tsd_interp 
    8693         WRITE(numout,*) '      damping of ocean T & S toward T &S input data        ln_tsd_tradmp = ', ln_tsd_tradmp 
    8794         WRITE(numout,*) 
     
    97104         ln_tsd_init = .FALSE. 
    98105      ENDIF 
     106      IF( ln_tsd_interp .AND. ln_tsd_tradmp ) THEN 
     107            CALL ctl_stop( 'dta_tsd_init: Tracer damping and vertical interpolation not yet configured' )   ;   RETURN 
     108      ENDIF 
     109      IF( ln_tsd_interp .AND. LEN(TRIM(sn_msk%wname)) > 0 ) THEN 
     110            CALL ctl_stop( 'dta_tsd_init: Using vertical interpolation and weights files not recommended' )   ;   RETURN 
     111      ENDIF 
    99112      ! 
    100113      !                             ! allocate the arrays (if necessary) 
    101114      IF(  ln_tsd_init .OR. ln_tsd_tradmp  ) THEN 
    102115         ! 
    103          ALLOCATE( sf_tsd(jpts), STAT=ierr0 ) 
     116         IF( ln_tsd_interp ) THEN 
     117           ALLOCATE( sf_tsd(jpts+2), STAT=ierr0 ) ! to carry the addtional depth information 
     118         ELSE 
     119           ALLOCATE( sf_tsd(jpts  ), STAT=ierr0 )  
     120         ENDIF  
    104121         IF( ierr0 > 0 ) THEN 
    105122            CALL ctl_stop( 'dta_tsd_init: unable to allocate sf_tsd structure' )   ;   RETURN 
    106123         ENDIF 
    107124         ! 
    108                                 ALLOCATE( sf_tsd(jp_tem)%fnow(jpi,jpj,jpk)   , STAT=ierr0 ) 
    109          IF( sn_tem%ln_tint )   ALLOCATE( sf_tsd(jp_tem)%fdta(jpi,jpj,jpk,2) , STAT=ierr1 ) 
    110                                 ALLOCATE( sf_tsd(jp_sal)%fnow(jpi,jpj,jpk)   , STAT=ierr2 ) 
    111          IF( sn_sal%ln_tint )   ALLOCATE( sf_tsd(jp_sal)%fdta(jpi,jpj,jpk,2) , STAT=ierr3 ) 
    112          ! 
    113          IF( ierr0 + ierr1 + ierr2 + ierr3 > 0 ) THEN 
     125         IF( ln_tsd_interp ) THEN 
     126            CALL iom_open ( trim(cn_dir) // trim(sn_dep%clname), inum_dta )  
     127            id = iom_varid( inum_dta, sn_dep%clvar, zdim ) 
     128            jpk_init = zdim(3) 
     129            IF(lwp) WRITE(numout,*) 'Dimension of veritcal coordinate in ICs: ', jpk_init 
     130            CALL iom_close( inum_dta )   ! Close the input file 
     131            ! 
     132                                 ALLOCATE( sf_tsd(jp_tem)%fnow(jpi,jpj,jpk_init  ) , STAT=ierr0 ) 
     133            IF( sn_tem%ln_tint ) ALLOCATE( sf_tsd(jp_tem)%fdta(jpi,jpj,jpk_init,2) , STAT=ierr1 ) 
     134                                 ALLOCATE( sf_tsd(jp_sal)%fnow(jpi,jpj,jpk_init  ) , STAT=ierr2 ) 
     135            IF( sn_sal%ln_tint ) ALLOCATE( sf_tsd(jp_sal)%fdta(jpi,jpj,jpk_init,2) , STAT=ierr3 )   
     136                                 ALLOCATE( sf_tsd(jp_dep)%fnow(jpi,jpj,jpk_init  ) , STAT=ierr4 ) 
     137                                 ALLOCATE( sf_tsd(jp_msk)%fnow(jpi,jpj,jpk_init  ) , STAT=ierr5 ) 
     138         ELSE 
     139                                 ALLOCATE( sf_tsd(jp_tem)%fnow(jpi,jpj,jpk)   , STAT=ierr0 ) 
     140            IF( sn_tem%ln_tint ) ALLOCATE( sf_tsd(jp_tem)%fdta(jpi,jpj,jpk,2) , STAT=ierr1 ) 
     141                                 ALLOCATE( sf_tsd(jp_sal)%fnow(jpi,jpj,jpk)   , STAT=ierr2 ) 
     142            IF( sn_sal%ln_tint ) ALLOCATE( sf_tsd(jp_sal)%fdta(jpi,jpj,jpk,2) , STAT=ierr3 )   
     143         ENDIF ! ln_tsd_interp 
     144 
     145         ! 
     146         IF( ierr0 + ierr1 + ierr2 + ierr3 + ierr4 + ierr5 > 0 ) THEN 
    114147            CALL ctl_stop( 'dta_tsd : unable to allocate T & S data arrays' )   ;   RETURN 
    115148         ENDIF 
    116149         !                         ! fill sf_tsd with sn_tem & sn_sal and control print 
    117150         slf_i(jp_tem) = sn_tem   ;   slf_i(jp_sal) = sn_sal 
     151         IF( ln_tsd_interp ) slf_i(jp_dep) = sn_dep   ;   slf_i(jp_msk) = sn_msk 
    118152         CALL fld_fill( sf_tsd, slf_i, cn_dir, 'dta_tsd', 'Temperature & Salinity data', 'namtsd', no_print ) 
    119153         ! 
     
    143177      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(  out) ::   ptsd   ! T & S data 
    144178      ! 
    145       INTEGER ::   ji, jj, jk, jl, jkk   ! dummy loop indicies 
    146       INTEGER ::   ik, il0, il1, ii0, ii1, ij0, ij1   ! local integers 
     179      INTEGER ::   ji, jj, jk, jl, jk_init   ! dummy loop indicies 
     180      INTEGER ::   ik, il0, il1, ii0, ii1, ij0, ij1        ! local integers 
    147181      REAL(wp)::   zl, zi 
    148       REAL(wp), POINTER, DIMENSION(:) ::  ztp, zsp   ! 1D workspace 
    149182      !!---------------------------------------------------------------------- 
    150183      ! 
     
    181214!!gm end 
    182215      ! 
    183       ptsd(:,:,:,jp_tem) = sf_tsd(jp_tem)%fnow(:,:,:)    ! NO mask 
    184       ptsd(:,:,:,jp_sal) = sf_tsd(jp_sal)%fnow(:,:,:)  
    185       ! 
    186       IF( ln_sco ) THEN                   !==   s- or mixed s-zps-coordinate   ==! 
    187          ! 
    188          CALL wrk_alloc( jpk, ztp, zsp ) 
    189          ! 
    190          IF( kt == nit000 .AND. lwp )THEN 
    191             WRITE(numout,*) 
    192             WRITE(numout,*) 'dta_tsd: interpolates T & S data onto the s- or mixed s-z-coordinate mesh' 
    193          ENDIF 
    194          ! 
    195          DO jj = 1, jpj                         ! vertical interpolation of T & S 
    196             DO ji = 1, jpi 
    197                DO jk = 1, jpk                        ! determines the intepolated T-S profiles at each (i,j) points 
     216      IF( kt == nit000 .AND. lwp )THEN 
     217         WRITE(numout,*) 
     218         WRITE(numout,*) 'dta_tsd: interpolates T & S data onto current mesh' 
     219      ENDIF 
     220      ! 
     221      IF( ln_tsd_interp ) THEN                 ! probably should use pointers in the following to make more readable 
     222      ! 
     223         DO jk = 1, jpk                        ! determines the intepolated T-S profiles at each (i,j) points 
     224            DO jj= 1, jpj 
     225               DO ji= 1, jpi 
    198226                  zl = gdept_0(ji,jj,jk) 
    199                   IF(     zl < gdept_1d(1  ) ) THEN          ! above the first level of data 
    200                      ztp(jk) =  ptsd(ji,jj,1    ,jp_tem) 
    201                      zsp(jk) =  ptsd(ji,jj,1    ,jp_sal) 
    202                   ELSEIF( zl > gdept_1d(jpk) ) THEN          ! below the last level of data 
    203                      ztp(jk) =  ptsd(ji,jj,jpkm1,jp_tem) 
    204                      zsp(jk) =  ptsd(ji,jj,jpkm1,jp_sal) 
    205                   ELSE                                      ! inbetween : vertical interpolation between jkk & jkk+1 
    206                      DO jkk = 1, jpkm1                                  ! when  gdept(jkk) < zl < gdept(jkk+1) 
    207                         IF( (zl-gdept_1d(jkk)) * (zl-gdept_1d(jkk+1)) <= 0._wp ) THEN 
    208                            zi = ( zl - gdept_1d(jkk) ) / (gdept_1d(jkk+1)-gdept_1d(jkk)) 
    209                            ztp(jk) = ptsd(ji,jj,jkk,jp_tem) + ( ptsd(ji,jj,jkk+1,jp_tem) - ptsd(ji,jj,jkk,jp_tem) ) * zi  
    210                            zsp(jk) = ptsd(ji,jj,jkk,jp_sal) + ( ptsd(ji,jj,jkk+1,jp_sal) - ptsd(ji,jj,jkk,jp_sal) ) * zi 
     227                  IF( zl < sf_tsd(jp_dep)%fnow(ji,jj,1) ) THEN                     ! above the first level of data 
     228                     ptsd(ji,jj,jk,jp_tem) = sf_tsd(jp_tem)%fnow(ji,jj,1)  
     229                     ptsd(ji,jj,jk,jp_sal) = sf_tsd(jp_sal)%fnow(ji,jj,1) 
     230                  ELSEIF( zl > sf_tsd(jp_dep)%fnow(ji,jj,jpk_init) ) THEN          ! below the last level of data 
     231                     ptsd(ji,jj,jk,jp_tem) = sf_tsd(jp_tem)%fnow(ji,jj,jpk_init) 
     232                     ptsd(ji,jj,jk,jp_sal) = sf_tsd(jp_sal)%fnow(ji,jj,jpk_init) 
     233                  ELSE                                                             ! inbetween : vertical interpolation between jk_init & jk_init+1 
     234                     DO jk_init = 1, jpk_init-1                                    ! when  gdept(jk_init) < zl < gdept(jk_init+1) 
     235                        IF( sf_tsd(jp_msk)%fnow(ji,jj,jk_init+1) == 0 ) THEN       ! if there is no data fill down 
     236                           sf_tsd(jp_tem)%fnow(ji,jj,jk_init+1) = sf_tsd(jp_tem)%fnow(ji,jj,jk_init) 
     237                           sf_tsd(jp_sal)%fnow(ji,jj,jk_init+1) = sf_tsd(jp_sal)%fnow(ji,jj,jk_init) 
     238                        ENDIF 
     239                        IF( (zl-sf_tsd(jp_dep)%fnow(ji,jj,jk_init)) * (zl-sf_tsd(jp_dep)%fnow(ji,jj,jk_init+1)) <= 0._wp ) THEN 
     240                           zi = ( zl - sf_tsd(jp_dep)%fnow(ji,jj,jk_init) ) / & 
     241                        &       (sf_tsd(jp_dep)%fnow(ji,jj,jk_init+1)-sf_tsd(jp_dep)%fnow(ji,jj,jk_init)) 
     242                           ptsd(ji,jj,jk,jp_tem) = sf_tsd(jp_tem)%fnow(ji,jj,jk_init) + & 
     243                        &                          (sf_tsd(jp_tem)%fnow(ji,jj,jk_init+1)-sf_tsd(jp_tem)%fnow(ji,jj,jk_init)) * zi 
     244                           ptsd(ji,jj,jk,jp_sal) = sf_tsd(jp_sal)%fnow(ji,jj,jk_init) + & 
     245                        &                          (sf_tsd(jp_sal)%fnow(ji,jj,jk_init+1)-sf_tsd(jp_sal)%fnow(ji,jj,jk_init)) * zi 
    211246                        ENDIF 
    212247                     END DO 
    213248                  ENDIF 
    214                END DO 
    215                DO jk = 1, jpkm1 
    216                   ptsd(ji,jj,jk,jp_tem) = ztp(jk) * tmask(ji,jj,jk)     ! mask required for mixed zps-s-coord 
    217                   ptsd(ji,jj,jk,jp_sal) = zsp(jk) * tmask(ji,jj,jk) 
    218                END DO 
    219                ptsd(ji,jj,jpk,jp_tem) = 0._wp 
    220                ptsd(ji,jj,jpk,jp_sal) = 0._wp 
    221             END DO 
     249               ENDDO 
     250            ENDDO 
    222251         END DO 
    223          !  
    224          CALL wrk_dealloc( jpk, ztp, zsp ) 
    225          !  
     252         ! 
     253         ptsd(:,:,:,jp_tem) = ptsd(:,:,:,jp_tem) *tmask(:,:,:) 
     254         ptsd(:,:,:,jp_sal) = ptsd(:,:,:,jp_sal) *tmask(:,:,:) 
    226255      ELSE                                !==   z- or zps- coordinate   ==! 
    227256         !                              
    228          ptsd(:,:,:,jp_tem) = ptsd(:,:,:,jp_tem) * tmask(:,:,:)    ! Mask 
    229          ptsd(:,:,:,jp_sal) = ptsd(:,:,:,jp_sal) * tmask(:,:,:) 
     257         ptsd(:,:,:,jp_tem) = sf_tsd(jp_tem)%fnow(:,:,:)  * tmask(:,:,:)  ! Mask 
     258         ptsd(:,:,:,jp_sal) = sf_tsd(jp_sal)%fnow(:,:,:) * tmask(:,:,:) 
    230259         ! 
    231260         IF( ln_zps ) THEN                      ! zps-coordinate (partial steps) interpolation at the last ocean level 
     
    257286                                        DEALLOCATE( sf_tsd(jp_sal)%fnow )     ! S arrays in the structure 
    258287         IF( sf_tsd(jp_sal)%ln_tint )   DEALLOCATE( sf_tsd(jp_sal)%fdta ) 
     288         IF( ln_tsd_interp )            DEALLOCATE( sf_tsd(jp_dep)%fnow )     ! T arrays in the structure 
     289         IF( ln_tsd_interp )            DEALLOCATE( sf_tsd(jp_msk)%fnow )     ! T arrays in the structure 
    259290                                        DEALLOCATE( sf_tsd              )     ! the structure itself 
    260291      ENDIF 
Note: See TracChangeset for help on using the changeset viewer.