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 10362 for NEMO/trunk/src/TOP/PISCES/SED/seddsr.F90 – NEMO

Ignore:
Timestamp:
2018-11-30T16:38:17+01:00 (5 years ago)
Author:
aumont
Message:

Various bug fixes and improvements

File:
1 edited

Legend:

Unmodified
Added
Removed

  • NEMO/trunk/src/TOP/PISCES/SED/seddsr.F90

    r10222 r10362  
    591591            zgamma = pwcp(ji,jk,jwalk) - 2.0 * zsedtra(4) 
    592592            zdelta = pwcp(ji,jk,jwpo4) - redfep * zsedtra(4) 
    593             zsedtra(4) = ( zsedtra(4) * zalpha ) / ( 0.25 * zsedtra(4) * ( exp( reac_fe2 * zalpha * dtsed2 / 2. ) - 1.0 )  & 
    594             &            + zalpha * exp( reac_fe2 * zalpha * dtsed2 / 2. ) + rtrn ) 
     593            IF ( zalpha == 0. ) THEN 
     594               zsedtra(4) = zsedtra(4) / ( 1.0 + zsedtra(4) * reac_fe2 * dtsed2 / 2.0 ) 
     595            ELSE 
     596               zsedtra(4) = ( zsedtra(4) * zalpha ) / ( 0.25 * zsedtra(4) * ( exp( reac_fe2 * zalpha * dtsed2 / 2. ) - 1.0 )  & 
     597               &            + zalpha * exp( reac_fe2 * zalpha * dtsed2 / 2. ) ) 
     598            ENDIF 
    595599            zsedtra(1) = zalpha + 0.25 * zsedtra(4)  
    596600            zsedtra(5) = zbeta  - zsedtra(4) 
    … …  
    601605            zbeta  = pwcp(ji,jk,jwso4) + zsedtra(2) 
    602606            zgamma = pwcp(ji,jk,jwalk) - 2.0 * zsedtra(2) 
    603             zsedtra(2) = ( zsedtra(2) * zalpha ) / ( 2.0 * zsedtra(2) * ( exp( reac_h2s * zalpha * dtsed2 / 2. ) - 1.0 )  & 
    604             &            + zalpha * exp( reac_h2s * zalpha * dtsed2 / 2. ) + rtrn ) 
     607            IF ( zalpha == 0. ) THEN 
     608               zsedtra(2) = zsedtra(2) / ( 1.0 + zsedtra(2) * reac_h2s * dtsed2 / 2.0 ) 
     609            ELSE 
     610               zsedtra(2) = ( zsedtra(2) * zalpha ) / ( 2.0 * zsedtra(2) * ( exp( reac_h2s * zalpha * dtsed2 / 2. ) - 1.0 )  & 
     611               &            + zalpha * exp( reac_h2s * zalpha * dtsed2 / 2. ) ) 
     612            ENDIF 
    605613            zsedtra(1) = zalpha + 2.0 * zsedtra(2) 
    606614            pwcp(ji,jk,jwalk) = zgamma + 2.0 * zsedtra(2) 
    … …  
    609617            zalpha = zsedtra(1) - 2.0 * zsedtra(3) 
    610618            zgamma = pwcp(ji,jk,jwalk) - 2.0 * zsedtra(3) 
    611             zsedtra(3) = ( zsedtra(3) * zalpha ) / ( 2.0 * zsedtra(3) * ( exp( reac_nh4 * zadsnh4 * zalpha * dtsed2 / 2. ) - 1.0 )  & 
    612             &            + zalpha * exp( reac_nh4 * zadsnh4 * zalpha * dtsed2 /2. ) + rtrn ) 
     619            IF ( zalpha == 0. ) THEN 
     620               zsedtra(3) = zsedtra(3) / ( 1.0 + zsedtra(3) * reac_nh4 * zadsnh4 * dtsed2 / 2.0 ) 
     621            ELSE 
     622               zsedtra(3) = ( zsedtra(3) * zalpha ) / ( 2.0 * zsedtra(3) * ( exp( reac_nh4 * zadsnh4 * zalpha * dtsed2 / 2. ) - 1.0 )  & 
     623               &            + zalpha * exp( reac_nh4 * zadsnh4 * zalpha * dtsed2 /2. ) ) 
     624            ENDIF 
    613625            zsedtra(1) = zalpha + 2.0 * zsedtra(3) 
    614626            pwcp(ji,jk,jwalk) = zgamma + 2.0 * zsedtra(3) 
    … …  
    617629            zbeta  = zsedtra(4) + zsedtra(6) 
    618630            zgamma = pwcp(ji,jk,jwso4) + zsedtra(6) 
    619             zsedtra(6) = ( zsedtra(6) * zalpha ) / ( 2.0 * zsedtra(6) * ( exp( reac_feso * zalpha * dtsed2 / 2. ) - 1.0 )  & 
    620             &            + zalpha * exp( reac_feso * zalpha * dtsed2 /2. ) + rtrn ) 
     631            IF ( zalpha == 0. ) THEN 
     632               zsedtra(6) = zsedtra(6) / ( 1.0 + zsedtra(6) * reac_feso * dtsed2 / 2.0 ) 
     633            ELSE 
     634               zsedtra(6) = ( zsedtra(6) * zalpha ) / ( 2.0 * zsedtra(6) * ( exp( reac_feso * zalpha * dtsed2 / 2. ) - 1.0 )  & 
     635               &            + zalpha * exp( reac_feso * zalpha * dtsed2 /2. ) ) 
     636            ENDIF 
    621637            zsedtra(1) = zalpha + 2.0 * zsedtra(6) 
    622638            zsedtra(4) = zbeta  - zsedtra(6) 
    … …  
    626642            zbeta  = zsedtra(4) + zsedtra(6) 
    627643            zgamma = pwcp(ji,jk,jwalk) - 2.0 * zsedtra(4) 
    628             zsedtra(4) = ( zsedtra(4) * zalpha ) / ( zsedtra(4) * ( exp( reac_fes * zalpha * dtsed2 / 2. ) - 1.0 )  & 
    629             &            + zalpha * exp( reac_fes * zalpha * dtsed2 /2. ) + rtrn ) 
     644            IF ( zalpha == 0. ) THEN 
     645               zsedtra(4) = zsedtra(4) / ( 1.0 + zsedtra(4) * reac_fes * dtsed2 / 2.0 ) 
     646            ELSE 
     647               zsedtra(4) = ( zsedtra(4) * zalpha ) / ( zsedtra(4) * ( exp( reac_fes * zalpha * dtsed2 / 2. ) - 1.0 )  & 
     648               &            + zalpha * exp( reac_fes * zalpha * dtsed2 /2. ) ) 
     649            ENDIF 
    630650            zsedtra(2) = zalpha + zsedtra(4) 
    631651            zsedtra(6) = zbeta  - zsedtra(4) 
    … …  
    637657            zdelta = pwcp(ji,jk,jwso4) + zsedtra(2) 
    638658            zepsi  = pwcp(ji,jk,jwpo4) + redfep * zsedtra(5) 
    639             zsedtra(2) = ( zsedtra(2) * zalpha ) / ( 2.0 * zsedtra(2) * ( exp( reac_feh2s * zalpha * dtsed2 ) - 1.0 )  & 
    640             &            + zalpha * exp( reac_feh2s * zalpha * dtsed2 ) + rtrn ) 
     659            IF ( zalpha == 0. ) THEN 
     660               zsedtra(2) = zsedtra(2) / ( 1.0 + zsedtra(2) * reac_feh2s * dtsed2 ) 
     661            ELSE 
     662               zsedtra(2) = ( zsedtra(2) * zalpha ) / ( 2.0 * zsedtra(2) * ( exp( reac_feh2s * zalpha * dtsed2 ) - 1.0 )  & 
     663               &            + zalpha * exp( reac_feh2s * zalpha * dtsed2 ) ) 
     664            ENDIF 
    641665            zsedtra(5) = zalpha + 2.0 * zsedtra(2) 
    642666            zsedtra(4) = zbeta  - zsedtra(5) 
    … …  
    648672            zbeta  = zsedtra(4) + zsedtra(6) 
    649673            zgamma = pwcp(ji,jk,jwalk) - 2.0 * zsedtra(4) 
    650             zsedtra(4) = ( zsedtra(4) * zalpha ) / ( zsedtra(4) * ( exp( reac_fes * zalpha * dtsed2 / 2. ) - 1.0 )  & 
    651             &            + zalpha * exp( reac_fes * zalpha * dtsed2 /2. ) + rtrn ) 
     674            IF ( zalpha == 0. ) THEN 
     675               zsedtra(4) = zsedtra(4) / ( 1.0 + zsedtra(4) * reac_fes * dtsed2 / 2.0 ) 
     676            ELSE 
     677               zsedtra(4) = ( zsedtra(4) * zalpha ) / ( zsedtra(4) * ( exp( reac_fes * zalpha * dtsed2 / 2. ) - 1.0 )  & 
     678               &            + zalpha * exp( reac_fes * zalpha * dtsed2 /2. ) ) 
     679            ENDIF 
    652680            zsedtra(2) = zalpha + zsedtra(4) 
    653681            zsedtra(6) = zbeta  - zsedtra(4) 
    … …  
    657685            zbeta  = zsedtra(4) + zsedtra(6) 
    658686            zgamma = pwcp(ji,jk,jwso4) + zsedtra(6) 
    659             zsedtra(6) = ( zsedtra(6) * zalpha ) / ( 2.0 * zsedtra(6) * ( exp( reac_feso * zalpha * dtsed2 / 2. ) - 1.0 )  & 
    660             &            + zalpha * exp( reac_feso * zalpha * dtsed2 /2. ) + rtrn ) 
     687            IF (zalpha == 0.) THEN 
     688               zsedtra(6) = zsedtra(6) / ( 1.0 + zsedtra(6) * reac_feso * dtsed2 / 2. ) 
     689            ELSE 
     690               zsedtra(6) = ( zsedtra(6) * zalpha ) / ( 2.0 * zsedtra(6) * ( exp( reac_feso * zalpha * dtsed2 / 2. ) - 1.0 )  & 
     691               &            + zalpha * exp( reac_feso * zalpha * dtsed2 /2. ) ) 
     692            ENDIF 
    661693            zsedtra(1) = zalpha + 2.0 * zsedtra(6) 
    662694            zsedtra(4) = zbeta  - zsedtra(6) 
    … …  
    665697            zalpha = zsedtra(1) - 2.0 * zsedtra(3) 
    666698            zgamma = pwcp(ji,jk,jwalk) - 2.0 * zsedtra(3) 
    667             zsedtra(3) = ( zsedtra(3) * zalpha ) / ( 2.0 * zsedtra(3) * ( exp( reac_nh4 * zadsnh4 * zalpha * dtsed2 / 2. ) - 1.0 )  & 
    668             &            + zalpha * exp( reac_nh4 * zadsnh4 * zalpha * dtsed2 /2. ) + rtrn ) 
     699            IF (zalpha == 0.) THEN 
     700               zsedtra(3) = zsedtra(3) / ( 1.0 + zsedtra(3) * reac_nh4 * zadsnh4 * dtsed2 / 2.0)  
     701            ELSE 
     702               zsedtra(3) = ( zsedtra(3) * zalpha ) / ( 2.0 * zsedtra(3) * ( exp( reac_nh4 * zadsnh4 * zalpha * dtsed2 / 2. ) - 1.0 )  & 
     703               &            + zalpha * exp( reac_nh4 * zadsnh4 * zalpha * dtsed2 /2. ) ) 
     704            ENDIF 
    669705            zsedtra(1) = zalpha + 2.0 * zsedtra(3) 
    670706            pwcp(ji,jk,jwalk) = zgamma + 2.0 * zsedtra(3) 
    … …  
    673709            zbeta  = pwcp(ji,jk,jwso4) + zsedtra(2) 
    674710            zgamma = pwcp(ji,jk,jwalk) - 2.0 * zsedtra(2) 
    675             zsedtra(2) = ( zsedtra(2) * zalpha ) / ( 2.0 * zsedtra(2) * ( exp( reac_h2s * zalpha * dtsed2 / 2. ) - 1.0 )  & 
    676             &            + zalpha * exp( reac_h2s * zalpha * dtsed2 / 2. ) + rtrn ) 
     711            IF ( zalpha == 0. ) THEN 
     712               zsedtra(2) = zsedtra(2) / ( 1.0 + zsedtra(2) * reac_h2s * dtsed2 / 2.0 ) 
     713            ELSE 
     714               zsedtra(2) = ( zsedtra(2) * zalpha ) / ( 2.0 * zsedtra(2) * ( exp( reac_h2s * zalpha * dtsed2 / 2. ) - 1.0 )  & 
     715               &            + zalpha * exp( reac_h2s * zalpha * dtsed2 / 2. ) ) 
     716            ENDIF 
    677717            zsedtra(1) = zalpha + 2.0 * zsedtra(2) 
    678718            pwcp(ji,jk,jwso4) = zbeta - zsedtra(2) 
    … …  
    683723            zgamma = pwcp(ji,jk,jwalk) - 2.0 * zsedtra(4) 
    684724            zdelta = pwcp(ji,jk,jwpo4) - redfep * zsedtra(4) 
    685             zsedtra(4) = ( zsedtra(4) * zalpha ) / ( 0.25 * zsedtra(4) * ( exp( reac_fe2 * zalpha * dtsed2 / 2. ) - 1.0 )  & 
    686             &            + zalpha * exp( reac_fe2 * zalpha * dtsed2 / 2. ) + rtrn ) 
     725            IF ( zalpha == 0. ) THEN 
     726               zsedtra(4) = zsedtra(4) / ( 1.0 + zsedtra(4) * reac_fe2 * dtsed2 / 2.0 ) 
     727            ELSE 
     728               zsedtra(4) = ( zsedtra(4) * zalpha ) / ( 0.25 * zsedtra(4) * ( exp( reac_fe2 * zalpha * dtsed2 / 2. ) - 1.0 )  & 
     729               &            + zalpha * exp( reac_fe2 * zalpha * dtsed2 / 2. ) ) 
     730            ENDIF 
    687731            zsedtra(1) = zalpha + 0.25 * zsedtra(4) 
    688732            zsedtra(5) = zbeta  - zsedtra(4) 
Note: See TracChangeset for help on using the changeset viewer.