! $Id: iniacademic.F90 1625 2012-05-09 13:14:48Z lguez $ SUBROUTINE iniacademic_loc(vcov, ucov, teta, q, masse, ps, phis, time_0) USE lmdz_filtreg, ONLY: inifilr USE lmdz_infotrac, ONLY: nqtot, niso, iqIsoPha, tracers, getKey, isoName, addPhase USE control_mod, ONLY: day_step, planet_type USE exner_hyb_m, ONLY: exner_hyb USE exner_milieu_m, ONLY: exner_milieu USE parallel_lmdz, ONLY: ijb_u, ije_u, ijb_v, ije_v USE IOIPSL, ONLY: getin USE lmdz_write_field USE comconst_mod, ONLY: cpp, kappa, g, daysec, dtvr, pi, im, jm USE logic_mod, ONLY: iflag_phys, read_start USE comvert_mod, ONLY: ap, bp, preff, pa, presnivs, pressure_exner USE temps_mod, ONLY: annee_ref, day_ini, day_ref USE ener_mod, ONLY: etot0, ptot0, ztot0, stot0, ang0 USE netcdf, ONLY: nf90_nowrite, nf90_open, nf90_noerr, nf90_inq_varid, nf90_close, nf90_get_var USE lmdz_ran1, ONLY: ran1 USE lmdz_iniprint, ONLY: lunout, prt_level USE lmdz_academic, ONLY: tetarappel, knewt_t, kfrict, knewt_g, clat4 USE lmdz_comgeom ! Author: Frederic Hourdin original: 15/01/93 ! The forcing defined here is from Held and Suarez, 1994, Bulletin ! of the American Meteorological Society, 75, 1825. USE lmdz_dimensions, ONLY: iim, jjm, llm, ndm USE lmdz_paramet IMPLICIT NONE ! Declararations: ! --------------- ! Arguments: ! ---------- REAL, INTENT(OUT) :: time_0 ! fields REAL, INTENT(OUT) :: vcov(ijb_v:ije_v, llm) ! meridional covariant wind REAL, INTENT(OUT) :: ucov(ijb_u:ije_u, llm) ! zonal covariant wind REAL, INTENT(OUT) :: teta(ijb_u:ije_u, llm) ! potential temperature (K) REAL, INTENT(OUT) :: q(ijb_u:ije_u, llm, nqtot) ! advected tracers (.../kg_of_air) REAL, INTENT(OUT) :: ps(ijb_u:ije_u) ! surface pressure (Pa) REAL, INTENT(OUT) :: masse(ijb_u:ije_u, llm) ! air mass in grid cell (kg) REAL, INTENT(OUT) :: phis(ijb_u:ije_u) ! surface geopotential ! Local: ! ------ REAL, ALLOCATABLE :: vcov_glo(:, :), ucov_glo(:, :), teta_glo(:, :) REAL, ALLOCATABLE :: q_glo(:, :), masse_glo(:, :), ps_glo(:) REAL, ALLOCATABLE :: phis_glo(:) REAL p (ip1jmp1, llmp1) ! pression aux interfac.des couches REAL pks(ip1jmp1) ! exner au sol REAL pk(ip1jmp1, llm) ! exner au milieu des couches REAL phi(ip1jmp1, llm) ! geopotentiel REAL ddsin, zsig, tetapv, w_pv ! variables auxiliaires REAL tetastrat ! potential temperature in the stratosphere, in K REAL tetajl(jjp1, llm) INTEGER i, j, l, lsup, ij, iq, iName, iPhase, iqParent INTEGER :: nid_relief, varid, ierr REAL, DIMENSION(iip1, jjp1) :: relief REAL teta0, ttp, delt_y, delt_z, eps ! Constantes pour profil de T REAL k_f, k_c_a, k_c_s ! Constantes de rappel LOGICAL ok_geost ! Initialisation vent geost. ou nul LOGICAL ok_pv ! Polar Vortex REAL phi_pv, dphi_pv, gam_pv, tetanoise ! Constantes pour polar vortex REAL zz INTEGER idum REAL zdtvr, tnat, alpha_ideal LOGICAL :: ltnat1 CHARACTER(LEN=*),parameter :: modname="iniacademic" CHARACTER(LEN=80) :: abort_message ! Sanity check: verify that options selected by user are not incompatible IF ((iflag_phys==1).AND. .NOT. read_start) THEN WRITE(lunout, *) trim(modname), " error: if read_start is set to ", & " false then iflag_phys should not be 1" WRITE(lunout, *) "You most likely want an aquaplanet initialisation", & " (iflag_phys >= 100)" CALL abort_gcm(modname, "incompatible iflag_phys==1 and read_start==.FALSE.", 1) ENDIF !----------------------------------------------------------------------- ! 1. Initializations for Earth-like case ! -------------------------------------- ! initialize planet radius, rotation rate,... CALL conf_planete time_0 = 0. day_ref = 1 ! annee_ref=0 im = iim jm = jjm day_ini = 1 dtvr = daysec / REAL(day_step) zdtvr = dtvr etot0 = 0. ptot0 = 0. ztot0 = 0. stot0 = 0. ang0 = 0. IF (llm == 1) THEN ! specific initializations for the shallow water case kappa = 1 ENDIF CALL iniconst CALL inigeom CALL inifilr ! Initialize pressure and mass field if read_start=.FALSE. IF (.NOT. read_start) THEN ! allocate global fields: ! allocate(vcov_glo(ip1jm,llm)) allocate(ucov_glo(ip1jmp1, llm)) allocate(teta_glo(ip1jmp1, llm)) allocate(ps_glo(ip1jmp1)) allocate(masse_glo(ip1jmp1, llm)) allocate(phis_glo(ip1jmp1)) ! surface pressure ps_glo(:) = preff !------------------------------------------------------------------ ! Lecture eventuelle d'un fichier de relief interpollee sur la grille ! du modele. ! On suppose que le fichier relief_in.nc est stoké sur une grille ! iim*jjp1 ! Facile a créer à partir de la commande ! cdo remapcon,fichier_output_phys.nc Relief.nc relief_in.nc !------------------------------------------------------------------ relief = 0. ierr = nf90_open ('relief_in.nc', nf90_nowrite, nid_relief) IF (ierr==nf90_noerr) THEN ierr = nf90_inq_varid(nid_relief, 'RELIEF', varid) IF (ierr==nf90_noerr) THEN ierr = nf90_get_var(nid_relief, varid, relief(1:iim, 1:jjp1)) relief(iip1, :) = relief(1, :) else CALL abort_gcm ('iniacademic', 'variable RELIEF pas la', 1) endif endif ierr = nf90_close (nid_relief) !------------------------------------------------------------------ ! Initialisation du geopotentiel au sol et de la pression !------------------------------------------------------------------ PRINT*, 'relief=', minval(relief), maxval(relief), 'g=', g DO j = 1, jjp1 DO i = 1, iip1 phis_glo((j - 1) * iip1 + i) = g * relief(i, j) enddo enddo PRINT*, 'phis=', minval(phis), maxval(phis), 'g=', g CALL pression (ip1jmp1, ap, bp, ps_glo, p) IF (pressure_exner) THEN CALL exner_hyb(ip1jmp1, ps_glo, p, pks, pk) else CALL exner_milieu(ip1jmp1, ps_glo, p, pks, pk) endif CALL massdair(p, masse_glo) ENDIF IF (llm == 1) THEN ! initialize fields for the shallow water case, if required IF (.NOT.read_start) THEN phis(ijb_u:ije_u) = 0. q(ijb_u:ije_u, 1:llm, 1:nqtot) = 0 CALL sw_case_williamson91_6_loc(vcov, ucov, teta, masse, ps) endif ENDIF academic_case: if (iflag_phys >= 2) THEN ! initializations ! 1. local parameters ! by convention, winter is in the southern hemisphere ! Geostrophic wind or no wind? ok_geost = .TRUE. CALL getin('ok_geost', ok_geost) ! Constants for Newtonian relaxation and friction k_f = 1. !friction CALL getin('k_j', k_f) k_f = 1. / (daysec * k_f) k_c_s = 4. !cooling surface CALL getin('k_c_s', k_c_s) k_c_s = 1. / (daysec * k_c_s) k_c_a = 40. !cooling free atm CALL getin('k_c_a', k_c_a) k_c_a = 1. / (daysec * k_c_a) ! Constants for Teta equilibrium profile teta0 = 315. ! mean Teta (S.H. 315K) CALL getin('teta0', teta0) ttp = 200. ! Tropopause temperature (S.H. 200K) CALL getin('ttp', ttp) eps = 0. ! Deviation to N-S symmetry(~0-20K) CALL getin('eps', eps) delt_y = 60. ! Merid Temp. Gradient (S.H. 60K) CALL getin('delt_y', delt_y) delt_z = 10. ! Vertical Gradient (S.H. 10K) CALL getin('delt_z', delt_z) ! Polar vortex ok_pv = .FALSE. CALL getin('ok_pv', ok_pv) phi_pv = -50. ! Latitude of edge of vortex CALL getin('phi_pv', phi_pv) phi_pv = phi_pv * pi / 180. dphi_pv = 5. ! Width of the edge CALL getin('dphi_pv', dphi_pv) dphi_pv = dphi_pv * pi / 180. gam_pv = 4. ! -dT/dz vortex (in K/km) CALL getin('gam_pv', gam_pv) tetanoise = 0.005 CALL getin('tetanoise', tetanoise) ! 2. Initialize fields towards which to relax ! Friction knewt_g = k_c_a DO l = 1, llm zsig = presnivs(l) / preff knewt_t(l) = (k_c_s - k_c_a) * MAX(0., (zsig - 0.7) / 0.3) kfrict(l) = k_f * MAX(0., (zsig - 0.7) / 0.3) ENDDO DO j = 1, jjp1 clat4((j - 1) * iip1 + 1:j * iip1) = cos(rlatu(j))**4 ENDDO ! Potential temperature DO l = 1, llm zsig = presnivs(l) / preff tetastrat = ttp * zsig**(-kappa) tetapv = tetastrat IF ((ok_pv).AND.(zsig<0.1)) THEN tetapv = tetastrat * (zsig * 10.)**(kappa * cpp * gam_pv / 1000. / g) ENDIF DO j = 1, jjp1 ! Troposphere ddsin = sin(rlatu(j)) tetajl(j, l) = teta0 - delt_y * ddsin * ddsin + eps * ddsin & - delt_z * (1. - ddsin * ddsin) * log(zsig) IF (planet_type=="giant") THEN tetajl(j, l) = teta0 + (delt_y * & ((sin(rlatu(j) * 3.14159 * eps + 0.0001))**2) & / ((rlatu(j) * 3.14159 * eps + 0.0001)**2)) & - delt_z * log(zsig) endif ! Profil stratospherique isotherme (+vortex) w_pv = (1. - tanh((rlatu(j) - phi_pv) / dphi_pv)) / 2. tetastrat = tetastrat * (1. - w_pv) + tetapv * w_pv tetajl(j, l) = MAX(tetajl(j, l), tetastrat) ENDDO ENDDO ! CALL writefield('theta_eq',tetajl) DO l = 1, llm DO j = 1, jjp1 DO i = 1, iip1 ij = (j - 1) * iip1 + i tetarappel(ij, l) = tetajl(j, l) enddo enddo enddo ! 3. Initialize fields (if necessary) IF (.NOT. read_start) THEN ! bulk initialization of temperature IF (iflag_phys>10000) THEN ! Particular case to impose a constant temperature T0=0.01*iflag_phys teta_glo(:, :) = 0.01 * iflag_phys / (pk(:, :) / cpp) ELSE teta_glo(:, :) = tetarappel(:, :) ENDIF ! geopotential CALL geopot(ip1jmp1, teta_glo, pk, pks, phis_glo, phi) ! winds IF (ok_geost) THEN CALL ugeostr(phi, ucov_glo) else ucov_glo(:, :) = 0. endif vcov(ijb_v:ije_v, 1:llm) = 0. ! bulk initialization of tracers IF (planet_type=="earth") THEN ltnat1 = .TRUE.; CALL getin('tnateq1', ltnat1)! Earth: first two tracers will be water DO iq = 1, nqtot q(ijb_u:ije_u, :, iq) = 0. IF(tracers(iq)%name == addPhase('H2O', 'g')) q(ijb_u:ije_u, :, iq) = 1.e-10 IF(tracers(iq)%name == addPhase('H2O', 'l')) q(ijb_u:ije_u, :, iq) = 1.e-15 ! CRisi: init des isotopes ! distill de Rayleigh très simplifiée iName = tracers(iq)%iso_iName if (niso <= 0 .OR. iName <= 0) CYCLE iPhase = tracers(iq)%iso_iPhase iqParent = tracers(iq)%iqParent IF(tracers(iq)%iso_iZone == 0) THEN IF(ltnat1) THEN tnat = 1.0 alpha_ideal = 1.0 WRITE(lunout, *) 'In '//TRIM(modname)//': !!! Beware: alpha_ideal put to 1 !!!' ELSE IF(getKey('tnat', tnat, isoName(iName)) .OR. getKey('alpha', alpha_ideal, isoName(iName))) & CALL abort_gcm(TRIM(modname), 'missing isotopic parameters', 1) END IF q(ijb_u:ije_u,:,iq) = q(ijb_u:ije_u,:,iqParent)*tnat*(q(ijb_u:ije_u,:,iqParent)/30.e-3)**(alpha_ideal-1.) ELSE !IF(tracers(iq)%iso_iZone == 0) THEN IF(tracers(iq)%iso_iZone == 1) THEN ! a verifier. ! correction le 14 mai 2024 pour que tous les traceurs soient de la couleur 1. ! Sinon, on va avoir des porblèmes de conservation de masse de traceurs. q(ijb_u:ije_u,:,iq) = q(ijb_u:ije_u,:,iqIsoPha(iName,iPhase)) else !IF(tracers(iq)%iso_iZone == 1) THEN q(ijb_u:ije_u,:,iq) = 0.0 endif !IF(tracers(iq)%iso_iZone == 1) THEN END IF !IF(tracers(iq)%iso_iZone == 0) THEN enddo else q(ijb_u:ije_u,:,:)=0 endif ! of if (planet_type=="earth") CALL check_isotopes(q, ijb_u, ije_u, 'iniacademic_loc') ! add random perturbation to temperature idum = -1 zz = ran1(idum) idum = 0 DO l = 1, llm DO ij = iip2, ip1jm teta_glo(ij, l) = teta_glo(ij, l) * (1. + tetanoise * ran1(idum)) enddo enddo ! maintain periodicity in longitude DO l = 1, llm DO ij = 1, ip1jmp1, iip1 teta_glo(ij + iim, l) = teta_glo(ij, l) enddo enddo ! copy data from global array to local array: teta(ijb_u:ije_u, :) = teta_glo(ijb_u:ije_u, :) ucov(ijb_u:ije_u, :) = ucov_glo(ijb_u:ije_u, :) ! vcov(ijb_v:ije_v,:)=vcov_glo(ijb_v:ije_v,:) masse(ijb_u:ije_u, :) = masse_glo(ijb_u:ije_u, :) ps(ijb_u:ije_u) = ps_glo(ijb_u:ije_u) phis(ijb_u:ije_u) = phis_glo(ijb_u:ije_u) deallocate(teta_glo) deallocate(ucov_glo) ! deallocate(vcov_glo) deallocate(masse_glo) deallocate(ps_glo) deallocate(phis_glo) ENDIF ! of IF (.NOT. read_start) ENDIF academic_case END SUBROUTINE iniacademic_loc