source: LMDZ6/branches/Amaury_dev/libf/phylmd/phyaqua_mod.F90 @ 5119

Last change on this file since 5119 was 5117, checked in by abarral, 4 months ago

rename modules properly lmdz_*
move some unused files to obsolete/
(lint) uppercase fortran keywords
  • Property copyright set to
    Name of program: LMDZ
    Creation date: 1984
    Version: LMDZ5
    License: CeCILL version 2
    Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
    See the license file in the root directory
  • Property svn:keywords set to Id
File size: 30.5 KB
Line 
1! $Id: phyaqua_mod.F90 5117 2024-07-24 14:23:34Z abarral $
2
3MODULE phyaqua_mod
4  ! Routines complementaires pour la physique planetaire.
5  USE lmdz_abort_physic, ONLY: abort_physic
6  IMPLICIT NONE
7
8CONTAINS
9
10  SUBROUTINE iniaqua(nlon, year_len, iflag_phys)
11
12    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
13    ! Creation d'un etat initial et de conditions aux limites
14    ! (resp startphy.nc et limit.nc) pour des configurations idealisees
15    ! du modele LMDZ dans sa version terrestre.
16    ! iflag_phys est un parametre qui controle
17    ! iflag_phys = N
18    ! de 100 a 199 : aqua planetes avec SST forcees
19    ! N-100 determine le type de SSTs
20    ! de 200 a 299 : terra planetes avec Ts calcule
21
22    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
23
24    USE dimphy, ONLY: klon
25    USE lmdz_geometry, ONLY: latitude
26    USE surface_data, ONLY: type_ocean, ok_veget
27    USE pbl_surface_mod, ONLY: pbl_surface_init
28    USE fonte_neige_mod, ONLY: fonte_neige_init
29    USE phys_state_var_mod
30    USE time_phylmdz_mod, ONLY: day_ref, ndays, pdtphys, &
31            day_ini, day_end
32    USE indice_sol_mod
33    USE lmdz_physical_constants, ONLY: pi
34    !    USE ioipsl
35    USE lmdz_phys_para, ONLY: is_master
36    USE lmdz_phys_transfert_para, ONLY: bcast
37    USE lmdz_grid_phy
38    USE lmdz_ioipsl_getin_p, ONLY: getin_p
39    USE phys_cal_mod, ONLY: calend, year_len_phy => year_len
40    IMPLICIT NONE
41
42    include "YOMCST.h"
43    include "clesphys.h"
44    include "dimsoil.h"
45
46    INTEGER, INTENT (IN) :: nlon, year_len, iflag_phys
47    ! IM ajout latfi, lonfi
48    !    REAL, INTENT (IN) :: lonfi(nlon), latfi(nlon)
49
50    INTEGER type_profil, type_aqua
51
52    ! Ajouts initialisation des surfaces
53    REAL :: run_off_lic_0(nlon)
54    REAL :: qsolsrf(nlon, nbsrf), snsrf(nlon, nbsrf)
55    REAL :: tsoil(nlon, nsoilmx, nbsrf)
56    REAL :: tslab(nlon), seaice(nlon)
57    REAL fder(nlon)
58
59
60
61    ! Arguments :
62    ! -----------
63
64    ! integer radpas
65    INTEGER it, unit, i, k, itap
66
67    REAL rugos, albedo
68    REAL tsurf
69    REAL time, timestep, day, day0
70    REAL qsol_f
71    REAL rugsrel(nlon)
72    LOGICAL alb_ocean
73
74    CHARACTER *80 ans, file_forctl, file_fordat, file_start
75    CHARACTER *100 file, var
76    CHARACTER *2 cnbl
77
78    REAL phy_nat(nlon, year_len)
79    REAL phy_alb(nlon, year_len)
80    REAL phy_sst(nlon, year_len)
81    REAL phy_bil(nlon, year_len)
82    REAL phy_rug(nlon, year_len)
83    REAL phy_ice(nlon, year_len)
84    REAL phy_fter(nlon, year_len)
85    REAL phy_foce(nlon, year_len)
86    REAL phy_fsic(nlon, year_len)
87    REAL phy_flic(nlon, year_len)
88
89    INTEGER, SAVE :: read_climoz = 0 ! read ozone climatology
90    !$OMP THREADPRIVATE(read_climoz)
91
92    ! -------------------------------------------------------------------------
93    ! declaration pour l'appel a phyredem
94    ! -------------------------------------------------------------------------
95
96    ! real pctsrf(nlon,nbsrf),ftsol(nlon,nbsrf)
97    REAL falbe(nlon, nbsrf), falblw(nlon, nbsrf)
98    ! real pbl_tke(nlon,llm,nbsrf)
99    ! real rain_fall(nlon),snow_fall(nlon)
100    ! real solsw(nlon), sollw(nlon),radsol(nlon)
101    ! real t_ancien(nlon,llm),q_ancien(nlon,llm),rnebcon(nlon,llm)
102    ! real ratqs(nlon,llm)
103    ! real clwcon(nlon,llm)
104
105    INTEGER longcles
106    PARAMETER (longcles = 20)
107    REAL clesphy0(longcles)
108
109
110    ! -----------------------------------------------------------------------
111    ! dynamial tendencies :
112    ! ---------------------
113
114    INTEGER l, ierr, aslun
115
116    REAL paire
117
118    ! Local
119    CHARACTER (LEN = 20) :: modname = 'phyaqua'
120    CHARACTER (LEN = 80) :: abort_message
121
122
123    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
124    ! INITIALISATIONS
125    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
126
127    ! -----------------------------------------------------------------------
128    ! Initialisations  des constantes
129    ! -------------------------------
130
131    !IF (calend .EQ. "earth_360d") Then
132    year_len_phy = year_len
133    !END IF
134
135    IF (year_len/=360) THEN
136      write (*, *) year_len
137      CALL abort_physic("iniaqua", 'iniaqua: 360 day calendar is required !', 1)
138    endif
139
140    type_aqua = iflag_phys / 100
141    type_profil = iflag_phys - type_aqua * 100
142    PRINT *, 'iniaqua:type_aqua, type_profil', type_aqua, type_profil
143
144    IF (klon/=nlon) THEN
145      WRITE (*, *) 'iniaqua: klon=', klon, ' nlon=', nlon
146      abort_message = 'probleme de dimensions dans iniaqua'
147      CALL abort_physic(modname, abort_message, 1)
148    END IF
149    CALL phys_state_var_init(read_climoz)
150
151    read_climoz = 0
152    day0 = 217.
153    day = day0
154    it = 0
155    time = 0.
156
157    ! -----------------------------------------------------------------------
158    ! initialisations de la physique
159    ! -----------------------------------------------------------------------
160
161    day_ini = day_ref
162    day_end = day_ini + ndays
163
164    nbapp_rad = 24
165    CALL getin_p('nbapp_rad', nbapp_rad)
166
167    ! ---------------------------------------------------------------------
168    ! Creation des conditions aux limites:
169    ! ------------------------------------
170    ! Initialisations des constantes
171    ! Ajouter les manquants dans planete.def... (albedo etc)
172    co2_ppm = 348.
173    CALL getin_p('co2_ppm', co2_ppm)
174
175    solaire = 1365.
176    CALL getin_p('solaire', solaire)
177
178    ! CALL getin('albedo',albedo) ! albedo is set below, depending on
179    ! type_aqua
180    alb_ocean = .TRUE.
181    CALL getin_p('alb_ocean', alb_ocean)
182
183    WRITE (*, *) 'iniaqua: co2_ppm=', co2_ppm
184    WRITE (*, *) 'iniaqua: solaire=', solaire
185    WRITE (*, *) 'iniaqua: alb_ocean=', alb_ocean
186
187    radsol = 0.
188    qsol_f = 10.
189
190    ! Conditions aux limites:
191    ! -----------------------
192
193    qsol(:) = qsol_f
194    rugsrel = 0.0 ! (rugsrel = rugoro)
195    rugoro = 0.0
196    u_ancien = 0.0
197    v_ancien = 0.0
198    agesno = 50.0
199    ! Relief plat
200    zmea = 0.
201    zstd = 0.
202    zsig = 0.
203    zgam = 0.
204    zthe = 0.
205    zpic = 0.
206    zval = 0.
207
208    ! Une seule surface
209    pctsrf = 0.
210    IF (type_aqua==1) THEN
211      rugos = 1.E-4
212      albedo = 0.19
213      pctsrf(:, is_oce) = 1.
214    ELSE IF (type_aqua==2) THEN
215      rugos = 0.03
216      albedo = 0.1
217      pctsrf(:, is_ter) = 1.
218    END IF
219
220    CALL getin_p('rugos', rugos)
221
222    WRITE (*, *) 'iniaqua: rugos=', rugos
223    zmasq(:) = pctsrf(:, is_ter)
224
225    ! pctsrf_pot(:,is_oce) = 1. - zmasq(:)
226    ! pctsrf_pot(:,is_sic) = 1. - zmasq(:)
227
228    ! Si alb_ocean on calcule un albedo oceanique moyen
229    ! if (alb_ocean) THEN
230    ! Voir pourquoi on avait ca.
231    ! CALL ini_alb_oce(phy_alb)
232    ! else
233    phy_alb(:, :) = albedo ! albedo land only (old value condsurf_jyg=0.3)
234    ! END IF !alb_ocean
235
236    DO i = 1, year_len
237      ! IM Terraplanete   phy_sst(:,i) = 260.+50.*cos(rlatd(:))**2
238      ! IM ajout calcul profil sst selon le cas considere (cf. FBr)
239
240      phy_nat(:, i) = 1.0 ! 0=ocean libre, 1=land, 2=glacier, 3=banquise
241      phy_bil(:, i) = 1.0 ! ne sert que pour les slab_ocean
242      phy_rug(:, i) = rugos ! longueur rugosite utilisee sur land only
243      phy_ice(:, i) = 0.0 ! fraction de glace (?)
244      phy_fter(:, i) = pctsrf(:, is_ter) ! fraction de glace (?)
245      phy_foce(:, i) = pctsrf(:, is_oce) ! fraction de glace (?)
246      phy_fsic(:, i) = pctsrf(:, is_sic) ! fraction de glace (?)
247      phy_flic(:, i) = pctsrf(:, is_lic) ! fraction de glace (?)
248    END DO
249    ! IM calcul profil sst
250    CALL profil_sst(nlon, latitude, type_profil, phy_sst)
251
252    IF (grid_type==unstructured) THEN
253      CALL writelim_unstruct(klon, phy_nat, phy_alb, phy_sst, phy_bil, phy_rug, phy_ice, &
254              phy_fter, phy_foce, phy_flic, phy_fsic)
255    ELSE
256
257      CALL writelim(klon, phy_nat, phy_alb, phy_sst, phy_bil, phy_rug, phy_ice, &
258              phy_fter, phy_foce, phy_flic, phy_fsic)
259    ENDIF
260
261    ! ---------------------------------------------------------------------
262    ! Ecriture de l'etat initial:
263    ! ---------------------------
264
265
266    ! Ecriture etat initial physique
267
268    timestep = pdtphys
269    radpas = nint(rday / timestep / float(nbapp_rad))
270
271    DO i = 1, longcles
272      clesphy0(i) = 0.
273    END DO
274    clesphy0(1) = float(iflag_con)
275    clesphy0(2) = float(nbapp_rad)
276    ! IF( cycle_diurne  ) clesphy0(3) =  1.
277    clesphy0(3) = 1. ! cycle_diurne
278    clesphy0(4) = 1. ! soil_model
279    clesphy0(5) = 1. ! new_oliq
280    clesphy0(6) = 0. ! ok_orodr
281    clesphy0(7) = 0. ! ok_orolf
282    clesphy0(8) = 0. ! ok_limitvrai
283
284
285    ! =======================================================================
286    ! Profils initiaux
287    ! =======================================================================
288
289    ! On initialise les temperatures de surfaces comme les sst
290    DO i = 1, nlon
291      ftsol(i, :) = phy_sst(i, 1)
292      tsoil(i, :, :) = phy_sst(i, 1)
293      tslab(i) = phy_sst(i, 1)
294    END DO
295
296    falbe(:, :) = albedo
297    falblw(:, :) = albedo
298    rain_fall(:) = 0.
299    snow_fall(:) = 0.
300    solsw(:) = 0.
301    sollw(:) = 0.
302    radsol(:) = 0.
303
304    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
305    ! intialisation bidon mais pas grave
306    t_ancien(:, :) = 0.
307    q_ancien(:, :) = 0.
308    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
309    rnebcon = 0.
310    ratqs = 0.
311    clwcon = 0.
312    pbl_tke = 1.E-8
313
314    ! variables supplementaires pour appel a plb_surface_init
315    fder(:) = 0.
316    seaice(:) = 0.
317    run_off_lic_0 = 0.
318    fevap = 0.
319
320
321    ! Initialisations necessaires avant phyredem
322    type_ocean = 'force'
323    CALL fonte_neige_init(run_off_lic_0)
324    qsolsrf(:, :) = qsol(1) ! humidite du sol des sous surface
325    snsrf(:, :) = 0. ! couverture de neige des sous surface
326    z0m(:, :) = rugos ! couverture de neige des sous surface
327    z0h = z0m
328
329    CALL pbl_surface_init(fder, snsrf, qsolsrf, tsoil)
330
331    PRINT *, 'iniaqua: before phyredem'
332
333    pbl_tke(:, :, :) = 1.e-8
334    falb1 = albedo
335    falb2 = albedo
336    zmax0 = 0.
337    f0 = 0.
338    sig1 = 0.
339    w01 = 0.
340    wake_deltat = 0.
341    wake_deltaq = 0.
342    wake_s = 0.
343    wake_dens = 0.
344    wake_cstar = 0.
345    wake_pe = 0.
346    wake_fip = 0.
347    fm_therm = 0.
348    entr_therm = 0.
349    detr_therm = 0.
350    ale_bl = 0.
351    ale_bl_trig = 0.
352    alp_bl = 0.
353    treedrg(:, :, :) = 0.
354
355    u10m = 0.
356    v10m = 0.
357
358    ql_ancien = 0.
359    qs_ancien = 0.
360    qbs_ancien = 0.
361    u_ancien = 0.
362    v_ancien = 0.
363    prw_ancien = 0.
364    prlw_ancien = 0.
365    prsw_ancien = 0.
366    prbsw_ancien = 0.
367
368    ale_wake = 0.
369    ale_bl_stat = 0.
370
371
372    !ym error : the sub surface dimension is the third not second : forgotten for iniaqua
373    !    falb_dir(:,is_ter,:)=0.08; falb_dir(:,is_lic,:)=0.6
374    !    falb_dir(:,is_oce,:)=0.5;  falb_dir(:,is_sic,:)=0.6
375    falb_dir(:, :, is_ter) = 0.08; falb_dir(:, :, is_lic) = 0.6
376    falb_dir(:, :, is_oce) = 0.5;  falb_dir(:, :, is_sic) = 0.6
377
378    !ym falb_dif has been forgotten, initialize with defaukt value found in phyetat0 or 0 ?
379    !ym probably the uninitialized value was 0 for standard (regular grid) case
380    falb_dif(:, :, :) = 0
381
382    CALL phyredem('startphy.nc')
383
384    PRINT *, 'iniaqua: after phyredem'
385    CALL phys_state_var_end
386
387  END SUBROUTINE iniaqua
388
389
390  ! ====================================================================
391  ! ====================================================================
392  SUBROUTINE zenang_an(cycle_diurne, gmtime, rlat, rlon, rmu0, fract)
393    USE dimphy
394    IMPLICIT NONE
395    ! ====================================================================
396    ! =============================================================
397    ! CALL zenang(cycle_diurne,gmtime,rlat,rlon,rmu0,fract)
398    ! Auteur : A. Campoy et F. Hourdin
399    ! Objet  : calculer les valeurs moyennes du cos de l'angle zenithal
400    ! et l'ensoleillement moyen entre gmtime1 et gmtime2
401    ! connaissant la declinaison, la latitude et la longitude.
402
403    ! Dans cette version particuliere, on calcule le rayonnement
404    ! moyen sur l'année à chaque latitude.
405    ! angle zenithal calculé pour obtenir un
406    ! Fit polynomial de  l'ensoleillement moyen au sommet de l'atmosphere
407    ! en moyenne annuelle.
408    ! Spécifique de la terre. Utilisé pour les aqua planetes.
409
410    ! Rque   : Different de la routine angle en ce sens que zenang
411    ! fournit des moyennes de pmu0 et non des valeurs
412    ! instantanees, du coup frac prend toutes les valeurs
413    ! entre 0 et 1.
414    ! Date   : premiere version le 13 decembre 1994
415    ! revu pour  GCM  le 30 septembre 1996
416    ! ===============================================================
417    ! longi----INPUT : la longitude vraie de la terre dans son plan
418    ! solaire a partir de l'equinoxe de printemps (degre)
419    ! gmtime---INPUT : temps universel en fraction de jour
420    ! pdtrad---INPUT : pas de temps du rayonnement (secondes)
421    ! lat------INPUT : latitude en degres
422    ! long-----INPUT : longitude en degres
423    ! pmu0-----OUTPUT: angle zenithal moyen entre gmtime et gmtime+pdtrad
424    ! frac-----OUTPUT: ensoleillement moyen entre gmtime et gmtime+pdtrad
425    ! ================================================================
426    include "YOMCST.h"
427    ! ================================================================
428    LOGICAL cycle_diurne
429    REAL gmtime
430    REAL rlat(klon), rlon(klon), rmu0(klon), fract(klon)
431    ! ================================================================
432    INTEGER i
433    REAL gmtime1, gmtime2
434    REAL pi_local
435
436    REAL rmu0m(klon), rmu0a(klon)
437
438    pi_local = 4.0 * atan(1.0)
439
440    ! ================================================================
441    ! Calcul de l'angle zenithal moyen sur la journee
442    ! ================================================================
443
444    DO i = 1, klon
445      fract(i) = 1.
446      ! Calcule du flux moyen
447      IF (abs(rlat(i))<=28.75) THEN
448        rmu0m(i) = (210.1924 + 206.6059 * cos(0.0174533 * rlat(i))**2) / 1365.
449      ELSE IF (abs(rlat(i))<=43.75) THEN
450        rmu0m(i) = (187.4562 + 236.1853 * cos(0.0174533 * rlat(i))**2) / 1365.
451      ELSE IF (abs(rlat(i))<=71.25) THEN
452        rmu0m(i) = (162.4439 + 284.1192 * cos(0.0174533 * rlat(i))**2) / 1365.
453      ELSE
454        rmu0m(i) = (172.8125 + 183.7673 * cos(0.0174533 * rlat(i))**2) / 1365.
455      END IF
456    END DO
457
458    ! ================================================================
459    ! Avec ou sans cycle diurne
460    ! ================================================================
461
462    IF (cycle_diurne) THEN
463
464      ! On redecompose flux  au sommet suivant un cycle diurne idealise
465      ! identique a toutes les latitudes.
466
467      DO i = 1, klon
468        rmu0a(i) = 2. * rmu0m(i) * sqrt(2.) * pi_local / (4. - pi_local)
469        rmu0(i) = rmu0a(i) * abs(sin(pi_local * gmtime + pi_local * rlon(i) / 360.)) - &
470                rmu0a(i) / sqrt(2.)
471      END DO
472
473      DO i = 1, klon
474        IF (rmu0(i)<=0.) THEN
475          rmu0(i) = 0.
476          fract(i) = 0.
477        ELSE
478          fract(i) = 1.
479        END IF
480      END DO
481
482      ! Affichage de l'angel zenitale
483      ! PRINT*,'************************************'
484      ! PRINT*,'************************************'
485      ! PRINT*,'************************************'
486      ! PRINT*,'latitude=',rlat(i),'longitude=',rlon(i)
487      ! PRINT*,'rmu0m=',rmu0m(i)
488      ! PRINT*,'rmu0a=',rmu0a(i)
489      ! PRINT*,'rmu0=',rmu0(i)
490
491    ELSE
492
493      DO i = 1, klon
494        fract(i) = 0.5
495        rmu0(i) = rmu0m(i) * 2.
496      END DO
497
498    END IF
499
500  END SUBROUTINE zenang_an
501
502  ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
503
504  SUBROUTINE writelim_unstruct(klon, phy_nat, phy_alb, phy_sst, phy_bil, phy_rug, &
505          phy_ice, phy_fter, phy_foce, phy_flic, phy_fsic)
506
507    USE lmdz_phys_para, ONLY: is_omp_master, klon_mpi
508    USE lmdz_phys_transfert_para, ONLY: gather_omp
509    USE lmdz_xios
510    IMPLICIT NONE
511
512    INTEGER, INTENT (IN) :: klon
513    REAL, INTENT (IN) :: phy_nat(klon, 360)
514    REAL, INTENT (IN) :: phy_alb(klon, 360)
515    REAL, INTENT (IN) :: phy_sst(klon, 360)
516    REAL, INTENT (IN) :: phy_bil(klon, 360)
517    REAL, INTENT (IN) :: phy_rug(klon, 360)
518    REAL, INTENT (IN) :: phy_ice(klon, 360)
519    REAL, INTENT (IN) :: phy_fter(klon, 360)
520    REAL, INTENT (IN) :: phy_foce(klon, 360)
521    REAL, INTENT (IN) :: phy_flic(klon, 360)
522    REAL, INTENT (IN) :: phy_fsic(klon, 360)
523
524    REAL :: phy_mpi(klon_mpi, 360) ! temporary variable, to store phy_***(:)
525    ! on the whole physics grid
526
527    IF (using_xios) THEN
528      PRINT *, 'writelim: Ecriture du fichier limit'
529
530      CALL gather_omp(phy_foce, phy_mpi)
531      IF (is_omp_master) CALL xios_send_field('foce_limout', phy_mpi)
532
533      CALL gather_omp(phy_fsic, phy_mpi)
534      IF (is_omp_master) CALL xios_send_field('fsic_limout', phy_mpi)
535
536      CALL gather_omp(phy_fter, phy_mpi)
537      IF (is_omp_master) CALL xios_send_field('fter_limout', phy_mpi)
538
539      CALL gather_omp(phy_flic, phy_mpi)
540      IF (is_omp_master) CALL xios_send_field('flic_limout', phy_mpi)
541
542      CALL gather_omp(phy_sst, phy_mpi)
543      IF (is_omp_master) CALL xios_send_field('sst_limout', phy_mpi)
544
545      CALL gather_omp(phy_bil, phy_mpi)
546      IF (is_omp_master) CALL xios_send_field('bils_limout', phy_mpi)
547
548      CALL gather_omp(phy_alb, phy_mpi)
549      IF (is_omp_master) CALL xios_send_field('alb_limout', phy_mpi)
550
551      CALL gather_omp(phy_rug, phy_mpi)
552      IF (is_omp_master) CALL xios_send_field('rug_limout', phy_mpi)
553    ENDIF
554  END SUBROUTINE writelim_unstruct
555
556
557  SUBROUTINE writelim(klon, phy_nat, phy_alb, phy_sst, phy_bil, phy_rug, &
558          phy_ice, phy_fter, phy_foce, phy_flic, phy_fsic)
559
560    USE lmdz_phys_para, ONLY: is_master
561    USE lmdz_grid_phy, ONLY: klon_glo
562    USE lmdz_phys_transfert_para, ONLY: gather
563    USE phys_cal_mod, ONLY: year_len
564    USE netcdf, ONLY: nf90_def_var, nf90_put_var, nf90_get_var, nf90_strerror, nf90_close, &
565            nf90_enddef, nf90_put_att, nf90_unlimited, nf90_noerr, nf90_global, nf90_clobber, &
566            nf90_64bit_offset, nf90_def_dim, nf90_create
567    USE lmdz_cppkeys_wrapper, ONLY: nf90_format
568    IMPLICIT NONE
569
570    INTEGER, INTENT (IN) :: klon
571    REAL, INTENT (IN) :: phy_nat(klon, year_len)
572    REAL, INTENT (IN) :: phy_alb(klon, year_len)
573    REAL, INTENT (IN) :: phy_sst(klon, year_len)
574    REAL, INTENT (IN) :: phy_bil(klon, year_len)
575    REAL, INTENT (IN) :: phy_rug(klon, year_len)
576    REAL, INTENT (IN) :: phy_ice(klon, year_len)
577    REAL, INTENT (IN) :: phy_fter(klon, year_len)
578    REAL, INTENT (IN) :: phy_foce(klon, year_len)
579    REAL, INTENT (IN) :: phy_flic(klon, year_len)
580    REAL, INTENT (IN) :: phy_fsic(klon, year_len)
581
582    REAL :: phy_glo(klon_glo, year_len) ! temporary variable, to store phy_***(:)
583    ! on the whole physics grid
584    INTEGER :: k
585    INTEGER ierr
586    INTEGER dimfirst(3)
587    INTEGER dimlast(3)
588
589    INTEGER nid, ndim, ntim
590    INTEGER dims(2), debut(2), epais(2)
591    INTEGER id_tim
592    INTEGER id_nat, id_sst, id_bils, id_rug, id_alb
593    INTEGER id_fter, id_foce, id_fsic, id_flic
594
595    IF (is_master) THEN
596
597      PRINT *, 'writelim: Ecriture du fichier limit'
598
599      ierr = nf90_create('limit.nc', IOR(nf90_clobber, nf90_64bit_offset), nid)
600
601      ierr = nf90_put_att(nid, nf90_global, 'title', &
602              'Fichier conditions aux limites')
603      !        ierr = nf90_def_dim (nid, "points_physiques", klon, ndim)
604      ierr = nf90_def_dim(nid, 'points_physiques', klon_glo, ndim)
605      ierr = nf90_def_dim(nid, 'time', nf90_unlimited, ntim)
606
607      dims(1) = ndim
608      dims(2) = ntim
609
610      ierr = nf90_def_var(nid, 'TEMPS', nf90_format, [ntim], id_tim)
611      ierr = nf90_put_att(nid, id_tim, 'title', 'Jour dans l annee')
612
613      ierr = nf90_def_var(nid, 'NAT', nf90_format, dims, id_nat)
614      ierr = nf90_put_att(nid, id_nat, 'title', &
615              'Nature du sol (0,1,2,3)')
616
617      ierr = nf90_def_var(nid, 'SST', nf90_format, dims, id_sst)
618      ierr = nf90_put_att(nid, id_sst, 'title', &
619              'Temperature superficielle de la mer')
620
621      ierr = nf90_def_var(nid, 'BILS', nf90_format, dims, id_bils)
622      ierr = nf90_put_att(nid, id_bils, 'title', &
623              'Reference flux de chaleur au sol')
624
625      ierr = nf90_def_var(nid, 'ALB', nf90_format, dims, id_alb)
626      ierr = nf90_put_att(nid, id_alb, 'title', 'Albedo a la surface')
627
628      ierr = nf90_def_var(nid, 'RUG', nf90_format, dims, id_rug)
629      ierr = nf90_put_att(nid, id_rug, 'title', 'Rugosite')
630
631      ierr = nf90_def_var(nid, 'FTER', nf90_format, dims, id_fter)
632      ierr = nf90_put_att(nid, id_fter, 'title', 'Frac. Land')
633      ierr = nf90_def_var(nid, 'FOCE', nf90_format, dims, id_foce)
634      ierr = nf90_put_att(nid, id_foce, 'title', 'Frac. Ocean')
635      ierr = nf90_def_var(nid, 'FSIC', nf90_format, dims, id_fsic)
636      ierr = nf90_put_att(nid, id_fsic, 'title', 'Frac. Sea Ice')
637      ierr = nf90_def_var(nid, 'FLIC', nf90_format, dims, id_flic)
638      ierr = nf90_put_att(nid, id_flic, 'title', 'Frac. Land Ice')
639
640      ierr = nf90_enddef(nid)
641      IF (ierr/=nf90_noerr) THEN
642        WRITE (*, *) 'writelim error: failed to end define mode'
643        WRITE (*, *) nf90_strerror(ierr)
644      END IF
645
646
647      ! write the 'times'
648      DO k = 1, year_len
649        ierr = nf90_put_var(nid, id_tim, k, [k])
650        IF (ierr/=nf90_noerr) THEN
651          WRITE (*, *) 'writelim error with temps(k),k=', k
652          WRITE (*, *) nf90_strerror(ierr)
653        END IF
654      END DO
655
656    END IF ! of if (is_master)
657
658    ! write the fields, after having collected them on master
659
660    CALL gather(phy_nat, phy_glo)
661    IF (is_master) THEN
662      ierr = nf90_put_var(nid, id_nat, phy_glo)
663      IF (ierr/=nf90_noerr) THEN
664        WRITE (*, *) 'writelim error with phy_nat'
665        WRITE (*, *) nf90_strerror(ierr)
666      END IF
667    END IF
668
669    CALL gather(phy_sst, phy_glo)
670    IF (is_master) THEN
671      ierr = nf90_put_var(nid, id_sst, phy_glo)
672      IF (ierr/=nf90_noerr) THEN
673        WRITE (*, *) 'writelim error with phy_sst'
674        WRITE (*, *) nf90_strerror(ierr)
675      END IF
676    END IF
677
678    CALL gather(phy_bil, phy_glo)
679    IF (is_master) THEN
680      ierr = nf90_put_var(nid, id_bils, phy_glo)
681      IF (ierr/=nf90_noerr) THEN
682        WRITE (*, *) 'writelim error with phy_bil'
683        WRITE (*, *) nf90_strerror(ierr)
684      END IF
685    END IF
686
687    CALL gather(phy_alb, phy_glo)
688    IF (is_master) THEN
689      ierr = nf90_put_var(nid, id_alb, phy_glo)
690      IF (ierr/=nf90_noerr) THEN
691        WRITE (*, *) 'writelim error with phy_alb'
692        WRITE (*, *) nf90_strerror(ierr)
693      END IF
694    END IF
695
696    CALL gather(phy_rug, phy_glo)
697    IF (is_master) THEN
698      ierr = nf90_put_var(nid, id_rug, phy_glo)
699      IF (ierr/=nf90_noerr) THEN
700        WRITE (*, *) 'writelim error with phy_rug'
701        WRITE (*, *) nf90_strerror(ierr)
702      END IF
703    END IF
704
705    CALL gather(phy_fter, phy_glo)
706    IF (is_master) THEN
707      ierr = nf90_put_var(nid, id_fter, phy_glo)
708      IF (ierr/=nf90_noerr) THEN
709        WRITE (*, *) 'writelim error with phy_fter'
710        WRITE (*, *) nf90_strerror(ierr)
711      END IF
712    END IF
713
714    CALL gather(phy_foce, phy_glo)
715    IF (is_master) THEN
716      ierr = nf90_put_var(nid, id_foce, phy_glo)
717      IF (ierr/=nf90_noerr) THEN
718        WRITE (*, *) 'writelim error with phy_foce'
719        WRITE (*, *) nf90_strerror(ierr)
720      END IF
721    END IF
722
723    CALL gather(phy_fsic, phy_glo)
724    IF (is_master) THEN
725      ierr = nf90_put_var(nid, id_fsic, phy_glo)
726      IF (ierr/=nf90_noerr) THEN
727        WRITE (*, *) 'writelim error with phy_fsic'
728        WRITE (*, *) nf90_strerror(ierr)
729      END IF
730    END IF
731
732    CALL gather(phy_flic, phy_glo)
733    IF (is_master) THEN
734      ierr = nf90_put_var(nid, id_flic, phy_glo)
735      IF (ierr/=nf90_noerr) THEN
736        WRITE (*, *) 'writelim error with phy_flic'
737        WRITE (*, *) nf90_strerror(ierr)
738      END IF
739    END IF
740
741    ! close file:
742    IF (is_master) THEN
743      ierr = nf90_close(nid)
744    END IF
745
746  END SUBROUTINE writelim
747
748  ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
749
750  SUBROUTINE profil_sst(nlon, rlatd, type_profil, phy_sst)
751    USE dimphy
752    USE phys_cal_mod, ONLY: year_len
753    IMPLICIT NONE
754
755    INTEGER nlon, type_profil, i, k, j
756    REAL :: rlatd(nlon), phy_sst(nlon, year_len)
757    INTEGER imn, imx, amn, amx, kmn, kmx
758    INTEGER p, pplus, nlat_max
759    PARAMETER (nlat_max = 72)
760    REAL x_anom_sst(nlat_max)
761    CHARACTER (LEN = 20) :: modname = 'profil_sst'
762    CHARACTER (LEN = 80) :: abort_message
763
764    IF (klon/=nlon) THEN
765      abort_message = 'probleme de dimensions dans profil_sst'
766      CALL abort_physic(modname, abort_message, 1)
767    ENDIF
768    WRITE (*, *) ' profil_sst: type_profil=', type_profil
769    DO i = 1, year_len
770      ! phy_sst(:,i) = 260.+50.*cos(rlatd(:))**2
771
772      ! Rajout fbrlmd
773
774      IF (type_profil==1) THEN
775        ! Méthode 1 "Control" faible plateau à l'Equateur
776        DO j = 1, klon
777          phy_sst(j, i) = 273. + 27. * (1 - sin(1.5 * rlatd(j))**2)
778          ! PI/3=1.047197551
779          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
780            phy_sst(j, i) = 273.
781          END IF
782        END DO
783      END IF
784      IF (type_profil==2) THEN
785        ! Méthode 2 "Flat" fort plateau à l'Equateur
786        DO j = 1, klon
787          phy_sst(j, i) = 273. + 27. * (1 - sin(1.5 * rlatd(j))**4)
788          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
789            phy_sst(j, i) = 273.
790          END IF
791        END DO
792      END IF
793
794      IF (type_profil==3) THEN
795        ! Méthode 3 "Qobs" plateau réel à l'Equateur
796        DO j = 1, klon
797          phy_sst(j, i) = 273. + 0.5 * 27. * (2 - sin(1.5 * rlatd(j))**2 - sin(1.5 * &
798                  rlatd(j))**4)
799          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
800            phy_sst(j, i) = 273.
801          END IF
802        END DO
803      END IF
804
805      IF (type_profil==4) THEN
806        ! Méthode 4 : Méthode 3 + SST+2 "Qobs" plateau réel à l'Equateur
807        DO j = 1, klon
808          phy_sst(j, i) = 273. + 0.5 * 29. * (2 - sin(1.5 * rlatd(j))**2 - sin(1.5 * &
809                  rlatd(j))**4)
810          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
811            phy_sst(j, i) = 273.
812          END IF
813        END DO
814      END IF
815
816      IF (type_profil==5) THEN
817        ! Méthode 5 : Méthode 3 + +2K "Qobs" plateau réel à l'Equateur
818        DO j = 1, klon
819          phy_sst(j, i) = 273. + 2. + 0.5 * 27. * (2 - sin(1.5 * rlatd(j))**2 - sin(1.5 &
820                  * rlatd(j))**4)
821          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
822            phy_sst(j, i) = 273. + 2.
823          END IF
824
825        END DO
826      END IF
827
828      IF (type_profil==6) THEN
829        ! Méthode 6 "cst" valeur constante de SST
830        DO j = 1, klon
831          phy_sst(j, i) = 288.
832        END DO
833      END IF
834
835      IF (type_profil==7) THEN
836        ! Méthode 7 "cst" valeur constante de SST +2
837        DO j = 1, klon
838          phy_sst(j, i) = 288. + 2.
839        END DO
840      END IF
841
842      p = 0
843      IF (type_profil==8) THEN
844        ! Méthode 8 profil anomalies SST du modèle couplé AR4
845        DO j = 1, klon
846          IF (rlatd(j)==rlatd(j - 1)) THEN
847            phy_sst(j, i) = 273. + x_anom_sst(pplus) + &
848                    0.5 * 27. * (2 - sin(1.5 * rlatd(j))**2 - sin(1.5 * rlatd(j))**4)
849            IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
850              phy_sst(j, i) = 273. + x_anom_sst(pplus)
851            END IF
852          ELSE
853            p = p + 1
854            pplus = 73 - p
855            phy_sst(j, i) = 273. + x_anom_sst(pplus) + &
856                    0.5 * 27. * (2 - sin(1.5 * rlatd(j))**2 - sin(1.5 * rlatd(j))**4)
857            IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
858              phy_sst(j, i) = 273. + x_anom_sst(pplus)
859            END IF
860            WRITE (*, *) rlatd(j), x_anom_sst(pplus), phy_sst(j, i)
861          END IF
862        END DO
863      END IF
864
865      IF (type_profil==9) THEN
866        ! Méthode 5 : Méthode 3 + -2K "Qobs" plateau réel à l'Equateur
867        DO j = 1, klon
868          phy_sst(j, i) = 273. - 2. + 0.5 * 27. * (2 - sin(1.5 * rlatd(j))**2 - sin(1.5 &
869                  * rlatd(j))**4)
870          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
871            phy_sst(j, i) = 273. - 2.
872          END IF
873        END DO
874      END IF
875
876      IF (type_profil==10) THEN
877        ! Méthode 10 : Méthode 3 + +4K "Qobs" plateau réel à l'Equateur
878        DO j = 1, klon
879          phy_sst(j, i) = 273. + 4. + 0.5 * 27. * (2 - sin(1.5 * rlatd(j))**2 - sin(1.5 &
880                  * rlatd(j))**4)
881          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
882            phy_sst(j, i) = 273. + 4.
883          END IF
884        END DO
885      END IF
886
887      IF (type_profil==11) THEN
888        ! Méthode 11 : Méthode 3 + 4CO2 "Qobs" plateau réel à l'Equateur
889        DO j = 1, klon
890          phy_sst(j, i) = 273. + 0.5 * 27. * (2 - sin(1.5 * rlatd(j))**2 - sin(1.5 * &
891                  rlatd(j))**4)
892          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
893            phy_sst(j, i) = 273.
894          END IF
895        END DO
896      END IF
897
898      IF (type_profil==12) THEN
899        ! Méthode 12 : Méthode 10 + 4CO2 "Qobs" plateau réel à l'Equateur
900        DO j = 1, klon
901          phy_sst(j, i) = 273. + 4. + 0.5 * 27. * (2 - sin(1.5 * rlatd(j))**2 - sin(1.5 &
902                  * rlatd(j))**4)
903          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
904            phy_sst(j, i) = 273. + 4.
905          END IF
906        END DO
907      END IF
908
909      IF (type_profil==13) THEN
910        ! Méthode 13 "Qmax" plateau réel à l'Equateur augmenté !
911        DO j = 1, klon
912          phy_sst(j, i) = 273. + 0.5 * 29. * (2 - sin(1.5 * rlatd(j))**2 - sin(1.5 * &
913                  rlatd(j))**4)
914          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
915            phy_sst(j, i) = 273.
916          END IF
917        END DO
918      END IF
919
920      IF (type_profil==14) THEN
921        ! Méthode 13 "Qmax2K" plateau réel à l'Equateur augmenté +2K !
922        DO j = 1, klon
923          phy_sst(j, i) = 273. + 2. + 0.5 * 29. * (2 - sin(1.5 * rlatd(j))**2 - sin(1.5 &
924                  * rlatd(j))**4)
925          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
926            phy_sst(j, i) = 273.
927          END IF
928        END DO
929      END IF
930
931      IF (type_profil==20) THEN
932        PRINT*, 'Profile SST 20'
933        !     Méthode 13 "Qmax2K" plateau réel é|  l'Equateur augmenté +2K
934
935        do j = 1, klon
936          phy_sst(j, i) = 248. + 55. * (1 - sin(rlatd(j))**2)
937        enddo
938      endif
939
940      IF (type_profil==21) THEN
941        PRINT*, 'Profile SST 21'
942        !     Méthode 13 "Qmax2K" plateau réel é|  l'Equateur augmenté +2K
943        do j = 1, klon
944          phy_sst(j, i) = 252. + 55. * (1 - sin(rlatd(j))**2)
945        enddo
946      endif
947
948    END DO
949
950    ! IM beg : verif profil SST: phy_sst
951    amn = min(phy_sst(1, 1), 1000.)
952    amx = max(phy_sst(1, 1), -1000.)
953    imn = 1
954    kmn = 1
955    imx = 1
956    kmx = 1
957    DO k = 1, year_len
958      DO i = 2, nlon
959        IF (phy_sst(i, k)<amn) THEN
960          amn = phy_sst(i, k)
961          imn = i
962          kmn = k
963        END IF
964        IF (phy_sst(i, k)>amx) THEN
965          amx = phy_sst(i, k)
966          imx = i
967          kmx = k
968        END IF
969      END DO
970    END DO
971
972    PRINT *, 'profil_sst: imn, kmn, phy_sst(imn,kmn) ', imn, kmn, amn
973    PRINT *, 'profil_sst: imx, kmx, phy_sst(imx,kmx) ', imx, kmx, amx
974    ! IM end : verif profil SST: phy_sst
975
976  END SUBROUTINE profil_sst
977
978END MODULE phyaqua_mod
Note: See TracBrowser for help on using the repository browser.